From b0168338761de48ad2d96ff4de1fa904e347a61b Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 4 May 2026 12:10:50 +0200 Subject: [PATCH] I think I have it, and I think that the averaring might not be required --- src/modules/init/moduleInput.f90 | 2 ++ src/modules/mesh/moduleMesh.f90 | 1 + src/modules/mesh/moduleMesh@boundaryEM.f90 | 35 ++++++++++++++-------- 3 files changed, 26 insertions(+), 12 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 6917835..25d1b88 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -367,6 +367,7 @@ submodule(moduleMesh) boundaryEM ! Update subroutine updateFreeCurrent(self) use moduleCaseParam, only: tauMin + use moduleRefParam, only: Vol_ref, n_ref implicit none class(boundaryEMGeneric), intent(inout):: self @@ -376,13 +377,19 @@ submodule(moduleMesh) boundaryEM real(8), allocatable:: den_nodes(:) class(meshNode), pointer:: node real(8):: mom_center, den_center, charge_center - real(8):: vel_avg_center + real(8):: l_c ! Characteristic length + integer:: every + real(8):: every_inv + + ! TODO: Make this an input parameter + every = 10 + every_inv = 1.0d0 / real(every) select type(self) type is(boundaryEMFreeCurrent) self%counter = self%counter + 1 - do e=1, self%nEdges + do e=1, self%nEdges mom_center = 0.0d0 den_center = 0.0d0 charge_center = 0.0d0 @@ -392,16 +399,14 @@ submodule(moduleMesh) boundaryEM nodes = edge%getNodes(edge%nNodes) allocate(mom_nodes(1:edge%nNodes)) allocate(den_nodes(1:edge%nNodes)) - mom_nodes = 0.0d0 - den_nodes = 0.0d0 do s = 1, nSpecies do n = 1, self%nNodes node => mesh%nodes(nodes(n))%obj ! Minus sign as we look at the values exiting the surface (opposite to the normal) - 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)/(node%v*Vol_ref*n_ref) + den_nodes(n) = node%output(s)%den/(node%v*Vol_ref*n_ref) end do @@ -412,20 +417,26 @@ submodule(moduleMesh) boundaryEM end do if (den_center > 1.0e-10) then - vel_avg_center = mom_center / den_center - else - vel_avg_center = 0.0d0 + self%deltaElectricField(e) = self%deltaElectricField(e) + mom_center * charge_center / den_center * tauMin end if + deallocate(nodes, mom_nodes, den_nodes) + end associate - self%electricField(e) = self%electricField(e) - vel_avg_center*charge_center*tauMin - - deallocate(nodes, mom_nodes, den_nodes) end do + if (self%counter == every) then + + self%electricField = self%electricField - self%deltaElectricField*every_inv + + ! reset for the next iteration + self%counter = 0 + self%deltaElectricField = 0.0d0 + + end if end select