From a174a50b81493966a4df2d6417b29de546314cc1 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 2 May 2026 21:17:48 +0200 Subject: [PATCH] Similar but different approach based on Gauss and an average velocity to keep it electrostatic --- src/modules/init/moduleInput.f90 | 2 - src/modules/mesh/moduleMesh.f90 | 1 - src/modules/mesh/moduleMesh@boundaryEM.f90 | 44 +++++++++++++--------- 3 files changed, 26 insertions(+), 21 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 56a698b..3bbaa3b 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -1351,8 +1351,6 @@ MODULE moduleInput 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 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 7d07754..5638ae7 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -1118,7 +1118,6 @@ MODULE moduleMesh 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 integer:: counter contains procedure, pass:: apply => applyFreeCurrent diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90 index d0c203b..6917835 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -373,51 +373,59 @@ submodule(moduleMesh) boundaryEM integer:: e, n, s integer, allocatable:: nodes(:) real(8), allocatable:: mom_nodes(:) + real(8), allocatable:: den_nodes(:) class(meshNode), pointer:: node - real(8):: mom_center, edgeDensityCurrent + real(8):: mom_center, den_center, charge_center + real(8):: vel_avg_center select type(self) type is(boundaryEMFreeCurrent) self%counter = self%counter + 1 do e=1, self%nEdges - edgeDensityCurrent = 0.0d0 + + mom_center = 0.0d0 + den_center = 0.0d0 + charge_center = 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)) + allocate(den_nodes(1:edge%nNodes)) + mom_nodes = 0.0d0 + den_nodes = 0.0d0 + 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 + ! Minus sign as we look at the values exiting the surface (opposite to the normal) mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal) + den_nodes(n) = node%output(s)%den end do - mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) - - edgeDensityCurrent = edgeDensityCurrent + qSpecies(s) * mom_center + mom_center = mom_center + edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) + den_center = den_center + edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes) + charge_center = charge_center + qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes) end do + if (den_center > 1.0e-10) then + vel_avg_center = mom_center / den_center + else + vel_avg_center = 0.0d0 + + end if + end associate - self%surfaceCharge(e) = self%surfaceCharge(e) + edgeDensityCurrent * tauMin + self%electricField(e) = self%electricField(e) - vel_avg_center*charge_center*tauMin - deallocate(nodes, mom_nodes) + deallocate(nodes, mom_nodes, den_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 @@ -438,7 +446,7 @@ submodule(moduleMesh) boundaryEM 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 + print*, self%edges(e)%obj%n, self%electricField(e)*EF_ref end do