JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Similar but different approach based on Gauss and an average velocity to keep it electrostatic

Browse source
This commit is contained in:
Jorge Gonzalez 2026-05-02 21:17:48 +02:00
commit a174a50b81
3 changed files with 26 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

44
src/modules/mesh/moduleMesh@boundaryEM.f90
View file

@ -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