From 21184e91d36579edead3f0f1944e3cede270987e Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 17 Jul 2023 12:02:24 +0200 Subject: [PATCH 01/11] Type for SEE Implementation of the type for Secondary Electron Emission (SEE) --- src/modules/moduleBoundary.f90 | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 83c815c..9fb212f 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -46,6 +46,14 @@ MODULE moduleBoundary END TYPE boundaryIonization + !Secondary electron emission (by ion impact) + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundarySEE + !Yield as a function of ion energy + TYPE(table1D):: yield + CONTAINS + + END TYPE boundarySEE + !Symmetry axis TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis CONTAINS @@ -137,4 +145,14 @@ MODULE moduleBoundary END SUBROUTINE initIonization + SUBROUTINE initSEE(boundary, tableFile) + IMPLICIT NONE + + CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary + CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile + + ALLOCATE(boundarySEE:: boundary) + + END SUBROUTINE initSEE + END MODULE moduleBoundary From e369bccf78ff5343c8df51a342daf8fe6f7ee92b Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 17 Jul 2023 13:58:57 +0200 Subject: [PATCH 02/11] Function to create electrons Still required to assign velocity: - In the direction normal to the surface - Which energy? --- src/modules/mesh/moduleMeshBoundary.f90 | 75 +++++++++++++++++++++++++ src/modules/moduleBoundary.f90 | 14 ++++- 2 files changed, 88 insertions(+), 1 deletion(-) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 4f72e10..34dc8a7 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -212,6 +212,79 @@ MODULE moduleMeshBoundary END SUBROUTINE symmetryAxis + SUBROUTINE secondaryEmission(edge, part) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + REAL(8):: vRel, eRel + INTEGER:: yield + INTEGER:: p + REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons + INTEGER:: cell + TYPE(particle), POINTER:: secondaryElectron + + SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj) + TYPE IS(boundarySEE) + !Get relative velocity + vRel = NORM2(part%v) + !Convert to relative energy + eRel = part%species%m*vRel**2*5.D-1 + + !Get number of secondary electrons macro-particles + yield = INT(part%weight*bound%yield%get(eRel) / bound%electron%weight) + + + !position of impacting ion + rElectron = edge%intersection(part%r) + XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron) + + !New cell of origin + IF (ASSOCIATED(edge%e1)) THEN + cell = edge%e1%n + + ELSEIF (ASSOCIATED(edge%e2)) THEN + cell = edge%e2%n + + END IF + + DO p = 1, yield + !Create new macro-particle + ALLOCATE(secondaryElectron) + + !Assign species to electron + secondaryElectron%species => bound%electron + + !Assign position to particle + secondaryElectron%r = rElectron + secondaryElectron%Xi = XiElectron + + !Assign cell to electron + secondaryElectron%cell = cell + + !Assign weight + secondaryElectron%weight = bound%electron%weight + + !Assume particle is inside the numerical domain + secondaryElectron%n_in = .TRUE. + + !Assign velocity + + !Add particle to list + CALL partSurfaces%setLock() + CALL partSurfaces%add(secondaryElectron) + CALL partSurfaces%unsetLock() + + END DO + + !Absorb impacting particle + CALL absorption(edge, part) + + END SELECT + + END SUBROUTINE secondaryEmission + !Points the boundary function to specific type SUBROUTINE pointBoundaryFunction(edge, s) USE moduleErrors @@ -239,6 +312,8 @@ MODULE moduleMeshBoundary TYPE IS(boundaryAxis) edge%fBoundary(s)%apply => symmetryAxis + TYPE IS(boundarySEE) + edge%fBoundary(s)%apply => secondaryEmission CLASS DEFAULT CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction') diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 9fb212f..bfdbd70 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -50,6 +50,7 @@ MODULE moduleBoundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundarySEE !Yield as a function of ion energy TYPE(table1D):: yield + CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission CONTAINS END TYPE boundarySEE @@ -145,13 +146,24 @@ MODULE moduleBoundary END SUBROUTINE initIonization - SUBROUTINE initSEE(boundary, tableFile) + SUBROUTINE initSEE(boundary, tableFile, speciesID) + USE moduleRefParam + USE moduleConstParam + USE moduleSpecies IMPLICIT NONE CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile + INTEGER:: speciesID ALLOCATE(boundarySEE:: boundary) + SELECT TYPE(boundary) + TYPE IS(boundarySEE) + CALL boundary%yield%init(tableFile) + CALL boundary%yield%convert(eV2J/(m_ref*v_ref**2), 1.D0) + boundary%electron => species(speciesID)%obj + + END SELECT END SUBROUTINE initSEE From 0f7d9919ec1ab769609b1f61c3e5397b4e8a83a0 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 17 Jul 2023 16:18:51 +0200 Subject: [PATCH 03/11] Ignore backups of bibliography Sorry that this change is in this branch, but I just noticed that the backups from jabref were not being ignored. --- doc/user-manual/.gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/user-manual/.gitignore b/doc/user-manual/.gitignore index b2ba4d5..61ad7e0 100644 --- a/doc/user-manual/.gitignore +++ b/doc/user-manual/.gitignore @@ -19,3 +19,4 @@ fpakc_UserManual-blx.bib *.gls *.ist fpakc_UserManual.run.xml +*.bib.sav From 38fa37c9957cc152439bbe9ab58f40b72167ef04 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 17 Jul 2023 16:29:52 +0200 Subject: [PATCH 04/11] Dummy velocity for testing I have set the velocity of the particle to be a dummy Maxwellian distribution to start testing things. Also I added a constant yield table for testing. Next step is to do the user input and run some tests. --- data/see/constant.dat | 4 ++++ src/modules/mesh/moduleMeshBoundary.f90 | 5 ++++- 2 files changed, 8 insertions(+), 1 deletion(-) create mode 100644 data/see/constant.dat diff --git a/data/see/constant.dat b/data/see/constant.dat new file mode 100644 index 0000000..060c8ad --- /dev/null +++ b/data/see/constant.dat @@ -0,0 +1,4 @@ +#Relative energy (eV) yield () +0.000 1.000E-01 +1.000 1.000E-01 + diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 34dc8a7..ffd1af9 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -193,7 +193,7 @@ MODULE moduleMeshBoundary END SELECT - !Removes ionizing electron regardless the number of pair created + !Removes ionizing electron regardless the number of pairs created part%n_in = .FALSE. END SUBROUTINE ionization @@ -214,6 +214,7 @@ MODULE moduleMeshBoundary SUBROUTINE secondaryEmission(edge, part) USE moduleSpecies + USE moduleRandom IMPLICIT NONE CLASS(meshEdge), INTENT(inout):: edge @@ -270,6 +271,7 @@ MODULE moduleMeshBoundary secondaryElectron%n_in = .TRUE. !Assign velocity + secondaryElectron%v = 2.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /) !Add particle to list CALL partSurfaces%setLock() @@ -314,6 +316,7 @@ MODULE moduleMeshBoundary TYPE IS(boundarySEE) edge%fBoundary(s)%apply => secondaryEmission + CLASS DEFAULT CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction') From cc3e28e5e7d86b9635b22a63dcc22d82dc69d1b3 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 17 Jul 2023 22:48:52 +0200 Subject: [PATCH 05/11] Input done The input is done and testing. WARNING: Velocity of secondary electrons is still a dummy one! --- src/modules/init/moduleInput.f90 | 11 ++++++++++- src/modules/mesh/moduleMeshBoundary.f90 | 21 +++++++++++++++++---- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 9891806..c867e23 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -806,7 +806,7 @@ MODULE moduleInput REAL(8):: effTime REAL(8):: eThreshold !Energy threshold INTEGER:: speciesID - CHARACTER(:), ALLOCATABLE:: speciesName, crossSection + CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield LOGICAL:: found INTEGER:: nTypes @@ -872,6 +872,15 @@ MODULE moduleInput CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw) + CASE('secondaryEmission') + CALL config%get(object // '.yield', yield, found) + IF (.NOT. found) CALL criticalError("missing parameter 'yield' for secondary emission", 'readBoundary') + CALL config%get(object // '.electron', speciesName, found) + IF (.NOT. found) CALL criticalError("missing parameter 'electron' for secondary emission", 'readBoundary') + speciesID = speciesName2Index(speciesName) + + CALL initSEE(boundary(i)%bTypes(s)%obj, yield, speciesID) + CASE('axis') ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index ffd1af9..90420c6 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -220,7 +220,9 @@ MODULE moduleMeshBoundary CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part REAL(8):: vRel, eRel - INTEGER:: yield + REAL(8):: yield + INTEGER:: nNewElectrons + REAL(8), ALLOCATABLE:: weight(:) INTEGER:: p REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons INTEGER:: cell @@ -234,9 +236,20 @@ MODULE moduleMeshBoundary eRel = part%species%m*vRel**2*5.D-1 !Get number of secondary electrons macro-particles - yield = INT(part%weight*bound%yield%get(eRel) / bound%electron%weight) + yield = part%weight*bound%yield%get(eRel) + nNewElectrons = FLOOR(yield / bound%electron%weight) + IF (REAL(nNewElectrons) * bound%electron%weight < yield) THEN + nNewElectrons = nNewElectrons + 1 + ALLOCATE(weight(1:nNewElectrons)) + weight(1:nNewElectrons-1) = bound%electron%weight + weight(nNewElectrons) = yield - SUM(weight(1:nNewElectrons-1)) + ELSE + ALLOCATE(weight(1:nNewElectrons)) + weight(1:nNewElectrons) = bound%electron%weight + END IF + !position of impacting ion rElectron = edge%intersection(part%r) XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron) @@ -250,7 +263,7 @@ MODULE moduleMeshBoundary END IF - DO p = 1, yield + DO p = 1, nNewElectrons !Create new macro-particle ALLOCATE(secondaryElectron) @@ -265,7 +278,7 @@ MODULE moduleMeshBoundary secondaryElectron%cell = cell !Assign weight - secondaryElectron%weight = bound%electron%weight + secondaryElectron%weight = weight(p) !Assume particle is inside the numerical domain secondaryElectron%n_in = .TRUE. From e1009a21df141e89f758eb2f57656886bf69ea1e Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Tue, 18 Jul 2023 15:35:16 +0200 Subject: [PATCH 06/11] Correction due to angle Moving forward to a complete model for SEE. Now the yield has a correction for the incident angle on the surface (still to check if this is correct) Velocity distribution is still a dummy one. --- src/modules/common/moduleTable.f90 | 3 + src/modules/mesh/moduleMeshBoundary.f90 | 126 ++++++++++++++---------- src/modules/moduleBoundary.f90 | 12 +++ 3 files changed, 89 insertions(+), 52 deletions(-) diff --git a/src/modules/common/moduleTable.f90 b/src/modules/common/moduleTable.f90 index e4da335..ac5cbc7 100644 --- a/src/modules/common/moduleTable.f90 +++ b/src/modules/common/moduleTable.f90 @@ -30,6 +30,8 @@ MODULE moduleTable READ(id, '(A)', iostat = stat) dummy !If EOF or error, exit file IF (stat /= 0) EXIT + !If empty line, exit file + IF (dummy == "") EXIT !Skip comment IF (INDEX(dummy,'#') /= 0) CYCLE !Add row @@ -55,6 +57,7 @@ MODULE moduleTable !TODO: Make this a function IF (INDEX(dummy,'#') /= 0) CYCLE IF (stat /= 0) EXIT + IF (dummy == "") EXIT !Add data !TODO: substitute with extracting information from dummy BACKSPACE(id) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 90420c6..d130a8a 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -212,19 +212,24 @@ MODULE moduleMeshBoundary END SUBROUTINE symmetryAxis + !Secondary emission of electrons by impacting particle. SUBROUTINE secondaryEmission(edge, part) USE moduleSpecies USE moduleRandom + USE moduleConstParam IMPLICIT NONE CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part REAL(8):: vRel, eRel + REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons (impacting particle) + REAL(8), DIMENSION(1:3):: rIncident !Vector from imapcting particle position to particle position + REAL(8):: theta !incident angle REAL(8):: yield + REAL(8):: energy !incident energy corrected by threshold and INTEGER:: nNewElectrons REAL(8), ALLOCATABLE:: weight(:) INTEGER:: p - REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons INTEGER:: cell TYPE(particle), POINTER:: secondaryElectron @@ -235,64 +240,81 @@ MODULE moduleMeshBoundary !Convert to relative energy eRel = part%species%m*vRel**2*5.D-1 - !Get number of secondary electrons macro-particles - yield = part%weight*bound%yield%get(eRel) - nNewElectrons = FLOOR(yield / bound%electron%weight) - IF (REAL(nNewElectrons) * bound%electron%weight < yield) THEN - nNewElectrons = nNewElectrons + 1 - ALLOCATE(weight(1:nNewElectrons)) - weight(1:nNewElectrons-1) = bound%electron%weight - weight(nNewElectrons) = yield - SUM(weight(1:nNewElectrons-1)) + !If energy is abound threshold calculate secondary electrons + IF (eRel >= bound%energyThreshold) THEN - ELSE - ALLOCATE(weight(1:nNewElectrons)) - weight(1:nNewElectrons) = bound%electron%weight - - END IF + !position of impacting ion + rElectron = edge%intersection(part%r) + XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron) - !position of impacting ion - rElectron = edge%intersection(part%r) - XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron) + !Calculate incident angle + rIncident = part%r - rElectron + theta = ACOS(DOT_PRODUCT(rIncident, edge%normal) / (NORM2(rIncident) * NORM2(edge%normal))) - !New cell of origin - IF (ASSOCIATED(edge%e1)) THEN - cell = edge%e1%n + !Calculate incident energy + energy = (eRel - bound%energyThreshold) * PI2 / (PI2 + theta**2) + bound%energyThreshold - ELSEIF (ASSOCIATED(edge%e2)) THEN - cell = edge%e2%n + !Get number of secondary electrons particles + yield = part%weight*bound%yield%get(eRel) * (1.D0 * theta**2 / PI2) !Check equation! + + !Convert the number to macro-particles + nNewElectrons = FLOOR(yield / bound%electron%weight) + + !If the weight of new macro-particles is below the yield, correct adding a particle + IF (REAL(nNewElectrons) * bound%electron%weight < yield) THEN + nNewElectrons = nNewElectrons + 1 + ALLOCATE(weight(1:nNewElectrons)) + weight(1:nNewElectrons-1) = bound%electron%weight + weight(nNewElectrons) = yield - SUM(weight(1:nNewElectrons-1)) + + ELSE + ALLOCATE(weight(1:nNewElectrons)) + weight(1:nNewElectrons) = bound%electron%weight + + END IF + + !New cell of origin + IF (ASSOCIATED(edge%e1)) THEN + cell = edge%e1%n + + ELSEIF (ASSOCIATED(edge%e2)) THEN + cell = edge%e2%n + + END IF + + !Create the new electron macro-particles + DO p = 1, nNewElectrons + !Create new macro-particle + ALLOCATE(secondaryElectron) + + !Assign species to electron + secondaryElectron%species => bound%electron + + !Assign position to particle + secondaryElectron%r = rElectron + secondaryElectron%Xi = XiElectron + + !Assign cell to electron + secondaryElectron%cell = cell + + !Assign weight + secondaryElectron%weight = weight(p) + + !Assume particle is inside the numerical domain + secondaryElectron%n_in = .TRUE. + + !Assign velocity + secondaryElectron%v = 2.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /) + + !Add particle to list + CALL partSurfaces%setLock() + CALL partSurfaces%add(secondaryElectron) + CALL partSurfaces%unsetLock() + + END DO END IF - DO p = 1, nNewElectrons - !Create new macro-particle - ALLOCATE(secondaryElectron) - - !Assign species to electron - secondaryElectron%species => bound%electron - - !Assign position to particle - secondaryElectron%r = rElectron - secondaryElectron%Xi = XiElectron - - !Assign cell to electron - secondaryElectron%cell = cell - - !Assign weight - secondaryElectron%weight = weight(p) - - !Assume particle is inside the numerical domain - secondaryElectron%n_in = .TRUE. - - !Assign velocity - secondaryElectron%v = 2.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /) - - !Add particle to list - CALL partSurfaces%setLock() - CALL partSurfaces%add(secondaryElectron) - CALL partSurfaces%unsetLock() - - END DO - !Absorb impacting particle CALL absorption(edge, part) diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index bfdbd70..b3c1d12 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -51,6 +51,7 @@ MODULE moduleBoundary !Yield as a function of ion energy TYPE(table1D):: yield CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission + REAL(8):: energyThreshold CONTAINS END TYPE boundarySEE @@ -155,6 +156,7 @@ MODULE moduleBoundary CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile INTEGER:: speciesID + INTEGER:: i ALLOCATE(boundarySEE:: boundary) SELECT TYPE(boundary) @@ -162,6 +164,16 @@ MODULE moduleBoundary CALL boundary%yield%init(tableFile) CALL boundary%yield%convert(eV2J/(m_ref*v_ref**2), 1.D0) boundary%electron => species(speciesID)%obj + !Search for the threshold energy in the table + DO i = 1, SIZE(boundary%yield%f) + IF (boundary%yield%f(i) > 0.D0) THEN + boundary%energyThreshold = boundary%yield%x(i) + + EXIT + + END IF + + END DO END SELECT From 541a6ff97a5d4ac8749c113c2d68ffca3ce4dac6 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 30 Jul 2023 12:55:52 +0200 Subject: [PATCH 07/11] Correction to injection velocity Small correction to injection velocity. Still not fully satisfied with it. I have to find a way to ensure that velocity is injected in the right direction and fulfills the distribution functions in each direction regardless the direction of injection, mean velocity or distribution function. --- src/modules/common/moduleRandom.f90 | 18 ++++++++++++++++++ src/modules/moduleInject.f90 | 19 +++++++++---------- 2 files changed, 27 insertions(+), 10 deletions(-) diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index cd553a8..657ca28 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -66,6 +66,24 @@ MODULE moduleRandom END FUNCTION randomMaxwellian + FUNCTION randomHalfMaxwellian() RESULT(rnd) + USE moduleConstParam, ONLY: PI + IMPLICIT NONE + + REAL(8):: rnd + REAL(8):: x, y + + rnd = 0.D0 + x = 0.D0 + DO WHILE (x == 0.D0) + CALL RANDOM_NUMBER(x) + END DO + y = random(-PI, PI) + + rnd = DSQRT(-2.D0*DLOG(x))*DCOS(y) + + 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 4e57083..644cfa7 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -254,10 +254,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 @@ -334,17 +331,19 @@ MODULE moduleInject direction = self%n - partInj(n)%v = 0.D0 - + !Sample initial velocity partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & 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 + !For each direction, velocities have to agree with the direction of injection + DO i = 1, 3 + DO WHILE (partInj(n)%v(i)*direction(i) < 0) + partInj(n)%v(i) = self%vMod*direction(i) + self%v(i)%obj%randomVel() - END IF + END DO + + END DO !Obtain natural coordinates of particle in cell partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) From db804bebc13a9d97e05f38949a1ac34c24fc8ac6 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 13 Jul 2024 23:34:43 +0200 Subject: [PATCH 08/11] Identifying issues I need a way to compute the energy (velocity) for the emitted particles. I think I also need a better way to calculate the velocity distribution of particles being injected or emitted from a surface. --- src/modules/mesh/moduleMeshBoundary.f90 | 8 +++++--- src/modules/moduleInject.f90 | 11 ++++------- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 5320aa0..37fa3cc 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -261,12 +261,12 @@ MODULE moduleMeshBoundary energy = (eRel - bound%energyThreshold) * PI2 / (PI2 + theta**2) + bound%energyThreshold !Get number of secondary electrons particles - yield = part%weight*bound%yield%get(eRel) * (1.D0 * theta**2 / PI2) !Check equation! + yield = part%weight*bound%yield%get(eRel) * (1.D0 * theta**2 / PI2) !Convert the number to macro-particles nNewElectrons = FLOOR(yield / bound%electron%weight) - !If the weight of new macro-particles is below the yield, correct adding a particle + !If the weight of new macro-particles is below the yield, add an additional particle with less weight IF (REAL(nNewElectrons) * bound%electron%weight < yield) THEN nNewElectrons = nNewElectrons + 1 ALLOCATE(weight(1:nNewElectrons)) @@ -310,7 +310,9 @@ MODULE moduleMeshBoundary secondaryElectron%n_in = .TRUE. !Assign velocity - secondaryElectron%v = 2.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /) + !TODO: Calculate energy and velocity + !TODO: Better way to inject particles from surface. For this and for general injection + secondaryElectron%v = 1.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /) !Add particle to list CALL partSurfaces%setLock() diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index f4c9c46..41f7626 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -367,14 +367,11 @@ MODULE moduleInject self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /) - !For each direction, velocities have to agree with the direction of injection - DO j = 1, 3 - DO WHILE (partInj(n)%v(i)*direction(i) < 0) - partInj(n)%v(i) = self%vMod*direction(i) + self%v(i)%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 DO - - END DO + END IF !Obtain natural coordinates of particle in cell partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) From b171f801fcb36f905e2b14af545d553af4725869 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 27 Jul 2024 21:04:48 +0200 Subject: [PATCH 09/11] Remove old code from moduleRandom This should not be here. --- src/modules/common/moduleRandom.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 0cbdd1d..9134815 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -105,7 +105,6 @@ MODULE moduleRandom END IF END DO - ! rnd = MINLOC(DABS(rnd0b - cumWeight), 1) END FUNCTION randomWeighted From 94ff12b1cbc9b7d6b012d07c6042dd1e0193d492 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 27 Jul 2024 21:06:59 +0200 Subject: [PATCH 10/11] New subroutines in table These allow to invert the x and f arrays and make the X axis cumulative sum. --- src/modules/common/moduleTable.f90 | 48 ++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/modules/common/moduleTable.f90 b/src/modules/common/moduleTable.f90 index 394fe37..a6b8c57 100644 --- a/src/modules/common/moduleTable.f90 +++ b/src/modules/common/moduleTable.f90 @@ -8,6 +8,8 @@ MODULE moduleTable PROCEDURE, PASS:: init => initTable1D PROCEDURE, PASS:: get => getValueTable1D PROCEDURE, PASS:: convert => convertUnits + PROCEDURE, PASS:: invertXF + PROCEDURE, PASS:: cumsumX END TYPE table1D @@ -126,4 +128,50 @@ MODULE moduleTable END SUBROUTINE convertUnits + ! Invert the arrays of x and f + SUBROUTINE invertXF(self) + IMPLICIT NONE + + CLASS(table1D), INTENT(inout):: self + REAL(8), ALLOCATABLE:: xTemp(:) + REAL(8):: xTempMin, xTempMax + + ! Store temporal x data + xTemp = self%x + xTempMin = self%xMin + xTempMax = self%xMax + + ! Replace x data with f data + self%x = self%f + self%xMin = self%fMin + self%xMax = self%fMax + + ! Put temporal x data in f data + self%f = xTemp + self%fMin = xTempMin + self%fMax = xTempMax + + END SUBROUTINE invertXF + + ! Makes the x axis its cumulative sum + SUBROUTINE cumSumX(self) + IMPLICIT NONE + + CLASS(table1D), INTENT(inout):: self + INTEGER:: nTotal, n + REAL(8), ALLOCATABLE:: cumX(:) + + nTotal = SIZE(self%x) + ALLOCATE(cumX(1:nTotal)) + cumX(1) = self%x(1) + DO n = 2, nTotal + cumX(n) = self%x(n) + cumX(n-1) + + END DO + + self%x = cumX / cumX(nTotal) + DEALLOCATE(cumX) + + END SUBROUTINE cumSumX + END MODULE moduleTable From 0f4e44459d7b6390feea6ee066d4e4b07cadae57 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 13 Oct 2024 17:22:47 +0200 Subject: [PATCH 11/11] Making some progress Things are better organized. --- data/see/constant_energy.dat | 3 +++ src/modules/init/moduleInput.f90 | 11 +++++++--- src/modules/mesh/moduleMeshBoundary.f90 | 4 ++-- src/modules/moduleBoundary.f90 | 27 +++++++++++-------------- 4 files changed, 25 insertions(+), 20 deletions(-) create mode 100644 data/see/constant_energy.dat diff --git a/data/see/constant_energy.dat b/data/see/constant_energy.dat new file mode 100644 index 0000000..d9de3ba --- /dev/null +++ b/data/see/constant_energy.dat @@ -0,0 +1,3 @@ +#Secondary energy (eV) Probability (a. u.) +0.1 0.5 +1.0 0.5 diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 618ba53..31d9959 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -95,7 +95,10 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(LEN=*), INTENT(in):: step - IF (config%failed()) CALL criticalError("Error reading the JSON input file", TRIM(step)) + IF (config%failed()) THEN + CALL criticalError("Error reading the JSON input file", TRIM(step)) + + END IF END SUBROUTINE checkStatus @@ -808,7 +811,7 @@ MODULE moduleInput REAL(8):: effTime REAL(8):: eThreshold !Energy threshold INTEGER:: speciesID, electronSecondaryID - CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield, electronSecondary + CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield, energySpectrum, electronSecondary LOGICAL:: found INTEGER:: nTypes @@ -886,11 +889,13 @@ MODULE moduleInput CASE('secondaryEmission') CALL config%get(object // '.yield', yield, found) IF (.NOT. found) CALL criticalError("missing parameter 'yield' for secondary emission", 'readBoundary') + CALL config%get(object // '.energySpectrum', energySpectrum, found) + IF (.NOT. found) CALL criticalError("missing parameter 'energySpectrum' for secondary emission", 'readBoundary') CALL config%get(object // '.electron', speciesName, found) IF (.NOT. found) CALL criticalError("missing parameter 'electron' for secondary emission", 'readBoundary') speciesID = speciesName2Index(speciesName) - CALL initSEE(boundary(i)%bTypes(s)%obj, yield, speciesID) + CALL initSEE(boundary(i)%bTypes(s)%obj, yield, energySpectrum, speciesID) CASE('axis') ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 37fa3cc..618c61e 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -232,7 +232,7 @@ MODULE moduleMeshBoundary REAL(8), DIMENSION(1:3):: rIncident !Vector from imapcting particle position to particle position REAL(8):: theta !incident angle REAL(8):: yield - REAL(8):: energy !incident energy corrected by threshold and + REAL(8):: energyIncident !incident energy corrected by threshold INTEGER:: nNewElectrons REAL(8), ALLOCATABLE:: weight(:) INTEGER:: p @@ -258,7 +258,7 @@ MODULE moduleMeshBoundary theta = ACOS(DOT_PRODUCT(rIncident, edge%normal) / (NORM2(rIncident) * NORM2(edge%normal))) !Calculate incident energy - energy = (eRel - bound%energyThreshold) * PI2 / (PI2 + theta**2) + bound%energyThreshold + energyIncident = (eRel - bound%energyThreshold) * PI2 / (PI2 + theta**2) + bound%energyThreshold !Get number of secondary electrons particles yield = part%weight*bound%yield%get(eRel) * (1.D0 * theta**2 / PI2) diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 0a95f5e..091cccf 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -49,8 +49,8 @@ MODULE moduleBoundary !Secondary electron emission (by ion impact) TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundarySEE - !Yield as a function of ion energy - TYPE(table1D):: yield + TYPE(table1D):: yield !Yield as a function of ion energy + TYPE(table1D):: energySpectrum !Spectrum of secondary particle emitted CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission REAL(8):: energyThreshold CONTAINS @@ -165,33 +165,30 @@ MODULE moduleBoundary END SUBROUTINE initIonization - SUBROUTINE initSEE(boundary, tableFile, speciesID) + SUBROUTINE initSEE(boundary, yieldFile, energySpectrumFile, speciesID) USE moduleRefParam USE moduleConstParam USE moduleSpecies IMPLICIT NONE CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary - CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile + CHARACTER(:), ALLOCATABLE, INTENT(in):: yieldFile + CHARACTER(:), ALLOCATABLE, INTENT(in):: energySpectrumFile INTEGER:: speciesID INTEGER:: i ALLOCATE(boundarySEE:: boundary) SELECT TYPE(boundary) TYPE IS(boundarySEE) - CALL boundary%yield%init(tableFile) + CALL boundary%yield%init(yieldFile) CALL boundary%yield%convert(eV2J/(m_ref*v_ref**2), 1.D0) boundary%electron => species(speciesID)%obj - !Search for the threshold energy in the table - DO i = 1, SIZE(boundary%yield%f) - IF (boundary%yield%f(i) > 0.D0) THEN - boundary%energyThreshold = boundary%yield%x(i) - - EXIT - - END IF - - END DO + ! Use first value from table as threshold + boundary%energyThreshold = boundary%yield%xMin + CALL boundary%energySpectrum%init(energySpectrumFile) + CALL boundary%energySpectrum%convert(eV2J/(m_ref*v_ref**2), 1.D0) + CALL boundary%energySpectrum%invertXF() + CALL boundary%energySpectrum%cumSumX() END SELECT