From ec128902ad44c1cf035446164a1634d461853318 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 28 Mar 2021 15:55:26 +0200 Subject: [PATCH] The integer part%sp that referenced the species index has been substituted for a pointer to the species. --- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 4 +-- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 4 +-- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 14 ++++---- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 14 ++++---- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 8 ++--- src/modules/mesh/moduleMesh.f90 | 4 +-- src/modules/mesh/moduleMeshBoundary.f90 | 16 ++++----- src/modules/moduleBoundary.f90 | 5 +-- src/modules/moduleCollisions.f90 | 36 ++++++++++---------- src/modules/moduleInject.f90 | 24 +++++++------ src/modules/moduleInput.f90 | 24 ++++++------- src/modules/moduleSolver.f90 | 4 +-- src/modules/moduleSpecies.f90 | 24 +++++++------ 13 files changed, 94 insertions(+), 87 deletions(-) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 020b971..46f5a00 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -395,12 +395,12 @@ MODULE moduleMesh1DCart w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%sp) + vertex => self%n1%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%sp) + vertex => self%n2%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 08609cd..12020c8 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -406,12 +406,12 @@ MODULE moduleMesh1DRad w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%sp) + vertex => self%n1%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%sp) + vertex => self%n2%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 10b12fe..be1f071 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -510,22 +510,22 @@ MODULE moduleMesh2DCart w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%sp) + vertex => self%n1%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%sp) + vertex => self%n2%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS - vertex => self%n3%output(part%sp) + vertex => self%n3%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(3) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS - vertex => self%n4%output(part%sp) + vertex => self%n4%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(4) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(4)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS @@ -869,17 +869,17 @@ MODULE moduleMesh2DCart w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%sp) + vertex => self%n1%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%sp) + vertex => self%n2%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS - vertex => self%n3%output(part%sp) + vertex => self%n3%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(3) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 888fce8..4318dff 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -531,22 +531,22 @@ MODULE moduleMesh2DCyl w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%sp) + vertex => self%n1%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%sp) + vertex => self%n2%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS - vertex => self%n3%output(part%sp) + vertex => self%n3%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(3) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS - vertex => self%n4%output(part%sp) + vertex => self%n4%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(4) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(4)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS @@ -898,17 +898,17 @@ MODULE moduleMesh2DCyl w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%sp) + vertex => self%n1%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%sp) + vertex => self%n2%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS - vertex => self%n3%output(part%sp) + vertex => self%n3%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(3) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index ec7f2b1..ce7f235 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -496,22 +496,22 @@ MODULE moduleMesh3DCart w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%sp) + vertex => self%n1%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%sp) + vertex => self%n2%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS - vertex => self%n3%output(part%sp) + vertex => self%n3%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(3) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS - vertex => self%n4%output(part%sp) + vertex => self%n4%output(part%species%n) vertex%den = vertex%den + part%weight*w_p(4) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(4)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 80b0e03..9aff9e5 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -375,7 +375,7 @@ MODULE moduleMesh CLASS IS (meshEdge) !Particle encountered a surface, apply boundary - CALL nextElement%fBoundary(part%sp)%apply(nextElement,part) + CALL nextElement%fBoundary(part%species%n)%apply(nextElement,part) !If particle is still inside the domain, call findCell IF (part%n_in) THEN @@ -439,7 +439,7 @@ MODULE moduleMesh part_i => partTemp(rnd)%part rnd = random(1, nPart) part_j => partTemp(rnd)%part - ij = interactionIndex(part_i%sp, part_j%sp) + ij = interactionIndex(part_i%species%n, part_j%species%n) sigmaVrelMaxNew = 0.D0 DO k = 1, interactionMatrix(ij)%amount CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 9a512b0..e89b8f3 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -91,7 +91,7 @@ MODULE moduleMeshBoundary INTEGER:: i !Modifies particle velocity according to wall temperature - SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) + SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj) TYPE IS(boundaryWallTemperature) DO i = 1, 3 part%v(i) = part%v(i) + bound%vTh*randomMaxwellian() @@ -123,9 +123,9 @@ MODULE moduleMeshBoundary TYPE(particle), POINTER:: newElectron TYPE(particle), POINTER:: newIon - SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) + SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj) TYPE IS(boundaryIonization) - mRel = (bound%m0*species(part%sp)%obj%m)*(bound%m0+species(part%sp)%obj%m) + mRel = (bound%m0*part%species%m)*(bound%m0+part%species%m) vRel = SUM(DABS(part%v-bound%v0)) eRel = mRel*vRel**2*5.D-1 @@ -133,7 +133,7 @@ MODULE moduleMeshBoundary ionizationRate = part%weight*bound%n0*bound%crossSection%get(eRel)*vRel !Rounds the number of particles up - ionizationPair = NINT(ionizationRate*bound%effectiveTime/species(bound%sp)%obj%weight) + ionizationPair = NINT(ionizationRate*bound%effectiveTime/bound%species%weight) !Create the new pair of particles DO p = 1, ionizationPair @@ -146,8 +146,8 @@ MODULE moduleMeshBoundary ALLOCATE(newElectron) ALLOCATE(newIon) - newElectron%sp = part%sp - newIon%sp = bound%sp + newElectron%species => part%species + newIon%species => bound%species newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0)) newIon%v = v0 @@ -162,13 +162,13 @@ MODULE moduleMeshBoundary newIon%xi = newElectron%xi newElectron%qm = part%qm - SELECT TYPE(spe => species(bound%sp)%obj) + SELECT TYPE(spe => bound%species) TYPE IS(speciesCharged) newIon%qm = spe%qm END SELECT - newElectron%weight = species(bound%sp)%obj%weight + newElectron%weight = bound%species%weight newIon%weight = newElectron%weight newElectron%n_in = .TRUE. diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 320bdad..2b61c67 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -1,5 +1,6 @@ MODULE moduleBoundary USE moduleTable + USE moduleSpecies !Generic type for boundaries TYPE, PUBLIC:: boundaryGeneric @@ -36,7 +37,7 @@ MODULE moduleBoundary !Ionization boundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals. - INTEGER:: sp !Ion species + CLASS(speciesGeneric), POINTER:: species !Ion species TYPE(table1D):: crossSection REAL(8):: effectiveTime REAL(8):: eThreshold @@ -123,7 +124,7 @@ MODULE moduleBoundary boundary%n0 = n0 * Vol_ref boundary%v0 = v0 / v_ref boundary%vTh = DSQRT(kb*T0/m0)/v_ref - boundary%sp = speciesID + boundary%species => species(speciesID)%obj boundary%effectiveTime = effTime / ti_ref CALL boundary%crossSection%init(crossSection) CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index a70bfb6..70bef90 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -199,8 +199,8 @@ MODULE moduleCollisions sigmaVrel = self%crossSec%get(eRel)*vRel sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN - m_i = species(part_i%sp)%obj%m - m_j = species(part_j%sp)%obj%m + m_i = part_i%species%m + m_j = part_j%species%m !Applies the collision vCM = velocityCM(m_i, part_i%v, m_j, part_j%v) vp = vRel*randomDirectionVHS() @@ -292,11 +292,11 @@ MODULE moduleCollisions sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN !Find which particle is the ionizing electron - IF (part_i%sp == self%electron%sp) THEN + IF (ASSOCIATED(part_i%species, self%electron)) THEN electron => part_i neutral => part_j - ELSEIF(part_j%sp == self%electron%sp) THEN + ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN electron => part_j neutral => part_i @@ -314,17 +314,17 @@ MODULE moduleCollisions !Creates a new electron from ionization ALLOCATE(newElectron) - newElectron%sp = electron%sp - newElectron%v = vp_n - newElectron%r = neutral%r - newElectron%xi = neutral%xi - newElectron%n_in = .TRUE. - newElectron%vol = neutral%vol - newElectron%weight = neutral%weight - newElectron%qm = electron%qm + newElectron%species => electron%species + newElectron%v = vp_n + newElectron%r = neutral%r + newElectron%xi = neutral%xi + newElectron%n_in = .TRUE. + newElectron%vol = neutral%vol + newElectron%weight = neutral%weight + newElectron%qm = electron%qm !Ionize neutral particle - SELECT TYPE(sp => species(neutral%sp)%obj) + SELECT TYPE(sp => neutral%species) TYPE IS(speciesNeutral) CALL sp%ionize(neutral) @@ -421,11 +421,11 @@ MODULE moduleCollisions sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN !Find which particle is the ionizing electron - IF (part_i%sp == self%electron%sp) THEN + IF (ASSOCIATED(part_i%species, self%electron)) THEN electron => part_i ion => part_j - ELSEIF(part_j%sp == self%electron%sp) THEN + ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN electron => part_j ion => part_i @@ -442,7 +442,7 @@ MODULE moduleCollisions electron%n_in = .FALSE. !Neutralize ion particle - SELECT TYPE(sp => species(ion%sp)%obj) + SELECT TYPE(sp => ion%species) TYPE IS(speciesCharged) CALL sp%neutralize(ion) @@ -501,7 +501,7 @@ MODULE moduleCollisions sigmaVrel = self%crossSec%get(eRel)*vRel sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN - SELECT TYPE(sp => species(part_i%sp)%obj) + SELECT TYPE(sp => part_i%species) TYPE IS (speciesNeutral) !Species i is neutral, ionize particle i CALL sp%ionize(part_i) @@ -512,7 +512,7 @@ MODULE moduleCollisions END SELECT - SELECT TYPE(sp => species(part_j%sp)%obj) + SELECT TYPE(sp => part_j%species) TYPE IS (speciesNeutral) !Species j is neutral, ionize particle j CALL sp%ionize(part_j) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 9a2659e..9660aa7 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -1,5 +1,6 @@ !injection of particles MODULE moduleInject + USE moduleSpecies !Generic type for velocity distribution function TYPE, ABSTRACT:: velDistGeneric @@ -51,7 +52,7 @@ MODULE moduleInject REAL(8):: T(1:3) !Temperature REAL(8):: n(1:3) !Direction of injection INTEGER:: nParticles !Number of particles to introduce each time step - INTEGER:: sp !Species of injection + CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges INTEGER, ALLOCATABLE:: edges(:) !Array with edges REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability @@ -86,11 +87,11 @@ MODULE moduleInject INTEGER:: e, et INTEGER:: phSurface(1:mesh%numEdges) - self%id = i - self%vMod = v/v_ref - self%n = n - self%T = T/T_ref - self%sp = sp + self%id = i + self%vMod = v/v_ref + self%n = n + self%T = T/T_ref + self%species => species(sp)%obj SELECT CASE(units) CASE ("sccm") !Standard cubic centimeter per minute @@ -224,14 +225,12 @@ MODULE moduleInject !$OMP SINGLE nMin = SUM(inject(1:(self%id-1))%nParticles) + 1 nMax = nMin + self%nParticles - 1 - !Assign particle type - partInj(nMin:nMax)%sp = self%sp !Assign weight to particle. - partInj(nMin:nMax)%weight = species(self%sp)%obj%weight + partInj(nMin:nMax)%weight = self%species%weight !Particle is considered to be outside the domain partInj(nMin:nMax)%n_in = .FALSE. !Assign charge/mass to charged particle. - SELECT TYPE(sp => species(self%sp)%obj) + SELECT TYPE(sp => self%species) TYPE IS(speciesCharged) partInj(nMin:nMax)%qm = sp%qm @@ -257,6 +256,9 @@ MODULE moduleInject END IF + !Assign particle type + partInj(n)%species => self%species + partInj(n)%v = (/ self%v(1)%obj%randomVel(), & self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /) @@ -264,7 +266,7 @@ MODULE moduleInject !Obtain natural coordinates of particle in cell partInj(n)%xi = mesh%vols(partInj(n)%vol)%obj%phy2log(partInj(n)%r) !Push new particle with the minimum time step - CALL solver%pusher(self%sp)%pushParticle(partInj(n), tauMin) + CALL solver%pusher(self%species%n)%pushParticle(partInj(n), tauMin) !Assign cell to new particle CALL solver%updateParticleCell(partInj(n)) diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 61c7117..99d031e 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -283,15 +283,15 @@ MODULE moduleInput !Allocate new particles DO p = 1, nNewPart ALLOCATE(partNew) - partNew%sp = sp - partNew%v(1) = velocity(1) + vTh*randomMaxwellian() - partNew%v(2) = velocity(2) + vTh*randomMaxwellian() - partNew%v(3) = velocity(3) + vTh*randomMaxwellian() - partNew%vol = e - partNew%r = mesh%vols(e)%obj%randPos() - partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r) - partNew%n_in = .TRUE. - partNew%weight = species(sp)%obj%weight + partNew%species => species(sp)%obj + partNew%v(1) = velocity(1) + vTh*randomMaxwellian() + partNew%v(2) = velocity(2) + vTh*randomMaxwellian() + partNew%v(3) = velocity(3) + vTh*randomMaxwellian() + partNew%vol = e + partNew%r = mesh%vols(e)%obj%randPos() + partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r) + partNew%n_in = .TRUE. + partNew%weight = species(sp)%obj%weight !If charged species, add qm to particle SELECT TYPE(sp => species(sp)%obj) TYPE IS (speciesCharged) @@ -442,8 +442,8 @@ MODULE moduleInput !Assign shared parameters for all species CALL config%get(object // '.name', species(i)%obj%name, found) CALL config%get(object // '.weight', species(i)%obj%weight, found) - species(i)%obj%sp = i - species(i)%obj%m = mass + species(i)%obj%n = i + species(i)%obj%m = mass END DO @@ -894,7 +894,7 @@ MODULE moduleInput REAL(8):: v, T, m !Reads species mass - m = species(inj%sp)%obj%m + m = inj%species%m !Reads distribution functions for velocity DO i = 1, 3 WRITE(istring, '(i2)') i diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 2927e21..b549104 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -154,7 +154,7 @@ MODULE moduleSolver !$OMP DO DO n=1, nPartOld !Select species type - sp = partOld(n)%sp + sp = partOld(n)%species%n !Checks if the species sp is update this iteration IF (solver%pusher(sp)%pushSpecies) THEN !Push particle @@ -722,7 +722,7 @@ MODULE moduleSolver part%weight = part%weight * fractionVolume - fractionWeight = part%weight / species(part%sp)%obj%weight + fractionWeight = part%weight / part%species%weight IF (fractionWeight >= 2.D0) THEN nSplit = FLOOR(fractionWeight) diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index f2904f9..7d89e3b 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -7,7 +7,7 @@ MODULE moduleSpecies TYPE, ABSTRACT:: speciesGeneric CHARACTER(:), ALLOCATABLE:: name REAL(8):: m=0.D0, weight=0.D0 - INTEGER:: sp=0 + INTEGER:: n=0 END TYPE speciesGeneric TYPE, EXTENDS(speciesGeneric):: speciesNeutral @@ -37,7 +37,7 @@ MODULE moduleSpecies TYPE particle REAL(8):: r(1:3) !Position REAL(8):: v(1:3) !Velocity - INTEGER:: sp !Particle species id + CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle INTEGER:: vol !Index of element in which the particle is located 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 @@ -48,6 +48,7 @@ MODULE moduleSpecies !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 @@ -60,15 +61,18 @@ MODULE moduleSpecies CHARACTER(:), ALLOCATABLE:: speciesName INTEGER:: sp - INTEGER:: n + INTEGER:: s sp = 0 - DO n = 1, nSpecies - IF (speciesName == species(n)%obj%name) THEN - sp = species(n)%obj%sp - EXIT + DO s = 1, nSpecies + IF (speciesName == species(s)%obj%name) THEN + sp = species(s)%obj%n + EXIT !If a species is found, exit the loop + END IF + END DO + !If no species is found, call a critical error IF (sp == 0) CALL criticalError('Species ' // speciesName // ' not found.', 'speciesName2Index') @@ -83,7 +87,7 @@ MODULE moduleSpecies TYPE(particle), INTENT(inout):: part IF (ASSOCIATED(self%ion)) THEN - part%sp = self%ion%sp + part%species => self%ion ELSE CALL criticalError('No ion defined for species' // self%name, 'ionizeNeutral') @@ -101,7 +105,7 @@ MODULE moduleSpecies TYPE(particle), INTENT(inout):: part IF (ASSOCIATED(self%ion)) THEN - part%sp = self%ion%sp + part%species => self%ion ELSE CALL criticalError('No ion defined for species' // self%name, 'ionizeCharged') @@ -119,7 +123,7 @@ MODULE moduleSpecies TYPE(particle), INTENT(inout):: part IF (ASSOCIATED(self%neutral)) THEN - part%sp = self%neutral%sp + part%species => self%neutral ELSE CALL criticalError('No neutral defined for species' // self%name, 'neutralizeCharged')