From c7db17ecd229c6e6ee5bc83ceaa181c1ce59e0d6 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 4 May 2026 19:09:45 +0200 Subject: [PATCH] Right way to accumulate the value. It is not an average, really --- src/modules/init/moduleInput.f90 | 2 ++ src/modules/mesh/moduleMesh.f90 | 1 + src/modules/mesh/moduleMesh@boundaryEM.f90 | 31 ++++++++++++++++------ 3 files changed, 26 insertions(+), 8 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 3bbaa3b..fe278eb 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -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 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 5638ae7..052de3c 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -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 diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90 index 51b6242..85780d2 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -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