From 8edfd803e7fbf81826b48bee3f15345431c7034e Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Thu, 7 May 2026 11:14:16 +0200 Subject: [PATCH] Finally staying in Ampere. The other one is not quite stable. Have to figure out a way to justify this in the paper (not too hard). I have rescaled all the nodes values by the node volume to avoid issues --- src/modules/mesh/moduleMesh@boundaryEM.f90 | 28 ++++++------------- .../mesh/moduleMesh@boundaryParticle.f90 | 8 ++++-- 2 files changed, 13 insertions(+), 23 deletions(-) diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90 index 0be427a..ded0181 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -239,6 +239,7 @@ submodule(moduleMesh) boundaryEM ! Update subroutine updateFloating(self) use moduleCaseParam, only: tauMin + use moduleRefParam, only: Vol_ref, n_ref implicit none class(boundaryEMGeneric), intent(inout):: self @@ -265,7 +266,7 @@ submodule(moduleMesh) boundaryEM node => mesh%nodes(nodes(n))%obj ! Minus sign to get the flux exiting the edge - mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal) + mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)/(node%v*Vol_ref*n_ref) end do @@ -387,51 +388,38 @@ submodule(moduleMesh) boundaryEM integer:: e, n, s integer, allocatable:: nodes(:) real(8), allocatable:: mom_nodes(:) - real(8), allocatable:: den_nodes(:) class(meshNode), pointer:: node - real(8):: mom_center, den_center, charge_center + real(8):: currentDensity_center select type(self) type is(boundaryEMFreeCurrent) self%counter = self%counter + 1 do e=1, self%nEdges - mom_center = 0.0d0 - den_center = 0.0d0 - charge_center = 0.0d0 + currentDensity_center = 0.0d0 associate(edge => self%edges(e)%obj) nodes = edge%getNodes(edge%nNodes) allocate(mom_nodes(1:edge%nNodes)) - allocate(den_nodes(1:edge%nNodes)) 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) - ! Density and momentum are reescaled as 'den' has no data about the node volume + ! Momentum is reescaled as value at the node has no data about the node volume 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 - ! 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) - - mom_center = mom_center + qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) + currentDensity_center = currentDensity_center + qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) end do - ! if (den_center > 1.0e-12) then - ! self%deltaElectricField(e) = self%deltaElectricField(e) + mom_center/den_center * charge_center * tauMin + self%deltaElectricField(e) = self%deltaElectricField(e) + currentDensity_center * tauMin - ! end if - self%deltaElectricField(e) = self%deltaElectricField(e) + mom_center * tauMin - - deallocate(nodes, mom_nodes, den_nodes) + deallocate(nodes, mom_nodes) end associate diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index 04116be..514e613 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -403,6 +403,7 @@ submodule(moduleMesh) boundaryParticle ! Update subroutine quasiNeutrality_update(self) + use moduleRefParam, only: Vol_ref, n_ref implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -429,7 +430,7 @@ submodule(moduleMesh) boundaryParticle do n = 1, edge%nNodes node => mesh%nodes(nodes(n))%obj - density_nodes(n) = node%output(s)%den + density_nodes(n) = node%output(s)%den/(node%v*Vol_ref*n_ref) end do @@ -541,6 +542,7 @@ submodule(moduleMesh) boundaryParticle ! Update subroutine outflowAdaptive_update(self) + use moduleRefParam, only: Vol_ref, n_ref implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -574,8 +576,8 @@ submodule(moduleMesh) boundaryParticle 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/(node%v*Vol_ref*n_ref) + mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)/(node%v*Vol_ref*n_ref) end do