JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'development' into feature/3DCart

Browse source
This commit is contained in:
Jorge Gonzalez 2021-03-20 13:08:55 +01:00
commit 3b01099aba
11 changed files with 247 additions and 105 deletions
Download patch file Download diff file Expand all files Collapse all files

62
src/modules/mesh/1DCart/moduleMesh1DCart.f90
View file

@ -7,6 +7,9 @@ MODULE moduleMesh1DCart
USE moduleMeshBoundary
IMPLICIT NONE
REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /)
REAL(8), PARAMETER:: wSeg(1:3) = (/ 5.D0/9.D0 , 8.D0/9.D0, 5.D0/9.D0 /)
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DCart
!Element coordinates
REAL(8):: x = 0.D0
@ -152,23 +155,7 @@ MODULE moduleMesh1DCart
ALLOCATE(self%fboundary(1:nSpecies))
!Assign functions to boundary
DO s = 1, nSpecies
SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
TYPE IS(boundaryAbsorption)
self%fBoundary(s)%apply => absorption
TYPE IS(boundaryReflection)
self%fBoundary(s)%apply => reflection
TYPE IS(boundaryTransparent)
self%fBoundary(s)%apply => transparent
TYPE IS(boundaryWallTemperature)
self%fBoundary(s)%apply => wallTemperature
CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'initEdge1DCart')
END SELECT
CALL pointBoundaryFunction(self, s)
END DO
@ -322,22 +309,29 @@ MODULE moduleMesh1DCart
END SUBROUTINE partialDerSegm
!Computes local stiffness matrix
PURE FUNCTION elemKSegm(self) RESULT(ke)
FUNCTION elemKSegm(self) RESULT(ke)
IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(in):: self
REAL(8):: ke(1:2,1:2)
REAL(8):: Xii(1:3)
REAL(8):: dPsi(1:1, 1:2)
REAL(8):: invJ
REAL(8):: invJ(1), detJ
INTEGER:: l
ke = 0.D0
Xii = (/ 0.D0, 0.D0, 0.D0 /)
dPsi = self%dPsi(Xii)
invJ = self%invJac(Xii, dPsi)
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
Xii = 0.D0
DO l = 1, 3
xii(1) = corSeg(l)
dPsi = self%dPsi(Xii)
detJ = self%detJac(Xii, dPsi)
invJ = self%invJac(Xii, dPsi)
ke = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
wSeg(l)/detJ
END DO
END FUNCTION elemKSegm
@ -348,14 +342,22 @@ MODULE moduleMesh1DCart
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:2)
REAL(8):: detJ
REAL(8):: detJ, f
REAL(8):: Xii(1:3)
INTEGER:: l
Xii = 0.D0
fPsi = self%fPsi(Xii)
detJ = self%detJac(Xii)
ALLOCATE(localF(1:2))
localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ
localF = 0.D0
Xii = 0.D0
DO l = 1, 3
xii(1) = corSeg(l)
detJ = self%detJac(Xii)
fPsi = self%fPsi(Xii)
f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*wSeg(l)*detJ
END DO
END FUNCTION elemFSegm

18
src/modules/mesh/1DRad/moduleMesh1DRad.f90
View file

@ -153,23 +153,7 @@ MODULE moduleMesh1DRad
ALLOCATE(self%fboundary(1:nSpecies))
!Assign functions to boundary
DO s = 1, nSpecies
SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
TYPE IS(boundaryAbsorption)
self%fBoundary(s)%apply => absorption
TYPE IS(boundaryReflection)
self%fBoundary(s)%apply => reflection
TYPE IS(boundaryTransparent)
self%fBoundary(s)%apply => transparent
TYPE IS(boundaryWallTemperature)
self%fBoundary(s)%apply => wallTemperature
CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'initEdge1DRad')
END SELECT
CALL pointBoundaryFunction(self, s)
END DO

18
src/modules/mesh/2DCart/moduleMesh2DCart.f90
View file

