From 776734db5805041241e3fb36252dea60b32dc3d2 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 9 Mar 2026 14:31:11 +0100 Subject: [PATCH 1/4] Calculates a reflection parameter (alpha) per edge based on the ratio of densities between the incident species and the rest --- src/modules/init/moduleInput.f90 | 8 ++- src/modules/mesh/moduleMesh.f90 | 5 +- .../mesh/moduleMesh@boundaryParticle.f90 | 55 +++++++++++++++++-- 3 files changed, 60 insertions(+), 8 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 9494517..492eaf5 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -826,6 +826,7 @@ MODULE moduleInput real(8):: eThreshold !Energy threshold integer:: speciesID, electronSecondaryID character(:), allocatable:: speciesName, crossSection, electronSecondary + integer:: s_incident ! Read models of particles object = 'boundaries.particles' @@ -900,7 +901,12 @@ MODULE moduleInput END IF case('quasiNeutrality') - call initQuasiNeutrality(self) + call config%get(object // '.incident', speciesName, found) + if (.not. found) call criticalError("Incident species name not found for quasiNeutrality boundary model", 'readBoundary') + + s_incident = speciesName2Index(speciesName) + + call initQuasiNeutrality(self, s_incident) CASE DEFAULT CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 9133cca..25d7175 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -675,8 +675,9 @@ MODULE moduleMesh end subroutine initIonization - module subroutine initQuasiNeutrality(boundary) + module subroutine initQuasiNeutrality(boundary, s_incident) class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + integer, intent(in):: s_incident end subroutine initQuasiNeutrality @@ -762,7 +763,7 @@ MODULE moduleMesh ! Ensures quasi-neutrality by changing the reflection coefficient type, public, extends(boundaryParticleGeneric):: boundaryQuasiNeutrality - real(8):: alpha ! Reflection parameter + real(8), allocatable:: alpha(:) ! Reflection parameter integer:: s_incident ! species index of the incident species type(meshEdgePointer), allocatable:: edges(:) !Array with edges contains diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index 040f633..ba01430 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -95,18 +95,23 @@ submodule(moduleMesh) boundaryParticle END SUBROUTINE initIonization - module subroutine initQuasiNeutrality(boundary) + module subroutine initQuasiNeutrality(boundary, s_incident) implicit none class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + integer, intent(in):: s_incident allocate(boundaryQuasiNeutrality:: boundary) select type(boundary) type is(boundaryQuasiNeutrality) - boundary%alpha = 0.d0 allocate(boundary%edges(0)) + boundary%s_incident = s_incident + + allocate(boundary%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here + boundary%alpha = 0.d0 + boundary%update => quasiNeutrality_update end select @@ -333,7 +338,7 @@ submodule(moduleMesh) boundaryParticle class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - if (random() <= self%alpha) then + if (random() <= self%alpha(edge%n)) then call genericReflection(edge, part) else @@ -347,12 +352,52 @@ submodule(moduleMesh) boundaryParticle implicit none class(boundaryParticleGeneric), intent(inout):: self + integer:: e, n, s + integer, allocatable:: nodes(:) + class(meshEdge), pointer:: edge + real(8), allocatable:: density_nodes(:) + class(meshNode), pointer:: node + real(8):: density_incident, density_rest + real(8):: alpha select type(self) type is(boundaryQuasiNeutrality) - self%alpha = 0.1d0 + do e = 1, size(self%edges) + edge => mesh%edges(e)%obj - print *, self%alpha + density_incident = 0.d0 + density_rest = 0.d0 + + nodes = edge%getNodes(edge%nNodes) + allocate(density_nodes(1:edge%nNodes)) + do s = 1, nSpecies + do n = 1, edge%nNodes + node => mesh%nodes(n)%obj + density_nodes(n) = node%output(s)%den + + end do + + if (s == self%s_incident) then + density_incident = edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes) + + else + density_rest = density_rest + edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes) + + end if + + end do + + alpha = 1.d0 - density_incident/density_rest + + ! Limit alpha between 0 and 1 + alpha = min(alpha, 1.d0) + alpha = max(alpha, 0.d0) + + self%alpha(edge%n) = alpha + + deallocate(density_nodes) + + end do end select From e7a47611d9c748223df26ae506378590281d537b Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 9 Mar 2026 18:07:06 +0100 Subject: [PATCH 2/4] This is now an iterative process, which makes things better --- src/modules/mesh/moduleMesh@boundaryParticle.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index ba01430..c51c26d 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -387,13 +387,15 @@ submodule(moduleMesh) boundaryParticle end do + ! Correction for this time step alpha = 1.d0 - density_incident/density_rest - ! Limit alpha between 0 and 1 - alpha = min(alpha, 1.d0) - alpha = max(alpha, 0.d0) + ! Apply correction with a factor of 0.1 to avoid fast changes + self%alpha(edge%n) = self%alpha(edge%n) + 0.1d0 * alpha - self%alpha(edge%n) = alpha + ! Limit alpha between 0 and 1 + self%alpha(edge%n) = min(self%alpha(edge%n), 1.d0) + self%alpha(edge%n) = max(self%alpha(edge%n), 0.d0) deallocate(density_nodes) From 4489da74d5017cea72c2329bd19d43b03a038193 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 9 Mar 2026 18:34:30 +0100 Subject: [PATCH 3/4] Fixing an issue with the average in text format --- src/modules/mesh/inout/text/moduleMeshOutputText.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 index e243c80..9ee551e 100644 --- a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 +++ b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 @@ -234,7 +234,7 @@ module moduleMeshOutputText fileIDDeviation = 67 do s = 1, nSpecies - fileNameMean = formatFileName('Average_mean', species(s)%obj%name, 'csv', timeStep) + fileNameMean = formatFileName('Average_mean', species(s)%obj%name, 'csv') write(*, "(6X,A15,A)") "Creating file: ", fileNameMean open (fileIDMean, file = path // folder // '/' // fileNameMean) @@ -243,7 +243,7 @@ module moduleMeshOutputText 'Velocity, mean (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', & 'Temperature, mean (K)' - fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv', timeStep) + fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv') write(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation open (fileIDDeviation, file = path // folder // '/' // fileNameDeviation) From d7d0d6b47f31542f090b40b0e8d0266e032b62d7 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 10 Mar 2026 12:32:46 +0100 Subject: [PATCH 4/4] Add quotes to headers --- .../mesh/inout/text/moduleMeshOutputText.f90 | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 index 9ee551e..0cd539b 100644 --- a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 +++ b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 @@ -137,10 +137,10 @@ module moduleMeshOutputText write(*, "(6X,A15,A)") "Creating file: ", fileName open (fileID, file = path // folder // '/' // fileName) - write(fileID, '(5(A,","),A)') 'Position (m)', & - 'Density (m^-3)', & - 'Velocity (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', & - 'Temperature (K)' + write(fileID, '(5(A,","),A)') '"Position (m)"', & + '"Density (m^-3)"', & + '"Velocity (m s^-1):0"', 'Velocity (m s^-1):1"', 'Velocity (m s^-1):2"', & + '"Temperature (K)"' call writeSpeciesOutput(self, fileID, s) @@ -206,10 +206,10 @@ module moduleMeshOutputText write(*, "(6X,A15,A)") "Creating file: ", fileName open (fileID, file = path // folder // '/' // fileName) - write(fileID, '(8(A,","),A)') 'Position (m)', & - 'Potential (V)', & - 'Electric Field (V m^-1):0', 'Electric Field (V m^-1):1', 'Electric Field (V m^-1):2', & - 'Magnetic Field (T):0', 'Magnetic Field (T):1', 'Magnetic Field (T):2' + write(fileID, '(8(A,","),A)') '"Position (m)"', & + '"Potential (V)"', & + '"Electric Field (V m^-1):0"', 'Electric Field (V m^-1):1"', 'Electric Field (V m^-1):2"', & + '"Magnetic Field (T):0"', 'Magnetic Field (T):1"', 'Magnetic Field (T):2"' call writeEMOutput(self, fileID) @@ -238,19 +238,19 @@ module moduleMeshOutputText write(*, "(6X,A15,A)") "Creating file: ", fileNameMean open (fileIDMean, file = path // folder // '/' // fileNameMean) - write(fileIDMean, '(5(A,","),A)') 'Position (m)', & - 'Density, mean (m^-3)', & - 'Velocity, mean (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', & - 'Temperature, mean (K)' + write(fileIDMean, '(5(A,","),A)') '"Position (m)"', & + '"Density, mean (m^-3)"', & + '"Velocity, mean (m s^-1):0"', '"Velocity (m s^-1):1"', '"Velocity (m s^-1):2"', & + '"Temperature, mean (K)"' fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv') write(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation open (fileIDDeviation, file = path // folder // '/' // fileNameDeviation) - write(fileIDDeviation, '(5(A,","),A)') 'Position (m)', & - 'Density, deviation (m^-3)', & - 'Velocity, deviation (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', & - 'Temperature, deviation (K)' + write(fileIDDeviation, '(5(A,","),A)') '"Position (m)"', & + '"Density, deviation (m^-3)"', & + '"Velocity, deviation (m s^-1):0"', 'Velocity (m s^-1):1"', 'Velocity (m s^-1):2"', & + '"Temperature, deviation (K)"' call writeAverage(self, fileIDMean, fileIDDeviation, s)