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

Charge screws up the momentum

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

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

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