diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index d45c7a0..366fa84 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -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 diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 1d2bdba..d0f8311 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -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 diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index ac22b28..40c218c 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -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 diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 8ef5c6d..e0f3a27 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -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 diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 1c713c9..e1b8907 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -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 diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 97602f4..2252b84 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -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 diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 0f54ee9..66a98df 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -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 diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 9a02cb1..becf0a6 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -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') diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 13d43e3..897ddb2 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -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') diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index 15c4b23..c08061a 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -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 diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 4b0124e..d75f073 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -437,7 +437,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 @@ -458,13 +458,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 @@ -522,6 +524,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