Implementation of charge exchange and structure for ionization
processes.
This commit is contained in:
parent
baf25c1157
commit
35936ea918
9 changed files with 1604 additions and 29 deletions
|
|
@ -125,7 +125,7 @@ MODULE moduleMesh
|
|||
!Volume index
|
||||
INTEGER:: n = 0
|
||||
!Maximum collision rate
|
||||
REAL(8):: sigmaVrelMax = 1.D-15
|
||||
REAL(8):: sigmaVrelMax = 1.D-17
|
||||
!Volume
|
||||
REAL(8):: volume = 0.D0
|
||||
!List of particles inside the volume
|
||||
|
|
|
|||
|
|
@ -1,11 +1,11 @@
|
|||
MODULE moduleCollisions
|
||||
USE moduleTable
|
||||
|
||||
!Abstract type for collision between two particles
|
||||
TYPE, ABSTRACT:: collisionBinary
|
||||
REAL(8):: rMass !reduced mass
|
||||
TYPE(table1D):: crossSec !cross section of collision
|
||||
CONTAINS
|
||||
PROCEDURE(initBinary_interface), PASS, DEFERRED:: init
|
||||
PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide
|
||||
|
||||
END TYPE collisionBinary
|
||||
|
|
@ -38,16 +38,31 @@ MODULE moduleCollisions
|
|||
|
||||
END TYPE collisionCont
|
||||
|
||||
!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:: init => initBinaryElastic
|
||||
PROCEDURE, PASS:: collide => collideBinaryElastic
|
||||
|
||||
END TYPE collisionBinaryElastic
|
||||
|
||||
!Ionization binary interaction
|
||||
TYPE, EXTENDS(collisionBinary):: collisionBinaryIonization
|
||||
REAL(8):: eThreshold !Minimum energy (non-dimensional units) required for ionization
|
||||
CONTAINS
|
||||
PROCEDURE, PASS:: collide => collideBinaryIonization
|
||||
|
||||
END TYPE collisionBinaryIonization
|
||||
|
||||
!Resonant charge-exchange
|
||||
TYPE, EXTENDS(collisionBinary):: collisionBinaryChargeExchange
|
||||
CONTAINS
|
||||
PROCEDURE, PASS:: collide => collideBinaryChargeExchange
|
||||
|
||||
END TYPE collisionBinaryChargeExchange
|
||||
|
||||
!Type for interaction matrix
|
||||
TYPE:: interactionsBinary
|
||||
INTEGER:: amount
|
||||
|
|
@ -58,11 +73,12 @@ MODULE moduleCollisions
|
|||
END TYPE interactionsBinary
|
||||
|
||||
!Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory
|
||||
TYPE(interactionsBinary), ALLOCATABLE:: interactionMatrix(:)
|
||||
TYPE(interactionsBinary), ALLOCATABLE, TARGET:: interactionMatrix(:)
|
||||
!Folder for collision cross section tables
|
||||
CHARACTER(:), ALLOCATABLE:: pathCollisions
|
||||
|
||||
CONTAINS
|
||||
!Inits the interaction matrix
|
||||
SUBROUTINE initInteractionMatrix(interactionMatrix)
|
||||
USE moduleSpecies
|
||||
IMPLICIT NONE
|
||||
|
|
@ -75,6 +91,7 @@ MODULE moduleCollisions
|
|||
|
||||
END SUBROUTINE initInteractionMatrix
|
||||
|
||||
!Gets the interaction index from the collision matrix from index i,j
|
||||
FUNCTION interactionIndex(i,j) RESULT(k)
|
||||
|
||||
INTEGER:: i, j
|
||||
|
|
@ -87,50 +104,49 @@ MODULE moduleCollisions
|
|||
|
||||
END FUNCTION interactionIndex
|
||||
|
||||
!Inits the binary interaction
|
||||
SUBROUTINE initInteractionBinary(self, amount)
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(interactionsBinary), INTENT(inout):: self
|
||||
INTEGER, INTENT(in):: amount
|
||||
INTEGER:: k
|
||||
|
||||
self%amount = amount
|
||||
|
||||
ALLOCATE(self%collisions(1:self%amount))
|
||||
DO k= 1, self%amount
|
||||
!TODO: make type dependent
|
||||
ALLOCATE(collisionBinaryElastic:: self%collisions(k)%obj)
|
||||
|
||||
END DO
|
||||
|
||||
END SUBROUTINE initInteractionBinary
|
||||
|
||||
SUBROUTINE initBinaryElastic(self, crossSectionFilename, mass_i, mass_j)
|
||||
!ELASTIC COLLISIONS
|
||||
!Inits binary elastic collision
|
||||
SUBROUTINE initBinaryElastic(collision, crossSectionFilename, mass_i, mass_j)
|
||||
USE moduleTable
|
||||
USE moduleRefParam
|
||||
USE moduleConstParam
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(collisionBinaryElastic), INTENT(inout):: self
|
||||
CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
|
||||
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
|
||||
REAL(8), INTENT(in):: mass_i, mass_j
|
||||
|
||||
ALLOCATE(collisionBinaryElastic:: collision)
|
||||
|
||||
!Reads data from file
|
||||
CALL self%crossSec%init(crossSectionFilename)
|
||||
CALL collision%crossSec%init(crossSectionFilename)
|
||||
|
||||
!Convert to no-dimensional units
|
||||
CALL self%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
|
||||
CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
|
||||
|
||||
!Calculates reduced mass
|
||||
self%rMass = (mass_i*mass_j)/(mass_i+mass_j)
|
||||
collision%rMass = (mass_i*mass_j)/(mass_i+mass_j)
|
||||
|
||||
END SUBROUTINE
|
||||
END SUBROUTINE initBinaryElastic
|
||||
|
||||
!Binary elastic process
|
||||
SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrelMaxNew, &
|
||||
part_i, part_j)
|
||||
USE moduleConstParam
|
||||
USE moduleSpecies
|
||||
USE moduleTable
|
||||
USE moduleConstParam
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(collisionBinaryElastic), INTENT(in):: self
|
||||
|
|
@ -145,7 +161,7 @@ MODULE moduleCollisions
|
|||
REAL(8):: alpha !random angle of scattering
|
||||
REAL(8):: rnd
|
||||
|
||||
!eRel (in units of [m][L]^2[s]^-2
|
||||
!eRel (in units of [m][L]^2[t]^-2
|
||||
vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1
|
||||
eRel = self%rMass*vRel**2
|
||||
sigmaVrel = self%crossSec%get(eRel)*vRel
|
||||
|
|
@ -169,7 +185,124 @@ MODULE moduleCollisions
|
|||
|
||||
END IF
|
||||
|
||||
END SUBROUTINE collideBinaryElastic
|
||||
|
||||
END SUBROUTINE
|
||||
!ELECTRON IMPACT IONIZATION
|
||||
!Inits electron impact ionization
|
||||
SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, mass_i, mass_j)
|
||||
USE moduleTable
|
||||
USE moduleRefParam
|
||||
USE moduleConstParam
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
|
||||
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
|
||||
REAL(8), INTENT(in):: energyThreshold
|
||||
REAL(8), INTENT(in):: mass_i, mass_j
|
||||
|
||||
ALLOCATE(collisionBinaryIonization:: collision)
|
||||
|
||||
!Reads data from file
|
||||
CALL collision%crossSec%init(crossSectionFilename)
|
||||
|
||||
!Convert to no-dimensional units
|
||||
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)
|
||||
|
||||
SELECT TYPE(collision)
|
||||
TYPE IS(collisionBinaryIonization)
|
||||
!Assign the energy threshold
|
||||
collision%eThreshold = energyThreshold/(m_ref*v_ref**2)
|
||||
|
||||
END SELECT
|
||||
|
||||
END SUBROUTINE initBinaryIonization
|
||||
|
||||
!Binary electron impact ionization process
|
||||
SUBROUTINE collideBinaryIonization(self, sigmaVrelMax, sigmaVrelMaxNew, &
|
||||
part_i, part_j)
|
||||
USE moduleSpecies
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(collisionBinaryIonization), INTENT(in):: self
|
||||
REAL(8), INTENT(in):: sigmaVrelMax
|
||||
REAL(8), INTENT(inout):: sigmaVrelMaxNew
|
||||
TYPE(particle), INTENT(inout):: part_i, part_j
|
||||
|
||||
END SUBROUTINE collideBinaryIonization
|
||||
|
||||
!RESONANT CHARGE EXCHANGE
|
||||
!Inits resonant charge exchange
|
||||
SUBROUTINE initBinaryChargeExchange(collision, crossSectionFilename, mass_i, mass_j)
|
||||
USE moduleTable
|
||||
USE moduleRefParam
|
||||
USE moduleConstParam
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
|
||||
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
|
||||
REAL(8), INTENT(in):: mass_i, mass_j
|
||||
|
||||
ALLOCATE(collisionBinaryChargeExchange:: collision)
|
||||
|
||||
!Reads data from file
|
||||
CALL collision%crossSec%init(crossSectionFilename)
|
||||
|
||||
!Convert to no-dimensional units
|
||||
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)
|
||||
|
||||
END SUBROUTINE initBinaryChargeExchange
|
||||
|
||||
SUBROUTINE collideBinaryChargeExchange(self, sigmaVrelMax, sigmaVrelMaxNew, &
|
||||
part_i, part_j)
|
||||
USE moduleSpecies
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(collisionBinaryChargeExchange), INTENT(in):: self
|
||||
REAL(8), INTENT(in):: sigmaVrelMax
|
||||
REAL(8), INTENT(inout):: sigmaVrelMaxNew
|
||||
TYPE(particle), INTENT(inout):: part_i, part_j
|
||||
REAL(8):: sigmaVrel
|
||||
REAL(8):: vRel !relative velocity
|
||||
REAL(8):: eRel !relative energy
|
||||
REAL(8):: rnd
|
||||
|
||||
!eRel (in units of [m][L]^2[t]^-2
|
||||
vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1
|
||||
eRel = self%rMass*vRel**2
|
||||
sigmaVrel = self%crossSec%get(eRel)*vRel
|
||||
sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
|
||||
CALL RANDOM_NUMBER(rnd)
|
||||
IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) THEN
|
||||
SELECT TYPE(sp => species(part_i%sp)%obj)
|
||||
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)
|
||||
|
||||
END SELECT
|
||||
|
||||
SELECT TYPE(sp => species(part_j%sp)%obj)
|
||||
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)
|
||||
|
||||
END SELECT
|
||||
|
||||
END IF
|
||||
|
||||
END SUBROUTINE collideBinaryChargeExchange
|
||||
|
||||
END MODULE moduleCollisions
|
||||
|
|
|
|||
|
|
@ -85,7 +85,6 @@ MODULE moduleInput
|
|||
|
||||
!Derived parameters
|
||||
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
|
||||
!TODO: Make this solver dependent
|
||||
IF (found_r) THEN
|
||||
sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section
|
||||
L_ref = 1.D0/(sigma_ref*n_ref) !mean free path
|
||||
|
|
@ -238,6 +237,8 @@ MODULE moduleInput
|
|||
REAL(8):: mass, charge
|
||||
LOGICAL:: found
|
||||
INTEGER:: i
|
||||
CHARACTER(:), ALLOCATABLE:: linkName
|
||||
INTEGER:: linkID
|
||||
|
||||
!Gets the number of species
|
||||
CALL config%info('species', found, n_children = nSpecies)
|
||||
|
|
@ -276,6 +277,39 @@ MODULE moduleInput
|
|||
|
||||
END DO
|
||||
|
||||
!Reads relations between species
|
||||
DO i = 1, nSpecies
|
||||
SELECT TYPE(sp => species(i)%obj)
|
||||
TYPE IS (speciesNeutral)
|
||||
!Gets species linked ion
|
||||
CALL config%get(object // '.ion', linkName, found)
|
||||
IF (found) THEN
|
||||
linkID = speciesName2Index(linkName)
|
||||
sp%ion => species(linkID)%obj
|
||||
|
||||
END IF
|
||||
|
||||
TYPE IS (speciesCharged)
|
||||
!Gets species linked neutral
|
||||
CALL config%get(object // '.neutral', linkName)
|
||||
IF (found) THEN
|
||||
linkID = speciesName2Index(linkName)
|
||||
sp%neutral => species(linkID)%obj
|
||||
|
||||
END IF
|
||||
|
||||
!Gets species linked ion
|
||||
CALL config%get(object // '.ion', linkName, found)
|
||||
IF (found) THEN
|
||||
linkID = speciesName2Index(linkName)
|
||||
sp%ion => species(linkID)%obj
|
||||
|
||||
END IF
|
||||
|
||||
END SELECT
|
||||
|
||||
END DO
|
||||
|
||||
!Set number of particles to 0 for init state
|
||||
!TODO: In a future, this should include the particles from init states
|
||||
nPartOld = 0
|
||||
|
|
@ -289,6 +323,7 @@ MODULE moduleInput
|
|||
SUBROUTINE readInteractions(config)
|
||||
USE moduleSpecies
|
||||
USE moduleCollisions
|
||||
USE moduleErrors
|
||||
USE json_module
|
||||
IMPLICIT NONE
|
||||
|
||||
|
|
@ -298,10 +333,12 @@ MODULE moduleInput
|
|||
CHARACTER(:), ALLOCATABLE:: species_i, species_j
|
||||
CHARACTER(:), ALLOCATABLE:: crossSecFile
|
||||
CHARACTER(:), ALLOCATABLE:: crossSecFilePath
|
||||
CHARACTER(:), ALLOCATABLE:: cType
|
||||
LOGICAL:: found
|
||||
INTEGER:: nInteractions, nCollisions
|
||||
INTEGER:: i, k, ij
|
||||
INTEGER:: pt_i, pt_j
|
||||
REAL(8):: energyThreshold
|
||||
|
||||
CALL initInteractionMatrix(interactionMatrix)
|
||||
|
||||
|
|
@ -315,15 +352,43 @@ MODULE moduleInput
|
|||
pt_i = speciesName2Index(species_i)
|
||||
CALL config%get(object // '.species_j', species_j, found)
|
||||
pt_j = speciesName2Index(species_j)
|
||||
CALL config%info(object // '.crossSections', found, n_children = nCollisions)
|
||||
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)
|
||||
|
||||
DO k = 1, nCollisions
|
||||
WRITE (kString, '(I2)') k
|
||||
CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile, found)
|
||||
object = 'interactions.collisions(' // TRIM(iString) // ').cTypes(' // TRIM(kString) // ')'
|
||||
!Reads the cross section file
|
||||
CALL config%get(object // '.crossSection', crossSecFile, found)
|
||||
crossSecFilePath = pathCollisions // crossSecFile
|
||||
CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
|
||||
IF (.NOT. found) CALL criticalError('crossSection not found for ' // object, 'readInteractions')
|
||||
!Reads the type of collision
|
||||
CALL config%get(object // '.type', cType, found)
|
||||
!Initialize collision type and reads required additional data
|
||||
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)
|
||||
|
||||
CASE ('chargeExchange')
|
||||
!Resonant charge exchange
|
||||
CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, &
|
||||
crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
|
||||
|
||||
CASE ('ionization')
|
||||
!Electorn impact ionization
|
||||
CALL config%get(object // '.energyThreshold', found)
|
||||
IF (.NOT. found) CALL criticalError('energyThreshold 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)
|
||||
|
||||
CASE DEFAULT
|
||||
CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions')
|
||||
|
||||
END SELECT
|
||||
|
||||
END DO
|
||||
|
||||
|
|
|
|||
|
|
@ -11,11 +11,18 @@ MODULE moduleSpecies
|
|||
END TYPE speciesGeneric
|
||||
|
||||
TYPE, EXTENDS(speciesGeneric):: speciesNeutral
|
||||
CLASS(speciesGeneric), POINTER:: ion => NULL()
|
||||
CONTAINS
|
||||
PROCEDURE, PASS:: ionize => ionizeNeutral
|
||||
|
||||
END TYPE speciesNeutral
|
||||
|
||||
TYPE, EXTENDS(speciesGeneric):: speciesCharged
|
||||
REAL(8):: q=0.D0, qm=0.D0
|
||||
CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL()
|
||||
CONTAINS
|
||||
PROCEDURE, PASS:: ionize => ionizeCharged
|
||||
PROCEDURE, PASS:: neutralize => neutralizeCharged
|
||||
|
||||
END TYPE speciesCharged
|
||||
|
||||
|
|
@ -25,7 +32,7 @@ MODULE moduleSpecies
|
|||
END TYPE
|
||||
|
||||
INTEGER:: nSpecies
|
||||
TYPE(speciesCont), ALLOCATABLE:: species(:)
|
||||
TYPE(speciesCont), ALLOCATABLE, TARGET:: species(:)
|
||||
|
||||
TYPE particle
|
||||
REAL(8):: r(1:3) !Position
|
||||
|
|
@ -68,4 +75,59 @@ MODULE moduleSpecies
|
|||
|
||||
END FUNCTION speciesName2Index
|
||||
|
||||
!Change particle type to corresponding ion (neutral species)
|
||||
SUBROUTINE ionizeNeutral(self, part)
|
||||
USE moduleErrors
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(speciesNeutral), INTENT(IN):: self
|
||||
TYPE(particle), INTENT(inout):: part
|
||||
|
||||
IF (ASSOCIATED(self%ion)) THEN
|
||||
part%sp = self%ion%sp
|
||||
|
||||
ELSE
|
||||
CALL criticalError('No ion defined for species' // self%name, 'ionizeNeutral')
|
||||
|
||||
END IF
|
||||
|
||||
END SUBROUTINE ionizeNeutral
|
||||
|
||||
!Change particle type to corresponding ion (charged species)
|
||||
SUBROUTINE ionizeCharged(self, part)
|
||||
USE moduleErrors
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(speciesCharged), INTENT(IN):: self
|
||||
TYPE(particle), INTENT(inout):: part
|
||||
|
||||
IF (ASSOCIATED(self%ion)) THEN
|
||||
part%sp = self%ion%sp
|
||||
|
||||
ELSE
|
||||
CALL criticalError('No ion defined for species' // self%name, 'ionizeCharged')
|
||||
|
||||
END IF
|
||||
|
||||
END SUBROUTINE ionizeCharged
|
||||
|
||||
!Change particle type to corresponding neutral
|
||||
SUBROUTINE neutralizeCharged(self, part)
|
||||
USE moduleErrors
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(speciesCharged), INTENT(in):: self
|
||||
TYPE(particle), INTENT(inout):: part
|
||||
|
||||
IF (ASSOCIATED(self%neutral)) THEN
|
||||
part%sp = self%neutral%sp
|
||||
|
||||
ELSE
|
||||
CALL criticalError('No neutral defined for species' // self%name, 'neutralizeCharged')
|
||||
|
||||
END IF
|
||||
|
||||
END SUBROUTINE neutralizeCharged
|
||||
|
||||
END MODULE moduleSpecies
|
||||
|
||||
|
|
|
|||
|
|
@ -28,10 +28,10 @@ MODULE moduleTable
|
|||
amount = 0
|
||||
DO
|
||||
READ(id, '(A)', iostat = stat) dummy
|
||||
!Skip comment
|
||||
IF (INDEX(dummy,'#') /= 0) CYCLE
|
||||
!If EOF or error, exit file
|
||||
IF (stat /= 0) EXIT
|
||||
!Skip comment
|
||||
IF (INDEX(dummy,'#') /= 0) CYCLE
|
||||
!Add row
|
||||
amount = amount + 1
|
||||
|
||||
|
|
@ -62,7 +62,7 @@ MODULE moduleTable
|
|||
|
||||
END DO
|
||||
|
||||
CLOSE(10)
|
||||
CLOSE(id)
|
||||
|
||||
self%xMin = self%x(1)
|
||||
self%xMax = self%x(amount)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue