diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index c3b66ab..6e8bb61 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 1b4e258..b1dd6dd 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -581,19 +581,6 @@ 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 8b7aaff..34f2466 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 wallTemperature_init(bound, Tw, cw) + CALL initWallTemperature(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 ionization_init(bound, me/m_ref, m0, n0, v0, T0, & + CALL initIonization(bound, me/m_ref, m0, n0, v0, T0, & speciesID, effTime, crossSection, eThreshold, electronSecondaryID) ELSE - CALL ionization_init(bound, me/m_ref, m0, n0, v0, T0, & + CALL initIonization(bound, me/m_ref, m0, n0, v0, T0, & speciesID, effTime, crossSection, eThreshold) END IF @@ -898,15 +898,10 @@ MODULE moduleInput s_incident = speciesName2Index(speciesName) - call quasiNeutrality_init(bound, s_incident) + call initQuasiNeutrality(bound, s_incident) CASE('outflowAdaptive') - 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) + ALLOCATE(boundaryOutflowAdaptive:: bound) CASE DEFAULT CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') @@ -1300,7 +1295,6 @@ 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 @@ -1310,24 +1304,6 @@ 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 7ce7565..71c9bce 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -664,15 +664,13 @@ MODULE moduleMesh end type boundaryParticleGeneric interface - ! Init interfaces - ! wallTemeprature - module subroutine wallTemperature_init(boundary, T, c) + module subroutine initWallTemperature(boundary, T, c) CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary REAL(8), INTENT(in):: T, c !Wall temperature and specific heat - end subroutine wallTemperature_init + end subroutine initWallTemperature - module subroutine ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) + module subroutine initIonization(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 @@ -682,23 +680,14 @@ MODULE moduleMesh real(8), intent(in):: eThreshold integer, optional, intent(in):: electronSecondary - end subroutine ionization_init + end subroutine initIonization - ! quasiNeutrality - module subroutine quasiNeutrality_init(boundary, s_incident) + module subroutine initQuasiNeutrality(boundary, s_incident) class(boundaryParticleGeneric), allocatable, intent(inout):: boundary integer, intent(in):: s_incident - end subroutine quasiNeutrality_init + end subroutine initQuasiNeutrality - ! 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 @@ -708,7 +697,6 @@ 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 @@ -727,7 +715,7 @@ MODULE moduleMesh end subroutine updateParticle_interface - ! Write the values of the particle boundary model + ! Update the values of the particle boundary model subroutine printParticle_interface(self, fileID) import boundaryParticleGeneric @@ -741,28 +729,28 @@ MODULE moduleMesh !Reflecting boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryReflection CONTAINS - procedure, pass:: apply => reflection_apply + procedure, pass:: apply => reflection END TYPE boundaryReflection !Absorption boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAbsorption CONTAINS - procedure, pass:: apply => absorption_apply + procedure, pass:: apply => absorption END TYPE boundaryAbsorption !Transparent boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryTransparent CONTAINS - procedure, pass:: apply => transparent_apply + procedure, pass:: apply => transparent END TYPE boundaryTransparent !Symmetry axis TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAxis CONTAINS - procedure, pass:: apply => symmetryAxis_apply + procedure, pass:: apply => symmetryAxis END TYPE boundaryAxis @@ -771,7 +759,7 @@ MODULE moduleMesh !Thermal velocity of the wall: square root(Wall temperature X specific heat) REAL(8):: vTh CONTAINS - procedure, pass:: apply => wallTemperature_apply + procedure, pass:: apply => wallTemperature END TYPE boundaryWallTemperature @@ -787,7 +775,7 @@ MODULE moduleMesh integer:: tally integer(KIND=OMP_LOCK_KIND):: tally_lock CONTAINS - procedure, pass:: apply => ionization_apply + procedure, pass:: apply => ionization END TYPE boundaryIonization @@ -797,92 +785,91 @@ MODULE moduleMesh integer:: s_incident ! species index of the incident species type(meshEdgePointer), allocatable:: edges(:) !Array with edges contains - procedure, pass:: apply => quasiNeutrality_apply + procedure, pass:: apply => quasiNeutrality end type boundaryQuasiNeutrality !Boundary for quasi-neutral outflow adjusting reflection coefficient type, public, extends(boundaryParticleGeneric):: boundaryOutflowAdaptive - real(8), allocatable:: velocity_shift(:) - integer:: s_incident ! species index of the incident species - type(meshEdgePointer), allocatable:: edges(:) !Array with edges + real(8):: outflowCurrent + real(8):: reflectionFraction contains - procedure, pass:: apply => outflowAdaptive_apply + procedure, pass:: apply => outflowAdaptive end type boundaryOutflowAdaptive interface - module subroutine reflection_apply(self, edge, part) + module subroutine reflection(self, edge, part) use moduleSpecies class(boundaryReflection), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine reflection_apply + end subroutine reflection - module subroutine absorption_apply(self, edge, part) + module subroutine absorption(self, edge, part) use moduleSpecies class(boundaryAbsorption), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine absorption_apply + end subroutine absorption - module subroutine transparent_apply(self, edge, part) + module subroutine transparent(self, edge, part) use moduleSpecies class(boundaryTransparent), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine transparent_apply + end subroutine transparent - module subroutine symmetryAxis_apply(self, edge, part) + module subroutine symmetryAxis(self, edge, part) use moduleSpecies class(boundaryAxis), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine symmetryAxis_apply + end subroutine symmetryAxis - module subroutine wallTemperature_apply(self, edge, part) + module subroutine wallTemperature(self, edge, part) use moduleSpecies class(boundaryWallTemperature), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine wallTemperature_apply + end subroutine wallTemperature - module subroutine ionization_apply(self, edge, part) + module subroutine ionization(self, edge, part) use moduleSpecies class(boundaryIonization), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine ionization_apply + end subroutine ionization - module subroutine quasiNeutrality_apply(self, edge, part) + module subroutine quasiNeutrality(self, edge, part) use moduleSpecies class(boundaryQuasiNeutrality), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine quasiNeutrality_apply + end subroutine quasiNeutrality - module subroutine outflowAdaptive_apply(self, edge, part) + module subroutine outflowAdaptive(self, edge, part) use moduleSpecies class(boundaryOutflowAdaptive), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - end subroutine outflowAdaptive_apply + end subroutine outflowAdaptive ! 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 0280dcf..963116b 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_apply(self, edge, part) + module SUBROUTINE reflection(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_apply + END SUBROUTINE reflection !Absoption in a surface - module SUBROUTINE absorption_apply(self, edge, part) + module SUBROUTINE absorption(self, edge, part) USE moduleCaseParam USE moduleSpecies IMPLICIT NONE @@ -87,10 +87,10 @@ submodule(moduleMesh) boundaryParticle END IF - END SUBROUTINE absorption_apply + END SUBROUTINE absorption !Transparent boundary condition - module SUBROUTINE transparent_apply(self, edge, part) + module SUBROUTINE transparent(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_apply + 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_apply(self, edge, part) + module SUBROUTINE symmetryAxis(self, edge, part) USE moduleSpecies IMPLICIT NONE @@ -116,12 +116,12 @@ submodule(moduleMesh) boundaryParticle CALL genericReflection(edge, part) - END SUBROUTINE symmetryAxis_apply + END SUBROUTINE symmetryAxis ! wallTemperature ! During the collision, the boundary transfers part of the energy to the incident particle ! Init - module SUBROUTINE wallTemperature_init(boundary, T, c) + module SUBROUTINE initWallTemperature(boundary, T, c) USE moduleRefParam IMPLICIT NONE @@ -138,10 +138,10 @@ submodule(moduleMesh) boundaryParticle end select - END SUBROUTINE wallTemperature_init + END SUBROUTINE initWallTemperature ! Apply - module SUBROUTINE wallTemperature_apply(self, edge, part) + module SUBROUTINE wallTemperature(self, edge, part) USE moduleSpecies USE moduleRandom IMPLICIT NONE @@ -159,12 +159,12 @@ submodule(moduleMesh) boundaryParticle CALL genericReflection(edge, part) - END SUBROUTINE wallTemperature_apply + 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 ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) + module SUBROUTINE initIonization(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_update + boundary%update => ionization_resetTally boundary%print => ionization_print boundary%tally = 0 @@ -220,10 +220,10 @@ submodule(moduleMesh) boundaryParticle END SELECT - END SUBROUTINE ionization_init + END SUBROUTINE initIonization ! Apply - module SUBROUTINE ionization_apply(self, edge, part) + module SUBROUTINE ionization(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_apply + END SUBROUTINE ionization ! Update - subroutine ionization_update(self) + subroutine ionization_resetTally(self) implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -338,7 +338,7 @@ submodule(moduleMesh) boundaryParticle end select - end subroutine ionization_update + end subroutine ionization_resetTally ! 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 quasiNeutrality_init(boundary, s_incident) + module subroutine initQuasiNeutrality(boundary, s_incident) implicit none class(boundaryParticleGeneric), allocatable, intent(inout):: boundary @@ -379,10 +379,10 @@ submodule(moduleMesh) boundaryParticle end select - end subroutine quasiNeutrality_init + end subroutine initQuasiNeutrality ! Apply - module subroutine quasiNeutrality_apply(self, edge, part) + module subroutine quasiNeutrality(self, edge, part) use moduleRandom implicit none @@ -398,7 +398,7 @@ submodule(moduleMesh) boundaryParticle end if - end subroutine quasiNeutrality_apply + end subroutine quasiNeutrality ! 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 => self%edges(e)%obj + edge => mesh%edges(e)%obj density_incident = 0.d0 density_rest = 0.d0 @@ -425,8 +425,7 @@ submodule(moduleMesh) boundaryParticle allocate(density_nodes(1:edge%nNodes)) do s = 1, nSpecies do n = 1, edge%nNodes - node => mesh%nodes(nodes(n))%obj - + node => mesh%nodes(n)%obj density_nodes(n) = node%output(s)%den end do @@ -480,7 +479,7 @@ submodule(moduleMesh) boundaryParticle type is(boundaryQuasiNeutrality) write(fileID, '(A,",",A)') '"Edge id"', '"alpha"' do e = 1, size(self%edges) - write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n) + write(fileID, '('//fmtColInt//','//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n) end do @@ -490,30 +489,7 @@ submodule(moduleMesh) boundaryParticle ! outflowAdaptive ! Adjust the reflection coefficient of the boundary to maintain a quasi-neutral outflow - ! 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) + module subroutine outflowAdaptive(self, edge, part) use moduleSpecies use moduleRefParam, only: v_ref implicit none @@ -521,106 +497,25 @@ submodule(moduleMesh) boundaryParticle class(boundaryOutflowAdaptive), intent(inout):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part + real(8):: v_cut - call genericReflection(edge, part) + v_cut = 2.d0 * 40.d3/v_ref !This will be the drag velocity of the ions in the future - part%v(1) = part%v(1) - self%velocity_shift(edge%n) - - if (dot_product(part%v,edge%normal) <= 0.d0) then + 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 * + else call genericTransparent(edge, part) end if - 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 - + end subroutine outflowAdaptive ! Generic boundary conditions for internal use ! reflection