diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index 40ce55e..d908dd5 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -63,6 +63,7 @@ MODULE moduleMesh0D !Inits dummy 0D volume SUBROUTINE initVol0D(self, n, p, nodes) USE moduleRefParam + USE moduleSpecies IMPLICIT NONE CLASS(meshVol0D), INTENT(out):: self @@ -80,6 +81,9 @@ MODULE moduleMesh0D CALL OMP_INIT_LOCK(self%lock) + ALLOCATE(self%listPart_in(1:nSpecies)) + ALLOCATE(self%totalWeight(1:nSpecies)) + END SUBROUTINE initVol0D PURE FUNCTION getNodes0D(self) RESULT(n) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index f3f8131..2ad9251 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -220,6 +220,9 @@ MODULE moduleMesh1DCart CALL OMP_INIT_LOCK(self%lock) + ALLOCATE(self%listPart_in(1:nSpecies)) + ALLOCATE(self%totalWeight(1:nSpecies)) + END SUBROUTINE initVol1DCartSegm !Calculates a random position in 1D volume diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index c03e91d..d8c0cb1 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -222,6 +222,9 @@ MODULE moduleMesh1DRad CALL OMP_INIT_LOCK(self%lock) + ALLOCATE(self%listPart_in(1:nSpecies)) + ALLOCATE(self%totalWeight(1:nSpecies)) + END SUBROUTINE initVol1DRadSegm !Calculates a random position in 1D volume diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 5c8a937..8be8473 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -311,6 +311,9 @@ MODULE moduleMesh2DCart CALL OMP_INIT_LOCK(self%lock) + ALLOCATE(self%listPart_in(1:nSpecies)) + ALLOCATE(self%totalWeight(1:nSpecies)) + END SUBROUTINE initVolQuad2DCart !Computes element area @@ -638,6 +641,9 @@ MODULE moduleMesh2DCart CALL OMP_INIT_LOCK(self%lock) + ALLOCATE(self%listPart_in(1:nSpecies)) + ALLOCATE(self%totalWeight(1:nSpecies)) + END SUBROUTINE initVolTria2DCart !Random position in quadrilateral volume diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index d4192f1..02ff673 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -299,6 +299,9 @@ MODULE moduleMesh2DCyl CALL OMP_INIT_LOCK(self%lock) + ALLOCATE(self%listPart_in(1:nSpecies)) + ALLOCATE(self%totalWeight(1:nSpecies)) + END SUBROUTINE initVolQuad2DCyl !Computes element area @@ -659,6 +662,9 @@ MODULE moduleMesh2DCyl CALL OMP_INIT_LOCK(self%lock) + ALLOCATE(self%listPart_in(1:nSpecies)) + ALLOCATE(self%totalWeight(1:nSpecies)) + END SUBROUTINE initVolTria2DCyl !Random position in quadrilateral volume diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 24454d2..23acdd2 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -285,6 +285,9 @@ MODULE moduleMesh3DCart CALL OMP_INIT_LOCK(self%lock) + ALLOCATE(self%listPart_in(1:nSpecies)) + ALLOCATE(self%totalWeight(1:nSpecies)) + END SUBROUTINE initVolTetra3DCart !Random position in volume tetrahedron diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 8cc2497..e7d98f3 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -152,13 +152,13 @@ MODULE moduleMesh !Volume REAL(8):: volume = 0.D0 !List of particles inside the volume - TYPE(listNode):: listPart_in + TYPE(listNode), ALLOCATABLE:: listPart_in(:) !Lock indicator for listPart_in INTEGER(KIND=OMP_LOCK_KIND):: lock !Number of collisions per volume INTEGER:: nColl = 0 !Total weight of particles inside cell - REAL(8):: totalWeight = 0.D0 + REAL(8), ALLOCATABLE:: totalWeight(:) CONTAINS PROCEDURE(initVol_interface), DEFERRED, PASS:: init PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes @@ -516,6 +516,7 @@ MODULE moduleMesh CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell REAL(8):: xi(1:3) CLASS(meshElement), POINTER:: nextElement + INTEGER:: sp xi = self%phy2log(part%r) !Checks if particle is inside 'self' cell @@ -525,8 +526,9 @@ MODULE moduleMesh part%n_in = .TRUE. !Assign particle to listPart_in CALL OMP_SET_LOCK(self%lock) - CALL self%listPart_in%add(part) - self%totalWeight = self%totalWeight + part%weight + sp = part%species%n + CALL self%listPart_in(sp)%add(part) + self%totalWeight(sp) = self%totalWeight(sp) + part%weight CALL OMP_UNSET_LOCK(self%lock) ELSE @@ -586,6 +588,7 @@ MODULE moduleMesh CLASS(meshVol), POINTER:: vol REAL(8), DIMENSION(1:3):: xii CLASS(meshElement), POINTER:: nextElement + INTEGER:: sp found = .FALSE. @@ -595,8 +598,9 @@ MODULE moduleMesh IF (vol%inside(xii)) THEN part%volColl = vol%n CALL OMP_SET_LOCK(vol%lock) - CALL vol%listPart_in%add(part) - vol%totalWeight = vol%totalWeight + part%weight + sp = part%species%n + CALL vol%listPart_in(sp)%add(part) + vol%totalWeight(sp) = vol%totalWeight(sp) + part%weight CALL OMP_UNSET_LOCK(vol%lock) found = .TRUE. @@ -660,13 +664,12 @@ MODULE moduleMesh INTEGER, INTENT(in):: t INTEGER:: e CLASS(meshVol), POINTER:: vol - INTEGER:: nPart !Number of particles inside the cell + INTEGER:: k, nPairs, i, j + INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell REAL(8):: pMax !Maximum probability of collision - TYPE(pointerArray), ALLOCATABLE:: partTemp(:) - INTEGER:: i, j !random particle indexes + TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:) TYPE(particle), POINTER:: part_i, part_j - INTEGER:: n !collision - INTEGER:: ij, k + INTEGER:: n, c REAL(8):: vRel, eRel REAL(8):: sigmaVrelTotal REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:) @@ -681,80 +684,111 @@ MODULE moduleMesh realCollisions = 0 vol => self%vols(e)%obj - nPart = vol%listPart_in%amount - !Resets the number of collisions - vol%nColl = 0 + !TODO: Simplify this, to many sublevels + !Iterate over the number of pairs + nPairs = SIZE(interactionMatrix) !TODO: This does not change, make a variable in a module + DO k = 1, nPairs + IF (interactionMatrix(k)%amount > 0) THEN + !Select the species for the collision pair + i = interactionMatrix(k)%sp_i%n + j = interactionMatrix(k)%sp_j%n - !Calculates number of collisions if there is more than one particle in the cell - IF (nPart > 1) THEN - !Probability of collision - pMax = vol%totalWeight*vol%sigmaVrelMax*tauColl/vol%volume + !Number of particles per species in the collision pair + nPart_i = vol%listPart_in(i)%amount + nPart_j = vol%listPart_in(j)%amount - !Number of collisions in the cell - vol%nColl = NINT(REAL(nPart)*pMax*0.5D0) + IF (nPart_i > 0 .AND. nPart_j > 0) THEN + !Total number of particles for the collision pair + nPart = nPart_i + nPart_j - IF (vol%nColl > 0) THEN - !Converts the list of particles to an array for easy access - partTemp = vol%listPart_in%convert2Array() + !Resets the number of collisions in the cell + vol%nColl = 0 - END IF + !Probability of collision for pair i-j + pMax = (vol%totalWeight(i) + vol%totalWeight(j))*vol%sigmaVrelMax*tauColl/vol%volume - DO n = 1, vol%nColl - !Select random different particles - i = 0 - j = 0 - DO WHILE (i == j) - i = random(1, nPart) - j = random(1, nPart) - END DO - part_i => partTemp(i)%part - part_j => partTemp(j)%part - !TODO: I think that from here forward it can be passed to a procedure in interactionMatrix - ij = interactionIndex(part_i%species%n, part_j%species%n) + !Number of collisions in the cell + vol%nColl = NINT(REAL(nPart)*pMax*0.5D0) - IF (interactionMatrix(ij)%amount > 0) THEN - !Obtain the cross sections for the different processes - vRel = NORM2(part_i%v-part_j%v) - eRel = interactionMatrix(ij)%rMass*vRel**2 - CALL interactionMatrix(ij)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) - - !Update maximum sigma*v_rel - IF (sigmaVrelTotal > vol%sigmaVrelMax) THEN - vol%sigmaVrelMax = sigmaVrelTotal + !Converts the list of particles to an array for easy access + IF (vol%nColl > 0) THEN + partTemp_i = vol%listPart_in(i)%convert2Array() + partTemp_j = vol%listPart_in(j)%convert2Array() END IF - ALLOCATE(probabilityColl(0:interactionMatrix(ij)%amount)) - probabilityColl(0) = 0.0 - probabilityColl(1:interactionMatrix(ij)%amount) = sigmaVrel/vol%sigmaVrelMax + DO n = 1, vol%nColl + !Select random particles + part_i => NULL() + part_j => NULL() + rnd = random(1, nPart_i) + part_i => partTemp_i(rnd)%part + rnd = random(1, nPart_j) + part_j => partTemp_j(rnd)%part + !If they are the same particle, skip + !TODO: Maybe try to improve this + IF (ASSOCIATED(part_i, part_j)) THEN + CYCLE - !Selects random number between 0 and 1 - rnd = random() + END IF - !If the random number is below the total probability of collision, do collisions - IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax) THEN + !If particles do not belong to the species, skip collision + !This can happen, for example, if particle has been previously ionized or removed + !TODO: Try to find a way to no lose these collisions. Check new 'k' and use that for the collision, maybe? + IF (part_i%species%n /= i .OR. & + part_j%species%n /= j) THEN + CYCLE - !Loop over collisions - DO k = 1, interactionMatrix(ij)%amount - IF (SUM(probabilityColl(0:k-1)) + rnd <= probabilityColl(k)) THEN - CALL interactionMatrix(ij)%collisions(k)%obj%collide(part_i, part_j, vRel) - realCollisions = realCollisions + 1 + END IF - END IF + !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) + eRel = interactionMatrix(k)%rMass*vRel**2 + CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) - END DO + !Update maximum sigma*v_rel + IF (sigmaVrelTotal > vol%sigmaVrelMax) THEN + vol%sigmaVrelMax = sigmaVrelTotal - END IF + END IF + + ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount)) + probabilityColl(0) = 0.0 + probabilityColl(1:interactionMatrix(k)%amount) = sigmaVrel/vol%sigmaVrelMax + + !Selects random number between 0 and 1 + rnd = random() + + !If the random number is below the total probability of collision, collide particles + IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax) THEN + + !Loop over collisions + DO c = 1, interactionMatrix(k)%amount + IF (SUM(probabilityColl(0:c-1)) + rnd <= probabilityColl(c)) THEN + CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel) + realCollisions = realCollisions + 1 + + !A collision has ocurred, exit the loop + EXIT + + END IF + + END DO + + END IF - !Deallocate arrays for next collision - DEALLOCATE(sigmaVrel, probabilityColl) + !Deallocate arrays for next collision + DEALLOCATE(sigmaVrel, probabilityColl) + + END DO END IF - END DO + END IF - END IF + END DO vol%nColl = realCollisions diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index b651236..01571c6 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -67,6 +67,8 @@ MODULE moduleCollisions !Type for interaction matrix TYPE:: interactionsBinary + CLASS(speciesGeneric), POINTER:: sp_i + CLASS(speciesGeneric), POINTER:: sp_j INTEGER:: amount REAL(8):: rMass !Reduced mass TYPE(collisionCont), ALLOCATABLE:: collisions(:) @@ -142,16 +144,24 @@ MODULE moduleCollisions END FUNCTION interactionIndex !Inits the binary interaction - SUBROUTINE initInteractionBinary(self, amount, mass_i, mass_j) + SUBROUTINE initInteractionBinary(self, amount, i, j) USE moduleMath + USE moduleSpecies IMPLICIT NONE CLASS(interactionsBinary), INTENT(inout):: self INTEGER, INTENT(in):: amount - REAL(8), INTENT(in):: mass_i, mass_j + INTEGER, INTENT(in):: i, j + REAL(8):: mass_i, mass_j + + self%sp_i => species(i)%obj + self%sp_j => species(j)%obj self%amount = amount + mass_i = species(i)%obj%m + mass_j = species(j)%obj%m + self%rMass = reducedMass(mass_i, mass_j) ALLOCATE(self%collisions(1:self%amount)) @@ -165,14 +175,14 @@ MODULE moduleCollisions REAL(8), INTENT(in):: vRel, eRel REAL(8), INTENT(out):: sigmaVrelTotal REAL(8), INTENT(out), ALLOCATABLE:: sigmaVrel(:) - INTEGER:: k + INTEGER:: c sigmaVrelTotal = 0.D0 ALLOCATE(sigmaVrel(1:self%amount)) - DO k = 1, self%amount - sigmaVrel(k) = self%collisions(k)%obj%crossSec%get(eRel)*vRel + DO c = 1, self%amount + sigmaVrel(c) = self%collisions(c)%obj%crossSec%get(eRel)*vRel END DO sigmaVrelTotal = SUM(sigmaVrel) @@ -222,8 +232,13 @@ MODULE moduleCollisions vp = vRel*randomDirectionVHS() !Assign velocities to particles + PRINT *, part_i%v part_i%v = vCM + m_j*vp / (m_i + m_j) + PRINT *, part_i%v + PRINT *, part_j%v part_j%v = vCM - m_i*vp / (m_i + m_j) + PRINT *, part_j%v + PRINT * END SUBROUTINE collideBinaryElastic diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 5d90f62..1e907ef 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -333,6 +333,7 @@ MODULE moduleInput !Mean velocity and temperature at particle position REAL(8):: velocityXi(1:3), temperatureXi INTEGER:: nNewPart = 0.D0 + CLASS(meshVol), POINTER:: vol TYPE(particle), POINTER:: partNew REAL(8):: vTh TYPE(lNode), POINTER:: partCurr, partNext @@ -356,7 +357,6 @@ MODULE moduleInput !Density at centroid of cell nodes = mesh%vols(e)%obj%getNodes() nNodes = SIZE(nodes) - !TODO: Procedure to obtain centroid from element (also for printing Electric Field) fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/)) ALLOCATE(source(1:nNodes)) DO j = 1, nNodes @@ -430,6 +430,13 @@ MODULE moduleInput !Assign particle to temporal list of particles CALL partInitial%add(partNew) + !Assign particle to list in volume + vol => meshforMCC%vols(partNew%volColl)%obj + CALL OMP_SET_LOCK(vol%lock) + CALL vol%listPart_in(sp)%add(partNew) + vol%totalWeight(sp) = vol%totalWeight(sp) + partNew%weight + CALL OMP_UNSET_LOCK(vol%lock) + END DO DEALLOCATE(source) @@ -685,7 +692,7 @@ MODULE moduleInput CALL config%info(object // '.cTypes', found, n_children = nCollisions) ij = interactionIndex(pt_i,pt_j) !Allocates the required number of collisions per each pair of species ij - CALL interactionMatrix(ij)%init(nCollisions, species(pt_i)%obj%m, species(pt_j)%obj%m) + CALL interactionMatrix(ij)%init(nCollisions, pt_i, pt_j) DO k = 1, nCollisions WRITE (kString, '(I2)') k diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index acdd1de..cf27388 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -71,6 +71,7 @@ MODULE moduleList DO n=1, self%amount !Point element in array to element in list partArray(n)%part => tempNode%part + !Go to next element tempNode => tempNode%next diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 48d3dbc..f038ee7 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -457,6 +457,7 @@ MODULE moduleSolver INTEGER, SAVE:: nPartNew INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:) + INTEGER:: s !$OMP SECTIONS !$OMP SECTION @@ -540,18 +541,30 @@ MODULE moduleSolver END DO !$OMP SECTION - !Erase the list of particles inside the cell - DO e = 1, mesh%numVols - mesh%vols(e)%obj%totalWeight = 0.D0 - CALL mesh%vols(e)%obj%listPart_in%erase() + !Erase the list of particles inside the cell if particles have been pushed + DO s = 1, nSpecies + DO e = 1, mesh%numVols + IF (solver%pusher(s)%pushSpecies) THEN + CALL mesh%vols(e)%obj%listPart_in(s)%erase() + mesh%vols(e)%obj%totalWeight(s) = 0.D0 + + END IF + + END DO END DO !$OMP SECTION !Erase the list of particles inside the cell in coll mesh - DO e = 1, meshColl%numVols - meshColl%vols(e)%obj%totalWeight = 0.D0 - CALL meshColl%vols(e)%obj%listPart_in%erase() + DO s = 1, nSpecies + DO e = 1, meshColl%numVols + IF (solver%pusher(s)%pushSpecies) THEN + CALL meshColl%vols(e)%obj%listPart_in(s)%erase() + meshColl%vols(e)%obj%totalWeight(s) = 0.D0 + + END IF + + END DO END DO @@ -685,6 +698,7 @@ MODULE moduleSolver REAL(8):: newWeight TYPE(particle), POINTER:: newPart INTEGER:: p + INTEGER:: sp newWeight = part%weight / nSplit @@ -697,12 +711,14 @@ MODULE moduleSolver ALLOCATE(newPart) !Copy data from original particle newPart = part + !Add particle to list of new particles from weighting scheme CALL OMP_SET_LOCK(lockWScheme) CALL partWScheme%add(newPart) CALL OMP_UNSET_LOCK(lockWScheme) !Add particle to cell list CALL OMP_SET_lock(vol%lock) - CALL vol%listPart_in%add(newPart) + sp = part%species%n + CALL vol%listPart_in(sp)%add(newPart) CALL OMP_UNSET_lock(vol%lock) END DO diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index 39c1b20..fd8494f 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -91,7 +91,7 @@ MODULE moduleSpecies part%species => self%ion ELSE - CALL criticalError('No ion defined for species' // self%name, 'ionizeNeutral') + CALL criticalError('No ion defined for species :' // self%name, 'ionizeNeutral') END IF @@ -109,7 +109,7 @@ MODULE moduleSpecies part%species => self%ion ELSE - CALL criticalError('No ion defined for species' // self%name, 'ionizeCharged') + CALL criticalError('No ion defined for species :' // self%name, 'ionizeCharged') END IF