Adjusted initial condition, dimensions of output, and output file header

doc/user-manual/fpakc_UserManual.tex

@@ -626,6 +626,15 @@ make
                     File must be located in \textbf{output.path}.
                     The first column is the time in $\unit{s}$.
                     The second column is the factor that will multiply the value set in \textbf{potential}.
+                    
+                 \item \textbf{floating}: The potential is adjusted over time to maintain a zero current on the surface, \textit{i. e.}, the surface is \textit{floating}.                                   
+                   \textbf{potential}: Real.
+                     Initial potential of the surface.
+                     The initial charge value for the simulation will be calculated by the potential times the capacitance.
+                   
+                   \textbf{capacitance}: Real.
+                     Capacitance of the surface in Farads.
+                     If none is provided, a default value of $1$ nF is used.
                      
               \end{itemize}
           

src/fpakc.f90

@@ -25,7 +25,6 @@ PROGRAM fpakc
                              doAverage, &
                              doInjects
   use moduleMesh,      only: boundariesEM_update,       &
-                             boundariesEM_update,       &
                              boundariesParticle_update, &
                              mesh,                      &
                              meshForMCC,                &

src/modules/init/moduleInput.f90

@@ -493,6 +493,7 @@ MODULE moduleInput
       CALL config%get(object // '.numColl', collOutput, found)
       CALL config%get(object // '.EMField', emOutput, found)
       call config%get(object // '.boundariesParticle', boundaryParticleOutput, found)
+      call config%get(object // '.boundariesEM', boundaryEMOutput, found)
       call config%get(object // '.injects', injectOutput, found)
 
       CALL config%get(object // '.triggerCPUTime', triggerCPUTime, found)
@@ -969,17 +970,20 @@ MODULE moduleInput
 
           select case(bType)
           case ("dirichlet")
-            ! Allocate boundary edge
             allocate(boundaryEMDirichlet:: self)
 
             CALL initDirichlet(self, config, object)
 
           case ("dirichlettime")
-            ! Allocate boundary edge
             allocate(boundaryEMDirichletTime:: self)
 
             call initDirichletTime(self, config, object)
 
+          case ("floating")
+            allocate(boundaryEMFloating:: self)
+
+            call initFloating(self, config, object)
+
           case default
             call criticalError('Boundary type ' // bType // ' not supported', 'readBoundaryEM')
 
@@ -1358,6 +1362,18 @@ MODULE moduleInput
 
             end do
 
+          type is(boundaryEMFloating)
+            ! Loop over all physical surfaces
+            do ps = 1, nPhysicalSurfaces
+              ! 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]
+
+                bound%edges = [bound%edges, physicalSurfaces(ps)%edges]
+              end if
+
+            end do
+
           end select
 
           bound%nNodes = size(bound%nodes)

src/modules/mesh/moduleMesh.f90

@@ -963,7 +963,8 @@ MODULE moduleMesh
     integer:: nNodes
     type(meshNodePointer), allocatable:: nodes(:)
 
-    procedure(updateEM_interface), pointer, pass:: update => null()
+    procedure(updateEM_interface),        pointer, pass:: update => null()
+    procedure(printEMboundary_interface), pointer, pass:: print  => null()
     contains
       procedure(applyEM_interface), deferred, pass:: apply
 
@@ -987,9 +988,19 @@ MODULE moduleMesh
 
     end subroutine updateEM_interface
 
+    ! Write the values of the particle boundary model
+    subroutine printEMboundary_interface(self, fileID)
+      import boundaryEMGeneric
+
+      class(boundaryEMGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+
+    end subroutine printEMboundary_interface
+
   end interface
 
   ! Extended types
+  ! Dirichlet: Constant potential value
   TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet
     REAL(8):: potential
     CONTAINS
@@ -1016,6 +1027,7 @@ MODULE moduleMesh
 
   end interface
 
+  ! Dirichlet time: time-dependent potential value
   TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
     real(8):: potential
     real(8):: timeFactor
@@ -1044,6 +1056,36 @@ MODULE moduleMesh
 
   end interface
 
+  ! Floating: Floating surface, ie, zero net current
+  type, extends(boundaryEMGeneric):: boundaryEMFloating
+    real(8):: potential
+    real(8):: invC ! Inverse of the capacitance of the surface
+    real(8):: charge ! Charge accumulated on surface
+    type(meshEdgePointer), allocatable:: edges(:) !Array with edges
+    contains
+      ! boundaryEMGeneric DEFERRED PROCEDURES
+      procedure, pass:: apply => applyFloating
+
+  end type boundaryEMFloating
+
+  interface
+    module subroutine initFloating(self, config, object)
+      use json_module
+
+      class(boundaryEMGeneric), allocatable, intent(inout):: self
+      type(json_file), intent(inout):: config
+      character(:), allocatable, intent(in):: object
+
+    end subroutine initFloating
+
+    module subroutine applyFloating(self, vectorF)
+      class(boundaryEMFloating), intent(in):: self
+      real(8), intent(inout):: vectorF(:)
+
+    end subroutine applyFloating
+
+  end interface
+
   ! Container for boundary conditions
   TYPE:: boundaryEMCont
     CLASS(boundaryEMGeneric), ALLOCATABLE:: obj
@@ -1068,6 +1110,10 @@ MODULE moduleMesh
 
     end subroutine boundariesEM_update
 
+    module subroutine boundariesEM_write()
+
+    end subroutine boundariesEM_write
+
   end interface
 
   ! PHYSICAL SURFACES LINKING TO MESH ELEMENTS

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -1,30 +1,7 @@
 submodule(moduleMesh) boundaryEM
   CONTAINS
-    module function boundaryEMName_to_Index(boundaryName) result(bp)
-      use moduleErrors
-      implicit none
-
-      character(:), allocatable:: boundaryName
-      integer:: bp
-      integer:: b
-
-      bp = 0
-      do b = 1, nBoundariesEM
-        if (boundaryName == boundariesEM(b)%obj%name) then
-          bp = boundariesEM(b)%obj%n
-
-        end if
-
-      end do
-
-      if (bp == 0) then
-        call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryEMName_to_Index')
-      
-      end if
-    
-    end function boundaryEMName_to_Index
-
-    ! Initialize Dirichlet boundary condition
+    ! Dirichlet boundary condition
+    ! Init
     module SUBROUTINE initDirichlet(self, config, object)
       use json_module
       USE moduleRefParam, ONLY: Volt_ref
@@ -51,7 +28,25 @@ submodule(moduleMesh) boundaryEM
 
     end subroutine initDirichlet
 
-    ! Initialize Dirichlet boundary condition
+    ! Apply
+    module SUBROUTINE applyDirichlet(self, vectorF)
+      USE moduleMesh
+      IMPLICIT NONE
+
+      CLASS(boundaryEMDirichlet), INTENT(in):: self
+      REAL(8), INTENT(inout):: vectorF(:)
+      integer:: n
+
+      DO n = 1, self%nNodes
+        self%nodes(n)%obj%emData%phi  = self%potential
+        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
+
+      END DO
+
+    END SUBROUTINE applyDirichlet
+
+    ! Dirichlet time-dependent boundary condition
+    ! Init
     module subroutine initDirichletTime(self, config, object)
       use json_module
       use moduleRefParam, ONLY: Volt_ref, ti_ref
@@ -96,24 +91,7 @@ submodule(moduleMesh) boundaryEM
 
     END SUBROUTINE initDirichletTime
 
-    !Apply Dirichlet boundary condition to the poisson equation
-    module SUBROUTINE applyDirichlet(self, vectorF)
-      USE moduleMesh
-      IMPLICIT NONE
-
-      CLASS(boundaryEMDirichlet), INTENT(in):: self
-      REAL(8), INTENT(inout):: vectorF(:)
-      integer:: n
-
-      DO n = 1, self%nNodes
-        self%nodes(n)%obj%emData%phi  = self%potential
-        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
-
-      END DO
-
-    END SUBROUTINE applyDirichlet
-
-    !Apply Dirichlet boundary condition with time temporal profile
+    ! Apply
     module SUBROUTINE applyDirichletTime(self, vectorF)
       USE moduleMesh
       USE moduleCaseParam, ONLY: timeStep, tauMin
@@ -131,7 +109,8 @@ submodule(moduleMesh) boundaryEM
 
     END SUBROUTINE applyDirichletTime
 
-    module subroutine updateDirichletTime(self)
+    ! Update
+    subroutine updateDirichletTime(self)
       implicit none
 
       class(boundaryEMGeneric), intent(inout):: self
@@ -144,6 +123,170 @@ submodule(moduleMesh) boundaryEM
 
     end subroutine updateDirichletTime
 
+    ! Floating surface
+    ! Init
+    module subroutine initFloating(self, config, object)
+      use json_module
+      use moduleRefParam, only: Volt_ref, L_ref
+      use moduleConstParam, only: eps_0
+      use moduleErrors
+      implicit none
+
+      class(boundaryEMGeneric), allocatable, intent(inout):: self
+      type(json_file), intent(inout):: config
+      character(:), allocatable, intent(in):: object
+      real(8):: potential, capacitance
+      logical:: found
+
+      select type(self)
+      type is(boundaryEMFloating)
+        call config%get(object // '.potential', potential, found)
+        if (.not. found) then
+          call criticalError('Required parameter "potential" for Floating boundary condition not found', &
+                             'readBoundaryEM')
+
+        end if
+
+        self%potential = potential / Volt_ref
+
+        call config%get(object // '.capacitance', capacitance, found)
+        if (.not. found) then
+          call warningError('Parameter "capacitance" for floating boundary condition not found. Using the default value of 1 nF.')
+          capacitance = 1.0d-9
+
+        end if
+    
+        ! Inverse of non-dimensional capacitance
+        self%invC = 1.0d0 / (capacitance / (eps_0 * L_ref))
+
+        self%charge = self%potential / self%invC
+
+        self%update => updateFloating
+
+        self%print  => writeFloating
+
+        allocate(self%edges(0))
+
+      end select
+
+    end subroutine initFloating
+    
+    ! Apply
+    module subroutine applyFloating(self, vectorF)
+      use moduleMesh
+      implicit none
+
+      class(boundaryEMFloating), intent(in):: self
+      real(8), intent(inout):: vectorF(:)
+      integer:: n
+
+      do n = 1, self%nNodes
+        self%nodes(n)%obj%emData%phi  = self%potential
+        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
+
+      end do
+
+    end subroutine applyFloating
+
+    ! Update
+    subroutine updateFloating(self)
+      use moduleMesh, only: qSpecies, meshNode
+      use moduleCaseParam, only: tauMin
+      implicit none
+
+      class(boundaryEMGeneric), intent(inout):: self
+      integer:: e, n, s
+      integer, allocatable:: nodes(:)
+      class(meshEdge), pointer:: edge
+      real(8), allocatable:: mom_nodes(:)
+      class(meshNode), pointer:: node
+      real(8):: mom_center, edgeDensityCurrent
+
+      select type(self)
+      type is(boundaryEMFloating)
+        do e = 1, size(self%edges)
+          edge => self%edges(e)%obj
+
+          edgeDensityCurrent = 0.0d0
+
+          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)
+
+            ! Only account for charge exiting the surface
+            mom_center = max(mom_center, 0.0d0)
+            
+            edgeDensityCurrent = edgeDensityCurrent + qSpecies(s) * mom_center
+
+          end do
+
+          ! Accumulate charge of speceis on surface
+          self%charge = self%charge + edgeDensityCurrent * edge%surface * tauMin
+
+        end do
+
+        self%potential = self%charge * self%invC
+
+      end select
+
+    end subroutine updateFloating
+
+    ! Write
+    subroutine writeFloating(self, fileID)
+      use moduleOutput, only: fmtColReal
+      use moduleConstParam, only: qe
+      use moduleRefParam, only: Volt_ref, v_ref, n_ref, L_ref, ti_ref
+      implicit none
+
+      class(boundaryEMGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+
+      write(fileID, '(A)') self%name
+      select type(self)
+      type is(boundaryEMFloating)
+        write(fileID, '(A,",",A)') 'Total charge (C)', 'Potential (V)'
+        write(fileID, '(*('//fmtColReal//'))') self%charge*qe*v_ref*n_ref*L_ref**2*ti_ref, self%potential * Volt_ref
+
+      end select
+
+    end subroutine writeFloating
+
+    ! Get the index of the boundary based on the name
+    module function boundaryEMName_to_Index(boundaryName) result(bp)
+      use moduleErrors
+      implicit none
+
+      character(:), allocatable:: boundaryName
+      integer:: bp
+      integer:: b
+
+      bp = 0
+      do b = 1, nBoundariesEM
+        if (boundaryName == boundariesEM(b)%obj%name) then
+          bp = boundariesEM(b)%obj%n
+
+        end if
+
+      end do
+
+      if (bp == 0) then
+        call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryEMName_to_Index')
+      
+      end if
+    
+    end function boundaryEMName_to_Index
+
+    ! Update all EM boundaries
     module subroutine boundariesEM_update()
       implicit none
 
@@ -151,7 +294,7 @@ submodule(moduleMesh) boundaryEM
 
       do b = 1, nBoundariesEM
         if (associated(boundariesEM(b)%obj%update)) then
-          call boundariesEM(b)%obj%update
+          call boundariesEM(b)%obj%update()
 
         end if
 
@@ -159,4 +302,37 @@ submodule(moduleMesh) boundaryEM
 
     end subroutine boundariesEM_update
 
+    ! Write all EM boundaries
+    module subroutine boundariesEM_write()
+      use moduleCaseparam, only: timeStep
+      use moduleOutput, only: fileID_boundaryEM, &
+                              formatFileName, &
+                              informFileCreation,&
+                              boundaryEMOutput,&
+                              generateFilePath, &
+                              prefix
+      implicit none
+
+      integer:: b
+      character(:), allocatable:: fileName
+
+      if (boundaryEMOutput) then
+        fileName = formatFileName(prefix, 'boundariesEM', 'csv', timeStep)
+        call informFileCreation(fileName)
+        open(fileID_boundaryEM, file = generateFilePath(fileName))
+
+        do b = 1, nBoundariesEM
+          if (associated(boundariesEM(b)%obj%update)) then
+            call boundariesEM(b)%obj%print(fileID_boundaryEM)
+
+          end if
+
+        end do
+
+        close(fileID_boundaryEM)
+
+      end if
+
+    end subroutine boundariesEM_write
+
 end submodule boundaryEM

src/modules/mesh/moduleMesh@elements.f90

@@ -176,6 +176,14 @@ submodule(moduleMesh) elements
 
           END DO
 
+        TYPE IS(boundaryEMFloating)
+          DO n = 1, boundary%nNodes
+            ni = boundary%nodes(n)%obj%n
+            self%K(ni, :)  = 0.D0
+            self%K(ni, ni) = 1.D0
+
+          END DO
+
         END SELECT
 
       END DO

src/modules/output/moduleOutput.f90

@@ -14,6 +14,7 @@ MODULE moduleOutput
   LOGICAL:: collOutput = .FALSE.
   LOGICAL:: emOutput   = .FALSE.
   logical:: boundaryParticleOutput = .false.
+  logical:: boundaryEMOutput = .false.
   logical:: injectOutput = .false.
 
   ! Prefix for iteration files
@@ -34,7 +35,7 @@ MODULE moduleOutput
   integer, parameter:: fileID_output = 20           ! Base id for species/collisions/EM output
   integer, parameter:: fileID_boundaryParticle = 30 ! Particle boundaries
   integer, parameter:: fileID_boundaryEM = 31       ! EM boundaries
-  integer, parameter:: fileID_inject = 32          ! Injects
+  integer, parameter:: fileID_inject = 32           ! Injects
   integer, parameter:: fileID_reference = 40        ! Reference values
   integer, parameter:: fileID_time =50              ! Computation time
 

src/modules/solver/moduleSolver.f90

@@ -548,6 +548,7 @@ MODULE moduleSolver
 
         ! Output of boundaries
         call boundariesParticle_write()
+        call boundariesEM_write()
 
         ! Output of inejcts
         call writeInjects()