diff --git a/data/collisions/EL_e-Ar.dat b/data/collisions/EL_e-Ar.dat new file mode 100644 index 0000000..097729e --- /dev/null +++ b/data/collisions/EL_e-Ar.dat @@ -0,0 +1,3 @@ +# Relative energy (eV) cross section (m^2) +0.0 1e-20 +1000.0 1e-20 diff --git a/data/collisions/IO_e-Li.dat b/data/collisions/IO_e-Li.dat new file mode 100644 index 0000000..5e646f4 --- /dev/null +++ b/data/collisions/IO_e-Li.dat @@ -0,0 +1,148 @@ +# Rusudan I., et al., Atoms 9(4):90 2021 +# Relative energy (eV) cross section (m^2) +5.393 1.061E-23 +6.000 5.885E-21 +7.000 1.352E-20 +8.000 1.927E-20 +9.000 2.364E-20 +10.000 2.698E-20 +11.000 2.956E-20 +12.000 3.156E-20 +13.000 3.309E-20 +14.000 3.427E-20 +15.000 3.517E-20 +16.000 3.585E-20 +17.000 3.635E-20 +18.000 3.670E-20 +19.000 3.694E-20 +20.000 3.708E-20 +21.000 3.714E-20 +22.000 3.714E-20 +23.000 3.709E-20 +24.000 3.700E-20 +25.000 3.686E-20 +26.000 3.671E-20 +27.000 3.652E-20 +28.000 3.632E-20 +29.000 3.610E-20 +30.000 3.588E-20 +31.000 3.564E-20 +32.000 3.539E-20 +33.000 3.514E-20 +34.000 3.488E-20 +35.000 3.462E-20 +36.000 3.435E-20 +37.000 3.409E-20 +38.000 3.383E-20 +39.000 3.356E-20 +40.000 3.330E-20 +41.000 3.304E-20 +42.000 3.277E-20 +43.000 3.252E-20 +44.000 3.226E-20 +45.000 3.201E-20 +46.000 3.175E-20 +47.000 3.151E-20 +48.000 3.126E-20 +49.000 3.102E-20 +50.000 3.078E-20 +51.000 3.054E-20 +52.000 3.031E-20 +53.000 3.008E-20 +54.000 2.985E-20 +55.000 2.963E-20 +56.000 2.941E-20 +57.000 2.919E-20 +58.000 2.898E-20 +59.000 2.877E-20 +60.000 2.856E-20 +61.000 2.835E-20 +62.000 2.815E-20 +63.000 2.795E-20 +64.000 2.776E-20 +65.000 2.756E-20 +66.000 2.737E-20 +67.000 2.719E-20 +68.000 2.700E-20 +69.000 2.682E-20 +70.000 2.664E-20 +71.000 2.646E-20 +72.000 2.629E-20 +73.000 2.612E-20 +74.000 2.595E-20 +75.000 2.578E-20 +76.000 2.562E-20 +77.000 2.546E-20 +78.000 2.530E-20 +79.000 2.514E-20 +80.000 2.499E-20 +81.000 2.483E-20 +82.000 2.468E-20 +83.000 2.453E-20 +84.000 2.439E-20 +85.000 2.424E-20 +86.000 2.410E-20 +87.000 2.396E-20 +88.000 2.382E-20 +89.000 2.369E-20 +90.000 2.355E-20 +91.000 2.342E-20 +92.000 2.329E-20 +93.000 2.316E-20 +94.000 2.303E-20 +95.000 2.290E-20 +96.000 2.278E-20 +97.000 2.266E-20 +98.000 2.253E-20 +99.000 2.241E-20 +100.000 2.230E-20 +101.000 2.218E-20 +102.000 2.206E-20 +103.000 2.195E-20 +104.000 2.184E-20 +105.000 2.173E-20 +106.000 2.162E-20 +107.000 2.151E-20 +108.000 2.140E-20 +109.000 2.129E-20 +110.000 2.119E-20 +111.000 2.109E-20 +112.000 2.098E-20 +113.000 2.088E-20 +114.000 2.078E-20 +115.000 2.068E-20 +116.000 2.059E-20 +117.000 2.049E-20 +118.000 2.039E-20 +119.000 2.030E-20 +120.000 2.021E-20 +121.000 2.011E-20 +122.000 2.002E-20 +123.000 1.993E-20 +124.000 1.984E-20 +125.000 1.976E-20 +126.000 1.967E-20 +127.000 1.958E-20 +128.000 1.950E-20 +129.000 1.941E-20 +130.000 1.933E-20 +131.000 1.924E-20 +132.000 1.916E-20 +133.000 1.908E-20 +134.000 1.900E-20 +135.000 1.892E-20 +136.000 1.884E-20 +137.000 1.877E-20 +138.000 1.869E-20 +139.000 1.861E-20 +140.000 1.854E-20 +141.000 1.846E-20 +142.000 1.839E-20 +143.000 1.831E-20 +144.000 1.824E-20 +145.000 1.817E-20 +146.000 1.810E-20 +147.000 1.803E-20 +148.000 1.796E-20 +149.000 1.789E-20 +150.000 1.782E-20 diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index d908dd5..8e8072d 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -77,8 +77,6 @@ MODULE moduleMesh0D self%volume = 1.D0 self%n1%v = 1.D0 - self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) - CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 2ad9251..0a846d9 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -216,8 +216,6 @@ MODULE moduleMesh1DCart self%n1%v = self%n1%v + self%arNodes(1) self%n2%v = self%n2%v + self%arNodes(2) - self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) - CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index d8c0cb1..57619ee 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -218,8 +218,6 @@ MODULE moduleMesh1DRad self%n1%v = self%n1%v + self%arNodes(1) self%n2%v = self%n2%v + self%arNodes(2) - self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) - CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 8be8473..d01c0e4 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -307,8 +307,6 @@ MODULE moduleMesh2DCart self%n3%v = self%n3%v + self%arNodes(3) self%n4%v = self%n4%v + self%arNodes(4) - self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) - CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) @@ -637,8 +635,6 @@ MODULE moduleMesh2DCart self%n2%v = self%n2%v + self%arNodes(2) self%n3%v = self%n3%v + self%arNodes(3) - self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) - CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 02ff673..0a6b7f4 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -295,8 +295,6 @@ MODULE moduleMesh2DCyl self%n3%v = self%n3%v + self%arNodes(3) self%n4%v = self%n4%v + self%arNodes(4) - self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) - CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) @@ -658,8 +656,6 @@ MODULE moduleMesh2DCyl self%n2%v = self%n2%v + self%arNodes(2) self%n3%v = self%n3%v + self%arNodes(3) - self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) - CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 23acdd2..0add18a 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -281,8 +281,6 @@ MODULE moduleMesh3DCart self%n3%v = self%n3%v + volNodes(3) self%n4%v = self%n4%v + volNodes(4) - self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) - CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index 13e25ff..15045b1 100644 --- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 @@ -45,6 +45,7 @@ MODULE moduleMeshOutput0D CLASS(meshGeneric), INTENT(inout):: self INTEGER, INTENT(in):: t CHARACTER(:), ALLOCATABLE:: fileName + INTEGER:: k fileName='OUTPUT_Collisions.dat' IF (t == 0) THEN @@ -56,7 +57,7 @@ MODULE moduleMeshOutput0D END IF OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE(20, "(ES20.6E3, I20)") REAL(t)*tauMin*ti_ref, self%vols(1)%obj%nColl + WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%vols(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) CLOSE(20) END SUBROUTINE printColl0D diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index e63a73c..40eda4f 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -98,6 +98,7 @@ MODULE moduleMeshOutputGmsh2 CLASS(meshGeneric), INTENT(inout):: self INTEGER, INTENT(in):: t INTEGER:: numEdges + INTEGER:: k, c INTEGER:: n REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName @@ -125,19 +126,23 @@ MODULE moduleMeshOutputGmsh2 WRITE(60, "(A)") '$MeshFormat' WRITE(60, "(A)") '2.2 0 8' WRITE(60, "(A)") '$EndMeshFormat' - WRITE(60, "(A)") '$ElementData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Collisions"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numVols - DO n=1, self%numVols - WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%nColl + DO k = 1, nCollPairs + DO c = 1, interactionMatrix(k)%amount + WRITE(60, "(A)") '$ElementData' + WRITE(60, "(A)") '1' + WRITE(60, "(5A,I2)") '"Pair ', interactionMatrix(k)%sp_i%name, '-', interactionMatrix(k)%sp_j%name, ' collision ', c + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numVols + DO n=1, self%numVols + WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%tallyColl(k)%tally(c) + END DO + WRITE(60, "(A)") '$EndElementData' + END DO END DO - WRITE(60, "(A)") '$EndElementData' CLOSE(60) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index e7d98f3..142b78d 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -3,6 +3,7 @@ MODULE moduleMesh USE moduleList USE moduleOutput USE moduleBoundary + USE moduleCollisions IMPLICIT NONE !Generic mesh element @@ -148,15 +149,15 @@ MODULE moduleMesh !Parent of Volume element TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshVol !Maximum collision rate - REAL(8):: sigmaVrelMax = 0.D0 + REAL(8), ALLOCATABLE:: sigmaVrelMax(:) + !Arrays for counting number of collisions + TYPE(tallyCollisions), ALLOCATABLE:: tallyColl(:) !Volume REAL(8):: volume = 0.D0 !List of particles inside the volume 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), ALLOCATABLE:: totalWeight(:) CONTAINS @@ -658,15 +659,17 @@ MODULE moduleMesh USE moduleList use moduleRefParam USE moduleRandom + USE moduleOutput IMPLICIT NONE CLASS(meshGeneric), INTENT(inout), TARGET:: self INTEGER, INTENT(in):: t INTEGER:: e CLASS(meshVol), POINTER:: vol - INTEGER:: k, nPairs, i, j + INTEGER:: k, i, j INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell REAL(8):: pMax !Maximum probability of collision + INTEGER:: nColl TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:) TYPE(particle), POINTER:: part_i, part_j INTEGER:: n, c @@ -674,21 +677,23 @@ MODULE moduleMesh REAL(8):: sigmaVrelTotal REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:) REAL(8):: rnd !Random number for collision - INTEGER:: realCollisions IF (MOD(t, everyColl) == 0) THEN !Collisions need to be performed in this iteration - !$OMP DO SCHEDULE(DYNAMIC) + !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) DO e=1, self%numVols - realCollisions = 0 - vol => self%vols(e)%obj !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 + DO k = 1, nCollPairs + !Reset tally of collisions + IF (collOutput) THEN + vol%tallyColl(k)%tally = 0 + + END IF + IF (interactionMatrix(k)%amount > 0) THEN !Select the species for the collision pair i = interactionMatrix(k)%sp_i%n @@ -703,22 +708,22 @@ MODULE moduleMesh nPart = nPart_i + nPart_j !Resets the number of collisions in the cell - vol%nColl = 0 + nColl = 0 !Probability of collision for pair i-j - pMax = (vol%totalWeight(i) + vol%totalWeight(j))*vol%sigmaVrelMax*tauColl/vol%volume + pMax = (vol%totalWeight(i) + vol%totalWeight(j))*vol%sigmaVrelMax(k)*tauColl/vol%volume !Number of collisions in the cell - vol%nColl = NINT(REAL(nPart)*pMax*0.5D0) + nColl = NINT(REAL(nPart)*pMax*0.5D0) !Converts the list of particles to an array for easy access - IF (vol%nColl > 0) THEN + IF (nColl > 0) THEN partTemp_i = vol%listPart_in(i)%convert2Array() partTemp_j = vol%listPart_in(j)%convert2Array() END IF - DO n = 1, vol%nColl + DO n = 1, nColl !Select random particles part_i => NULL() part_j => NULL() @@ -735,13 +740,12 @@ MODULE moduleMesh !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? + !TODO: Try to find a way to no lose these collisions. Maybe check new 'k' and use that for the collision, maybe? IF (part_i%species%n /= i .OR. & part_j%species%n /= j) THEN CYCLE 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) @@ -749,26 +753,34 @@ MODULE moduleMesh CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) !Update maximum sigma*v_rel - IF (sigmaVrelTotal > vol%sigmaVrelMax) THEN - vol%sigmaVrelMax = sigmaVrelTotal + IF (sigmaVrelTotal > vol%sigmaVrelMax(k)) THEN + vol%sigmaVrelMax(k) = sigmaVrelTotal END IF ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount)) - probabilityColl(0) = 0.0 - probabilityColl(1:interactionMatrix(k)%amount) = sigmaVrel/vol%sigmaVrelMax + probabilityColl = 0.0 + DO c = 1, interactionMatrix(k)%amount + probabilityColl(c) = sigmaVrel(c)/vol%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1)) + + END DO !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 + IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax(k)) THEN !Loop over collisions DO c = 1, interactionMatrix(k)%amount - IF (SUM(probabilityColl(0:c-1)) + rnd <= probabilityColl(c)) THEN + IF (rnd <= probabilityColl(c)) THEN CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel) - realCollisions = realCollisions + 1 + + !If collisions are gonna be output, count the collision + IF (collOutput) THEN + vol%tallyColl(k)%tally(c) = vol%tallyColl(k)%tally(c) + 1 + + END IF !A collision has ocurred, exit the loop EXIT @@ -782,16 +794,17 @@ MODULE moduleMesh !Deallocate arrays for next collision DEALLOCATE(sigmaVrel, probabilityColl) + !End loop collisions in cell END DO END IF END IF + !End loop collision pairs END DO - vol%nColl = realCollisions - + !End loop volumes END DO !$OMP END DO diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 4f54e92..5f3bc26 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -159,13 +159,6 @@ MODULE moduleMeshBoundary newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r) newIon%xi = newElectron%xi - newElectron%qm = part%qm - SELECT TYPE(spe => bound%species) - TYPE IS(speciesCharged) - newIon%qm = spe%qm - - END SELECT - newElectron%weight = bound%species%weight newIon%weight = newElectron%weight diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 01571c6..47cace6 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -78,6 +78,14 @@ MODULE moduleCollisions END TYPE interactionsBinary + !Type to count number of collisions + TYPE:: tallyCollisions + INTEGER, ALLOCATABLE:: tally(:) + + END TYPE + + !Number of collision pairs (nSpecies*(nSpecies+1)/2) + INTEGER:: nCollPairs = 0 !Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory TYPE(interactionsBinary), ALLOCATABLE, TARGET:: interactionMatrix(:) !Folder for collision cross section tables @@ -120,10 +128,9 @@ MODULE moduleCollisions IMPLICIT NONE TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:) - INTEGER:: nInteractions - nInteractions = (nSpecies*(nSpecies+1))/2 - ALLOCATE(interactionMatrix(1:nInteractions)) + nCollPairs = (nSpecies*(nSpecies+1))/2 + ALLOCATE(interactionMatrix(1:nCollPairs)) interactionMatrix(:)%amount = 0 interactionMatrix(:)%rMass = 0.D0 @@ -222,8 +229,7 @@ MODULE moduleCollisions TYPE(particle), INTENT(inout), TARGET:: part_i, part_j REAL(8), INTENT(in):: vRel REAL(8):: m_i, m_j - REAL(8), DIMENSION(1:3):: vCM - REAL(8):: vp(1:3) + REAL(8), DIMENSION(1:3):: vCM, vp m_i = part_i%species%m m_j = part_j%species%m @@ -232,13 +238,8 @@ 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 @@ -273,16 +274,21 @@ MODULE moduleCollisions !Input energy is in eV. Convert to J with ev2J and then to !non-dimensional units. collision%eThreshold = energyThreshold*eV2J/(m_ref*v_ref**2) + !species for impacting electron electronIndex = speciesName2Index(electron) SELECT TYPE(sp => species(electronIndex)%obj) TYPE IS(speciesCharged) collision%electron => sp CLASS DEFAULT - CALL criticalError("Species " // sp%name // " chosen for ionization is not a charged species", 'initBinaryIonization') + CALL criticalError("Species " // sp%name // " chosen for " // & + "secondary electron is not a charged species", 'initBinaryIonization') END SELECT + !momentum change per ionization process + collision%deltaV = sqrt(collision%eThreshold / collision%electron%m) + END SELECT END SUBROUTINE initBinaryIonization @@ -302,8 +308,8 @@ MODULE moduleCollisions REAL(8), INTENT(in):: vRel REAL(8):: rMass, eRel TYPE(particle), POINTER:: electron => NULL(), neutral => NULL() + REAL(8), DIMENSION(1:3):: vChange TYPE(particle), POINTER:: newElectron - REAL(8), DIMENSION(1:3):: vp_e, vp_n rMass = reducedMass(part_i%species%m, part_j%species%m) eRel = rMass*vRel**2 @@ -322,35 +328,35 @@ MODULE moduleCollisions END IF - !Exchange energy between - vp_e = electron%v*(1.D0 - self%deltaV/NORM2(electron%v)) - vp_n = neutral%v* (1.D0 + self%deltaV/NORM2(neutral%v) ) - - !Changes velocity of impacting electron - electron%v = vp_e - - !Creates a new electron from ionization - ALLOCATE(newElectron) - 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%volColl = neutral%volColl - newElectron%weight = neutral%weight - newElectron%qm = electron%qm - !Ionize neutral particle SELECT TYPE(sp => neutral%species) TYPE IS(speciesNeutral) CALL sp%ionize(neutral) CLASS DEFAULT - CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization') + ! CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization') + RETURN END SELECT + !Exchange of velocity between particles + vChange = self%deltaV*randomDirectionVHS() + + !Energy is loss by the primary electron + electron%v = electron%v - vChange + + !Creates a new electron from ionization + ALLOCATE(newElectron) + newElectron%species => electron%species + !Secondary electorn gains energy from ionization + newElectron%v = vChange + newElectron%r = neutral%r + newElectron%xi = neutral%xi + newElectron%n_in = .TRUE. + newElectron%vol = neutral%vol + newElectron%volColl = neutral%volColl + newElectron%weight = neutral%weight + !Adds new electron to list of new particles from collisions CALL OMP_SET_LOCK(lockCollisions) CALL partCollisions%add(newElectron) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 9bfb35f..26f5ce1 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -82,6 +82,7 @@ MODULE moduleInject INTEGER, INTENT(in):: i REAL(8), INTENT(in):: v, n(1:3), T(1:3) INTEGER, INTENT(in):: sp, physicalSurface + REAL(8):: tauInject REAL(8), INTENT(in):: flow CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e, et @@ -93,18 +94,19 @@ MODULE moduleInject self%n = n / NORM2(n) self%T = T / T_ref self%species => species(sp)%obj + tauInject = tau(self%species%n) SELECT CASE(units) CASE ("sccm") !Standard cubic centimeter per minute - self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight) + self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight) CASE ("A") !Input current in Ampers - self%nParticles = INT(flow*tauMin*ti_ref/(qe*species(sp)%obj%weight)) + self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight)) CASE ("part/s") !Input current in Ampers - self%nParticles = INT(flow*tauMin*ti_ref/species(sp)%obj%weight) + self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight) CASE DEFAULT CALL criticalError("No support for units: " // units, 'initInject') @@ -161,8 +163,6 @@ MODULE moduleInject END DO self%sumWeight = self%cumWeight(self%nEdges) - nPartInj = nPartInj + self%nParticles - END SUBROUTINE initInject !Injection of particles @@ -174,12 +174,24 @@ MODULE moduleInject INTEGER:: i !$OMP SINGLE + nPartInj = 0 + DO i = 1, nInject + IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN + nPartInj = nPartInj + inject(i)%nParticles + + END IF + + END DO + IF (ALLOCATED(partInj)) DEALLOCATE(partInj) ALLOCATE(partInj(1:nPartInj)) !$OMP END SINGLE DO i=1, nInject - CALL inject(i)%addParticles() + IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN + CALL inject(i)%addParticles() + + END IF END DO END SUBROUTINE doInjects @@ -240,7 +252,7 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(in):: self INTEGER:: randomX INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array - INTEGER:: n + INTEGER:: n, sp CLASS(meshEdge), POINTER:: randomEdge !Insert particles @@ -251,12 +263,6 @@ MODULE moduleInject 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 => self%species) - TYPE IS(speciesCharged) - partInj(nMin:nMax)%qm = sp%qm - - END SELECT !$OMP END SINGLE !$OMP DO @@ -278,6 +284,7 @@ MODULE moduleInject END IF partInj(n)%volColl = randomEdge%eColl%n + sp = self%species%n !Assign particle type partInj(n)%species => self%species @@ -289,7 +296,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%species%n)%pushParticle(partInj(n), tauMin) + CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) !Assign cell to new particle CALL solver%updateParticleCell(partInj(n)) diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 1e907ef..db9ea74 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -132,6 +132,15 @@ MODULE moduleInput CALL config%get(object // '.temperature', T_ref, found) IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference') + !Derived parameters + L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length + v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity + ti_ref = L_ref/v_ref !reference time + Vol_ref = L_ref**3 !reference volume + EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field + Volt_ref = EF_ref*L_ref !reference voltage + B_ref = m_ref / (ti_ref * qe) !reference magnetic field + !If a reference cross section is given, it is used CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found) @@ -147,15 +156,6 @@ MODULE moduleInput END IF END IF - - !Derived parameters - L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length - v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity - ti_ref = L_ref/v_ref !reference time - Vol_ref = L_ref**3 !reference volume - EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field - Volt_ref = EF_ref*L_ref !reference voltage - B_ref = m_ref / (ti_ref * qe) !reference magnetic field END SUBROUTINE readReference @@ -417,15 +417,6 @@ MODULE moduleInput END IF 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) - partNew%qm = sp%qm - - CLASS DEFAULT - partNew%qm = 0.D0 - - END SELECT !Assign particle to temporal list of particles CALL partInitial%add(partNew) @@ -642,11 +633,13 @@ MODULE moduleInput CHARACTER(:), ALLOCATABLE:: crossSecFilePath CHARACTER(:), ALLOCATABLE:: cType LOGICAL:: found - INTEGER:: nInteractions, nCollisions + INTEGER:: nPairs, nCollisions INTEGER:: i, k, ij INTEGER:: pt_i, pt_j REAL(8):: energyThreshold, energyBinding CHARACTER(:), ALLOCATABLE:: electron + INTEGER:: e + CLASS(meshVol), POINTER:: vol !Firstly, checks if the object 'interactions' exists CALL config%info('interactions', found) @@ -681,8 +674,8 @@ MODULE moduleInput !Inits lock for list of particles CALL OMP_INIT_LOCK(lockCollisions) - CALL config%info('interactions.collisions', found, n_children = nInteractions) - DO i = 1, nInteractions + CALL config%info('interactions.collisions', found, n_children = nPairs) + DO i = 1, nPairs WRITE(iString, '(I2)') i object = 'interactions.collisions(' // TRIM(iString) // ')' CALL config%get(object // '.species_i', species_i, found) @@ -691,6 +684,7 @@ MODULE moduleInput pt_j = speciesName2Index(species_j) 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, pt_i, pt_j) @@ -732,7 +726,7 @@ MODULE moduleInput crossSecFilePath, energyBinding, electron) CASE DEFAULT - CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions') + CALL criticalError('Collision type' // cType // 'not defined', 'readInteractions') END SELECT @@ -740,6 +734,26 @@ MODULE moduleInput END DO + !Init the required arrays in each volume to account for MCC. + DO e = 1, meshForMCC%numVols + vol => meshForMCC%vols(e)%obj + + !Allocate Maximum cross section per collision pair and assign the initial collision rate + ALLOCATE(vol%sigmaVrelMax(1:nCollPairs)) + vol%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) + + IF (collOutput) THEN + ALLOCATE(vol%tallyColl(1:nCollPairs)) + DO k = 1, nCollPairs + ALLOCATE(vol%tallyColl(k)%tally(1:interactionmatrix(k)%amount)) + vol%tallyColl(k)%tally = 0 + + END DO + + END IF + + END DO + END IF END IF diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index f038ee7..017a321 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -165,7 +165,7 @@ MODULE moduleSolver INTEGER:: sp !$OMP DO - DO n=1, nPartOld + DO n = 1, nPartOld !Select species type sp = partOld(n)%species%n !Checks if the species sp is update this iteration @@ -204,7 +204,7 @@ MODULE moduleSolver REAL(8):: qmEFt(1:3) !Get the electric field at particle position - qmEFt = part%qm*gatherElecField(part)*tauIn + qmEFt = part%species%qm*gatherElecField(part)*tauIn !Update velocity part%v = part%v + qmEFt @@ -230,7 +230,7 @@ MODULE moduleSolver tauInHalf = tauIn *0.5D0 !Half of the force o f the electric field - qmEFt = part%qm*gatherElecField(part)*tauInHalf + qmEFt = part%species%qm*gatherElecField(part)*tauInHalf !Half step for electrostatic v_minus = part%v + qmEFt @@ -239,7 +239,7 @@ MODULE moduleSolver B = gatherMagnField(part) BNorm = NORM2(B) IF (BNorm > 0.D0) THEN - fn = DTAN(part%qm * tauInHalf*BNorm) / Bnorm + fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm 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) @@ -305,7 +305,7 @@ MODULE moduleSolver part_temp = part !Get electric field at particle position - qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn + qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauIn !z part_temp%v(1) = part%v(1) + qmEFt(1) part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn @@ -381,7 +381,7 @@ MODULE moduleSolver part_temp = part !Get electric field at particle position - qmEFt = part_temp%qm*gatherElecField(part_temp)*tauMin + qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauMin !r,theta v_p_oh_star(1) = part%v(1) + qmEFt(1) x_new = part%r(1) + v_p_oh_star(1)*tauIn @@ -586,7 +586,7 @@ MODULE moduleSolver !Loops over the particles to scatter them !$OMP DO - DO n=1, nPartOld + DO n = 1, nPartOld CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n)) END DO diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index fd8494f..d19ff28 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -6,7 +6,7 @@ MODULE moduleSpecies TYPE, ABSTRACT:: speciesGeneric CHARACTER(:), ALLOCATABLE:: name - REAL(8):: m=0.D0, weight=0.D0 + REAL(8):: m=0.D0, weight=0.D0, qm=0.D0 INTEGER:: n=0 END TYPE speciesGeneric @@ -18,7 +18,7 @@ MODULE moduleSpecies END TYPE speciesNeutral TYPE, EXTENDS(speciesGeneric):: speciesCharged - REAL(8):: q=0.D0, qm=0.D0 + REAL(8):: q=0.D0 CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL() CONTAINS PROCEDURE, PASS:: ionize => ionizeCharged @@ -43,7 +43,6 @@ MODULE moduleSpecies 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 - REAL(8):: qm = 0.D0 !charge over mass END TYPE particle