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

Had to remove charge, it is working now and it only accounts for charged species

Browse source
This commit is contained in:
Jorge Gonzalez 2026-03-16 12:42:49 +01:00
commit 086ce4f589
1 changed files with 27 additions and 19 deletions
Download patch file Download diff file Expand all files Collapse all files

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

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