From 286e858d66f6b18222af29faa2b17c357fe4b374 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 22 Sep 2023 18:18:46 +0200 Subject: [PATCH] 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