From 4ba08e74af5648c55d40f19779394cb88217f901 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 2 Jan 2021 12:50:22 +0100 Subject: [PATCH] Method to divide collisions from a collisional iteration into multiple pushing iterations. --- runs/Argon_Expansion/CX_case.json | 11 ++-- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 8 +-- src/modules/mesh/2DCyl/moduleMeshCyl.f90 | 22 +++---- src/modules/mesh/moduleMesh.f90 | 42 +++++++++---- src/modules/moduleCollisions.f90 | 72 +++++++++++++++------- src/modules/moduleConstParam.f90 | 4 +- src/modules/moduleSolver.f90 | 15 ++--- 7 files changed, 110 insertions(+), 64 deletions(-) diff --git a/runs/Argon_Expansion/CX_case.json b/runs/Argon_Expansion/CX_case.json index 94fc92f..54b2161 100644 --- a/runs/Argon_Expansion/CX_case.json +++ b/runs/Argon_Expansion/CX_case.json @@ -3,7 +3,7 @@ "path": "./runs/Argon_Expansion/", "triggerOutput": 10, "cpuTime": false, - "numColl": false + "numColl": true }, "geometry": { "type": "2DCyl", @@ -12,7 +12,7 @@ }, "species": [ {"name": "Argon", "type": "neutral", "mass": 6.633e-26, "weight": 1.0e8, "ion": "Argon+"}, - {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "weight": 1.0e8, "neutral": "Argon"} + {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "weight": 1.0e8, "charge": 1.0, "neutral": "Argon"} ], "boundary": [ {"name": "Injection", "physicalSurface": 1, "bTypes": [ @@ -29,9 +29,9 @@ ]} ], "inject": [ - {"name": "Exhausts Ar", "species": "Argon", "flow": 1.7, "units": "sccm", "v": 300.0, "T": [300.0, 300.0, 300.0], + {"name": "Exhausts Ar", "species": "Argon", "flow": 0.7, "units": "sccm", "v": 300.0, "T": [300.0, 300.0, 300.0], "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1}, - {"name": "Exhausts Ar+", "species": "Argon+", "flow": 1.3, "units": "sccm", "v": 40000.0, "T": [300.0, 300.0, 300.0], + {"name": "Exhausts Ar+", "species": "Argon+", "flow": 0.3, "units": "sccm", "v": 40000.0, "T": [300.0, 300.0, 300.0], "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1} ], "reference": { @@ -41,13 +41,14 @@ "radius": 1.88e-10 }, "case": { - "tau": [1.0e-4, 1.0e-6], + "tau": [1.0e-5, 1.0e-6], "time": 4.0e-3, "pusher": ["2DCylNeutral", "2DCylNeutral"], "WeightingScheme": "Volume" }, "interactions": { "folderCollisions": "./data/collisions/", + "tauCollisions": 1e-4, "collisions": [ {"species_i": "Argon", "species_j": "Argon", "cTypes": [ diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index eb78eb3..c05bba6 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -352,7 +352,7 @@ MODULE moduleMesh1DRad !Computes local stiffness matrix PURE FUNCTION elemKRad(self) RESULT(ke) - USE moduleConstParam, ONLY: PI + USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVol1DRadSegm), INTENT(in):: self @@ -371,12 +371,12 @@ MODULE moduleMesh1DRad ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) ke = 2.D0*ke*invJ - ke = ke*r*2.D0*PI + ke = ke*r*PI2 END FUNCTION elemKRad PURE FUNCTION elemFRad(self, source) RESULT(localF) - USE moduleConstParam, ONLY: PI + USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVol1DRadSegm), INTENT(in):: self @@ -393,7 +393,7 @@ MODULE moduleMesh1DRad r = DOT_PRODUCT(fPsi,self%r) ALLOCATE(localF(1:2)) localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ - localF = localF*r*2.D0*PI + localF = localF*r*PI2 END FUNCTION elemFRad diff --git a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 b/src/modules/mesh/2DCyl/moduleMeshCyl.f90 index 251573e..721361b 100644 --- a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMeshCyl.f90 @@ -369,7 +369,7 @@ MODULE moduleMeshCyl !Computes element area PURE SUBROUTINE areaQuad(self) - USE moduleConstParam + USE moduleConstParam, ONLY: PI8 IMPLICIT NONE CLASS(meshVolCylQuad), INTENT(inout):: self @@ -381,7 +381,7 @@ MODULE moduleMeshCyl self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral xi = 0.D0 - detJ = self%detJac(xi)*8.D0*PI !4*2*pi + detJ = self%detJac(xi)*PI8 !4*2*pi fPsi = self%fPsi(xi) r = DOT_PRODUCT(fPsi,self%r) self%volume = r*detJ @@ -466,7 +466,7 @@ MODULE moduleMeshCyl !Computes element local stiffness matrix PURE FUNCTION elemKQuad(self) RESULT(ke) - USE moduleConstParam, ONLY: PI + USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVolCylQuad), INTENT(in):: self @@ -493,13 +493,13 @@ MODULE moduleMeshCyl END DO END DO - ke = ke*2.D0*PI + ke = ke*PI2 END FUNCTION elemKQuad !Computes the local source vector for a force f PURE FUNCTION elemFQuad(self, source) RESULT(localF) - USE moduleConstParam + USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVolCylQuad), INTENT(in):: self @@ -525,7 +525,7 @@ MODULE moduleMeshCyl END DO END DO - localF = localF*2.D0*PI + localF = localF*PI2 END FUNCTION elemFQuad @@ -760,7 +760,7 @@ MODULE moduleMeshCyl !Calculates area for triangular element PURE SUBROUTINE areaTria(self) - USE moduleConstParam + USE moduleConstParam, ONLY: PI IMPLICIT NONE CLASS(meshVolCylTria), INTENT(inout):: self @@ -851,7 +851,7 @@ MODULE moduleMeshCyl !Computes element local stiffness matrix PURE FUNCTION elemKTria(self) RESULT(ke) - USE moduleConstParam + USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVolCylTria), INTENT(in):: self @@ -875,13 +875,13 @@ MODULE moduleMeshCyl ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ END DO - ke = ke*2.D0*PI + ke = ke*PI2 END FUNCTION elemKTria !Computes element local source vector PURE FUNCTION elemFTria(self, source) RESULT(localF) - USE moduleConstParam + USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVolCylTria), INTENT(in):: self @@ -906,7 +906,7 @@ MODULE moduleMeshCyl localF = localF + r*f*fPsi*wTria(l)*detJ END DO - localF = localF*2.D0*PI + localF = localF*PI2 END FUNCTION elemFTria diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index b0b2bd2..4446f89 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -379,7 +379,7 @@ MODULE moduleMesh END SUBROUTINE findCell !Computes collisions in element - SUBROUTINE collision(self) + SUBROUTINE collision(self, t) USE moduleCollisions USE moduleSpecies USE moduleList @@ -388,8 +388,12 @@ MODULE moduleMesh IMPLICIT NONE CLASS(meshVol), INTENT(inout):: self + INTEGER, INTENT(in):: t + INTEGER:: modCollisions !Remain of current iteration and everyCollisions + INTEGER:: iterToCollisions !Number of iterations from current to next collision INTEGER:: nPart !Number of particles inside the cell REAL(8):: pMax !Maximum probability of collision + INTEGER:: nCollIter !Number of collisions to be computed in this iteration INTEGER:: rnd !random index TYPE(particle), POINTER:: part_i, part_j INTEGER:: n !collision @@ -397,19 +401,36 @@ MODULE moduleMesh REAL(8):: sigmaVrelMaxNew TYPE(pointerArray), ALLOCATABLE:: partTemp(:) - self%nColl = 0 + modCollisions = MOD(t, everyCollisions) + iterToCollisions = everyCollisions - modCollisions + nPart = self%listPart_in%amount IF (nPart > 1) THEN - pMax = self%totalWeight*self%sigmaVrelMax*tauCollisions/self%volume - self%nColl = INT(REAL(nPart)*pMax*0.5D0) - - !Converts the list of particles to an array for easy access - IF (self%nColl > 0) THEN - partTemp = self%listPart_in%convert2Array() + IF (modCollisions == 0) THEN + !Collisional iteration, computes the number of iterations + pMax = self%totalWeight*self%sigmaVrelMax*tauCollisions/self%volume + self%nColl = INT(REAL(nPart)*pMax*0.5D0) END IF - DO n = 1, self%nColl + IF (self%nColl > iterToCollisions) THEN + nCollIter = self%nColl / iterToCollisions + + ELSE + nCollIter = self%nColl + + END IF + + IF (nCollIter > 0) THEN + !Converts the list of particles to an array for easy access + partTemp = self%listPart_in%convert2Array() + + !Removes collisions from this iteration form the total in the cell + self%nColl = self%nColl - nCollIter + + END IF + + DO n = 1, nCollIter !Select random numbers rnd = random(1, nPart) part_i => partTemp(rnd)%part @@ -429,11 +450,8 @@ MODULE moduleMesh END IF END DO - END IF - self%totalWeight = 0.D0 - END SUBROUTINE collision SUBROUTINE printOutputGmsh(self, t) diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 4ade202..e446bb7 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -4,7 +4,8 @@ MODULE moduleCollisions !Abstract type for collision between two particles TYPE, ABSTRACT:: collisionBinary - REAL(8):: rMass !reduced mass + REAL(8):: rMass !Reduced mass + REAL(8):: sMass !Summed mass TYPE(table1D):: crossSec !cross section of collision CONTAINS PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide @@ -33,9 +34,6 @@ MODULE moduleCollisions !Binary elastic interaction TYPE, EXTENDS(collisionBinary):: collisionBinaryElastic - !Weight distribution for Maxwellian function - REAL(8):: w_i = (1.D0+DSQRT(3.D0))/2.D0 - REAL(8):: w_j = (DSQRT(3.D0)-1.D0)/2.D0 CONTAINS PROCEDURE, PASS:: collide => collideBinaryElastic @@ -86,6 +84,35 @@ MODULE moduleCollisions INTEGER:: everyCollisions CONTAINS + !Velocity of center of mass of two particles + PURE FUNCTION velocityCM(m_i, v_i, m_j, v_j) RESULT(vCM) + IMPLICIT NONE + + REAL(8), INTENT(in):: m_i, m_j + REAL(8), INTENT(in), DIMENSION(1:3):: v_i, v_j + REAL(8):: vCM(1:3) + + vCM = (m_i*v_i + m_j*v_j)/(m_i + m_j) + + END FUNCTION velocityCM + + !Random direction for hard sphere collisions + FUNCTION randomDirectionVHS() RESULT(n) + USE moduleConstParam + USE moduleRandom + IMPLICIT NONE + + REAL(8):: n(1:3) + REAL(8):: cosXii, sinXii, eps + + cosXii = random(-1.D0, 1.D0) + sinXii = DSQRT(1.D0 - cosXii**2) + eps = random(0.D0, PI2) + + n = (/ cosXii, sinXii*DCOS(eps), sinXii*DSIN(eps) /) + + END FUNCTION randomDirectionVHS + !Inits the interaction matrix SUBROUTINE initInteractionMatrix(interactionMatrix) USE moduleSpecies @@ -146,7 +173,8 @@ MODULE moduleCollisions CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) !Calculates reduced mass - collision%rMass = (mass_i*mass_j)/(mass_i+mass_j) + collision%sMass = mass_i+mass_j + collision%rMass = (mass_i*mass_j)/collision%sMass END SUBROUTINE initBinaryElastic @@ -165,9 +193,9 @@ MODULE moduleCollisions REAL(8):: sigmaVrel REAL(8):: vRel !relative velocity REAL(8):: eRel !relative energy - REAL(8):: vp_i, vp_j, v_i, v_j !post and pre-collision velocities - REAL(8):: v_ij !sum of velocities modules - REAL(8):: alpha !random angle of scattering + REAL(8):: m_i, m_j + REAL(8), DIMENSION(1:3):: vCM + REAL(8):: vp(1:3) !eRel (in units of [m][L]^2[t]^-2 vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1 @@ -175,18 +203,15 @@ 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 !Applies the collision - v_i = NORM2(part_i%v) - v_j = NORM2(part_j%v) - v_ij = v_i+v_j - vp_j = v_ij*self%w_i - vp_i = v_ij*self%w_j - alpha = PI*random() - part_i%v(1) = vp_i*DCOS(alpha) - part_i%v(2) = vp_i*DSIN(alpha) - alpha = PI*random() - part_j%v(1) = vp_j*DCOS(alpha) - part_j%v(2) = vp_j*DSIN(alpha) + vCM = velocityCM(m_i, part_i%v, m_j, part_j%v) + vp = vRel*randomDirectionVHS() + + !Assign velocities to particles + part_i%v = vCM + m_j*vp/self%sMass + part_j%v = vCM - m_i*vp/self%sMass END IF @@ -218,7 +243,8 @@ MODULE moduleCollisions CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) !Calculates reduced mass - collision%rMass = (mass_i*mass_j)/(mass_i+mass_j) + collision%sMass = mass_i+mass_j + collision%rMass = (mass_i*mass_j)/collision%sMass !Specific parameters for ionization collision SELECT TYPE(collision) @@ -348,7 +374,8 @@ MODULE moduleCollisions CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) !Calculates reduced mass - collision%rMass = (mass_i*mass_j)/(mass_i+mass_j) + collision%sMass = mass_i+mass_j + collision%rMass = (mass_i*mass_j)/collision%sMass !Specific parameters for ionization collision SELECT TYPE(collision) @@ -453,7 +480,8 @@ MODULE moduleCollisions CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) !Calculates reduced mass - collision%rMass = (mass_i*mass_j)/(mass_i+mass_j) + collision%sMass = mass_i+mass_j + collision%rMass = (mass_i*mass_j)/collision%sMass END SUBROUTINE initBinaryChargeExchange diff --git a/src/modules/moduleConstParam.f90 b/src/modules/moduleConstParam.f90 index 634f845..cb12a23 100644 --- a/src/modules/moduleConstParam.f90 +++ b/src/modules/moduleConstParam.f90 @@ -4,7 +4,9 @@ MODULE moduleConstParam PUBLIC - REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi + REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi + REAL(8), PARAMETER:: PI2 = 2.D0*PI !2*pi + REAL(8), PARAMETER:: PI8 = 8.D0*PI !2*pi REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1 REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 801c8b3..3d39273 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -297,20 +297,16 @@ MODULE moduleSolver !Do the collisions in all the cells SUBROUTINE doCollisions(t) USE moduleMesh - USE moduleCollisions IMPLICIT NONE INTEGER, INTENT(in):: t INTEGER:: e - IF (MOD(t, everyCollisions) == 0) THEN - !$OMP DO SCHEDULE(DYNAMIC) - DO e=1, mesh%numVols - CALL mesh%vols(e)%obj%collision() - END DO - !$OMP END DO - - END IF + !$OMP DO SCHEDULE(DYNAMIC) + DO e=1, mesh%numVols + CALL mesh%vols(e)%obj%collision(t) + END DO + !$OMP END DO END SUBROUTINE doCollisions @@ -417,6 +413,7 @@ MODULE moduleSolver !$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() END DO