This Ampere's Law should not work. The field is electrostatic and we should treat it as such

src/modules/mesh/moduleMesh.f90

@@ -1119,6 +1119,7 @@ MODULE moduleMesh
     type(meshEdgePointer), allocatable:: edges(:) ! Edges included in the boundary
     real(8), allocatable:: electricField(:) ! Electric field normal to the edge that must be applied with a Neumann boundary condition
     real(8), allocatable:: surfaceCharge(:) ! Surface charge density
+    integer:: counter
     contains
       procedure, pass:: apply => applyFreeCurrent
 

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -322,6 +322,8 @@ submodule(moduleMesh) boundaryEM
       type is(boundaryEMFreeCurrent)
         allocate(self%edges(0))
 
+        self%counter = 0
+
         self%update => updateFreeCurrent
         self%print  => writeFreeCurrent
 
@@ -376,6 +378,7 @@ submodule(moduleMesh) boundaryEM
 
       select type(self)
       type is(boundaryEMFreeCurrent)
+        self%counter = self%counter + 1
         do e=1, self%nEdges
           edgeDensityCurrent = 0.0d0
 
@@ -402,12 +405,20 @@ submodule(moduleMesh) boundaryEM
 
           end associate
 
-          self%surfaceCharge(e) =   self%surfaceCharge(e) + edgeDensityCurrent * tauMin
+          self%surfaceCharge(e) = self%surfaceCharge(e) + edgeDensityCurrent * tauMin
 
-          self%electricField(e) = - self%surfaceCharge(e)
+          deallocate(nodes, mom_nodes)
 
         end do
 
+        if (self%counter == 10) then
+          self%counter = 0
+
+          self%electricField = self%electricField - self%surfaceCharge / 10.d0
+          self%surfaceCharge = 0.0d0
+
+        end if
+
       end select
 
     end subroutine updateFreeCurrent