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

Right way to accumulate the value. It is not an average, really

Browse source
This commit is contained in:
Jorge Gonzalez 2026-05-04 19:09:45 +02:00
commit c7db17ecd2
3 changed files with 26 additions and 8 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/modules/init/moduleInput.f90
View file

@ -1351,6 +1351,8 @@ MODULE moduleInput
bound%nEdges = size(bound%edges)
allocate(bound%electricField(1:bound%nEdges))
bound%electricField = 0.0d0
allocate(bound%deltaElectricField(1:bound%nEdges))
bound%deltaElectricField = 0.0d0
end if

1
src/modules/mesh/moduleMesh.f90
View file

@ -1118,6 +1118,7 @@ MODULE moduleMesh
integer:: nEdges
type(meshEdgePointer), allocatable:: edges(:) ! Edges included in the boundary
real(8), allocatable:: electricField(:) ! Electric field normal to the edge that must be applied with a Neumann boundary condition
real(8), allocatable:: deltaElectricField(:) ! Accumulated change in E
integer:: counter
contains
procedure, pass:: apply => applyFreeCurrent

31
src/modules/mesh/moduleMesh@boundaryEM.f90
View file

@ -378,9 +378,16 @@ submodule(moduleMesh) boundaryEM
class(meshNode), pointer:: node
real(8):: mom_center, den_center, charge_center
real(8):: l_c ! Characteristic length
integer:: every
! TODO: Make an input parameter
every = 10
select type(self)
type is(boundaryEMFreeCurrent)
self%counter = self%counter + 1
do e=1, self%nEdges
mom_center = 0.0d0
den_center = 0.0d0
@ -397,8 +404,13 @@ 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 needed to be reescaled as 'den' has no data about the node volume
! Density is 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)
end do
@ -409,16 +421,11 @@ submodule(moduleMesh) boundaryEM
end do
if (den_center > 1.0e-10) then
l_c = mom_center / den_center * tauMin
else
l_c = 0.0d0
if (den_center > 1.0e-12) then
self%deltaElectricField(e) = self%deltaElectricField(e) + mom_center / den_center * charge_center * tauMin
end if
self%electricField(e) = self%electricField(e) - l_c * charge_center
deallocate(nodes, mom_nodes, den_nodes)
end associate
@ -426,6 +433,14 @@ submodule(moduleMesh) boundaryEM
end do
if (self%counter == every) then
self%electricField = self%electricField - self%deltaElectricField!*every_inv
self%counter = 0
self%deltaElectricField = 0.0d0
end if
end select
end subroutine updateFreeCurrent