Using average flux over the control volume

This commit is contained in:
Jorge Gonzalez 2026-05-05 12:15:21 +02:00
commit 770277a569

View file

@ -403,27 +403,22 @@ submodule(moduleMesh) boundaryEM
node => mesh%nodes(nodes(n))%obj
! Minus sign as we look at the values exiting the surface (opposite to the normal)
! Density is reescaled as 'den' has no data about the node volume
! Density and momentum are reescaled as 'den' has no data about the node volume
mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)/(node%v*Vol_ref*n_ref)
! Avoid the contribution of returning particles close to the edge
! if (mom_nodes(n) < 0.0d0) then
! mom_nodes(n) = 0.0d0
! end if
! den_nodes(n) = node%output(s)%den/(node%v*Vol_ref*n_ref)
den_nodes(n) = node%output(s)%den/(node%v*Vol_ref*n_ref)
end do
mom_center = mom_center + qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
! den_center = den_center + edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
! charge_center = charge_center + qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
mom_center = mom_center + edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
den_center = den_center + edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
charge_center = charge_center + qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
end do
! if (den_center > 1.0e-12) then
self%deltaElectricField(e) = self%deltaElectricField(e) + mom_center * tauMin
if (den_center > 1.0e-12) then
self%deltaElectricField(e) = self%deltaElectricField(e) + mom_center/den_center * charge_center * tauMin
! end if
end if
deallocate(nodes, mom_nodes, den_nodes)