@ -200,23 +200,7 @@ MODULE moduleMesh2DCart
ALLOCATE(self%fboundary(1:nSpecies))
!Assign functions to boundary
DO s = 1, nSpecies
SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
TYPE IS(boundaryAbsorption)
self%fBoundary(s)%apply => absorption
TYPE IS(boundaryReflection)
self%fBoundary(s)%apply => reflection
TYPE IS(boundaryTransparent)
self%fBoundary(s)%apply => transparent
TYPE IS(boundaryWallTemperature)
self%fBoundary(s)%apply => wallTemperature
CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'initEdge2DCart')
END SELECT
CALL pointBoundaryFunction(self, s)
END DO

25
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -200,26 +200,7 @@ MODULE moduleMesh2DCyl
ALLOCATE(self%fboundary(1:nSpecies))
!Assign functions to boundary
DO s = 1, nSpecies
SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
TYPE IS(boundaryAbsorption)
self%fBoundary(s)%apply => absorption
TYPE IS(boundaryReflection)
self%fBoundary(s)%apply => reflection
TYPE IS(boundaryTransparent)
self%fBoundary(s)%apply => transparent
TYPE IS(boundaryWallTemperature)
self%fBoundary(s)%apply => wallTemperature
TYPE IS(boundaryAxis)
self%fBoundary(s)%apply => symmetryAxis
CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'initEdge2DCyl')
END SELECT
CALL pointBoundaryFunction(self, s)
END DO
@ -459,7 +440,9 @@ MODULE moduleMesh2DCyl
detJ = self%detJac(xi,dPsi)
invJ = self%invJac(xi,dPsi)
r = DOT_PRODUCT(fPsi,self%r)
ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wQuad(l)*wQuad(m)/detJ
ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
MATMUL(invJ,dPsi))* &
r*wQuad(l)*wQuad(m)/detJ
END DO
END DO

18
src/modules/mesh/3DCart/moduleMesh3DCart.f90
View file

@ -166,23 +166,7 @@ MODULE moduleMesh3DCart
ALLOCATE(self%fBoundary(1:nSpecies))
!Assign functions to boundary
DO s = 1, nSpecies
SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
TYPE IS(boundaryAbsorption)
self%fBoundary(s)%apply => absorption
TYPE IS(boundaryReflection)
self%fBoundary(s)%apply => reflection
TYPE IS(boundaryTransparent)
self%fBoundary(s)%apply => transparent
TYPE IS(boundaryWallTemperature)
self%fBoundary(s)%apply => wallTemperature
CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'initEdge3DCart')
END SELECT
CALL pointBoundaryFunction(self, s)
END DO

113
src/modules/mesh/moduleMeshBoundary.f90
View file

