From f4dab6aa71752c65e1c16be24ca6448f864f9ba9 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 15 Mar 2026 15:07:26 +0100 Subject: [PATCH] Including charge in boundaries. Testing now --- .../mesh/moduleMesh@boundaryParticle.f90 | 38 ++++++++----------- 1 file changed, 15 insertions(+), 23 deletions(-) diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index 0280dcf..d2fa524 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -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