diff --git a/runs/ALPHIE_Grid/inputIonization_0.10mA.json b/runs/ALPHIE_Grid/inputIonization_0.10mA.json index 2a44da4..11a12f3 100644 --- a/runs/ALPHIE_Grid/inputIonization_0.10mA.json +++ b/runs/ALPHIE_Grid/inputIonization_0.10mA.json @@ -5,6 +5,7 @@ "cpuTime": false, "numColl": false, "EMField": true, + "boundariesParticle": true, "folder": "ionization_0.10mA" }, "geometry": { diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 2e587a1..d1a6920 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -694,9 +694,9 @@ MODULE moduleMesh use moduleSpecies import boundaryParticleGeneric, meshEdge - class(boundaryParticleGeneric), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(boundaryParticleGeneric), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine @@ -712,7 +712,7 @@ MODULE moduleMesh subroutine printParticle_interface(self, fileID) import boundaryParticleGeneric - class(boundaryParticleGeneric), intent(in):: self + class(boundaryParticleGeneric), intent(inout):: self integer, intent(in):: fileID end subroutine printParticle_interface @@ -765,6 +765,8 @@ MODULE moduleMesh REAL(8):: effectiveTime REAL(8):: eThreshold REAL(8):: deltaV + integer:: tally + integer(KIND=OMP_LOCK_KIND):: tally_lock CONTAINS procedure, pass:: apply => ionization @@ -785,63 +787,63 @@ MODULE moduleMesh module subroutine reflection(self, edge, part) use moduleSpecies - class(boundaryReflection), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(boundaryReflection), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine reflection module subroutine absorption(self, edge, part) use moduleSpecies - class(boundaryAbsorption), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(boundaryAbsorption), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine absorption module subroutine transparent(self, edge, part) use moduleSpecies - class(boundaryTransparent), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(boundaryTransparent), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine transparent module subroutine symmetryAxis(self, edge, part) use moduleSpecies - class(boundaryAxis), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(boundaryAxis), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine symmetryAxis module subroutine wallTemperature(self, edge, part) use moduleSpecies - class(boundaryWallTemperature), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(boundaryWallTemperature), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine wallTemperature module subroutine ionization(self, edge, part) use moduleSpecies - class(boundaryIonization), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(boundaryIonization), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine ionization module subroutine quasiNeutrality(self, edge, part) use moduleSpecies - class(boundaryQuasiNeutrality), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(boundaryQuasiNeutrality), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine quasiNeutrality diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index 14a35fd..a940dab 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -1,6 +1,7 @@ !moduleMeshBoundary: Boundary functions for the mesh edges submodule(moduleMesh) boundaryParticle contains + ! Returns the array index of the boundary based on the name module function boundaryParticleName_to_Index(boundaryName) result(bp) use moduleErrors implicit none @@ -25,6 +26,101 @@ submodule(moduleMesh) boundaryParticle end function boundaryParticleName_to_Index + module SUBROUTINE reflection(self, edge, part) + USE moduleCaseParam + USE moduleSpecies + IMPLICIT NONE + + class(boundaryReflection), intent(inout):: self + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + !rp = intersection between particle and edge + !rpp = final position of particle + !vpp = final velocity of particle + REAL(8), DIMENSION(1:3):: rp, vpp + + !Reflect particle velocity + vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal + part%v = vpp + + rp = edge%intersection(part%r) + + part%r = 2.D0*(rp - part%r) + part%r + + !particle is assumed to be inside + part%n_in = .TRUE. + + END SUBROUTINE reflection + + !Absoption in a surface + module SUBROUTINE absorption(self, edge, part) + USE moduleCaseParam + USE moduleSpecies + IMPLICIT NONE + + class(boundaryAbsorption), intent(inout):: self + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + REAL(8):: rpp(1:3) !Position of particle projected to the edge + REAL(8):: d !Distance from particle to edge + + rpp = edge%intersection(part%r) + + d = NORM2(rpp - part%r) + + IF (d >= 0.D0) THEN + part%weight = part%weight/d + + END IF + + !Assign new position to particle + part%r = rpp + !Remove particle from the domain + part%n_in = .FALSE. + + !Scatter particle in associated volume + IF (ASSOCIATED(edge%e1)) THEN + CALL edge%e1%scatter(edge%e1%nNodes, part) + + ELSE + CALL edge%e2%scatter(edge%e2%nNodes, part) + + END IF + + END SUBROUTINE absorption + + !Transparent boundary condition + module SUBROUTINE transparent(self, edge, part) + USE moduleSpecies + IMPLICIT NONE + + class(boundaryTransparent), intent(inout):: self + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + + !Removes particle from domain + part%n_in = .FALSE. + + END SUBROUTINE transparent + + !Symmetry axis. Reflects particles. + !Although this function should never be called, it is set as a reflective boundary + !to properly deal with possible particles reaching a corner and selecting this boundary. + module SUBROUTINE symmetryAxis(self, edge, part) + USE moduleSpecies + IMPLICIT NONE + + class(boundaryAxis), intent(inout):: self + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + + CALL genericReflection(edge, part) + + END SUBROUTINE symmetryAxis + + ! wallTemperature + ! During the collision, the boundary transfers part of the energy to the incident particle + ! Init module SUBROUTINE initWallTemperature(boundary, T, c) USE moduleRefParam IMPLICIT NONE @@ -44,12 +140,37 @@ submodule(moduleMesh) boundaryParticle END SUBROUTINE initWallTemperature + ! Apply + module SUBROUTINE wallTemperature(self, edge, part) + USE moduleSpecies + USE moduleRandom + IMPLICIT NONE + + class(boundaryWallTemperature), intent(inout):: self + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + INTEGER:: i + + !Modifies particle velocity according to wall temperature + DO i = 1, 3 + part%v(i) = part%v(i) + self%vTh*randomMaxwellian() + + END DO + + CALL genericReflection(edge, part) + + END SUBROUTINE wallTemperature + + ! ionization + ! An electron will pass through the surface and create a series of ion-electron pairs based on a neutral background + !init module SUBROUTINE initIonization(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) use moduleRefParam use moduleSpecies use moduleCaseParam use moduleConstParam use moduleErrors, only: criticalError + use omp_lib, only: omp_init_lock implicit none class(boundaryParticleGeneric), allocatable, intent(inout):: boundary @@ -91,147 +212,17 @@ submodule(moduleMesh) boundaryParticle boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2) boundary%deltaV = DSQRT(boundary%eThreshold/mImpact) + boundary%update => ionization_resetTally + boundary%print => ionization_print + + boundary%tally = 0 + call omp_init_lock(boundary%tally_lock) + END SELECT END SUBROUTINE initIonization - module subroutine initQuasiNeutrality(boundary, s_incident) - implicit none - - class(boundaryParticleGeneric), allocatable, intent(inout):: boundary - integer, intent(in):: s_incident - - allocate(boundaryQuasiNeutrality:: boundary) - - select type(boundary) - type is(boundaryQuasiNeutrality) - allocate(boundary%edges(0)) - - boundary%s_incident = s_incident - - - boundary%update => quasiNeutrality_update - boundary%print => quasiNeutrality_print - - end select - - end subroutine initQuasiNeutrality - - module SUBROUTINE reflection(self, edge, part) - USE moduleCaseParam - USE moduleSpecies - IMPLICIT NONE - - class(boundaryReflection), intent(in):: self - CLASS(meshEdge), INTENT(inout):: edge - CLASS(particle), INTENT(inout):: part - !rp = intersection between particle and edge - !rpp = final position of particle - !vpp = final velocity of particle - REAL(8), DIMENSION(1:3):: rp, vpp - - !Reflect particle velocity - vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal - part%v = vpp - - rp = edge%intersection(part%r) - - part%r = 2.D0*(rp - part%r) + part%r - - !particle is assumed to be inside - part%n_in = .TRUE. - - END SUBROUTINE reflection - - !Absoption in a surface - module SUBROUTINE absorption(self, edge, part) - USE moduleCaseParam - USE moduleSpecies - IMPLICIT NONE - - class(boundaryAbsorption), intent(in):: self - CLASS(meshEdge), INTENT(inout):: edge - CLASS(particle), INTENT(inout):: part - REAL(8):: rpp(1:3) !Position of particle projected to the edge - REAL(8):: d !Distance from particle to edge - - rpp = edge%intersection(part%r) - - d = NORM2(rpp - part%r) - - IF (d >= 0.D0) THEN - part%weight = part%weight/d - - END IF - - !Assign new position to particle - part%r = rpp - !Remove particle from the domain - part%n_in = .FALSE. - - !Scatter particle in associated volume - IF (ASSOCIATED(edge%e1)) THEN - CALL edge%e1%scatter(edge%e1%nNodes, part) - - ELSE - CALL edge%e2%scatter(edge%e2%nNodes, part) - - END IF - - END SUBROUTINE absorption - - !Transparent boundary condition - module SUBROUTINE transparent(self, edge, part) - USE moduleSpecies - IMPLICIT NONE - - class(boundaryTransparent), intent(in):: self - CLASS(meshEdge), INTENT(inout):: edge - CLASS(particle), INTENT(inout):: part - - !Removes particle from domain - part%n_in = .FALSE. - - END SUBROUTINE transparent - - !Symmetry axis. Reflects particles. - !Although this function should never be called, it is set as a reflective boundary - !to properly deal with possible particles reaching a corner and selecting this boundary. - module SUBROUTINE symmetryAxis(self, edge, part) - USE moduleSpecies - IMPLICIT NONE - - class(boundaryAxis), intent(in):: self - CLASS(meshEdge), INTENT(inout):: edge - CLASS(particle), INTENT(inout):: part - - CALL genericReflection(edge, part) - - END SUBROUTINE symmetryAxis - - !Wall with temperature - module SUBROUTINE wallTemperature(self, edge, part) - USE moduleSpecies - USE moduleRandom - IMPLICIT NONE - - class(boundaryWallTemperature), intent(in):: self - CLASS(meshEdge), INTENT(inout):: edge - CLASS(particle), INTENT(inout):: part - INTEGER:: i - - !Modifies particle velocity according to wall temperature - DO i = 1, 3 - part%v(i) = part%v(i) + self%vTh*randomMaxwellian() - - END DO - - CALL genericReflection(edge, part) - - END SUBROUTINE wallTemperature - - !Ionization surface: an electron will pass through the surface - ! and create an ion-electron pair based on a neutral background + ! Apply module SUBROUTINE ionization(self, edge, part) USE moduleList USE moduleSpecies @@ -239,9 +230,10 @@ submodule(moduleMesh) boundaryParticle USE moduleRefParam USE moduleRandom USE moduleMath + use omp_lib, only: omp_set_lock, omp_unset_lock IMPLICIT NONE - class(boundaryIonization), intent(in):: self + class(boundaryIonization), intent(inout):: self CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass @@ -320,6 +312,11 @@ submodule(moduleMesh) boundaryParticle !Reduce number of possible ionizations nIonizations = nIonizations - 1 + ! Increase the tally + call omp_set_lock(self%tally_lock) + self%tally = self%tally + 1 + call omp_unset_lock(self%tally_lock) + END IF END DO @@ -329,11 +326,68 @@ submodule(moduleMesh) boundaryParticle END SUBROUTINE ionization + ! Update + subroutine ionization_resetTally(self) + implicit none + + class(boundaryParticleGeneric), intent(inout):: self + + select type(self) + type is(boundaryIonization) + self%tally = 0 + + end select + + end subroutine ionization_resetTally + + ! Print + subroutine ionization_print(self, fileID) + implicit none + + class(boundaryParticleGeneric), intent(inout):: self + integer, intent(in):: fileID + integer:: e + + write(fileID, '(A)') self%name + select type(self) + type is(boundaryIonization) + write(fileID, '(A)') 'Number of successful ionization events per iteration' + write(fileID, '('//fmtInt//')') self%tally + + end select + + end subroutine ionization_print + + ! quasiNeutrality + ! Maintains n_incident = n_rest by modifying the reflection parameter of the edge. + ! Init + module subroutine initQuasiNeutrality(boundary, s_incident) + implicit none + + class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + integer, intent(in):: s_incident + + allocate(boundaryQuasiNeutrality:: boundary) + + select type(boundary) + type is(boundaryQuasiNeutrality) + allocate(boundary%edges(0)) + + boundary%s_incident = s_incident + + boundary%update => quasiNeutrality_update + boundary%print => quasiNeutrality_print + + end select + + end subroutine initQuasiNeutrality + + ! Apply module subroutine quasiNeutrality(self, edge, part) use moduleRandom implicit none - class(boundaryQuasiNeutrality), intent(in):: self + class(boundaryQuasiNeutrality), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part @@ -347,6 +401,7 @@ submodule(moduleMesh) boundaryParticle end subroutine quasiNeutrality + ! Update subroutine quasiNeutrality_update(self) implicit none @@ -404,10 +459,11 @@ submodule(moduleMesh) boundaryParticle end subroutine quasiNeutrality_update + ! Print output subroutine quasiNeutrality_print(self, fileID) implicit none - class(boundaryParticleGeneric), intent(in):: self + class(boundaryParticleGeneric), intent(inout):: self integer, intent(in):: fileID integer:: e @@ -416,7 +472,7 @@ submodule(moduleMesh) boundaryParticle type is(boundaryQuasiNeutrality) write(fileID, '(A,",",A)') '"Edge id"', '"alpha"' do e = 1, size(self%edges) - write(fileID, '(I0,",",ES0.6E3)') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n) + write(fileID, '('//fmtInt//',",",'//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n) end do @@ -425,6 +481,7 @@ submodule(moduleMesh) boundaryParticle end subroutine quasiNeutrality_print ! Generic boundary conditions for internal use + ! reflection module subroutine genericReflection(edge, part) use moduleCaseParam use moduleSpecies @@ -450,6 +507,7 @@ submodule(moduleMesh) boundaryParticle end subroutine genericReflection + ! transparent module subroutine genericTransparent(edge, part) use moduleSpecies implicit none @@ -462,6 +520,7 @@ submodule(moduleMesh) boundaryParticle end subroutine genericTransparent + ! Updates the boundary condition values when %update procedure is associated module subroutine boundariesParticle_update() implicit none @@ -477,6 +536,7 @@ submodule(moduleMesh) boundaryParticle end subroutine boundariesParticle_update + ! Writes the output into the Step_XXXXX_boundaryParticles file when the %print procedure is associated module subroutine boundariesParticle_write() use moduleCaseparam, only: timeStep use moduleOutput, only:fileID_boundaryParticle, formatFileName, informFileCreation