Implementation done, but values of E are too high. Too much charge being accumulated

src/modules/init/moduleInput.f90

@@ -972,6 +972,11 @@ MODULE moduleInput
 
             call initFloating(self, config, object)
 
+          case ("freeCurrent")
+            allocate(boundaryEMFreeCurrent:: self)
+
+            call initFreeCurrent(self, config, object)
+
           case default
             call criticalError('Boundary type ' // bType // ' not supported', 'readBoundaryEM')
 
@@ -1334,8 +1339,25 @@ MODULE moduleInput
             ! If the boundary for the species is linked to the one analysing, add the edges
             if (associated(physicalSurfaces(ps)%EM, bound)) then
               bound%nodes = [bound%nodes, physicalSurfaces(ps)%nodes]
+
             end if
 
+            ! Specific assignments per type
+            select type(bound)
+            type is(boundaryEMFreeCurrent)
+              if (associated(physicalSurfaces(ps)%EM, bound)) then
+                bound%edges = [bound%edges, physicalSurfaces(ps)%edges]
+
+                bound%nEdges = size(bound%edges)
+                allocate(bound%electricField(1:bound%nEdges))
+                bound%electricField = 0.0d0
+                allocate(bound%surfaceCharge(1:bound%nEdges))
+                bound%surfaceCharge = 0.0d0
+
+              end if
+
+            end select
+
           end do
 
           bound%nNodes = size(bound%nodes)

src/modules/mesh/moduleMesh.f90

@@ -1115,7 +1115,10 @@ MODULE moduleMesh
 
   ! Calculates the electric field normal to the surface based on Ampere's law (without magnetic field)
   type, extends(boundaryEMGeneric):: boundaryEMFreeCurrent
-    type(meshEdgePointer), allocatable:: electricField(:) ! Electric field in each edge
+    integer:: nEdges
+    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
     contains
       procedure, pass:: apply => applyFreeCurrent
 

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -238,7 +238,6 @@ submodule(moduleMesh) boundaryEM
 
     ! Update
     subroutine updateFloating(self)
-      use moduleMesh, only: qSpecies, meshNode
       use moduleCaseParam, only: tauMin
       implicit none
 
@@ -309,6 +308,133 @@ submodule(moduleMesh) boundaryEM
 
     end subroutine writeFloating
 
+    ! Free current
+    ! Init
+    module subroutine initFreeCurrent(self, config, object)
+      use json_module
+      implicit none
+
+      class(boundaryEMGeneric), allocatable, intent(inout):: self
+      type(json_file), intent(inout):: config
+      character(:), allocatable, intent(in):: object
+
+      select type(self)
+      type is(boundaryEMFreeCurrent)
+        allocate(self%edges(0))
+
+        self%update => updateFreeCurrent
+        self%print  => writeFreeCurrent
+
+      end select
+
+    end subroutine initFreeCurrent
+
+    ! Apply
+    module subroutine applyFreeCurrent(self, vectorF)
+      implicit none
+
+      class(boundaryEMFreeCurrent), intent(in):: self
+      real(8), intent(inout):: vectorF(:)
+      integer:: e, n
+      integer, allocatable:: nodes(:)
+      real(8), allocatable:: fPsi(:)
+
+      do e = 1, self%nEdges
+        associate(edge => self%edges(e)%obj)
+          fPsi = edge%fPsi(edge%centerXi(), edge%nNodes)
+
+          nodes = edge%getNodes(edge%nNodes)
+
+          do n = 1, edge%nNodes
+            associate(node => mesh%nodes(nodes(n))%obj)
+              ! Assign to each node the corresponding weight of the Neumann BC
+              vectorF(node%n) = vectorF(node%n) - fPsi(n) * self%electricField(e)/node%v
+
+            end associate
+
+          end do
+
+        end associate
+
+        deallocate(fPsi, nodes)
+
+      end do
+
+    end subroutine applyFreeCurrent
+
+    ! Update
+    subroutine updateFreeCurrent(self)
+      use moduleCaseParam, only: tauMin
+      implicit none
+
+      class(boundaryEMGeneric), intent(inout):: self
+      integer:: e, n, s
+      integer, allocatable:: nodes(:)
+      real(8), allocatable:: mom_nodes(:)
+      class(meshNode), pointer:: node
+      real(8):: mom_center, edgeDensityCurrent
+
+      select type(self)
+      type is(boundaryEMFreeCurrent)
+        do e=1, self%nEdges
+          edgeDensityCurrent = 0.0d0
+
+          associate(edge => self%edges(e)%obj)
+
+            ! Gather the current density at the edge
+            nodes = edge%getNodes(edge%nNodes)
+            allocate(mom_nodes(1:edge%nNodes))
+            do s = 1, nSpecies
+              mom_center = 0.0d0
+              do n = 1, self%nNodes
+                node => mesh%nodes(nodes(n))%obj
+
+                ! Minus sign to get the flux exiting the edge
+                mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)
+
+              end do
+
+              mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
+
+              edgeDensityCurrent = edgeDensityCurrent + qSpecies(s) * mom_center
+
+            end do
+
+          end associate
+
+          self%surfaceCharge(e) =   self%surfaceCharge(e) + edgeDensityCurrent * tauMin
+
+          self%electricField(e) = - self%surfaceCharge(e)
+
+        end do
+
+      end select
+
+    end subroutine updateFreeCurrent
+
+    ! Write
+    subroutine writeFreeCurrent(self, fileID)
+      ! use moduleOutput, only: fmtColReal
+      use moduleConstParam, only: qe, eps_0
+      use moduleRefParam, only: EF_ref, L_ref, n_ref, v_ref, ti_ref
+      implicit none
+
+      class(boundaryEMGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+      integer:: e
+
+      write(fileID, '(A)') self%name
+      select type(self)
+      type is(boundaryEMFreeCurrent)
+        do e = 1, self%nEdges
+          print*, self%edges(e)%obj%n, self%electricField(e)*EF_ref, self%surfaceCharge(e)*qe*n_ref*v_ref*ti_ref
+
+        end do
+
+      end select
+
+    end subroutine writeFreeCurrent
+
     ! Get the index of the boundary based on the name
     module function boundaryEMName_to_Index(boundaryName) result(bp)
       use moduleErrors