fpakc/src/modules/moduleInject.f90
Jorge Gonzalez 2a843547b8 Ionization boundary condition fully tested.
Documentation updated properly.

3D Cartesian geometry also tested.
Documentation updated properly.

Added weighting probability in the injection of particles.
2021-03-27 11:38:18 +01:00

276 lines
7.7 KiB
Fortran

!injection of particles
MODULE moduleInject
!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
INTEGER:: sp !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), INTENT(in):: flow
CHARACTER(:), ALLOCATABLE, INTENT(in):: units
INTEGER:: e, et
INTEGER:: phSurface(1:mesh%numEdges)
self%id = i
self%vMod = v/v_ref
self%n = n
self%T = T/T_ref
self%sp = sp
SELECT CASE(units)
CASE ("sccm")
!Standard cubic centimeter per minute
self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight)
CASE ("A")
!Input current in Ampers
self%nParticles = INT(flow*tauMin*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)
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
END IF
END DO
!Calculates cumulative probability
ALLOCATE(self%cumWeight(1:self%nEdges))
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)
nPartInj = nPartInj + self%nParticles
END SUBROUTINE initInject
!Injection of particles
SUBROUTINE doInjects()
USE moduleSpecies
USE moduleSolver
IMPLICIT NONE
INTEGER:: i
!$OMP SINGLE
IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
ALLOCATE(partInj(1:nPartInj))
!$OMP END SINGLE
DO i=1, nInject
CALL inject(i)%addParticles()
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:: n
CLASS(meshEdge), POINTER:: randomEdge
!Insert particles
!$OMP SINGLE
nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
nMax = nMin + self%nParticles - 1
!Assign particle type
partInj(nMin:nMax)%sp = self%sp
!Assign weight to particle.
partInj(nMin:nMax)%weight = species(self%sp)%obj%weight
!Particle is considered to be outside the domain
partInj(nMin:nMax)%n_in = .FALSE.
!Assign charge/mass to charged particle.
SELECT TYPE(sp => species(self%sp)%obj)
TYPE IS(speciesCharged)
partInj(nMin:nMax)%qm = sp%qm
END SELECT
!$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)%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(self%sp)%pushParticle(partInj(n), tauMin)
!Assign cell to new particle
CALL solver%updateParticleCell(partInj(n))
END DO
!$OMP END DO
END SUBROUTINE addParticles
END MODULE moduleInject