fpakc/src/modules/moduleInject.f90

319 lines
8.8 KiB
Fortran

!injection of particles
MODULE moduleInject
USE moduleSpecies
!Generic type for velocity distribution function
TYPE, ABSTRACT:: velDistGeneric
CONTAINS
!Returns random velocity from distribution function
PROCEDURE(randomVel_interface), DEFERRED, PASS:: randomVel
END TYPE velDistGeneric
ABSTRACT INTERFACE
FUNCTION randomVel_interface(self) RESULT(v)
IMPORT velDistGeneric
CLASS(velDistGeneric), INTENT(in):: self
REAL(8):: v
END FUNCTION randomVel_interface
END INTERFACE
!Container for velocity distributions
TYPE:: velDistCont
CLASS(velDistGeneric), ALLOCATABLE:: obj
END TYPE velDistCont
!Maxwellian distribution function
TYPE, EXTENDS(velDistGeneric):: velDistMaxwellian
REAL(8):: v !Velocity
REAL(8):: vTh !Thermal Velocity
CONTAINS
PROCEDURE, PASS:: randomVel => randomVelMaxwellian
END TYPE velDistMaxwellian
!Dirac's delta distribution function
TYPE, EXTENDS(velDistGeneric):: velDistDelta
REAL(8):: v !Velocity
CONTAINS
PROCEDURE, PASS:: randomVel => randomVelDelta
END TYPE velDistDelta
!Generic injection of particles
TYPE:: injectGeneric
INTEGER:: id
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: vMod !Velocity (module)
REAL(8):: T(1:3) !Temperature
REAL(8):: n(1:3) !Direction of injection
INTEGER:: nParticles !Number of particles to introduce each time step
CLASS(speciesGeneric), POINTER:: species !Species of injection
INTEGER:: nEdges
INTEGER, ALLOCATABLE:: edges(:) !Array with edges
REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability
REAL(8):: sumWeight
TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
CONTAINS
PROCEDURE, PASS:: init => initInject
PROCEDURE, PASS:: addParticles
END TYPE injectGeneric
INTEGER:: nInject
TYPE(injectGeneric), ALLOCATABLE:: inject(:)
CONTAINS
!Initialize an injection of particles
SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
USE moduleMesh
USE moduleRefParam
USE moduleConstParam
USE moduleSpecies
USE moduleSolver
USE moduleErrors
IMPLICIT NONE
CLASS(injectGeneric), INTENT(inout):: self
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
INTEGER:: phSurface(1:mesh%numEdges)
INTEGER:: nVolColl
self%id = i
self%vMod = v / v_ref
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*tauInject*ti_ref/species(sp)%obj%weight)
CASE ("A")
!Input current in Ampers
self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight))
CASE ("part/s")
!Input current in Ampers
self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight)
CASE DEFAULT
CALL criticalError("No support for units: " // units, 'initInject')
END SELECT
!Scale particles for different species steps
IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
!Gets the edge elements from which particles are injected
DO e = 1, mesh%numEdges
phSurface(e) = mesh%edges(e)%obj%physicalSurface
END DO
self%nEdges = COUNT(phSurface == physicalSurface)
ALLOCATE(inject(i)%edges(1:self%nEdges))
et = 0
DO e=1, mesh%numEdges
IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
et = et + 1
self%edges(et) = mesh%edges(e)%obj%n
!Assign connectivity between injection edge and meshColl volume
IF (doubleMesh) THEN
nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos())
IF (nVolColl > 0) THEN
mesh%edges(e)%obj%eColl => meshColl%vols(nVolColl)%obj
ELSE
CALL criticalError("No connection between edge and meshColl", "initInject")
END IF
ELSE
IF (ASSOCIATED(mesh%edges(e)%obj%e1)) THEN
mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e1
ELSE
mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e2
END IF
END IF
END IF
END DO
!Calculates cumulative probability
ALLOCATE(self%cumWeight(1:self%nEdges))
et = 1
self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight
DO et = 2, self%nEdges
self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1)
END DO
self%sumWeight = self%cumWeight(self%nEdges)
END SUBROUTINE initInject
!Injection of particles
SUBROUTINE doInjects()
USE moduleSpecies
USE moduleSolver
IMPLICIT NONE
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
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
CALL inject(i)%addParticles()
END IF
END DO
END SUBROUTINE doInjects
SUBROUTINE initVelDistMaxwellian(velDist, v, T, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: v, T, m
velDist = velDistMaxwellian(v = v, vTh = DSQRT(T/m))
END SUBROUTINE initVelDistMaxwellian
SUBROUTINE initVelDistDelta(velDist, v)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: v
velDist = velDistDelta(v = v)
END SUBROUTINE initVelDistDelta
!Random velocity from Maxwellian distribution
FUNCTION randomVelMaxwellian(self) RESULT (v)
USE moduleRandom
IMPLICIT NONE
CLASS(velDistMaxwellian), INTENT(in):: self
REAL(8):: v
v = 0.D0
v = self%v + self%vTh*randomMaxwellian()
END FUNCTION randomVelMaxwellian
!Random velocity from Dirac's delta distribution
PURE FUNCTION randomVelDelta(self) RESULT(v)
IMPLICIT NONE
CLASS(velDistDelta), INTENT(in):: self
REAL(8):: v
v = self%v
END FUNCTION randomVelDelta
!Add particles for the injection
SUBROUTINE addParticles(self)
USE moduleSpecies
USE moduleSolver
USE moduleMesh
USE moduleRandom
USE moduleErrors
IMPLICIT NONE
CLASS(injectGeneric), INTENT(in):: self
INTEGER:: randomX
INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
INTEGER:: i
INTEGER:: n, sp
CLASS(meshEdge), POINTER:: randomEdge
!Insert particles
!$OMP SINGLE
nMin = 0
DO i = 1, self%id -1
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
nMin = nMin + inject(i)%nParticles
END IF
END DO
nMin = nMin + 1
nMax = nMin + self%nParticles - 1
!Assign weight to particle.
partInj(nMin:nMax)%weight = self%species%weight
!Particle is considered to be outside the domain
partInj(nMin:nMax)%n_in = .FALSE.
!$OMP END SINGLE
!$OMP DO
DO n = nMin, nMax
randomX = randomWeighted(self%cumWeight, self%sumWeight)
randomEdge => mesh%edges(self%edges(randomX))%obj
!Random position in edge
partInj(n)%r = randomEdge%randPos()
!Volume associated to the edge:
IF (ASSOCIATED(randomEdge%e1)) THEN
partInj(n)%vol = randomEdge%e1%n
ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
partInj(n)%vol = randomEdge%e2%n
ELSE
CALL criticalError("No Volume associated to edge", 'addParticles')
END IF
partInj(n)%volColl = randomEdge%eColl%n
sp = self%species%n
!Assign particle type
partInj(n)%species => self%species
partInj(n)%v = (/ self%v(1)%obj%randomVel(), &
self%v(2)%obj%randomVel(), &
self%v(3)%obj%randomVel() /)
!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(sp)%pushParticle(partInj(n), tau(sp))
!Assign cell to new particle
CALL solver%updateParticleCell(partInj(n))
END DO
!$OMP END DO
END SUBROUTINE addParticles
END MODULE moduleInject