diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 34f2466..6f89660 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -898,10 +898,10 @@ MODULE moduleInput s_incident = speciesName2Index(speciesName) - call initQuasiNeutrality(bound, s_incident) + call quasiNeutrality_init(bound, s_incident) CASE('outflowAdaptive') - ALLOCATE(boundaryOutflowAdaptive:: bound) + call outflowAdaptive_init(bound) CASE DEFAULT CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') @@ -1295,6 +1295,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 +1305,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..b8ae18d 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -682,11 +682,17 @@ MODULE moduleMesh end subroutine initIonization - module subroutine initQuasiNeutrality(boundary, s_incident) + 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 + + module subroutine outflowAdaptive_init(boundary)!, s_incident) + class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + ! integer, intent(in):: s_incident + + end subroutine outflowAdaptive_init module function boundaryParticleName_to_Index(boundaryName) result(bp) character(:), allocatable:: boundaryName @@ -791,8 +797,9 @@ MODULE moduleMesh !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 diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index d5c4f62..011c9b8 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -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,7 +379,7 @@ submodule(moduleMesh) boundaryParticle end select - end subroutine initQuasiNeutrality + end subroutine quasiNeutrality_init ! Apply module subroutine quasiNeutrality(self, edge, part) @@ -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 @@ -489,6 +490,28 @@ submodule(moduleMesh) boundaryParticle ! outflowAdaptive ! Adjust the reflection coefficient of the boundary to maintain a quasi-neutral outflow + ! Init + subroutine outflowAdaptive_init(boundary) + implicit none + + class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + + allocate(boundaryOutflowAdaptive:: boundary) + + select type(boundary) + type is(boundaryOutflowAdaptive) + allocate(boundary%edges(0)) + + boundary%s_incident = 1 !TODO Make this an input parameter + + boundary%update => outflowAdaptive_update + boundary%print => outflowAdaptive_print + + end select + + end subroutine outflowAdaptive_init + + ! Apply module subroutine outflowAdaptive(self, edge, part) use moduleSpecies use moduleRefParam, only: v_ref @@ -497,13 +520,10 @@ 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) - part%v(1) = part%v(1) + v_cut + part%v(1) = part%v(1) - self%velocity_shift(edge%n) if (dot_product(part%v,edge%normal) <= 0.d0) then @@ -513,6 +533,94 @@ submodule(moduleMesh) boundaryParticle end subroutine outflowAdaptive + ! 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 module subroutine genericReflection(edge, part)