From a3bdf8230a50fd21f0dd23e55c831df5d8f7cc06 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 19 May 2024 10:55:20 +0200 Subject: [PATCH 01/17] Implementation of Boltzmann electrons Still not working, just saving code. --- src/modules/init/moduleInput.f90 | 4 + .../solver/electromagnetic/moduleEM.f90 | 95 +++++++++++++++++-- src/modules/solver/moduleSolver.f90 | 2 + 3 files changed, 94 insertions(+), 7 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 46b627b..7f6c036 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -273,6 +273,10 @@ MODULE moduleInput !Read BC CALL readEMBoundary(config) + CASE("ElectrostaticBoltzmann") + !Read BC + CALL readEMBoundary(config) + CASE("ConstantB") !Read BC CALL readEMBoundary(config) diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index bdf6b03..c7c74ee 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -49,7 +49,7 @@ MODULE moduleEM END SUBROUTINE !Assemble the source vector based on the charge density to solve Poisson's equation - SUBROUTINE assembleSourceVector(vectorF) + SUBROUTINE assembleSourceVector(vectorF, n_e) USE moduleMesh USE moduleRefParam IMPLICIT NONE @@ -58,15 +58,16 @@ MODULE moduleEM REAL(8), ALLOCATABLE:: localF(:) INTEGER, ALLOCATABLE:: nodes(:) REAL(8), ALLOCATABLE:: rho(:) + REAL(8), INTENT(in), OPTIONAL:: n_e(1:mesh%numNodes) INTEGER:: nNodes INTEGER:: e, i, ni CLASS(meshNode), POINTER:: node - !$OMP SINGLE + ! !$OMP SINGLE vectorF = 0.D0 - !$OMP END SINGLE + ! !$OMP END SINGLE - !$OMP DO REDUCTION(+:vectorF) + ! !$OMP DO REDUCTION(+:vectorF) DO e = 1, mesh%numCells nNodes = mesh%cells(e)%obj%nNodes nodes = mesh%cells(e)%obj%getNodes(nNodes) @@ -77,6 +78,10 @@ MODULE moduleEM ni = nodes(i) node => mesh%nodes(ni)%obj rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref)) + IF (PRESENT(n_e)) THEN + rho(i) = rho(i) - n_e(i) + + END IF END DO @@ -94,10 +99,10 @@ MODULE moduleEM DEALLOCATE(nodes, rho) END DO - !$OMP END DO + ! !$OMP END DO !Apply boundary conditions - !$OMP DO + ! !$OMP DO DO i = 1, mesh%numNodes node => mesh%nodes(i)%obj @@ -108,7 +113,7 @@ MODULE moduleEM END SELECT END DO - !$OMP END DO + ! !$OMP END DO END SUBROUTINE assembleSourceVector @@ -156,4 +161,80 @@ MODULE moduleEM END SUBROUTINE solveElecField + FUNCTION BoltzmannElectron(phi, n) RESULT(n_e) + USE moduleRefParam + USE moduleConstParam + IMPLICIT NONE + + INTEGER, INTENT(in):: n + REAL(8), INTENT(in):: phi(1:n) + REAL(8):: n_e(1:n) + REAL(8):: n_e0 = 1.0D16, phi_0 = -520.0D0, T_e = 11604.0 + + n_e = n_e0 / n_ref * EXP(-qe * (phi*Volt_ref - phi_0) / (kb * T_e)) + + RETURN + + END FUNCTION BoltzmannElectron + + SUBROUTINE solveElecFieldBoltzmann + USE moduleMesh + USE moduleErrors + IMPLICIT NONE + + INTEGER, SAVE:: INFO + INTEGER:: n + REAL(8), ALLOCATABLE, SAVE:: tempF(:) + REAL(8), ALLOCATABLE:: n_e(:), phi_old(:) + INTEGER:: k + EXTERNAL:: dgetrs + + !$OMP SINGLE + ALLOCATE(tempF(1:mesh%numNodes)) + ALLOCATE(n_e(1:mesh%numNodes)) + ALLOCATE(phi_old(1:mesh%numNodes)) + DO n = 1, mesh%numNodes + tempF(n) = mesh%nodes(n)%obj%emData%phi + + END DO + n_e = BoltzmannElectron(tempF, mesh%numNodes) + !$OMP END SINGLE + + CALL assembleSourceVector(tempF, n_e) + + !$OMP SINGLE + DO k = 1, 5 + phi_old = tempF + CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, & + mesh%IPIV, tempF, mesh%numNodes, info) + + PRINT*, k, "diff = ", MAXVAL(ABS(tempF - phi_old)) + n_e = BoltzmannElectron(tempF, mesh%numNodes) + CALL assembleSourceVector(tempF, n_e) + + END DO + !$OMP END SINGLE + + IF (info == 0) THEN + !Suscessful resolution of Poission equation + !$OMP DO + DO n = 1, mesh%numNodes + mesh%nodes(n)%obj%emData%phi = tempF(n) + + END DO + !$OMP END DO + + ELSE + !$OMP SINGLE + CALL criticalError('Poisson equation failed', 'solveElecFieldBoltzmann') + !$OMP END SINGLE + + END IF + + !$OMP SINGLE + DEALLOCATE(tempF) + !$OMP END SINGLE + + END SUBROUTINE solveElecFieldBoltzmann + END MODULE moduleEM diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index e539fed..c7a0785 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -138,6 +138,8 @@ MODULE moduleSolver CASE('Electrostatic','ConstantB') self%solveEM => solveElecField + CASE('ElectrostaticBoltzmann') + self%solveEM => solveElecFieldBoltzmann END SELECT END SUBROUTINE initEM From e4f7987f9078a44dec5484df643402c7f0e20280 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 19 May 2024 16:45:03 +0200 Subject: [PATCH 02/17] Trying to solve Still I don't understand this basic thing... --- .../solver/electromagnetic/moduleEM.f90 | 30 ++++++++++++------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index c7c74ee..0f02bed 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -170,8 +170,18 @@ MODULE moduleEM REAL(8), INTENT(in):: phi(1:n) REAL(8):: n_e(1:n) REAL(8):: n_e0 = 1.0D16, phi_0 = -520.0D0, T_e = 11604.0 + INTEGER:: i - n_e = n_e0 / n_ref * EXP(-qe * (phi*Volt_ref - phi_0) / (kb * T_e)) + DO i =1, n + IF (phi(i)*Volt_ref >= phi_0) THEN + n_e(i) = n_e0 / n_ref * EXP(-qe * (phi(i)*Volt_ref - phi_0) / (kb * T_e)) + + ELSE + n_e(i) = 0.D0 + + END IF + + END DO RETURN @@ -185,7 +195,7 @@ MODULE moduleEM INTEGER, SAVE:: INFO INTEGER:: n REAL(8), ALLOCATABLE, SAVE:: tempF(:) - REAL(8), ALLOCATABLE:: n_e(:), phi_old(:) + REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:) INTEGER:: k EXTERNAL:: dgetrs @@ -193,21 +203,19 @@ MODULE moduleEM ALLOCATE(tempF(1:mesh%numNodes)) ALLOCATE(n_e(1:mesh%numNodes)) ALLOCATE(phi_old(1:mesh%numNodes)) - DO n = 1, mesh%numNodes - tempF(n) = mesh%nodes(n)%obj%emData%phi - - END DO - n_e = BoltzmannElectron(tempF, mesh%numNodes) - !$OMP END SINGLE - + n_e = 0.D0 CALL assembleSourceVector(tempF, n_e) + !$OMP END SINGLE + !$OMP SINGLE - DO k = 1, 5 + DO k = 1, 10 phi_old = tempF CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, & mesh%IPIV, tempF, mesh%numNodes, info) + PRINT *, MAXVAL(n_e), MINVAL(n_e) + PRINT *, MAXVAL(tempF), MINVAL(tempF) PRINT*, k, "diff = ", MAXVAL(ABS(tempF - phi_old)) n_e = BoltzmannElectron(tempF, mesh%numNodes) CALL assembleSourceVector(tempF, n_e) @@ -232,7 +240,7 @@ MODULE moduleEM END IF !$OMP SINGLE - DEALLOCATE(tempF) + DEALLOCATE(tempF, n_e, phi_old) !$OMP END SINGLE END SUBROUTINE solveElecFieldBoltzmann From 98ee3e9c9c38ec6d461c50268a44a99912c8eb70 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Mon, 30 Sep 2024 17:06:25 +0200 Subject: [PATCH 03/17] Still not working, but will be fixed I have the solution in the plasma expansion code, but I need to do other stuff. --- src/modules/moduleInject.f90 | 3 +- .../solver/electromagnetic/moduleEM.f90 | 50 +++++++++---------- 2 files changed, 26 insertions(+), 27 deletions(-) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 4e57083..496ea6a 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -108,7 +108,8 @@ MODULE moduleInject CASE ("A") !Input current in Ampers - self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight)) + ! TODO: Make this only available for charge species + self%nParticles = INT(flow*tauInject*ti_ref/(qe*abs(species(sp)%obj%qm*species(sp)%obj%m)*species(sp)%obj%weight)) CASE ("part/s") !Input current in Ampers diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index 0f02bed..a8a5323 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -46,7 +46,7 @@ MODULE moduleEM END DO - END SUBROUTINE + END SUBROUTINE apply !Assemble the source vector based on the charge density to solve Poisson's equation SUBROUTINE assembleSourceVector(vectorF, n_e) @@ -169,25 +169,16 @@ MODULE moduleEM INTEGER, INTENT(in):: n REAL(8), INTENT(in):: phi(1:n) REAL(8):: n_e(1:n) - REAL(8):: n_e0 = 1.0D16, phi_0 = -520.0D0, T_e = 11604.0 + REAL(8):: n_e0 = 1.0D16, phi_0 = -500.0D0, T_e = 11604.0 INTEGER:: i - DO i =1, n - IF (phi(i)*Volt_ref >= phi_0) THEN - n_e(i) = n_e0 / n_ref * EXP(-qe * (phi(i)*Volt_ref - phi_0) / (kb * T_e)) - - ELSE - n_e(i) = 0.D0 - - END IF - - END DO + n_e = n_e0 / n_ref * exp(qe * (phi*Volt_ref - phi_0) / (kb * T_e)) RETURN END FUNCTION BoltzmannElectron - SUBROUTINE solveElecFieldBoltzmann + SUBROUTINE solveElecFieldBoltzmann() USE moduleMesh USE moduleErrors IMPLICIT NONE @@ -195,7 +186,7 @@ MODULE moduleEM INTEGER, SAVE:: INFO INTEGER:: n REAL(8), ALLOCATABLE, SAVE:: tempF(:) - REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:) + REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:), phi(:) INTEGER:: k EXTERNAL:: dgetrs @@ -203,22 +194,29 @@ MODULE moduleEM ALLOCATE(tempF(1:mesh%numNodes)) ALLOCATE(n_e(1:mesh%numNodes)) ALLOCATE(phi_old(1:mesh%numNodes)) - n_e = 0.D0 - CALL assembleSourceVector(tempF, n_e) - + ALLOCATE(phi(1:mesh%numNodes)) !$OMP END SINGLE + !$OMP DO + DO n = 1, mesh%numNodes + phi_old(n) = mesh%nodes(n)%obj%emData%phi + + END DO + !$OMP END DO + !$OMP SINGLE - DO k = 1, 10 - phi_old = tempF + DO k = 1, 100 + n_e = BoltzmannElectron(phi_old, mesh%numNodes) + CALL assembleSourceVector(tempF, n_e) + CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, & mesh%IPIV, tempF, mesh%numNodes, info) + phi = tempF - PRINT *, MAXVAL(n_e), MINVAL(n_e) - PRINT *, MAXVAL(tempF), MINVAL(tempF) - PRINT*, k, "diff = ", MAXVAL(ABS(tempF - phi_old)) - n_e = BoltzmannElectron(tempF, mesh%numNodes) - CALL assembleSourceVector(tempF, n_e) + PRINT *, MAXVAL(n_e), MINVAL(n_e) + PRINT *, MAXVAL(phi), MINVAL(phi) + PRINT*, k, "diff = ", MAXVAL(ABS(phi - phi_old)) + phi_old = phi END DO !$OMP END SINGLE @@ -227,7 +225,7 @@ MODULE moduleEM !Suscessful resolution of Poission equation !$OMP DO DO n = 1, mesh%numNodes - mesh%nodes(n)%obj%emData%phi = tempF(n) + mesh%nodes(n)%obj%emData%phi = phi_old(n) END DO !$OMP END DO @@ -240,7 +238,7 @@ MODULE moduleEM END IF !$OMP SINGLE - DEALLOCATE(tempF, n_e, phi_old) + DEALLOCATE(tempF, n_e, phi_old, phi) !$OMP END SINGLE END SUBROUTINE solveElecFieldBoltzmann From d28dd16c2e14c42c0ce712b05f33ecbf6a1afd01 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Thu, 17 Jul 2025 18:34:11 +0200 Subject: [PATCH 04/17] Average fix and data for Xe --- data/collisions/IO_e-Xe.dat | 52 ++++++++++++++++++++++++++++++++ src/modules/init/moduleInput.f90 | 6 +++- 2 files changed, 57 insertions(+), 1 deletion(-) create mode 100644 data/collisions/IO_e-Xe.dat diff --git a/data/collisions/IO_e-Xe.dat b/data/collisions/IO_e-Xe.dat new file mode 100644 index 0000000..3067578 --- /dev/null +++ b/data/collisions/IO_e-Xe.dat @@ -0,0 +1,52 @@ +# EL cross sections extracted from PROGRAM MAGBOLTZ, VERSION 7.1 JUNE 2004 www.lxcat.net/Biagi-v7.1 +# Relative energy (eV) cross section (m^2) +1.21E+01 0 +1.41E+01 3.923E-21 +1.64E+01 1.194E-20 +1.91E+01 2.1E-20 +2.22E+01 2.946E-20 +2.58E+01 3.65E-20 +3.00E+01 4.185E-20 +3.49E+01 4.552E-20 +4.06E+01 4.766E-20 +4.72E+01 4.85E-20 +5.49E+01 4.828E-20 +6.39E+01 5.031E-20 +7.43E+01 5.1E-20 +8.64E+01 5.1E-20 +1.01E+02 5.032E-20 +1.17E+02 4.906E-20 +1.36E+02 4.732E-20 +1.58E+02 4.521E-20 +1.84E+02 4.283E-20 +2.14E+02 4.029E-20 +2.49E+02 3.764E-20 +2.90E+02 3.497E-20 +3.37E+02 3.233E-20 +3.92E+02 2.975E-20 +4.56E+02 2.726E-20 +5.31E+02 2.489E-20 +6.17E+02 2.266E-20 +7.18E+02 2.056E-20 +8.35E+02 1.861E-20 +9.72E+02 1.68E-20 +1.13E+03 1.514E-20 +1.32E+03 1.361E-20 +1.53E+03 1.221E-20 +1.78E+03 1.094E-20 +2.07E+03 9.781E-21 +2.41E+03 8.735E-21 +2.80E+03 7.789E-21 +3.26E+03 6.938E-21 +3.79E+03 6.171E-21 +4.41E+03 5.484E-21 +5.13E+03 4.868E-21 +5.97E+03 4.316E-21 +6.94E+03 3.824E-21 +8.07E+03 3.385E-21 +9.39E+03 2.994E-21 +1.09E+04 2.646E-21 +1.27E+04 2.336E-21 +1.48E+04 2.062E-21 +1.72E+04 1.818E-21 +2.00E+04 1.602E-21 diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 5090942..18c6e65 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -1331,6 +1331,7 @@ MODULE moduleInput USE moduleCaseParam, ONLY: tauMin USE moduleMesh, ONLY: mesh USE moduleSpecies, ONLY: nSpecies + USE moduleRefParam, ONLY: ti_ref IMPLICIT NONE TYPE(json_file), INTENT(inout):: config @@ -1344,8 +1345,11 @@ MODULE moduleInput CALL config%get('average.startTime', tStart, found) IF (found) THEN - tAverageStart = INT(tStart / tauMin) + tAverageStart = INT(tStart / ti_ref / tauMin) + ELSE + tAverageStart = 0 + END IF ALLOCATE(averageScheme(1:mesh%numNodes)) From a2f9957f32219725485e77a03a7475854d96a213 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 18 Jul 2025 16:31:52 +0200 Subject: [PATCH 05/17] I am dumb The Poisson equation was not working because I didn't finish implementing the new type of BCs. Dirichlet is probably untested. I should stop doing shitty developments and no testing. --- src/modules/output/moduleOutput.f90 | 1 - .../solver/electromagnetic/moduleEM.f90 | 33 ++++++++----------- 2 files changed, 13 insertions(+), 21 deletions(-) diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index e6dc91f..2af3c70 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -22,7 +22,6 @@ MODULE moduleOutput !Type for EM data in node TYPE emNode - CHARACTER(:), ALLOCATABLE:: type REAL(8):: phi REAL(8):: B(1:3) diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index 740fd14..7f6891c 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -11,7 +11,6 @@ MODULE moduleEM CONTAINS PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply - !PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better. END TYPE boundaryEMGeneric @@ -168,8 +167,8 @@ MODULE moduleEM INTEGER:: n, ni DO n = 1, self%nNodes - self%nodes(n)%obj%emData%phi = self%potential - vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi + self%nodes(n)%obj%emData%phi = self%potential + vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi END DO @@ -189,8 +188,8 @@ MODULE moduleEM timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin) DO n = 1, self%nNodes - self%nodes(n)%obj%emData%phi = self%potential * timeFactor - vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi + self%nodes(n)%obj%emData%phi = self%potential * timeFactor + vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi END DO @@ -211,11 +210,11 @@ MODULE moduleEM INTEGER:: e, i, ni, b CLASS(meshNode), POINTER:: node - ! !$OMP SINGLE + !$OMP SINGLE vectorF = 0.D0 - ! !$OMP END SINGLE + !$OMP END SINGLE - ! !$OMP DO REDUCTION(+:vectorF) + !$OMP DO REDUCTION(+:vectorF) DO e = 1, mesh%numCells nNodes = mesh%cells(e)%obj%nNodes nodes = mesh%cells(e)%obj%getNodes(nNodes) @@ -247,21 +246,15 @@ MODULE moduleEM DEALLOCATE(nodes, rho) END DO - ! !$OMP END DO + !$OMP END DO !Apply boundary conditions - ! !$OMP DO - DO i = 1, mesh%numNodes - node => mesh%nodes(i)%obj + !$OMP SINGLE + do b = 1, nBoundaryEM + call boundaryEM(b)%obj%apply(vectorF) - SELECT CASE(node%emData%type) - CASE ("dirichlet") - vectorF(i) = node%emData%phi - - END SELECT - - END DO - ! !$OMP END DO + end do + !$OMP END SINGLE END SUBROUTINE assembleSourceVector From 69215ef66d52557f84232c8ffedd906d0cb2bc90 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez - Galactica Date: Tue, 22 Jul 2025 19:52:39 +0200 Subject: [PATCH 06/17] Change in calculating ionization I don't know why I normalizing density n_0 by Vol_ref and not n_ref --- src/modules/moduleBoundary.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 0b76105..b5bda2f 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -126,7 +126,7 @@ MODULE moduleBoundary SELECT TYPE(boundary) TYPE IS(boundaryIonization) boundary%m0 = m0 / m_ref - boundary%n0 = n0 * Vol_ref + boundary%n0 = n0 / n_ref boundary%v0 = v0 / v_ref boundary%vTh = DSQRT(kb*T0/m0)/v_ref boundary%species => species(ion)%obj From 7f73b69dc2f42933344d75e887a63a76697e9c89 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 27 Jul 2025 17:14:38 +0200 Subject: [PATCH 07/17] Fix injection Half-Maxwellian distribution should inject particles correctly --- src/modules/common/moduleRandom.f90 | 17 +++++++++++++++++ src/modules/moduleInject.f90 | 5 +---- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index ae5c548..7d9cec6 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -66,6 +66,23 @@ MODULE moduleRandom END FUNCTION randomMaxwellian + !Returns a random number in a Maxwellian distribution of mean 0 and width 1 + FUNCTION randomHalfMaxwellian() RESULT(rnd) + IMPLICIT NONE + + REAL(8):: rnd + REAL(8):: x + + rnd = 0.D0 + x = 0.D0 + DO WHILE (x == 0.D0) + CALL RANDOM_NUMBER(x) + END DO + + rnd = DSQRT(-DLOG(x)) + + END FUNCTION randomHalfMaxwellian + !Returns a random number weighted with the cumWeight array FUNCTION randomWeighted(cumWeight, sumWeight) RESULT(rnd) IMPLICIT NONE diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index b4e0ed2..b6fc64b 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -298,10 +298,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - DO WHILE (v <= 0.D0) - v = self%vTh*randomMaxwellian() - - END DO + v = self%vTh*randomHalfMaxwellian() END FUNCTION randomVelHalfMaxwellian From 8e531ede08bea0d68863279b17d459fa155b1aaa Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 27 Jul 2025 17:16:57 +0200 Subject: [PATCH 08/17] Vol_ref was the right answer --- src/modules/moduleBoundary.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index b5bda2f..0b76105 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -126,7 +126,7 @@ MODULE moduleBoundary SELECT TYPE(boundary) TYPE IS(boundaryIonization) boundary%m0 = m0 / m_ref - boundary%n0 = n0 / n_ref + boundary%n0 = n0 * Vol_ref boundary%v0 = v0 / v_ref boundary%vTh = DSQRT(kb*T0/m0)/v_ref boundary%species => species(ion)%obj From d1e73297eb9f51a6586c78e07bb53dc39ed8c4ea Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 27 Jul 2025 18:17:01 +0200 Subject: [PATCH 09/17] Adjust flux when no particlesPerEdge is used This does not affects the cases of the IEPC, but I am also doing other stuff. --- src/modules/moduleInject.f90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index b6fc64b..000823e 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -199,17 +199,21 @@ MODULE moduleInject END DO + self%nParticles = SUM(self%particlesPerEdge) + ELSE ! No particles assigned per edge, use the species weight self%weightPerEdge = self%species%weight DO et = 1, self%nEdges - self%particlesPerEdge(et) = FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface /self%species%weight) - + self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight)) END DO - END IF + self%nParticles = SUM(self%particlesPerEdge) - self%nParticles = SUM(self%particlesPerEdge) + !Rescale weight to match flux + self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles)) + + END IF !Scale particles for different species steps IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') From 78cb9a245335ce68720c5f3e13cf97c7ccc0139c Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez - Galactica Date: Sat, 2 Aug 2025 10:31:06 +0200 Subject: [PATCH 10/17] No reflection of particles at injection, that should be a boundary condition --- src/modules/moduleInject.f90 | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 000823e..4fb855a 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -385,12 +385,6 @@ MODULE moduleInject self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /) - !If velocity is not in the right direction, invert it - IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN - partInj(n)%v = - partInj(n)%v - - END IF - !Obtain natural coordinates of particle in cell partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) !Push new particle with the minimum time step From 4b040c35c39064c0832dbc6631dc99b2385ed263 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 2 Aug 2025 13:25:48 +0200 Subject: [PATCH 11/17] Fixes to random variables After reading some works and reviewing what I had, I've done some corrections to how the randomb velicities in Maxwellian distributions are calculated. These should be correct now. --- src/modules/common/moduleRandom.f90 | 19 ++++++++++--------- src/modules/moduleInject.f90 | 4 ++-- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 7d9cec6..156f5e4 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -47,22 +47,23 @@ MODULE moduleRandom END FUNCTION randomIntAB - !Returns a random number in a Maxwellian distribution of mean 0 and width 1 + !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method FUNCTION randomMaxwellian() RESULT(rnd) USE moduleConstParam, ONLY: PI IMPLICIT NONE REAL(8):: rnd - REAL(8):: x, y + REAL(8):: v1, v2, Rsquare - rnd = 0.D0 - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - END DO - CALL RANDOM_NUMBER(y) + Rsquare = 1.D0 + do while (Rsquare >= 1.D0 .and. Rsquare > 0.D0) + v1 = 2.D0 * random() - 1.D0 + v2 = 2.D0 * random() - 1.D0 + Rsquare = v1**2 + v2**2 - rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + end do + + rnd = v2 * sqrt(-2.D0 * log(Rsquare) / Rsquare) END FUNCTION randomMaxwellian diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 4fb855a..5e72835 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -289,7 +289,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian() + v = sqrt(2.0)*self%vTh*randomMaxwellian() END FUNCTION randomVelMaxwellian @@ -302,7 +302,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomHalfMaxwellian() + v = sqrt(2.0)*self%vTh*randomHalfMaxwellian() END FUNCTION randomVelHalfMaxwellian From d86c8480f3cc772bd7bac19c4eb3e921a86fe32d Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez - Galactica Date: Sat, 2 Aug 2025 16:51:00 +0200 Subject: [PATCH 12/17] Fixed issue with some velocities. Still, at some point I need to rething all the injection thing. --- src/modules/moduleInject.f90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 5e72835..2a2c393 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -385,6 +385,13 @@ MODULE moduleInject self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /) + !If injecting a no-drift distribution and velocity is negative, reflect + if ((self%vMod == 0.D0) .and. & + (dot_product(direction, partInj(n)%v) < 0.D0)) then + partInj(n)%v = - partInj(n)%v + + end if + !Obtain natural coordinates of particle in cell partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) !Push new particle with the minimum time step From 7e193b9fa8133a105326c4e96c5f1b32779ebcb8 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 3 Aug 2025 15:32:55 +0200 Subject: [PATCH 13/17] Minor changes, no improvement made yet --- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 2 ++ src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 4 ++-- src/modules/solver/moduleSolver.f90 | 1 + 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index dbc8b25..c7308bf 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -701,6 +701,8 @@ MODULE moduleMesh2DCart pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%y(1:3)), & DOT_PRODUCT(dPsi(2,1:3),self%y(1:3)) /) + pDer(3,3) = 1.D0 + END FUNCTION partialDerTria !Gather electric field at position Xi diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index ae1eb92..bb94440 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -705,8 +705,8 @@ MODULE moduleMesh2DCyl dPsi = 0.D0 - dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index f85d812..df8d5e8 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -140,6 +140,7 @@ MODULE moduleSolver CASE('ElectrostaticBoltzmann') self%solveEM => solveElecFieldBoltzmann + END SELECT END SUBROUTINE initEM From 102fd013f3d30f39b425806bd1968d9d233de7b7 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 3 Aug 2025 20:46:12 +0200 Subject: [PATCH 14/17] Issue with triangles fixed Now they give the right electric field. I have to change 2DCyl. However, there was some insonsistency between the change of coordinates in phy2log and the Jacobian for the K matrix. I fixed it putting a transpose() in phy2log, but I don't like that solution. I need to review the basic procedure of phy2log. --- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index c7308bf..eb834e5 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -520,7 +520,7 @@ MODULE moduleMesh2DCart fPsi = self%fPsi(XiO, 4) f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), & DOT_PRODUCT(fPsi,self%y) /) - r(1:2) - Xi = XiO - MATMUL(invJ, f)/detJ + Xi = XiO - MATMUL(transpose(invJ), f)/detJ conv = MAXVAL(DABS(Xi-XiO),1) XiO = Xi @@ -831,14 +831,14 @@ MODULE moduleMesh2DCart REAL(8):: deltaR(1:3) REAL(8):: dPsi(1:3, 1:3) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ + REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates Xi = 0.D0 deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /) dPsi = self%dPsi(Xi, 3) pDer = self%partialDer(3, dPsi) - invJ = self%invJac(pDer) + invJ = transpose(self%invJac(pDer)) detJ = self%detJac(pDer) Xi = MATMUL(invJ,deltaR)/detJ @@ -915,8 +915,8 @@ MODULE moduleMesh2DCart invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) - invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) + invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /) invJ(3, 3) = 1.D0 END FUNCTION invJ2DCart From 3d7b1ce476a804d45fce35fbcb7a20d64c6d3186 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 3 Aug 2025 22:14:19 +0200 Subject: [PATCH 15/17] Fixed! So it seems that rectangles and triangles are now working properly. I have also checked the phy2log routine, now it seems a bit more complicated, but it is much clearer. Maybe in the future is worth rethinking to improve speed (specially for quad elements) --- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 58 ++++++++++---------- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 56 ++++++++++--------- 2 files changed, 58 insertions(+), 56 deletions(-) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index eb834e5..126d40c 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -495,34 +495,36 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Transform physical coordinates to element coordinates + !Transform physical coordinates to element coordinates with a Taylor series PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv = 1.D0 - XiO = 0.D0 + conv = 1.D0 + Xi0 = 0.D0 + Xi(3) = 0.D0 - f(3) = 0.D0 DO WHILE(conv > 1.D-4) - dPsi = self%dPsi(XiO, 4) - pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer) - invJ = self%invJac(pDer) - fPsi = self%fPsi(XiO, 4) - f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), & - DOT_PRODUCT(fPsi,self%y) /) - r(1:2) - Xi = XiO - MATMUL(transpose(invJ), f)/detJ - conv = MAXVAL(DABS(Xi-XiO),1) - XiO = Xi + fPsi = self%fPsi(Xi0, 4) + x0(1) = dot_product(fPsi, self%x) + x0(2) = dot_product(fPsi, self%y) + deltaR = r(1:2) - x0 + dPsi = self%dPsi(Xi0, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ + conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1) + Xi0(1:2) = Xi(1:2) END DO @@ -680,8 +682,8 @@ MODULE moduleMesh2DCart dPsi = 0.D0 - dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria @@ -701,8 +703,6 @@ MODULE moduleMesh2DCart pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%y(1:3)), & DOT_PRODUCT(dPsi(2,1:3),self%y(1:3)) /) - pDer(3,3) = 1.D0 - END FUNCTION partialDerTria !Gather electric field at position Xi @@ -828,19 +828,19 @@ MODULE moduleMesh2DCart CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:3) + REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) + REAL(8):: dPsi(1:3,1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - Xi = 0.D0 - deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /) - dPsi = self%dPsi(Xi, 3) - pDer = self%partialDer(3, dPsi) - invJ = transpose(self%invJac(pDer)) - detJ = self%detJac(pDer) - Xi = MATMUL(invJ,deltaR)/detJ + Xi(3) = 0.D0 + deltaR = (/ r(1) - self%x(1), r(2) - self%y(1) /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ END FUNCTION phy2logTria diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index bb94440..1f5f33c 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -510,34 +510,36 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad - !Transform physical coordinates to element coordinates + !Transform physical coordinates to element coordinates with a Taylor series PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv = 1.D0 - XiO = 0.D0 + conv = 1.D0 + Xi0 = 0.D0 + Xi(3) = 0.D0 - f(3) = 0.D0 DO WHILE(conv > 1.D-4) - dPsi = self%dPsi(XiO, 4) - pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer) - invJ = self%invJac(pDer) - fPsi = self%fPsi(XiO, 4) - f(1:2) = (/ DOT_PRODUCT(fPsi,self%z), & - DOT_PRODUCT(fPsi,self%r) /) - r(1:2) - Xi = XiO - MATMUL(invJ, f)/detJ - conv = MAXVAL(DABS(Xi-XiO),1) - XiO = Xi + fPsi = self%fPsi(Xi0, 4) + x0(1) = dot_product(fPsi, self%z) + x0(2) = dot_product(fPsi, self%r) + deltaR = r(1:2) - x0 + dPsi = self%dPsi(Xi0, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ + conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1) + Xi0(1:2) = Xi(1:2) END DO @@ -858,19 +860,19 @@ MODULE moduleMesh2DCyl CLASS(meshCell2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:3) + REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) + REAL(8):: dPsi(1:3,1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - Xi = 0.D0 - deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /) - dPsi = self%dPsi(Xi, 3) - pDer = self%partialDer(3, dPsi) - invJ = self%invJac(pDer) - detJ = self%detJac(pDer) - Xi = MATMUL(invJ,deltaR)/detJ + Xi(3) = 0.D0 + deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ END FUNCTION phy2logTria @@ -948,8 +950,8 @@ MODULE moduleMesh2DCyl invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) - invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) + invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /) invJ(3, 3) = 1.D0 END FUNCTION invJ2DCyl From 5166a650d28a65524800cbfc81a0127a20b248de Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 6 Aug 2025 10:59:03 +0200 Subject: [PATCH 16/17] Data for Kr, just to test --- data/collisions/IO_e-Kr.dat | 52 +++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 data/collisions/IO_e-Kr.dat diff --git a/data/collisions/IO_e-Kr.dat b/data/collisions/IO_e-Kr.dat new file mode 100644 index 0000000..532ca4d --- /dev/null +++ b/data/collisions/IO_e-Kr.dat @@ -0,0 +1,52 @@ +# D108525 "refs": {"B56": {"note": "CLM-R294 (1989)"}} +# Relative energy (eV) cross section (m^2) +1.40E+01 0 +1.62E+01 7.249E-21 +1.88E+01 1.199E-20 +2.18E+01 1.644E-20 +2.53E+01 2.1E-20 +2.94E+01 2.542E-20 +3.41E+01 2.937E-20 +3.95E+01 3.26E-20 +4.58E+01 3.499E-20 +5.32E+01 3.653E-20 +6.17E+01 3.726E-20 +7.15E+01 3.728E-20 +8.29E+01 3.671E-20 +9.62E+01 3.566E-20 +1.12E+02 3.426E-20 +1.29E+02 3.259E-20 +1.50E+02 3.075E-20 +1.74E+02 2.881E-20 +2.02E+02 2.682E-20 +2.34E+02 2.484E-20 +2.72E+02 2.289E-20 +3.15E+02 2.101E-20 +3.65E+02 1.922E-20 +4.24E+02 1.751E-20 +4.91E+02 1.592E-20 +5.70E+02 1.443E-20 +6.61E+02 1.305E-20 +7.67E+02 1.177E-20 +8.89E+02 1.06E-20 +1.03E+03 9.526E-21 +1.20E+03 8.547E-21 +1.39E+03 7.658E-21 +1.61E+03 6.851E-21 +1.87E+03 6.121E-21 +2.16E+03 5.462E-21 +2.51E+03 4.868E-21 +2.91E+03 4.334E-21 +3.38E+03 3.855E-21 +3.92E+03 3.426E-21 +4.54E+03 3.041E-21 +5.27E+03 2.698E-21 +6.11E+03 2.391E-21 +7.09E+03 2.118E-21 +8.22E+03 1.875E-21 +9.53E+03 1.658E-21 +1.11E+04 1.466E-21 +1.28E+04 1.295E-21 +1.49E+04 1.143E-21 +1.72E+04 1.009E-21 +2.00E+04 8.898E-22 From dff9a87f0d5c39b74076f81537f6313457d9207a Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 8 Aug 2025 19:27:27 +0200 Subject: [PATCH 17/17] Implemenint injecting particles without direction I was almost sure this was implemented in the past, but it was not working. Now, if n = 0 or if n is not provided, particles are injected with the normal to the surface. --- src/modules/moduleInject.f90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 2a2c393..528cb91 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -377,7 +377,13 @@ MODULE moduleInject !Assign particle type partInj(n)%species => self%species - direction = self%n + if (all(self%n == 0.D0)) then + direction = randomEdge%normal + + else + direction = self%n + + end if partInj(n)%v = 0.D0