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

Including charge in boundaries. Testing now

Browse source
This commit is contained in:
Jorge Gonzalez 2026-03-15 15:07:26 +01:00
commit f4dab6aa71
1 changed files with 15 additions and 23 deletions
Download patch file Download diff file Expand all files Collapse all files

38
src/modules/mesh/moduleMesh@boundaryParticle.f90
View file

@ -410,6 +410,7 @@ submodule(moduleMesh) boundaryParticle
class(meshEdge), pointer:: edge
real(8), allocatable:: density_nodes(:)
class(meshNode), pointer:: node
real(8):: den_center
real(8):: density_incident, density_rest
real(8):: alpha
@ -431,11 +432,14 @@ submodule(moduleMesh) boundaryParticle
end do
! Gather the density at the edge center and multiply by the species charge
den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
if (s == self%s_incident) then
density_incident = edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
density_incident = den_center
else
density_rest = density_rest + edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
density_rest = density_rest + den_center
end if
@ -444,7 +448,7 @@ submodule(moduleMesh) boundaryParticle
! Correction for this time step
if (density_rest > 1.0d-10) then
! If there is a rest population, correct
alpha = 1.d0 - density_incident/density_rest
alpha = 1.d0 + density_incident/density_rest
else
! If not, alpha is assumed unchaged
@ -542,11 +546,10 @@ submodule(moduleMesh) boundaryParticle
integer:: e, n, s
integer, allocatable:: nodes(:)
class(meshEdge), pointer:: edge
real(8), allocatable:: den_nodes(:)
real(8), allocatable:: mom_nodes(:)
class(meshNode), pointer:: node
real(8):: den_center, mom_center, vel_center
real(8):: velocity_incident, velocity_rest
real(8):: mom_center
real(8):: mom_incident, mom_rest
real(8):: alpha
select type(self)
@ -554,45 +557,34 @@ submodule(moduleMesh) boundaryParticle
do e = 1, size(self%edges)
edge => self%edges(e)%obj
velocity_incident = 0.d0
velocity_rest = 0.d0
mom_incident = 0.d0
mom_rest = 0.d0
nodes = edge%getNodes(edge%nNodes)
allocate(den_nodes(1:edge%nNodes))
allocate(mom_nodes(1:edge%nNodes))
do s = 1, nSpecies
do n = 1, edge%nNodes
node => mesh%nodes(nodes(n))%obj
den_nodes(n) = node%output(s)%den
mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
end do
den_center = edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
if (den_center > 1.0d-10) then
vel_center = mom_center / den_center
else
vel_center = 0.0d0
end if
mom_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
if (s == self%s_incident) then
velocity_incident = vel_center
mom_incident = mom_center
else
velocity_rest = velocity_rest + vel_center
mom_rest = mom_rest + mom_center
end if
end do
alpha = velocity_rest - velocity_incident
alpha = mom_rest + mom_incident
! Apply correction with a factor of 0.1 to avoid fast changes
self%velocity_shift(edge%n) = self%velocity_shift(edge%n) + 1.0d-2 * alpha
deallocate(den_nodes)
deallocate(mom_nodes)
end do