diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index 93fd4bd..40ce55e 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -76,7 +76,7 @@ MODULE moduleMesh0D self%volume = 1.D0 self%n1%v = 1.D0 - self%sigmaVrelMax = sigma_ref/L_ref**2 + self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index fe755fa..f3f8131 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -216,7 +216,7 @@ MODULE moduleMesh1DCart self%n1%v = self%n1%v + self%arNodes(1) self%n2%v = self%n2%v + self%arNodes(2) - self%sigmaVrelMax = sigma_ref/L_ref**2 + self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 6fc3c66..c03e91d 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -218,7 +218,7 @@ MODULE moduleMesh1DRad self%n1%v = self%n1%v + self%arNodes(1) self%n2%v = self%n2%v + self%arNodes(2) - self%sigmaVrelMax = sigma_ref/L_ref**2 + self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 4b23ff1..5c8a937 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -307,7 +307,7 @@ MODULE moduleMesh2DCart self%n3%v = self%n3%v + self%arNodes(3) self%n4%v = self%n4%v + self%arNodes(4) - self%sigmaVrelMax = sigma_ref/L_ref**2 + self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) CALL OMP_INIT_LOCK(self%lock) @@ -634,7 +634,7 @@ MODULE moduleMesh2DCart self%n2%v = self%n2%v + self%arNodes(2) self%n3%v = self%n3%v + self%arNodes(3) - self%sigmaVrelMax = sigma_ref/L_ref**2 + self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 1e4e580..d4192f1 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -295,7 +295,7 @@ MODULE moduleMesh2DCyl self%n3%v = self%n3%v + self%arNodes(3) self%n4%v = self%n4%v + self%arNodes(4) - self%sigmaVrelMax = sigma_ref/L_ref**2 + self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) CALL OMP_INIT_LOCK(self%lock) @@ -655,7 +655,7 @@ MODULE moduleMesh2DCyl self%n2%v = self%n2%v + self%arNodes(2) self%n3%v = self%n3%v + self%arNodes(3) - self%sigmaVrelMax = sigma_ref/L_ref**2 + self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index d6a6391..24454d2 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -281,7 +281,7 @@ MODULE moduleMesh3DCart self%n3%v = self%n3%v + volNodes(3) self%n4%v = self%n4%v + volNodes(4) - self%sigmaVrelMax = sigma_ref/L_ref**2 + self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 9547304..8cc2497 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -662,17 +662,24 @@ MODULE moduleMesh CLASS(meshVol), POINTER:: vol INTEGER:: nPart !Number of particles inside the cell REAL(8):: pMax !Maximum probability of collision - INTEGER:: rnd !random index + TYPE(pointerArray), ALLOCATABLE:: partTemp(:) + INTEGER:: i, j !random particle indexes TYPE(particle), POINTER:: part_i, part_j INTEGER:: n !collision INTEGER:: ij, k - REAL(8):: sigmaVrelMaxNew - TYPE(pointerArray), ALLOCATABLE:: partTemp(:) + REAL(8):: vRel, eRel + 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) DO e=1, self%numVols + + realCollisions = 0 + vol => self%vols(e)%obj nPart = vol%listPart_in%amount @@ -694,21 +701,54 @@ MODULE moduleMesh END IF 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) - + !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) - !Update maximum cross section*v_rel per each collision - IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN - vol%sigmaVrelMax = sigmaVrelMaxNew + 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 + + END IF + + ALLOCATE(probabilityColl(0:interactionMatrix(ij)%amount)) + probabilityColl(0) = 0.0 + probabilityColl(1:interactionMatrix(ij)%amount) = sigmaVrel/vol%sigmaVrelMax + + !Selects random number between 0 and 1 + rnd = random() + + !If the random number is below the total probability of collision, do collisions + IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax) THEN + + !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 DO + + END IF + + !Deallocate arrays for next collision + DEALLOCATE(sigmaVrel, probabilityColl) END IF @@ -716,6 +756,8 @@ MODULE moduleMesh END IF + vol%nColl = realCollisions + END DO !$OMP END DO diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 2f22204..b651236 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -7,8 +7,6 @@ MODULE moduleCollisions !Abstract type for collision between two particles TYPE, ABSTRACT:: collisionBinary - REAL(8):: rMass !Reduced mass - REAL(8):: sMassInv !Summed mass TYPE(table1D):: crossSec !cross section of collision CONTAINS PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide @@ -16,14 +14,13 @@ MODULE moduleCollisions END TYPE collisionBinary ABSTRACT INTERFACE - SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrelMaxNew, part_i, part_j) + SUBROUTINE collideBinary_interface(self, part_i, part_j, vRel) USE moduleSpecies IMPORT:: collisionBinary CLASS(collisionBinary), INTENT(in):: self - REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(inout):: sigmaVrelMaxNew TYPE(particle), INTENT(inout), TARGET:: part_i, part_j + REAL(8), INTENT(in):: vRel END SUBROUTINE @@ -71,9 +68,11 @@ MODULE moduleCollisions !Type for interaction matrix TYPE:: interactionsBinary INTEGER:: amount + REAL(8):: rMass !Reduced mass TYPE(collisionCont), ALLOCATABLE:: collisions(:) CONTAINS PROCEDURE, PASS:: init => initInteractionBinary + PROCEDURE, PASS:: getSigmaVrel => getSigmaVrelBinary END TYPE interactionsBinary @@ -85,13 +84,14 @@ MODULE moduleCollisions 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) + vCM = (m_i*v_i + m_j*v_j) / (m_i + m_j) END FUNCTION velocityCM @@ -123,6 +123,9 @@ MODULE moduleCollisions nInteractions = (nSpecies*(nSpecies+1))/2 ALLOCATE(interactionMatrix(1:nInteractions)) + interactionMatrix(:)%amount = 0 + interactionMatrix(:)%rMass = 0.D0 + END SUBROUTINE initInteractionMatrix !Gets the interaction index from the collision matrix from index i,j @@ -139,21 +142,46 @@ MODULE moduleCollisions END FUNCTION interactionIndex !Inits the binary interaction - SUBROUTINE initInteractionBinary(self, amount) + SUBROUTINE initInteractionBinary(self, amount, mass_i, mass_j) + USE moduleMath IMPLICIT NONE CLASS(interactionsBinary), INTENT(inout):: self INTEGER, INTENT(in):: amount + REAL(8), INTENT(in):: mass_i, mass_j self%amount = amount + self%rMass = reducedMass(mass_i, mass_j) + ALLOCATE(self%collisions(1:self%amount)) END SUBROUTINE initInteractionBinary + SUBROUTINE getSigmaVrelBinary (self, vRel, eRel, sigmaVrelTotal, sigmaVrel) + IMPLICIT NONE + + CLASS(interactionsBinary), INTENT(in):: self + REAL(8), INTENT(in):: vRel, eRel + REAL(8), INTENT(out):: sigmaVrelTotal + REAL(8), INTENT(out), ALLOCATABLE:: sigmaVrel(:) + INTEGER:: k + + sigmaVrelTotal = 0.D0 + + ALLOCATE(sigmaVrel(1:self%amount)) + + DO k = 1, self%amount + sigmaVrel(k) = self%collisions(k)%obj%crossSec%get(eRel)*vRel + + END DO + sigmaVrelTotal = SUM(sigmaVrel) + + END SUBROUTINE getSigmaVrelBinary + !ELASTIC COLLISIONS !Inits binary elastic collision - SUBROUTINE initBinaryElastic(collision, crossSectionFilename, mass_i, mass_j) + SUBROUTINE initBinaryElastic(collision, crossSectionFilename) USE moduleTable USE moduleRefParam USE moduleConstParam @@ -161,7 +189,6 @@ MODULE moduleCollisions CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename - REAL(8), INTENT(in):: mass_i, mass_j ALLOCATE(collisionBinaryElastic:: collision) @@ -171,15 +198,10 @@ MODULE moduleCollisions !Convert to no-dimensional units CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) - !Calculates reduced mass - collision%sMassInv = 1.D0/(mass_i+mass_j) - collision%rMass = (mass_i*mass_j)*collision%sMassInv - END SUBROUTINE initBinaryElastic !Binary elastic process - SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrelMaxNew, & - part_i, part_j) + SUBROUTINE collideBinaryElastic(self, part_i, part_j, vRel) USE moduleSpecies USE moduleConstParam USE moduleRandom @@ -187,38 +209,27 @@ MODULE moduleCollisions IMPLICIT NONE CLASS(collisionBinaryElastic), INTENT(in):: self - REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(inout):: sigmaVrelMaxNew TYPE(particle), INTENT(inout), TARGET:: part_i, part_j - REAL(8):: sigmaVrel - REAL(8):: vRel !relative velocity - REAL(8):: eRel !relative energy + REAL(8), INTENT(in):: vRel REAL(8):: m_i, m_j REAL(8), DIMENSION(1:3):: vCM REAL(8):: vp(1:3) - vRel = NORM2(part_i%v-part_j%v) - eRel = self%rMass*vRel**2 - sigmaVrel = self%crossSec%get(eRel)*vRel - sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel - IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN - 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() + 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() - !Assign velocities to particles - part_i%v = vCM + m_j*vp*self%sMassInv - part_j%v = vCM - m_i*vp*self%sMassInv - - END IF + !Assign velocities to particles + part_i%v = vCM + m_j*vp / (m_i + m_j) + part_j%v = vCM - m_i*vp / (m_i + m_j) END SUBROUTINE collideBinaryElastic !ELECTRON IMPACT IONIZATION !Inits electron impact ionization - SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, mass_i, mass_j, electron) + SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron) USE moduleTable USE moduleRefParam USE moduleConstParam @@ -229,7 +240,6 @@ MODULE moduleCollisions CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename REAL(8), INTENT(in):: energyThreshold - REAL(8), INTENT(in):: mass_i, mass_j CHARACTER(:), ALLOCATABLE, INTENT(in):: electron INTEGER:: electronIndex @@ -241,10 +251,6 @@ MODULE moduleCollisions !Convert to no-dimensional units CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) - !Calculates reduced mass - collision%sMassInv = 1.D0/(mass_i+mass_j) - collision%rMass = (mass_i*mass_j)*collision%sMassInv - !Specific parameters for ionization collision SELECT TYPE(collision) TYPE IS(collisionBinaryIonization) @@ -267,8 +273,7 @@ MODULE moduleCollisions END SUBROUTINE initBinaryIonization !Binary electron impact ionization process - SUBROUTINE collideBinaryIonization(self, sigmaVrelMax, sigmaVrelMaxNew, & - part_i, part_j) + SUBROUTINE collideBinaryIonization(self, part_i, part_j, vRel) USE moduleSpecies USE moduleErrors USE moduleList @@ -278,80 +283,71 @@ MODULE moduleCollisions IMPLICIT NONE CLASS(collisionBinaryIonization), INTENT(in):: self - REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(inout):: sigmaVrelMaxNew TYPE(particle), INTENT(inout), TARGET:: part_i, part_j + REAL(8), INTENT(in):: vRel + REAL(8):: rMass, eRel TYPE(particle), POINTER:: electron => NULL(), neutral => NULL() TYPE(particle), POINTER:: newElectron - REAL(8):: vRel, eRel - REAL(8):: sigmaVrel REAL(8), DIMENSION(1:3):: vp_e, vp_n - !eRel (in units of [m][L]^2[t]^-2 - vRel = NORM2(part_i%v-part_j%v) - eRel = self%rMass*vRel**2 + rMass = reducedMass(part_i%species%m, part_j%species%m) + eRel = rMass*vRel**2 !Relative energy must be higher than threshold IF (eRel > self%eThreshold) THEN - sigmaVrel = self%crossSec%get(eRel)*vRel - sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel - IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN - !Find which particle is the ionizing electron - IF (ASSOCIATED(part_i%species, self%electron)) THEN - electron => part_i - neutral => part_j + IF (ASSOCIATED(part_i%species, self%electron)) THEN + electron => part_i + neutral => part_j - ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN - electron => part_j - neutral => part_i + ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN + electron => part_j + neutral => part_i - ELSE - CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization') - - 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') - - END SELECT - - !Adds new electron to list of new particles from collisions - CALL OMP_SET_LOCK(lockCollisions) - CALL partCollisions%add(newElectron) - CALL OMP_UNSET_LOCK(lockCollisions) + ELSE + CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization') 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') + + END SELECT + + !Adds new electron to list of new particles from collisions + CALL OMP_SET_LOCK(lockCollisions) + CALL partCollisions%add(newElectron) + CALL OMP_UNSET_LOCK(lockCollisions) + END IF END SUBROUTINE collideBinaryIonization !ELECTRON ION RESONANT RECOMBINATION !Inits electron ion recombination - SUBROUTINE initBinaryRecombination(collision, crossSectionFilename, energyBinding, mass_i, mass_j, electron) + SUBROUTINE initBinaryRecombination(collision, crossSectionFilename, energyBinding, electron) USE moduleTable USE moduleRefParam USE moduleConstParam @@ -362,7 +358,6 @@ MODULE moduleCollisions CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename REAL(8), INTENT(in):: energyBinding - REAL(8), INTENT(in):: mass_i, mass_j CHARACTER(:), ALLOCATABLE, INTENT(in):: electron INTEGER:: electronIndex @@ -374,10 +369,6 @@ MODULE moduleCollisions !Convert to no-dimensional units CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) - !Calculates reduced mass - collision%sMassInv = 1.D0/(mass_i+mass_j) - collision%rMass = (mass_i*mass_j)*collision%sMassInv - !Specific parameters for ionization collision SELECT TYPE(collision) TYPE IS(collisionBinaryRecombination) @@ -399,9 +390,8 @@ MODULE moduleCollisions END SUBROUTINE initBinaryRecombination - !Binary electron impact ionization process - SUBROUTINE collideBinaryRecombination(self, sigmaVrelMax, sigmaVrelMaxNew, & - part_i, part_j) + !Binary recombination + SUBROUTINE collideBinaryRecombination(self, part_i, part_j, vRel) USE moduleSpecies USE moduleErrors USE moduleList @@ -410,59 +400,47 @@ MODULE moduleCollisions IMPLICIT NONE CLASS(collisionBinaryRecombination), INTENT(in):: self - REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(inout):: sigmaVrelMaxNew + REAL(8), INTENT(in):: vRel TYPE(particle), INTENT(inout), TARGET:: part_i, part_j TYPE(particle), POINTER:: electron => NULL(), ion => NULL() - REAL(8):: vRel, eRel REAL(8):: sigmaVrel REAL(8), DIMENSION(1:3):: vp_i - !eRel (in units of [m][L]^2[t]^-2 - vRel = NORM2(part_i%v-part_j%v) - eRel = self%rMass*vRel**2 - !Relative energy must be higher than threshold - sigmaVrel = self%crossSec%get(eRel)*vRel - sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel - IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN - !Find which particle is the ionizing electron - IF (ASSOCIATED(part_i%species, self%electron)) THEN - electron => part_i - ion => part_j + IF (ASSOCIATED(part_i%species, self%electron)) THEN + electron => part_i + ion => part_j - ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN - electron => part_j - ion => part_i + ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN + electron => part_j + ion => part_i - ELSE - CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization') - - END IF - - !Excess energy - !TODO: This energy should be transformed into photons - vp_i = ion%v* (1.D0 - (vRel + self%deltaV)/NORM2(ion%v)) - - !Remove electron from simulation - electron%n_in = .FALSE. - - !Neutralize ion particle - SELECT TYPE(sp => ion%species) - TYPE IS(speciesCharged) - CALL sp%neutralize(ion) - - CLASS DEFAULT - CALL criticalError(sp%name // " is not a charge", 'collideBinaryRecombination') - - END SELECT + ELSE + CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization') END IF + !Excess energy + !TODO: This energy should be transformed into photons + vp_i = ion%v* (1.D0 - (vRel + self%deltaV)/NORM2(ion%v)) + + !Remove electron from simulation + electron%n_in = .FALSE. + + !Neutralize ion particle + SELECT TYPE(sp => ion%species) + TYPE IS(speciesCharged) + CALL sp%neutralize(ion) + + CLASS DEFAULT + CALL criticalError(sp%name // " is not a charge", 'collideBinaryRecombination') + + END SELECT + END SUBROUTINE collideBinaryRecombination !RESONANT CHARGE EXCHANGE !Inits resonant charge exchange - SUBROUTINE initBinaryChargeExchange(collision, crossSectionFilename, mass_i, mass_j) + SUBROUTINE initBinaryChargeExchange(collision, crossSectionFilename) USE moduleTable USE moduleRefParam USE moduleConstParam @@ -470,7 +448,6 @@ MODULE moduleCollisions CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename - REAL(8), INTENT(in):: mass_i, mass_j ALLOCATE(collisionBinaryChargeExchange:: collision) @@ -480,56 +457,39 @@ MODULE moduleCollisions !Convert to no-dimensional units CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) - !Calculates reduced mass - collision%sMassInv = 1.D0/(mass_i+mass_j) - collision%rMass = (mass_i*mass_j)/collision%sMassInv - END SUBROUTINE initBinaryChargeExchange - SUBROUTINE collideBinaryChargeExchange(self, sigmaVrelMax, sigmaVrelMaxNew, & - part_i, part_j) + SUBROUTINE collideBinaryChargeExchange(self, part_i, part_j, vRel) USE moduleSpecies USE moduleRandom USE moduleMath IMPLICIT NONE CLASS(collisionBinaryChargeExchange), INTENT(in):: self - REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(inout):: sigmaVrelMaxNew + REAL(8), INTENT(in):: vRel TYPE(particle), INTENT(inout), TARGET:: part_i, part_j - REAL(8):: sigmaVrel - REAL(8):: vRel !relative velocity - REAL(8):: eRel !relative energy - !eRel (in units of [m][L]^2[t]^-2 - vRel = NORM2(part_i%v-part_j%v) - eRel = self%rMass*vRel**2 - sigmaVrel = self%crossSec%get(eRel)*vRel - sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel - IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN - SELECT TYPE(sp => part_i%species) - TYPE IS (speciesNeutral) - !Species i is neutral, ionize particle i - CALL sp%ionize(part_i) + SELECT TYPE(sp => part_i%species) + TYPE IS (speciesNeutral) + !Species i is neutral, ionize particle i + CALL sp%ionize(part_i) - TYPE IS (speciesCharged) - !Species i is charged, neutralize particle i - CALL sp%neutralize(part_i) + TYPE IS (speciesCharged) + !Species i is charged, neutralize particle i + CALL sp%neutralize(part_i) - END SELECT + END SELECT - SELECT TYPE(sp => part_j%species) - TYPE IS (speciesNeutral) - !Species j is neutral, ionize particle j - CALL sp%ionize(part_j) + SELECT TYPE(sp => part_j%species) + TYPE IS (speciesNeutral) + !Species j is neutral, ionize particle j + CALL sp%ionize(part_j) - TYPE IS (speciesCharged) - !Species j is charged, neutralize particle j - CALL sp%neutralize(part_j) + TYPE IS (speciesCharged) + !Species j is charged, neutralize particle j + CALL sp%neutralize(part_j) - END SELECT - - END IF + END SELECT END SUBROUTINE collideBinaryChargeExchange diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 31e2588..5d90f62 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -133,16 +133,16 @@ MODULE moduleInput IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference') !If a reference cross section is given, it is used - CALL config%get(object // '.crossSection', sigma_ref, found) + CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found) !If not, the reference radius is searched IF (.NOT. found) THEN CALL config%get(object // '.radius', r_ref, found) IF (found) THEN - sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section + sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity ELSE - sigma_ref = 0.D0 !Assume no collisions + sigmaVrel_ref = 0.D0 !Assume no collisions END IF @@ -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, species(pt_i)%obj%m, species(pt_j)%obj%m) DO k = 1, nCollisions WRITE (kString, '(I2)') k @@ -700,13 +700,11 @@ MODULE moduleInput SELECT CASE(cType) CASE ('elastic') !Elastic collision - CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m) + CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath) CASE ('chargeExchange') !Resonant charge exchange - CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m) + CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath) CASE ('ionization') !Electorn impact ionization @@ -715,7 +713,7 @@ MODULE moduleInput CALL config%get(object // '.electron', electron, found) IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions') CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron) + crossSecFilePath, energyThreshold, electron) CASE ('recombination') !Electorn impact ionization @@ -724,7 +722,7 @@ MODULE moduleInput CALL config%get(object // '.electron', electron, found) IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions') CALL initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, energyBinding, species(pt_i)%obj%m, species(pt_j)%obj%m, electron) + crossSecFilePath, energyBinding, electron) CASE DEFAULT CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions') diff --git a/src/modules/moduleMath.f90 b/src/modules/moduleMath.f90 index 9bc8453..ca8d780 100644 --- a/src/modules/moduleMath.f90 +++ b/src/modules/moduleMath.f90 @@ -49,4 +49,14 @@ MODULE moduleMath END FUNCTION norm1 + PURE FUNCTION reducedMass(m_i, m_j) RESULT(rMass) + IMPLICIT NONE + + REAL(8), INTENT(in):: m_i, m_j + REAL(8):: rMass + + rMass = (m_i * m_j) / (m_i + m_j) + + END FUNCTION + END MODULE moduleMath diff --git a/src/modules/moduleRefParam.f90 b/src/modules/moduleRefParam.f90 index 2dacc21..051190b 100644 --- a/src/modules/moduleRefParam.f90 +++ b/src/modules/moduleRefParam.f90 @@ -1,7 +1,7 @@ !Reference parameters MODULE moduleRefParam !Parameters that define the problem (inputs) - REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref + REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigmaVrel_ref !Reference parameters for non-dimensional problem REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 4142549..48d3dbc 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