@ -108,6 +108,85 @@ MODULE moduleMeshBoundary
END SUBROUTINE wallTemperature
!Ionization surface: an electron will pass through the surface
! and create an ion-electron pair based on a neutral background
SUBROUTINE ionization(edge, part)
USE moduleList
USE moduleSpecies
USE moduleMesh
USE moduleRefParam
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
REAL(8):: ionizationRate
INTEGER:: ionizationPair, p
TYPE(particle), POINTER:: newElectron
TYPE(particle), POINTER:: newIon
SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
TYPE IS(boundaryIonization)
mRel = (bound%m0*species(part%sp)%obj%m)*(bound%m0+species(part%sp)%obj%m)
vRel = SUM(DABS(part%v-bound%v0))
eRel = mRel*vRel**2*5.D-1
IF (eRel > bound%eThreshold) THEN
!TODO: Check units
ionizationRate = bound%n0*bound%crossSection%get(eRel)
ionizationPair = INT(ionizationRate*tauMin*ti_ref/species(bound%sp)%obj%weight)
!Create the new pair of particles
DO p = 1, ionizationPair
ALLOCATE(newElectron)
ALLOCATE(newIon)
newElectron%sp = part%sp
newIon%sp = bound%sp
newElectron%v = 10.D0*bound%v0 + (1.D0 + bound%deltaV/NORM2(bound%v0))
newIon%v = bound%v0
newElectron%r = edge%randPos()
newIon%r = newElectron%r
newElectron%vol = part%vol
newIon%vol = part%vol
newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r)
newIon%xi = newElectron%xi
newElectron%qm = part%qm
SELECT TYPE(spe => species(bound%sp)%obj)
TYPE IS(speciesCharged)
newIon%qm = spe%qm
END SELECT
newElectron%weight = species(bound%sp)%obj%weight
newIon%weight = newElectron%weight
newElectron%n_in = .TRUE.
newIon%n_in = .TRUE.
!Add particles to list
CALL OMP_SET_LOCK(lockSurfaces)
CALL partSurfaces%add(newElectron)
CALL partSurfaces%add(newIon)
CALL OMP_UNSET_LOCK(lockSurfaces)
END DO
!Removes ionizing electron
part%n_in = .FALSE.
END IF
END SELECT
END SUBROUTINE ionization
!Symmetry axis. Dummy function
SUBROUTINE symmetryAxis(edge, part)
USE moduleSpecies
@ -118,4 +197,38 @@ MODULE moduleMeshBoundary
END SUBROUTINE symmetryAxis
!Points the boundary function to specific type
SUBROUTINE pointBoundaryFunction(edge, s)
USE moduleErrors
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
INTEGER, INTENT(in):: s !Species index
SELECT TYPE(obj => edge%boundary%bTypes(s)%obj)
TYPE IS(boundaryAbsorption)
edge%fBoundary(s)%apply => absorption
TYPE IS(boundaryReflection)
edge%fBoundary(s)%apply => reflection
TYPE IS(boundaryTransparent)
edge%fBoundary(s)%apply => transparent
TYPE IS(boundaryWallTemperature)
edge%fBoundary(s)%apply => wallTemperature
TYPE IS(boundaryIonization)
edge%fBoundary(s)%apply => ionization
TYPE IS(boundaryAxis)
edge%fBoundary(s)%apply => symmetryAxis
CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')
END SELECT
END SUBROUTINE pointBoundaryFunction
END MODULE moduleMeshBoundary

45
src/modules/moduleBoundary.f90
View file

@ -1,4 +1,5 @@
MODULE moduleBoundary
USE moduleTable
!Generic type for boundaries
TYPE, PUBLIC:: boundaryGeneric
@ -24,7 +25,7 @@ MODULE moduleBoundary
END TYPE boundaryTransparent
!Transparent boundary
!Wall Temperature boundary
TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature
!Thermal velocity of the wall: square root(Wall temperature X specific heat)
REAL(8):: vTh
@ -32,6 +33,18 @@ MODULE moduleBoundary
END TYPE boundaryWallTemperature
!Ionization boundary
TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
REAL(8):: m0, n0, v0(1:3), T0 !Properties of background neutrals.
INTEGER:: sp !Ion species
TYPE(table1D):: crossSection
REAL(8):: eThreshold
REAL(8):: deltaV
CONTAINS
END TYPE boundaryIonization
!Symmetry axis
TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
CONTAINS
@ -86,4 +99,34 @@ MODULE moduleBoundary
END SUBROUTINE initWallTemperature
SUBROUTINE initIonization(boundary, m0, n0, v0, T0, speciesID, crossSection, eThreshold)
USE moduleRefParam
USE moduleSpecies
USE moduleCaseParam
USE moduleConstParam
IMPLICIT NONE
CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
REAL(8), INTENT(in):: m0, n0, v0(1:3), T0
INTEGER:: speciesID
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
REAL(8), INTENT(in):: eThreshold
ALLOCATE(boundaryIonization:: boundary)
SELECT TYPE(boundary)
TYPE IS(boundaryIonization)
boundary%m0 = m0 / m_ref
boundary%n0 = n0
boundary%v0 = v0 / v_ref
boundary%T0 = T0 / T_ref
boundary%sp = speciesID
CALL boundary%crossSection%init(crossSection)
CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2)
END SELECT
END SUBROUTINE initIonization
END MODULE moduleBoundary

