From e1dda5b7e092d89ff6f7ede8f759f6aa39d60c0a Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 15 Mar 2026 21:23:41 +0100 Subject: [PATCH 1/2] Charge screws up the momentum --- .../mesh/moduleMesh@boundaryParticle.f90 | 27 +++++++++++++------ 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index d2fa524..334ce7d 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -546,10 +546,11 @@ 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):: mom_center - real(8):: mom_incident, mom_rest + real(8):: den_center, vel_center + real(8):: vel_incident, vel_rest real(8):: alpha select type(self) @@ -557,34 +558,44 @@ submodule(moduleMesh) boundaryParticle do e = 1, size(self%edges) edge => self%edges(e)%obj - mom_incident = 0.d0 - mom_rest = 0.d0 + vel_incident = 0.d0 + vel_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 - mom_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) + den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes) + if (den_center > 1.0d-10) then + vel_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) / den_center + + else + vel_center = 0.0d0 + + end if if (s == self%s_incident) then - mom_incident = mom_center + vel_incident = vel_center else - mom_rest = mom_rest + mom_center + vel_rest = vel_rest + vel_center end if end do - alpha = mom_rest + mom_incident + alpha = vel_rest + vel_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 From 086ce4f589aa05fae4b4866facf0223b5561024f Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 16 Mar 2026 12:42:49 +0100 Subject: [PATCH 2/2] Had to remove charge, it is working now and it only accounts for charged species --- .../mesh/moduleMesh@boundaryParticle.f90 | 46 +++++++++++-------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index 334ce7d..a20ffc2 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -549,7 +549,7 @@ submodule(moduleMesh) boundaryParticle real(8), allocatable:: den_nodes(:) real(8), allocatable:: mom_nodes(:) class(meshNode), pointer:: node - real(8):: den_center, vel_center + real(8):: den_center, mom_center, vel_center real(8):: vel_incident, vel_rest real(8):: alpha @@ -565,32 +565,40 @@ submodule(moduleMesh) boundaryParticle 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_center = 0.0d0 + mom_center = 0.0d0 + vel_center = 0.0d0 + select type(sp => species(s)%obj) + type is(speciesCharged) + 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) + den_nodes(n) = node%output(s)%den + mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal) - end do + end do - den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes) - if (den_center > 1.0d-10) then - vel_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) / den_center + den_center = edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes) + if (den_center > 1.0d-10) then + mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) + vel_center = mom_center / den_center - else - vel_center = 0.0d0 - - end if + else + vel_center = 0.0d0 + + end if - if (s == self%s_incident) then - vel_incident = vel_center - else - vel_rest = vel_rest + vel_center - end if + if (s == self%s_incident) then + vel_incident = vel_center + else + vel_rest = vel_rest + vel_center + end if + + end select end do - alpha = vel_rest + vel_incident + alpha = vel_rest - vel_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