diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index 93fd4bd..6bae02f 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 fe755fa..9c8b039 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 6fc3c66..20c732e 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 4b23ff1..e0765da 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 1e4e580..a12aec8 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 d6a6391..b310457 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 9547304..a4b88cf 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,61 +664,83 @@ 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_i(:), partTemp_j(:) INTEGER:: rnd !random index TYPE(particle), POINTER:: part_i, part_j - INTEGER:: n !collision - INTEGER:: ij, k + INTEGER:: n, c REAL(8):: sigmaVrelMaxNew - TYPE(pointerArray), ALLOCATABLE:: partTemp(:) IF (MOD(t, everyColl) == 0) THEN !Collisions need to be performed in this iteration !$OMP DO SCHEDULE(DYNAMIC) DO e=1, self%numVols 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 numbers - rnd = random(1, nPart) - part_i => partTemp(rnd)%part - rnd = random(1, nPart) - part_j => partTemp(rnd)%part - 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(vol%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) + !Number of collisions in the cell + vol%nColl = NINT(REAL(nPart)*pMax*0.5D0) - END DO + !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() - !Update maximum cross section*v_rel per each collision - IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN - vol%sigmaVrelMax = sigmaVrelMaxNew + END IF + + !Do collisions + DO n = 1, vol%nColl + !Select random particles + rnd = random(1, nPart_i) + part_i => partTemp_i(rnd)%part + rnd = random(1, nPart_j) + part_j => partTemp_j(rnd)%part + + !Check all possible collisions in the pair + !TODO: Stop when collision occurs (and count it) + sigmaVrelMaxNew = 0.D0 + DO c = 1, interactionMatrix(k)%amount + CALL interactionMatrix(k)%collisions(c)%obj%collide(vol%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) + + END DO + + !Update maximum cross section*v_rel per each collision + IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN + vol%sigmaVrelMax = sigmaVrelMaxNew + + END IF + + END DO END IF - END DO + END IF - END IF + END DO END DO !$OMP END DO diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 2f22204..394b3cb 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -70,6 +70,8 @@ MODULE moduleCollisions !Type for interaction matrix TYPE:: interactionsBinary + CLASS(speciesGeneric), POINTER:: sp_i + CLASS(speciesGeneric), POINTER:: sp_j INTEGER:: amount TYPE(collisionCont), ALLOCATABLE:: collisions(:) CONTAINS @@ -139,11 +141,15 @@ MODULE moduleCollisions END FUNCTION interactionIndex !Inits the binary interaction - SUBROUTINE initInteractionBinary(self, amount) + SUBROUTINE initInteractionBinary(self, amount, i, j) IMPLICIT NONE CLASS(interactionsBinary), INTENT(inout):: self INTEGER, INTENT(in):: amount + INTEGER, INTENT(in):: i, j + + self%sp_i => species(i)%obj + self%sp_j => species(j)%obj self%amount = amount diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 31e2588..4ceca04 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -685,7 +685,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) + CALL interactionMatrix(ij)%init(nCollisions, pt_i, pt_j) DO k = 1, nCollisions WRITE (kString, '(I2)') k diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 4142549..b89c7e8 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -243,8 +243,6 @@ MODULE moduleSolver v_prime = v_minus + fn * crossProduct(v_minus, B) v_plus = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B) - PRINT *, v_minus, v_plus - END IF !Half step for electrostatic @@ -459,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 @@ -545,7 +544,10 @@ MODULE moduleSolver !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() + DO s = 1, nSpecies + CALL mesh%vols(e)%obj%listPart_in(s)%erase() + + END DO END DO @@ -553,7 +555,10 @@ MODULE moduleSolver !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 + CALL meshColl%vols(e)%obj%listPart_in(s)%erase() + + END DO END DO @@ -687,6 +692,7 @@ MODULE moduleSolver REAL(8):: newWeight TYPE(particle), POINTER:: newPart INTEGER:: p + INTEGER:: sp newWeight = part%weight / nSplit @@ -699,12 +705,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