4
src/modules/moduleInject.f90
View file

@ -98,6 +98,10 @@ MODULE moduleInject
!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')

26
src/modules/moduleInput.f90
View file

@ -599,6 +599,13 @@ MODULE moduleInput
CHARACTER(2):: istring, sString
CHARACTER(:), ALLOCATABLE:: object, bType
REAL(8):: Tw, cw !Wall temperature and specific heat
!Neutral Properties
REAL(8):: m0, n0
REAL(8), DIMENSION(:), ALLOCATABLE:: v0
REAL(8):: T0
REAL(8):: eThreshold !Energy threshold
INTEGER:: speciesID
CHARACTER(:), ALLOCATABLE:: speciesName, crossSection
LOGICAL:: found
INTEGER:: nTypes
@ -628,6 +635,25 @@ MODULE moduleInput
CASE('transparent')
ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj)
CASE('ionization')
CALL config%get(object // '.neutral.name', speciesName, found)
IF (.NOT. found) CALL criticalError("missing parameter 'name' for neutrals in ionization", 'readBoundary')
speciesID = speciesName2Index(speciesName)
CALL config%get(object // '.neutral.mass', m0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'mass' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.neutral.density', n0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.neutral.velocity', v0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.neutral.temperature', T0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.energyThreshold', eThreshold, found)
IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
CALL config%get(object // '.crossSection', crossSection, found)
IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
CALL initIonization(boundary(i)%bTypes(s)%obj, m0, n0, v0, T0, speciesID, crossSection, eThreshold)
CASE('wallTemperature')
CALL config%get(object // '.temperature', Tw, found)
IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')

2
src/modules/moduleList.f90
View file

@ -24,6 +24,8 @@ MODULE moduleList
INTEGER(KIND=OMP_LOCK_KIND):: lockWScheme !Lock for the NA list of particles
TYPE(listNode):: partCollisions !Particles created in collisional process
INTEGER(KIND=OMP_LOCK_KIND):: lockCollisions !Lock for the NA list of particles
TYPE(listNode):: partSurfaces !Particles created in surface interactions
INTEGER(KIND=OMP_LOCK_KIND):: lockSurfaces !Lock for the NA list of particles
TYPE(listNode):: partInitial !Initial distribution of particles
TYPE pointerArray

21
src/modules/moduleSolver.f90
View file

@ -506,7 +506,7 @@ MODULE moduleSolver
INTEGER:: nn, n, e
INTEGER, SAVE:: nPartNew
INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions
INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces
TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
TYPE(lNode), POINTER:: partCurr, partNext
@ -527,13 +527,15 @@ MODULE moduleSolver
nWScheme = partWScheme%amount
!$OMP SECTION
nCollisions = partCollisions%amount
!$OMP SECTION
nSurfaces = partSurfaces%amount
!$OMP END SECTIONS
!$OMP BARRIER
!$OMP SINGLE
CALL MOVE_ALLOC(partOld, partTemp)
nPartNew = nInjIn + nOldIn + nWScheme + nCollisions
nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces
ALLOCATE(partOld(1:nPartNew))
!$OMP END SINGLE
@ -591,6 +593,21 @@ MODULE moduleSolver
IF (ASSOCIATED(partCollisions%tail)) NULLIFY(partCollisions%tail)
partCollisions%amount = 0
!$OMP SECTION
!Reset particles from surface process
nn = nInjIn + nOldIn + nWScheme + nCollisions
partCurr => partSurfaces%head
DO n = 1, nSurfaces
partNext => partCurr%next
partOld(nn+n) = partCurr%part
DEALLOCATE(partCurr)
partCurr => partNext
END DO
IF (ASSOCIATED(partSurfaces%head)) NULLIFY(partSurfaces%head)
IF (ASSOCIATED(partSurfaces%tail)) NULLIFY(partSurfaces%tail)
partSurfaces%amount = 0
!$OMP SECTION
!Reset output in nodes
DO e = 1, mesh%numNodes