diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index 6e8bb61..c3b66ab 100644 Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index b1dd6dd..1b4e258 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -581,6 +581,19 @@ make \textbf{crossSection}: Character. Complete path to the cross-section data for the ionization process. + + \item \textbf{quasiNeutrality}: Maintains equal density between the species defined by the \textbf{incident} parameter and the others in the case by modifying the reflection coefficient and reflecting incident particles of the species. + The required input is: + + \textbf{incident}: Character. + Name of the incident species that will be reflected. + + \item \textbf{outflowAdaptive}: Changes the velocity of the reflected particles in the species defined by the \textbf{incident} parameter to maintain the mean velocity equal to the average mean velocity of the others in the case. + The required input is: + + \textbf{incident}: Character. + Name of the incident species that will be reflected. + \end{itemize} \end{itemize} \item \textbf{EM}. Array of objects.determines the boundary conditions for the electromagnetic field. diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 34f2466..8b7aaff 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -864,7 +864,7 @@ MODULE moduleInput CALL config%get(object // '.specificHeat', cw, found) IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary') - CALL initWallTemperature(bound, Tw, cw) + CALL wallTemperature_init(bound, Tw, cw) CASE('ionization') !Neutral parameters @@ -882,10 +882,10 @@ MODULE moduleInput CALL config%get(object // '.electronSecondary', electronSecondary, found) IF (found) THEN electronSecondaryID = speciesName2Index(electronSecondary) - CALL initIonization(bound, me/m_ref, m0, n0, v0, T0, & + CALL ionization_init(bound, me/m_ref, m0, n0, v0, T0, & speciesID, effTime, crossSection, eThreshold, electronSecondaryID) ELSE - CALL initIonization(bound, me/m_ref, m0, n0, v0, T0, & + CALL ionization_init(bound, me/m_ref, m0, n0, v0, T0, & speciesID, effTime, crossSection, eThreshold) END IF @@ -898,10 +898,15 @@ MODULE moduleInput s_incident = speciesName2Index(speciesName) - call initQuasiNeutrality(bound, s_incident) + call quasiNeutrality_init(bound, s_incident) CASE('outflowAdaptive') - ALLOCATE(boundaryOutflowAdaptive:: bound) + call config%get(object // '.incident', speciesName, found) + if (.not. found) call criticalError("Incident species name not found for quasiNeutrality boundary model", 'readBoundary') + + s_incident = speciesName2Index(speciesName) + + call outflowAdaptive_init(bound, s_incident) CASE DEFAULT CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') @@ -1295,6 +1300,7 @@ MODULE moduleInput ! If the boundary for the species is linked to the one analysing, add the edges if (associated(physicalSurfaces(ps)%particles(s)%obj, bound)) then bound%edges = [bound%edges, physicalSurfaces(ps)%edges] + end if end do @@ -1304,6 +1310,24 @@ MODULE moduleInput allocate(bound%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here bound%alpha = 0.d0 + type is(boundaryOutflowAdaptive) + ! Loop over all physical surfaces + do ps = 1, nPhysicalSurfaces + ! Loop over all species + do s = 1, nSpecies + ! If the boundary for the species is linked to the one analysing, add the edges + if (associated(physicalSurfaces(ps)%particles(s)%obj, bound)) then + bound%edges = [bound%edges, physicalSurfaces(ps)%edges] + + end if + + end do + + end do + + allocate(bound%velocity_shift(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here + bound%velocity_shift = 0.d0 + end select end do diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 71c9bce..7ce7565 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -664,13 +664,15 @@ MODULE moduleMesh end type boundaryParticleGeneric interface - module subroutine initWallTemperature(boundary, T, c) + ! Init interfaces + ! wallTemeprature + module subroutine wallTemperature_init(boundary, T, c) CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary REAL(8), INTENT(in):: T, c !Wall temperature and specific heat - end subroutine initWallTemperature + end subroutine wallTemperature_init - module subroutine initIonization(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) + module subroutine ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) class(boundaryParticleGeneric), allocatable, intent(inout):: boundary real(8), intent(in):: mImpact real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties @@ -680,14 +682,23 @@ MODULE moduleMesh real(8), intent(in):: eThreshold integer, optional, intent(in):: electronSecondary - end subroutine initIonization + end subroutine ionization_init - module subroutine initQuasiNeutrality(boundary, s_incident) + ! quasiNeutrality + module subroutine quasiNeutrality_init(boundary, s_incident) class(boundaryParticleGeneric), allocatable, intent(inout):: boundary integer, intent(in):: s_incident - end subroutine initQuasiNeutrality + end subroutine quasiNeutrality_init + ! outflowAdaptive + module subroutine outflowAdaptive_init(boundary, s_incident) + class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + integer, intent(in):: s_incident + + end subroutine outflowAdaptive_init + + ! Get the index of the species from its name module function boundaryParticleName_to_Index(boundaryName) result(bp) character(:), allocatable:: boundaryName integer:: bp @@ -697,6 +708,7 @@ MODULE moduleMesh end interface abstract interface + ! Apply the effects of the boundary to the particle subroutine applyParticle_interface(self, edge, part) use moduleSpecies import boundaryParticleGeneric, meshEdge @@ -715,7 +727,7 @@ MODULE moduleMesh end subroutine updateParticle_interface - ! Update the values of the particle boundary model + ! Write the values of the particle boundary model subroutine printParticle_interface(self, fileID) import boundaryParticleGeneric @@ -729,28 +741,28 @@ MODULE moduleMesh !Reflecting boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryReflection CONTAINS - procedure, pass:: apply => reflection + procedure, pass:: apply => reflection_apply END TYPE boundaryReflection !Absorption boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAbsorption CONTAINS - procedure, pass:: apply => absorption + procedure, pass:: apply => absorption_apply END TYPE boundaryAbsorption !Transparent boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryTransparent CONTAINS - procedure, pass:: apply => transparent + procedure, pass:: apply => transparent_apply END TYPE boundaryTransparent !Symmetry axis TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAxis CONTAINS - procedure, pass:: apply => symmetryAxis + procedure, pass:: apply => symmetryAxis_apply END TYPE boundaryAxis @@ -759,7 +771,7 @@ MODULE moduleMesh !Thermal velocity of the wall: square root(Wall temperature X specific heat) REAL(8):: vTh CONTAINS - procedure, pass:: apply => wallTemperature + procedure, pass:: apply => wallTemperature_apply END TYPE boundaryWallTemperature @@ -775,7 +787,7 @@ MODULE moduleMesh integer:: tally integer(KIND=OMP_LOCK_KIND):: tally_lock CONTAINS - procedure, pass:: apply => ionization + procedure, pass:: apply => ionization_apply END TYPE boundaryIonization @@ -785,91 +797,92 @@ MODULE moduleMesh integer:: s_incident ! species index of the incident species type(meshEdgePointer), allocatable:: edges(:) !Array with edges contains - procedure, pass:: apply => quasiNeutrality + procedure, pass:: apply => quasiNeutrality_apply end type boundaryQuasiNeutrality !Boundary for quasi-neutral outflow adjusting reflection coefficient type, public, extends(boundaryParticleGeneric):: boundaryOutflowAdaptive - real(8):: outflowCurrent - real(8):: reflectionFraction + real(8), allocatable:: velocity_shift(:) + integer:: s_incident ! species index of the incident species + type(meshEdgePointer), allocatable:: edges(:) !Array with edges contains - procedure, pass:: apply => outflowAdaptive + procedure, pass:: apply => outflowAdaptive_apply end type boundaryOutflowAdaptive interface - module subroutine reflection(self, edge, part) + module subroutine reflection_apply(self, edge, part) use moduleSpecies class(boundaryReflection), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine reflection + end subroutine reflection_apply - module subroutine absorption(self, edge, part) + module subroutine absorption_apply(self, edge, part) use moduleSpecies class(boundaryAbsorption), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine absorption + end subroutine absorption_apply - module subroutine transparent(self, edge, part) + module subroutine transparent_apply(self, edge, part) use moduleSpecies class(boundaryTransparent), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine transparent + end subroutine transparent_apply - module subroutine symmetryAxis(self, edge, part) + module subroutine symmetryAxis_apply(self, edge, part) use moduleSpecies class(boundaryAxis), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine symmetryAxis + end subroutine symmetryAxis_apply - module subroutine wallTemperature(self, edge, part) + module subroutine wallTemperature_apply(self, edge, part) use moduleSpecies class(boundaryWallTemperature), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine wallTemperature + end subroutine wallTemperature_apply - module subroutine ionization(self, edge, part) + module subroutine ionization_apply(self, edge, part) use moduleSpecies class(boundaryIonization), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine ionization + end subroutine ionization_apply - module subroutine quasiNeutrality(self, edge, part) + module subroutine quasiNeutrality_apply(self, edge, part) use moduleSpecies class(boundaryQuasiNeutrality), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine quasiNeutrality + end subroutine quasiNeutrality_apply - module subroutine outflowAdaptive(self, edge, part) + module subroutine outflowAdaptive_apply(self, edge, part) use moduleSpecies class(boundaryOutflowAdaptive), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine outflowAdaptive + end subroutine outflowAdaptive_apply ! Generic basic boundary conditions to use internally in the code module subroutine genericReflection(edge, part) diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index 963116b..0280dcf 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -26,7 +26,7 @@ submodule(moduleMesh) boundaryParticle end function boundaryParticleName_to_Index - module SUBROUTINE reflection(self, edge, part) + module SUBROUTINE reflection_apply(self, edge, part) USE moduleCaseParam USE moduleSpecies IMPLICIT NONE @@ -50,10 +50,10 @@ submodule(moduleMesh) boundaryParticle !particle is assumed to be inside part%n_in = .TRUE. - END SUBROUTINE reflection + END SUBROUTINE reflection_apply !Absoption in a surface - module SUBROUTINE absorption(self, edge, part) + module SUBROUTINE absorption_apply(self, edge, part) USE moduleCaseParam USE moduleSpecies IMPLICIT NONE @@ -87,10 +87,10 @@ submodule(moduleMesh) boundaryParticle END IF - END SUBROUTINE absorption + END SUBROUTINE absorption_apply !Transparent boundary condition - module SUBROUTINE transparent(self, edge, part) + module SUBROUTINE transparent_apply(self, edge, part) USE moduleSpecies IMPLICIT NONE @@ -101,12 +101,12 @@ submodule(moduleMesh) boundaryParticle !Removes particle from domain part%n_in = .FALSE. - END SUBROUTINE transparent + END SUBROUTINE transparent_apply !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) + module SUBROUTINE symmetryAxis_apply(self, edge, part) USE moduleSpecies IMPLICIT NONE @@ -116,12 +116,12 @@ submodule(moduleMesh) boundaryParticle CALL genericReflection(edge, part) - END SUBROUTINE symmetryAxis + END SUBROUTINE symmetryAxis_apply ! wallTemperature ! During the collision, the boundary transfers part of the energy to the incident particle ! Init - module SUBROUTINE initWallTemperature(boundary, T, c) + module SUBROUTINE wallTemperature_init(boundary, T, c) USE moduleRefParam IMPLICIT NONE @@ -138,10 +138,10 @@ submodule(moduleMesh) boundaryParticle end select - END SUBROUTINE initWallTemperature + END SUBROUTINE wallTemperature_init ! Apply - module SUBROUTINE wallTemperature(self, edge, part) + module SUBROUTINE wallTemperature_apply(self, edge, part) USE moduleSpecies USE moduleRandom IMPLICIT NONE @@ -159,12 +159,12 @@ submodule(moduleMesh) boundaryParticle CALL genericReflection(edge, part) - END SUBROUTINE wallTemperature + END SUBROUTINE wallTemperature_apply ! 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) + module SUBROUTINE ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) use moduleRefParam use moduleSpecies use moduleCaseParam @@ -212,7 +212,7 @@ submodule(moduleMesh) boundaryParticle boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2) boundary%deltaV = DSQRT(boundary%eThreshold/mImpact) - boundary%update => ionization_resetTally + boundary%update => ionization_update boundary%print => ionization_print boundary%tally = 0 @@ -220,10 +220,10 @@ submodule(moduleMesh) boundaryParticle END SELECT - END SUBROUTINE initIonization + END SUBROUTINE ionization_init ! Apply - module SUBROUTINE ionization(self, edge, part) + module SUBROUTINE ionization_apply(self, edge, part) USE moduleList USE moduleSpecies USE moduleMesh @@ -324,10 +324,10 @@ submodule(moduleMesh) boundaryParticle !Removes ionizing electron regardless the number of pair created part%n_in = .FALSE. - END SUBROUTINE ionization + END SUBROUTINE ionization_apply ! Update - subroutine ionization_resetTally(self) + subroutine ionization_update(self) implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -338,7 +338,7 @@ submodule(moduleMesh) boundaryParticle end select - end subroutine ionization_resetTally + end subroutine ionization_update ! Print subroutine ionization_print(self, fileID) @@ -360,7 +360,7 @@ submodule(moduleMesh) boundaryParticle ! quasiNeutrality ! Maintains n_incident = n_rest by modifying the reflection parameter of the edge. ! Init - module subroutine initQuasiNeutrality(boundary, s_incident) + module subroutine quasiNeutrality_init(boundary, s_incident) implicit none class(boundaryParticleGeneric), allocatable, intent(inout):: boundary @@ -379,10 +379,10 @@ submodule(moduleMesh) boundaryParticle end select - end subroutine initQuasiNeutrality + end subroutine quasiNeutrality_init ! Apply - module subroutine quasiNeutrality(self, edge, part) + module subroutine quasiNeutrality_apply(self, edge, part) use moduleRandom implicit none @@ -398,7 +398,7 @@ submodule(moduleMesh) boundaryParticle end if - end subroutine quasiNeutrality + end subroutine quasiNeutrality_apply ! Update subroutine quasiNeutrality_update(self) @@ -416,7 +416,7 @@ submodule(moduleMesh) boundaryParticle select type(self) type is(boundaryQuasiNeutrality) do e = 1, size(self%edges) - edge => mesh%edges(e)%obj + edge => self%edges(e)%obj density_incident = 0.d0 density_rest = 0.d0 @@ -425,7 +425,8 @@ submodule(moduleMesh) boundaryParticle allocate(density_nodes(1:edge%nNodes)) do s = 1, nSpecies do n = 1, edge%nNodes - node => mesh%nodes(n)%obj + node => mesh%nodes(nodes(n))%obj + density_nodes(n) = node%output(s)%den end do @@ -479,7 +480,7 @@ submodule(moduleMesh) boundaryParticle type is(boundaryQuasiNeutrality) write(fileID, '(A,",",A)') '"Edge id"', '"alpha"' do e = 1, size(self%edges) - write(fileID, '('//fmtColInt//','//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n) + write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n) end do @@ -489,7 +490,30 @@ submodule(moduleMesh) boundaryParticle ! outflowAdaptive ! Adjust the reflection coefficient of the boundary to maintain a quasi-neutral outflow - module subroutine outflowAdaptive(self, edge, part) + ! Init + subroutine outflowAdaptive_init(boundary, s_incident) + implicit none + + class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + integer, intent(in):: s_incident + + allocate(boundaryOutflowAdaptive:: boundary) + + select type(boundary) + type is(boundaryOutflowAdaptive) + allocate(boundary%edges(0)) + + boundary%s_incident = s_incident + + boundary%update => outflowAdaptive_update + boundary%print => outflowAdaptive_print + + end select + + end subroutine outflowAdaptive_init + + ! Apply + module subroutine outflowAdaptive_apply(self, edge, part) use moduleSpecies use moduleRefParam, only: v_ref implicit none @@ -497,25 +521,106 @@ submodule(moduleMesh) boundaryParticle class(boundaryOutflowAdaptive), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - real(8):: v_cut - v_cut = 2.d0 * 40.d3/v_ref !This will be the drag velocity of the ions in the future + call genericReflection(edge, part) - if (dot_product(part%v,-edge%normal) > v_cut) then - ! print *,part%v(1), v_cut - ! part%v = part%v + v_cut*edge%normal - part%v(1) = part%v(1) - v_cut - ! print *,part%v(1) - call genericReflection(edge, part) - ! print *,part%v(1) - ! print * + part%v(1) = part%v(1) - self%velocity_shift(edge%n) + + if (dot_product(part%v,edge%normal) <= 0.d0) then - else call genericTransparent(edge, part) end if - end subroutine outflowAdaptive + end subroutine outflowAdaptive_apply + + ! Update + subroutine outflowAdaptive_update(self) + implicit none + + class(boundaryParticleGeneric), intent(inout):: self + integer:: e, n, s + integer, allocatable:: nodes(:) + class(meshEdge), pointer:: edge + real(8), allocatable:: den_nodes(:) + real(8), allocatable:: mom_nodes(:) + class(meshNode), pointer:: node + real(8):: den_center, mom_center, vel_center + real(8):: velocity_incident, velocity_rest + real(8):: alpha + + select type(self) + type is(boundaryOutflowAdaptive) + do e = 1, size(self%edges) + edge => self%edges(e)%obj + + velocity_incident = 0.d0 + velocity_rest = 0.d0 + + nodes = edge%getNodes(edge%nNodes) + allocate(den_nodes(1:edge%nNodes)) + allocate(mom_nodes(1:edge%nNodes)) + do s = 1, nSpecies + do n = 1, edge%nNodes + node => mesh%nodes(nodes(n))%obj + + den_nodes(n) = node%output(s)%den + mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal) + + end do + + den_center = edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes) + mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) + + if (den_center > 1.0d-10) then + vel_center = mom_center / den_center + + else + vel_center = 0.0d0 + end if + + if (s == self%s_incident) then + velocity_incident = vel_center + else + velocity_rest = velocity_rest + vel_center + end if + + end do + + alpha = velocity_rest - velocity_incident + + ! Apply correction with a factor of 0.1 to avoid fast changes + self%velocity_shift(edge%n) = self%velocity_shift(edge%n) + 1.0d-2 * alpha + + deallocate(den_nodes) + deallocate(mom_nodes) + + end do + + end select + + end subroutine outflowAdaptive_update + + ! Print + subroutine outflowAdaptive_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(boundaryOutflowAdaptive) + write(fileID, '(A,",",A)') '"Edge id"', '"vel shift non-dimensional"' + do e = 1, size(self%edges) + write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%velocity_shift(self%edges(e)%obj%n) + end do + + end select + + end subroutine outflowAdaptive_print + ! Generic boundary conditions for internal use ! reflection