From 21184e91d36579edead3f0f1944e3cede270987e Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 17 Jul 2023 12:02:24 +0200 Subject: [PATCH 01/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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 b73c8531e852552ede1c421cc90d7068b109969b Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 6 Aug 2023 19:45:48 +0200 Subject: [PATCH 08/10] Fix an issue with different time steps Sometimes the division between number of steps and time step was not giving the right results. Nevertheless, this just indicates that the species have to be in separated arrays so that the assignment of particles in list, pushing and scattering can be dealt independently. Thus, this is the first step in creating separate arrays of particles per species. --- data/see/{constant.dat => constant_yield.dat} | 0 src/modules/mesh/moduleMesh.f90 | 3 +-- src/modules/solver/moduleSolver.f90 | 7 ++++++- 3 files changed, 7 insertions(+), 3 deletions(-) rename data/see/{constant.dat => constant_yield.dat} (100%) diff --git a/data/see/constant.dat b/data/see/constant_yield.dat similarity index 100% rename from data/see/constant.dat rename to data/see/constant_yield.dat diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 7e30326..e8e38f7 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -816,9 +816,8 @@ MODULE moduleMesh IF (MOD(t, everyColl) == 0) THEN !Collisions need to be performed in this iteration - !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) + !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j, cell) DO e=1, self%numCells - cell => self%cells(e)%obj !TODO: Simplify this, to many sublevels diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index e539fed..de544ba 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -123,7 +123,12 @@ MODULE moduleSolver END SELECT self%pushSpecies = .FALSE. - self%every = INT(tauSp/tau) + self%every = FLOOR(tauSp/tau) + !Correct value if not fulfilled + IF (tau*REAL(self%every) < tauSp) THEN + self%every = self%every + 1 + + END IF END SUBROUTINE initPusher From 286e858d66f6b18222af29faa2b17c357fe4b374 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 22 Sep 2023 18:18:46 +0200 Subject: [PATCH 09/10] Intermediate push I still have to figure out how to do this properly, but I am tired of working on my laptop. --- src/modules/init/moduleInput.f90 | 40 ++----- src/modules/mesh/moduleMesh.f90 | 6 +- src/modules/mesh/moduleMeshBoundary.f90 | 4 +- src/modules/moduleCollisions.f90 | 12 +- src/modules/moduleCoulomb.f90 | 10 +- src/modules/moduleInject.f90 | 20 ++-- src/modules/moduleSpecies.f90 | 74 ++++++++---- src/modules/output/moduleOutput.f90 | 4 +- src/modules/solver/moduleSolver.f90 | 133 ++++++++++++--------- src/modules/solver/pusher/modulePusher.f90 | 2 +- 10 files changed, 166 insertions(+), 139 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index c867e23..9f16195 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -339,10 +339,9 @@ MODULE moduleInput !Mean velocity and temperature at particle position REAL(8):: velocityXi(1:3), temperatureXi INTEGER:: nNewPart = 0.D0 - CLASS(meshCell), POINTER:: cell TYPE(particle), POINTER:: partNew + CLASS(meshCell), POINTER:: cell REAL(8):: vTh - TYPE(lNode), POINTER:: partCurr, partNext CALL config%info('solver.initial', found, n_children = nInitial) @@ -371,9 +370,12 @@ MODULE moduleInput !Calculate number of particles nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) - !Allocate new particles + !Allocate array of particles for this species + ALLOCATE(partOld(sp)%p(1:nNewPart)) + + !Create new particles DO p = 1, nNewPart - ALLOCATE(partNew) + partNew = partOld(sp)%p(p) partNew%species => species(sp)%obj partNew%r = mesh%cells(e)%obj%randPos() partNew%Xi = mesh%cells(e)%obj%phy2log(partNew%r) @@ -408,9 +410,6 @@ MODULE moduleInput partNew%weight = species(sp)%obj%weight - !Assign particle to temporal list of particles - CALL partInitial%add(partNew) - !Assign particle to list in volume IF (listInCells) THEN cell => meshforMCC%cells(partNew%cellColl)%obj @@ -430,30 +429,10 @@ MODULE moduleInput END DO + nPartOld(sp) = SIZE(partOld(sp)%p) + END DO - !Convert temporal list of particles into initial partOld array - !Deallocate the list of initial particles - nNewPart = partInitial%amount - IF (nNewPart > 0) THEN - ALLOCATE(partOld(1:nNewPart)) - partCurr => partInitial%head - DO p = 1, nNewPart - partNext => partCurr%next - partOld(p) = partCurr%part - DEALLOCATE(partCurr) - partCurr => partNext - - END DO - - IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head) - IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail) - partInitial%amount = 0 - - END IF - - nPartOld = SIZE(partOld) - END IF END SUBROUTINE readInitial @@ -600,7 +579,10 @@ MODULE moduleInput END DO + !Allocate the wrapper array that contains particles + ALLOCATE(partOld(1:nSpecies)) !Set number of particles to 0 for init state + ALLOCATE(nPartOld(1:nSpecies)) nPartOld = 0 !Initialize the lock for the non-analogue (NA) list of particles diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index e8e38f7..203047f 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -884,7 +884,7 @@ MODULE moduleMesh !Obtain the cross sections for the different processes !TODO: From here it might be a procedure in interactionMatrix vRel = NORM2(part_i%v-part_j%v) - rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m) + rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass) eRel = rMass*vRel**2 CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) @@ -1079,7 +1079,7 @@ MODULE moduleMesh !Compute changes in velocity for each particle deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v - mass_ij(p) = pair%sp_i%m*partTemp%part%weight + mass_ij(p) = pair%sp_i%mass*partTemp%part%weight p_ij(p,1:3) = mass_ij(p)*partTemp%part%v !Move to the next particle in the list @@ -1162,7 +1162,7 @@ MODULE moduleMesh !Compute changes in velocity for each particle deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v - mass_ji(p) = pair%sp_j%m*partTemp%part%weight + mass_ji(p) = pair%sp_j%mass*partTemp%part%weight p_ji(p,1:3) = mass_ji(p)*partTemp%part%v !Move to the next particle in the list diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index d130a8a..e7e4d65 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -125,7 +125,7 @@ MODULE moduleMeshBoundary SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj) TYPE IS(boundaryIonization) - mRel = reducedMass(bound%m0, part%species%m) + mRel = reducedMass(bound%m0, part%species%mass) vRel = SUM(DABS(part%v-bound%v0)) eRel = mRel*vRel**2*5.D-1 @@ -238,7 +238,7 @@ MODULE moduleMeshBoundary !Get relative velocity vRel = NORM2(part%v) !Convert to relative energy - eRel = part%species%m*vRel**2*5.D-1 + eRel = part%species%mass*vRel**2*5.D-1 !If energy is abound threshold calculate secondary electrons IF (eRel >= bound%energyThreshold) THEN diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index ccca930..f01c20b 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -164,8 +164,8 @@ MODULE moduleCollisions self%amount = amount - mass_i = species(i)%obj%m - mass_j = species(j)%obj%m + mass_i = species(i)%obj%mass + mass_j = species(j)%obj%mass ALLOCATE(self%collisions(1:self%amount)) @@ -227,8 +227,8 @@ MODULE moduleCollisions REAL(8):: m_i, m_j REAL(8), DIMENSION(1:3):: vCM, vp - m_i = part_i%species%m - m_j = part_j%species%m + m_i = part_i%species%mass + m_j = part_j%species%mass !Applies the collision vCM = velocityCM(part_i%weight*m_i, part_i%v, part_j%weight*m_j, part_j%v) vp = vRel*randomDirectionVHS() @@ -283,7 +283,7 @@ MODULE moduleCollisions END SELECT !momentum change per ionization process - collision%deltaV = sqrt(collision%eThreshold / collision%electron%m) + collision%deltaV = sqrt(collision%eThreshold / collision%electron%mass) END SELECT @@ -307,7 +307,7 @@ MODULE moduleCollisions REAL(8), DIMENSION(1:3):: vChange TYPE(particle), POINTER:: newElectron => NULL(), remainingNeutral => NULL() - rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m) + rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass) eRel = rMass*vRel**2 !Relative energy must be higher than threshold IF (eRel > self%eThreshold) THEN diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index 166b6b2..da03ad8 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -61,7 +61,7 @@ MODULE moduleCoulomb self%sp_i => species(i)%obj self%sp_j => species(j)%obj - self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m + self%one_plus_massRatio_ij = 1.D0 + self%sp_i%mass/self%sp_j%mass Z_i = 0.D0 Z_j = 0.D0 @@ -87,11 +87,11 @@ MODULE moduleCoulomb scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3) - self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent - self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent + self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%mass**2) * scaleFactor !Missing density because it's cell dependent + self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%mass**2) * scaleFactor !Missing density because it's cell dependent - self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent - self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent + self%l2_j = self%sp_j%mass / 2.D0 !Missing temperature because it's cell dependent + self%l2_i = self%sp_i%mass / 2.D0 !Missing temperature because it's cell dependent END SUBROUTINE initInteractionCoulomb diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 644cfa7..b5babfb 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -179,16 +179,21 @@ MODULE moduleInject IMPLICIT NONE INTEGER:: i + INTEGER, DIMENSION(1:nInject):: nMin, nMax !$OMP SINGLE nPartInj = 0 DO i = 1, nInject IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN + nMin(i) = nPartInj + 1 nPartInj = nPartInj + inject(i)%nParticles + nMax(i) = nPartInj END IF END DO + PRINT *, nMin + PRINT *, nMax IF (ALLOCATED(partInj)) DEALLOCATE(partInj) ALLOCATE(partInj(1:nPartInj)) @@ -196,7 +201,7 @@ MODULE moduleInject DO i=1, nInject IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - CALL inject(i)%addParticles() + CALL inject(i)%addParticles(nMin(i), nMax(i)) END IF END DO @@ -279,8 +284,8 @@ MODULE moduleInject IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self + INTEGER, INTENT(in):: nMin, nMax !Min and Max index in partInj array INTEGER:: randomX - INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array INTEGER:: i INTEGER:: n, sp CLASS(meshEdge), POINTER:: randomEdge @@ -288,17 +293,6 @@ MODULE moduleInject !Insert particles !$OMP SINGLE - nMin = 0 - DO i = 1, self%id -1 - IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - nMin = nMin + inject(i)%nParticles - - END IF - - END DO - - nMin = nMin + 1 - nMax = nMin + self%nParticles - 1 !Assign weight to particle. partInj(nMin:nMax)%weight = self%species%weight !Particle is considered to be outside the domain diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index ab08f08..e1b552b 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -4,12 +4,56 @@ MODULE moduleSpecies USE OMP_LIB IMPLICIT NONE + !Basic type that defines a macro-particle + TYPE:: particle + REAL(8):: r(1:3) !Position + REAL(8):: v(1:3) !Velocity + CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle + INTEGER:: cell !Index of element in which the particle is located. TODO: Make these pointers + INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh + REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p. + LOGICAL:: n_in !Flag that indicates if a particle is in the domain + REAL(8):: weight=0.D0 !weight of particle + + END TYPE particle + + !Wrapper to store the particles per species + TYPE:: particleArray + TYPE(particle), ALLOCATABLE, DIMENSION(:):: p + + END TYPE particleArray + + !Array of pointers for the species to be pushed + TYPE:: particleArray_pointer + TYPE(particle), POINTER, DIMENSION(:):: p + + END TYPE particleArray_pointer + + !Number of old particles + INTEGER, ALLOCATABLE, DIMENSION(:):: nPartOld + INTEGER:: nPartOldTotal + !Number of injected particles + INTEGER:: nPartInj + !Arrays that contain the particles + TYPE(particle), ALLOCATABLE, TARGET, DIMENSION(:):: partInj !array of inject particles + TYPE(particleArray), ALLOCATABLE, TARGET, DIMENSION(:):: partOld !array of particles from previous iteration + TYPE(particleArray_pointer), ALLOCATABLE, DIMENSION(:):: particlesToPush !particles pushed in each iteration + + INTEGER:: nSpeciesToPush + INTEGER, ALLOCATABLE, DIMENSION(:):: nPartOldToPush + + !Generic species type TYPE, ABSTRACT:: speciesGeneric - CHARACTER(:), ALLOCATABLE:: name - REAL(8):: m=0.D0, weight=0.D0, qm=0.D0 - INTEGER:: n=0 + INTEGER:: n=0 !Index of species + CHARACTER(:), ALLOCATABLE:: name !Name of species + !Mass, default weight of species and charge over mass + REAL(8):: mass=0.D0, weight=0.D0, qm=0.D0 + INTEGER:: every !How many interations between advancing the species + LOGICAL:: pushSpecies !Boolean to indicate if the species is moved in the iteration + END TYPE speciesGeneric + !Neutral species TYPE, EXTENDS(speciesGeneric):: speciesNeutral CLASS(speciesGeneric), POINTER:: ion => NULL() CONTAINS @@ -17,6 +61,7 @@ MODULE moduleSpecies END TYPE speciesNeutral + !Charged species TYPE, EXTENDS(speciesGeneric):: speciesCharged REAL(8):: q=0.D0 CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL() @@ -26,34 +71,17 @@ MODULE moduleSpecies END TYPE speciesCharged + !Wrapper for species TYPE:: speciesCont CLASS(speciesGeneric), ALLOCATABLE:: obj END TYPE + !Number of species INTEGER:: nSpecies + !Array for species TYPE(speciesCont), ALLOCATABLE, TARGET:: species(:) - TYPE particle - REAL(8):: r(1:3) !Position - REAL(8):: v(1:3) !Velocity - CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle - INTEGER:: cell !Index of element in which the particle is located - INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh - REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p. - LOGICAL:: n_in !Flag that indicates if a particle is in the domain - REAL(8):: weight=0.D0 !weight of particle - - END TYPE particle - - !Number of old particles - INTEGER:: nPartOld - !Number of injected particles - INTEGER:: nPartInj - !Arrays that contain the particles - TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration - TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles - CONTAINS FUNCTION speciesName2Index(speciesName) RESULT(sp) USE moduleErrors diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index 18fbb7f..2cb94f6 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -148,8 +148,8 @@ MODULE moduleOutput formatValues%density = rawValues%den*tempVol formatValues%velocity(:) = tempVel IF (tensorTrace(tensorTemp) > 0.D0) THEN - formatValues%pressure = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0 - formatValues%temperature = formatValues%pressure/(formatValues%density*kb) + formatValues%pressure = speciesIn%mass*tensorTrace(tensorTemp)*tempVol/3.D0 + formatValues%temperature = formatValues%pressure/(formatValues%density*kb) END IF END IF diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index de544ba..96e3f78 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -3,7 +3,7 @@ MODULE moduleSolver !Generic type for pusher of particles TYPE, PUBLIC:: pusherGeneric - PROCEDURE(push_interafece), POINTER, NOPASS:: pushParticle => NULL() + PROCEDURE(push_interface), POINTER, NOPASS:: pushParticle => NULL() !Boolean to indicate if the species is moved in the iteration LOGICAL:: pushSpecies !How many interations between advancing the species @@ -29,13 +29,13 @@ MODULE moduleSolver INTERFACE !Push a particle - PURE SUBROUTINE push_interafece(part, tauIn) + PURE SUBROUTINE push_interface(part, tauIn) USE moduleSpecies TYPE(particle), INTENT(inout):: part REAL(8), INTENT(in):: tauIn - END SUBROUTINE push_interafece + END SUBROUTINE push_interface !Solve the electromagnetic field SUBROUTINE solveEM_interface() @@ -169,24 +169,75 @@ MODULE moduleSolver USE moduleMesh IMPLICIT NONE - INTEGER:: n - INTEGER:: sp + INTEGER:: p + INTEGER:: s, sp + INTEGER:: e - !$OMP DO - DO n = 1, nPartOld - !Select species type - sp = partOld(n)%species%n + !$OMP SINGLE + !Allocate the array of particles to push + nSpeciesToPush = COUNT(solver%pusher(:)%pushSpecies) + ALLOCATE(particlesToPush(1:nSpeciesToPush)) + ALLOCATE(nPartOldToPush(1:nSpeciesToPush)) + !Point the arrays to the location of the particles + sp = 0 + DO s = 1, nSpecies !Checks if the species sp is update this iteration IF (solver%pusher(sp)%pushSpecies) THEN - !Push particle - CALL solver%pusher(sp)%pushParticle(partOld(n), tau(sp)) - !Find cell in wich particle reside - CALL solver%updateParticleCell(partOld(n)) + sp = sp + 1 + particlesToPush(sp)%p = partOld(s)%p + nPartOldToPush(sp) = nPartOld(s) END IF END DO - !$OMP END DO + + !Delete list of particles in cells for collisions if the particle is pushed + IF (listInCells) THEN + DO e = 1, mesh%numCells + DO s = 1, nSpecies + IF (solver%pusher(s)%pushSpecies) THEN + CALL mesh%cells(e)%obj%listPart_in(s)%erase() + mesh%cells(e)%obj%totalWeight(s) = 0.D0 + + END IF + + END DO + + END DO + + END IF + + !Erase the list of particles inside the cell in coll mesh if the particle is pushed + IF (doubleMesh) THEN + DO e = 1, meshColl%numCells + DO s = 1, nSpecies + IF (solver%pusher(s)%pushSpecies) THEN + CALL meshColl%cells(e)%obj%listPart_in(s)%erase() + meshColl%cells(e)%obj%totalWeight(s) = 0.D0 + + END IF + + END DO + + END DO + + END IF + !$OMP END SINGLE + + !Now, push particles + !$OMP DO + DO sp = 1, nSpeciesToPush + DO p = 1, nPartOldToPush(sp) + !Push particle + CALL solver%pusher(sp)%pushParticle(particlesToPush(sp)%p(p), tau(sp)) + !Find cell in which particle reside + CALL solver%updateParticleCell(particlesToPush(sp)%p(p)) + + END DO + + END DO + !$END OMP DO + END SUBROUTINE doPushes @@ -247,7 +298,10 @@ MODULE moduleSolver !$OMP SECTION nOldIn = 0 IF (ALLOCATED(partOld)) THEN - nOldIn = COUNT(partOld%n_in) + DO s = 1, nSpecies + nOldIn = nOldin + COUNT(partOld(s)%p(:)%n_in) + + END DO END IF !$OMP SECTION @@ -324,40 +378,6 @@ MODULE moduleSolver END DO - !$OMP SECTION - !Erase the list of particles inside the cell if particles have been pushed - IF (listInCells) THEN - DO s = 1, nSpecies - DO e = 1, mesh%numCells - IF (solver%pusher(s)%pushSpecies) THEN - CALL mesh%cells(e)%obj%listPart_in(s)%erase() - mesh%cells(e)%obj%totalWeight(s) = 0.D0 - - END IF - - END DO - - END DO - - END IF - - !$OMP SECTION - !Erase the list of particles inside the cell in coll mesh - IF (doubleMesh) THEN - DO s = 1, nSpecies - DO e = 1, meshColl%numCells - IF (solver%pusher(s)%pushSpecies) THEN - CALL meshColl%cells(e)%obj%listPart_in(s)%erase() - meshColl%cells(e)%obj%totalWeight(s) = 0.D0 - - END IF - - END DO - - END DO - - END IF - !$OMP END SECTIONS !$OMP SINGLE @@ -372,14 +392,17 @@ MODULE moduleSolver USE moduleMesh IMPLICIT NONE - INTEGER:: n + INTEGER:: n, sp CLASS(meshCell), POINTER:: cell !Loops over the particles to scatter them !$OMP DO - DO n = 1, nPartOld - cell => mesh%cells(partOld(n)%cell)%obj - CALL cell%scatter(cell%nNodes, partOld(n)) + DO sp = 1, nSpeciesToPush + DO n = 1, nPartOldToPush(sp) + cell => mesh%cells(particlesToPush(sp)%p(n)%cell)%obj + CALL cell%scatter(cell%nNodes, particlesToPush(sp)%p(n)) + + END DO END DO !$OMP END DO @@ -549,8 +572,8 @@ MODULE moduleSolver END IF - IF (nPartOld > 0) THEN - WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns" + IF (nPartOldTotal > 0) THEN + WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOldTotal), "ns" END IF WRITE(*,*) diff --git a/src/modules/solver/pusher/modulePusher.f90 b/src/modules/solver/pusher/modulePusher.f90 index f75c733..4c7e372 100644 --- a/src/modules/solver/pusher/modulePusher.f90 +++ b/src/modules/solver/pusher/modulePusher.f90 @@ -14,7 +14,7 @@ MODULE modulePusher END SUBROUTINE pushCartNeutral PURE SUBROUTINE pushCartElectrostatic(part, tauIn) - USE moduleSPecies + USE moduleSpecies USE moduleMesh IMPLICIT NONE From 0c4972384851536ca9e5bb5709b7ae1ed9bc6607 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Tue, 21 Nov 2023 09:23:18 +0100 Subject: [PATCH 10/10] NOT WORKING I need to fix something in another place and make some test, so this commit does not comile. --- src/modules/moduleSpecies.f90 | 34 ++++++++++++++--------------- src/modules/output/moduleOutput.f90 | 2 +- 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index e1b552b..d039247 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -8,10 +8,10 @@ MODULE moduleSpecies TYPE:: particle REAL(8):: r(1:3) !Position REAL(8):: v(1:3) !Velocity - CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle - INTEGER:: cell !Index of element in which the particle is located. TODO: Make these pointers + CLASS(speciesGeneric), POINTER:: species !Pointer to particle's species + INTEGER:: cell !Index of element in which the particle is located TODO: Make these pointers INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh - REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p. + REAL(8):: Xi(1:3) !Logical coordinates of particle in cell LOGICAL:: n_in !Flag that indicates if a particle is in the domain REAL(8):: weight=0.D0 !weight of particle @@ -23,24 +23,22 @@ MODULE moduleSpecies END TYPE particleArray - !Array of pointers for the species to be pushed - TYPE:: particleArray_pointer - TYPE(particle), POINTER, DIMENSION(:):: p + !Array of pointers for the particles to be pushed + TYPE:: particlePointer + TYPE(particle), POINTER:: p - END TYPE particleArray_pointer + END TYPE particlePointer - !Number of old particles - INTEGER, ALLOCATABLE, DIMENSION(:):: nPartOld - INTEGER:: nPartOldTotal - !Number of injected particles - INTEGER:: nPartInj !Arrays that contain the particles - TYPE(particle), ALLOCATABLE, TARGET, DIMENSION(:):: partInj !array of inject particles - TYPE(particleArray), ALLOCATABLE, TARGET, DIMENSION(:):: partOld !array of particles from previous iteration - TYPE(particleArray_pointer), ALLOCATABLE, DIMENSION(:):: particlesToPush !particles pushed in each iteration - - INTEGER:: nSpeciesToPush - INTEGER, ALLOCATABLE, DIMENSION(:):: nPartOldToPush + TYPE(particleArray), ALLOCATABLE, TARGET, DIMENSION(:):: particles !array of particles in the domain + TYPE(particle), ALLOCATABLE, TARGET, DIMENSION(:):: particlesInjection !array of inject particles + TYPE(particlePointer), ALLOCATABLE, DIMENSION(:):: particlesToPush !particles pushed in each iteration + !Integers to track number of particles in domain + INTEGER, ALLOCATABLE, DIMENSION(:):: nParticles + INTEGER:: nParticlesTotal + !Number of injected particles + INTEGER, ALLOCATABLE, DIMENSION(:):: nParticlesInject + INTEGER:: nPariclesInjectTotal !Generic species type TYPE, ABSTRACT:: speciesGeneric diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index 2cb94f6..16ed1b4 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -187,7 +187,7 @@ MODULE moduleOutput OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField + WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nParticlesTotal, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField CLOSE(20)