! Copyright (c) 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2019 Jana Gericke <jana.gericke@uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 2021-2022 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ************************************************************************** ! !> author: Kannan Masilamani !! Average Interpolation of flow quantities between different grid levels !! ! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de> ! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2014-2015, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, this ! list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE ! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE ! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE ! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR ! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, ! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Make sure loglvl is defined to an integer value. ! Usually this should be defined on the command line with -Dloglvl= ! to some value. module mus_interpolate_average_module use iso_c_binding, only: c_loc, c_ptr, c_f_pointer ! include treelm modules use env_module, only: rk use tem_aux_module, only: tem_abort use tem_element_module, only: eT_GhostFromCoarser, & & eT_ghostFromFiner use tem_param_module, only: cs2inv, cs2, PI, div1_2, div1_9, div4_9,& & div1_36, div2_3, div2_9, div1_18, rho0, rho0Inv use tem_comm_env_module, only: tem_comm_env_type use tem_debug_module, only: dbgUnit use tem_construction_module, only: tem_levelDesc_type use tem_matrix_module, only: tem_matrix_type use tem_stencil_module, only: tem_stencilHeader_type use tem_logging_module, only: logUnit use tem_varSys_module, only: tem_varSys_type use tem_time_module, only: tem_time_type ! include musubi modules use mus_interpolate_debug_module, only: TGV_2D use mus_scheme_layout_module, only: mus_scheme_layout_type use mus_physics_module, only: mus_physics_type use mus_interpolate_header_module, only: mus_interpolation_method_type use mus_field_prop_module, only: mus_field_prop_type use mus_fluid_module, only: mus_fluid_type use mus_relaxationParam_module, only: mus_calcOmegaFromVisc use mus_derVarPos_module, only: mus_derVarPos_type use mus_derivedQuantities_module2, only: getNonEqFac_intp_fine_to_coarse, & & getNonEqFac_intp_coarse_to_fine use mus_varSys_module, only: mus_varSys_data_type implicit none private public :: fillArbiMyGhostsFromFiner_avg public :: fillMyGhostsFromFiner_avg_feq_fneq public :: fillMyGhostsFromFiner_avgLES_feq_fneq public :: fillMyGhostsFromFiner_avgIncomp_feq_fneq public :: fillMyGhostsFromFiner_avgIncompLES_feq_fneq public :: fillMyGhostsFromFiner_avg2D_feq_fneq public :: fillMyGhostsFromFiner_avg2DIncomp_feq_fneq public :: fillArbiFinerGhostsFromMe_weighAvg public :: fillFinerGhostsFromMe_weighAvg_feq_fneq public :: fillFinerGhostsFromMe_weighAvgLES_feq_fneq public :: fillFinerGhostsFromMe_weighAvgIncomp_feq_fneq public :: fillFinerGhostsFromMe_weighAvgIncompLES_feq_fneq public :: fillFinerGhostsFromMe_weighAvg2D_feq_fneq public :: fillFinerGhostsFromMe_weighAvg2DIncomp_feq_fneq contains ! **************************************************************************** ! !> Interpolate auxiliary field from fine source to coarse target !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine_arbitraryVal]] in intp/[[mus_interpolate_header_module]].f90 !! in order to be callable via !! [[mus_interpolation_method_type:do_intpArbiVal]] function pointer. subroutine fillArbiMyGhostsFromFiner_avg( method, tLevelDesc, level, & & stencil, sVal, tVal, nTargets, & & targetList, nScalars ) ! ------------------------------------------------------------------ ! class(mus_interpolation_method_type), intent(inout) :: method !> my refinement level integer, intent(in) :: level !> stencil header type(tem_stencilHeader_type), intent(in) :: stencil !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sVal(:) !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tVal(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> number of scalars to interpolate integer, intent(in) :: nScalars ! ------------------------------------------------------------------ ! ! -------------------------------------------------------------------------- integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: inv_nSourceElems real(kind=rk) :: tArbi(nScalars) ! target auxField real(kind=rk) :: sArbi(nScalars) ! temp source ArbiField ! -------------------------------------------------------------------------- sourceLevel = level + 1 targetLevel = level ! Treat all coarse target elements do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner) ! Get how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromFiner( iElem )%elem%nVals inv_nSourceElems = 1.0_rk / dble(nSourceElems) tArbi = 0.0_rk ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element sourceElem = tLevelDesc%depFromFiner( iElem ) & & %elem%val( iSourceElem ) ! Get souce auxilary variables sArbi(1:nScalars) = sVal((sourceElem-1)*nScalars+1:sourceElem*nScalars) ! target aux tArbi = sArbi + tArbi end do ! iSourceElem ! interpolate all auxiliary variables by average tVal((targetElem-1)*nScalars+1: targetElem*nScalars) & & = tArbi(:) * inv_nSourceElems enddo end subroutine fillArbiMyGhostsFromFiner_avg ! **************************************************************************** ! ! **************************************************************************** ! !> [Average interpolation](../page/features/intp_methods.html) of ghostFromFiner !! The interpolation procedure used in this routine is:\n !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillMyGhostsFromFiner_avg_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f( layout%fStencil%QQ ), f_eq( layout%fStencil%QQ ) integer :: QQ real(kind=rk) :: inv_nSourceElems real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation type(mus_fluid_type), pointer :: fluid real(kind=rk) :: cVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level + 1 targetLevel = level dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all coarse target elements do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner) ! Get how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromFiner( iElem )%elem%nVals inv_nSourceElems = 1.0_rk / dble(nSourceElems) t_f_eq = 0._rk t_f_neq = 0._rk ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element sourceElem = tLevelDesc%depFromFiner( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) vel(3) = sAuxField(elemOff + vel_pos(3)) do iDir = 1, QQ ! calculate equilibrium density f_eq(idir) = layout%weight( idir ) & & * rho & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) t_f_eq(iDir) = t_f_eq(iDir) + f_eq(iDir) t_f_neq(iDir) = t_f_neq(iDir) + (f(iDir) - f_eq(iDir)) enddo end do ! iSourceElem ! get normalized kinematic viscosity on target (coarse) element cVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_f = 2 v^s_c ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(2_rk*cVisc) cOmegaKine = mus_calcOmegaFromVisc(cVisc) !evaluate scaling factor fac = getNonEqFac_intp_fine_to_coarse( cOmegaKine, fOmegaKine) ! interpolate all pdfs by average and scale accordingly do iDir = 1, QQ t_f_eq(iDir) = t_f_eq(iDir) * inv_nSourceElems t_f_neq(iDir) = t_f_neq(iDir) * inv_nSourceElems * fac f(iDir) = t_f_eq(iDir) + t_f_neq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = & & f(iDir) enddo enddo end subroutine fillMyGhostsFromFiner_avg_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> [Average interpolation](../page/features/intp_methods.html) of ghostFromFiner !! The interpolation procedure used in this routine is:\n !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillMyGhostsFromFiner_avgLES_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f( layout%fStencil%QQ ), f_eq( layout%fStencil%QQ ) integer :: QQ real(kind=rk) :: inv_nSourceElems real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation type(mus_fluid_type), pointer :: fluid real(kind=rk) :: cVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac real(kind=rk) :: tTurbVisc integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level + 1 targetLevel = level dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all coarse target elements do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner) ! Get how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromFiner( iElem )%elem%nVals inv_nSourceElems = 1.0_rk / dble(nSourceElems) t_f_eq = 0._rk t_f_neq = 0._rk tTurbVisc = 0.0_rk ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element sourceElem = tLevelDesc%depFromFiner( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) vel(3) = sAuxField(elemOff + vel_pos(3)) do iDir = 1, QQ ! calculate equilibrium density f_eq(idir) = layout%weight( idir ) & & * rho & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) t_f_eq(iDir) = t_f_eq(iDir) + f_eq(iDir) t_f_neq(iDir) = t_f_neq(iDir) + (f(iDir) - f_eq(iDir)) enddo tTurbVisc = tTurbVisc + fluid%turbulence%dataOnLvl(sourceLevel) & & %visc(sourceElem) end do ! iSourceElem ! interpolate turbulent viscosity on target element tTurbVisc = tTurbVisc * inv_nSourceElems ! scale interpolated turbulent viscosity to target element fluid%turbulence%dataOnLvl(targetLevel)%visc(targetElem) & & = fluid%turbulence%fac_f2c*tTurbVisc ! get normalized kinematic viscosity on target (coarse) element cVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_f = 2 v^s_c ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(2_rk*cVisc + tTurbVisc) cOmegaKine = mus_calcOmegaFromVisc(cVisc + fluid%turbulence%fac_f2c*tTurbVisc) !evaluate scaling factor fac = getNonEqFac_intp_fine_to_coarse(cOmegaKine, fOmegaKine) ! interpolate all pdfs by average and scale accordingly do iDir = 1, QQ t_f_eq(iDir) = t_f_eq(iDir) * inv_nSourceElems t_f_neq(iDir) = t_f_neq(iDir) * inv_nSourceElems * fac f(iDir) = t_f_eq(iDir) + t_f_neq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo enddo end subroutine fillMyGhostsFromFiner_avgLES_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> [Average interpolation](../page/features/intp_methods.html) of ghostFromFiner !! The interpolation procedure used in this routine is:\n !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillMyGhostsFromFiner_avgIncomp_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f( layout%fStencil%QQ ), f_eq( layout%fStencil%QQ ) integer :: QQ real(kind=rk) :: inv_nSourceElems real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation type(mus_fluid_type), pointer :: fluid real(kind=rk) :: cVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level + 1 targetLevel = level dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all coarse target elements do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner) ! Get how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromFiner( iElem )%elem%nVals inv_nSourceElems = 1.0_rk / dble(nSourceElems) t_f_eq = 0._rk t_f_neq = 0._rk ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element sourceElem = tLevelDesc%depFromFiner( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) vel(3) = sAuxField(elemOff + vel_pos(3)) do iDir = 1, QQ f_eq(idir) = layout%weight( idir ) & & * ( rho + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:) * vel(:)) * 0.5_rk * cs2inv ) ) t_f_eq(iDir) = t_f_eq(iDir) + f_eq(iDir) t_f_neq(iDir) = t_f_neq(iDir) + (f(iDir) - f_eq(iDir)) enddo end do ! iSourceElem ! get normalized kinematic viscosity on target (coarse) element cVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_f = 2 v^s_c ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(2_rk*cVisc) cOmegaKine = mus_calcOmegaFromVisc(cVisc) !evaluate scaling factor fac = getNonEqFac_intp_fine_to_coarse( cOmegaKine, fOmegaKine) ! interpolate all pdfs by average and scale accordingly do iDir = 1, QQ t_f_eq(iDir) = t_f_eq(iDir) * inv_nSourceElems t_f_neq(iDir) = t_f_neq(iDir) * inv_nSourceElems * fac f(iDir) = t_f_eq(iDir) + t_f_neq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo enddo end subroutine fillMyGhostsFromFiner_avgIncomp_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> [Average interpolation](../page/features/intp_methods.html) of ghostFromFiner !! The interpolation procedure used in this routine is:\n !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillMyGhostsFromFiner_avgIncompLES_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f( layout%fStencil%QQ ), f_eq( layout%fStencil%QQ ) integer :: QQ real(kind=rk) :: inv_nSourceElems real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation type(mus_fluid_type), pointer :: fluid real(kind=rk) :: cVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac real(kind=rk) :: tTurbVisc integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level + 1 targetLevel = level dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all coarse target elements do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner) ! Get how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromFiner( iElem )%elem%nVals inv_nSourceElems = 1.0_rk / dble(nSourceElems) t_f_eq = 0._rk t_f_neq = 0._rk tTurbVisc = 0.0_rk ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element sourceElem = tLevelDesc%depFromFiner( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) vel(3) = sAuxField(elemOff + vel_pos(3)) do iDir = 1, QQ f_eq(idir) = layout%weight( idir ) & & * ( rho + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:) * vel(:)) * 0.5_rk * cs2inv ) ) t_f_eq(iDir) = t_f_eq(iDir) + f_eq(iDir) t_f_neq(iDir) = t_f_neq(iDir) + (f(iDir) - f_eq(iDir)) enddo !interpolate turbulent viscosity to target element tTurbVisc = tTurbVisc + fluid%turbulence%dataOnLvl(sourceLevel) & & %visc(sourceElem) end do ! iSourceElem ! interpolate turbulent viscosity on target element tTurbVisc = tTurbVisc * inv_nSourceElems ! scale interpolated turbulent viscosity to target element fluid%turbulence%dataOnLvl(targetLevel)%visc(targetElem) & & = fluid%turbulence%fac_f2c*tTurbVisc ! get normalized kinematic viscosity on target (coarse) element cVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_f = 2 v^s_c ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(2_rk*cVisc + tTurbVisc) cOmegaKine = mus_calcOmegaFromVisc(cVisc + fluid%turbulence%fac_f2c*tTurbVisc) !evaluate scaling factor fac = getNonEqFac_intp_fine_to_coarse( cOmegaKine, fOmegaKine) ! interpolate all pdfs by average and scale accordingly do iDir = 1, QQ t_f_eq(iDir) = t_f_eq(iDir) * inv_nSourceElems t_f_neq(iDir) = t_f_neq(iDir) * inv_nSourceElems * fac f(iDir) = t_f_eq(iDir) + t_f_neq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo enddo end subroutine fillMyGhostsFromFiner_avgIncompLES_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> Fill coarse target ghost from fine source fluid by average interpolation. !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillMyGhostsFromFiner_avg2D_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f( layout%fStencil%QQ ), f_eq( layout%fStencil%QQ ) integer :: QQ real(kind=rk) :: inv_nSourceElems ! source elements' pdf real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation type(mus_fluid_type), pointer :: fluid real(kind=rk) :: cVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- ! for initialization do something else vel = 0._rk fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level + 1 targetLevel = level dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all coarse target elements do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element treeId targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner) ! Find out how many fine source elements we have for interpolation. nSourceELems = tLevelDesc%depFromFiner( iElem )%elem%nVals inv_nSourceElems = 1.0_rk / dble(nSourceElems) ! Now loop over all fine source elements for this target: t_f_eq = 0._rk t_f_neq = 0._rk do iSourceElem = 1, nSourceElems ! Get the source element position sourceElem = tLevelDesc%depFromFiner( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-velocity vel(1) = sum( f * layout%fStencil%cxDirRK(1, :) ) / rho vel(2) = sum( f * layout%fStencil%cxDirRK(2, :) ) / rho do iDir = 1, QQ ! calculate equilibrium density f_eq(idir) = layout%weight( idir ) & & * rho & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) t_f_eq(iDir) = t_f_eq(iDir) + f_eq(iDir) t_f_neq(iDir) = t_f_neq(iDir) + (f(iDir) - f_eq(iDir)) enddo end do ! iSourceElem ! get normalized kinematic viscosity on target (coarse) element cVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_f = 2 v^s_c ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(2_rk*cVisc) cOmegaKine = mus_calcOmegaFromVisc(cVisc) !evaluate scaling factor fac = getNonEqFac_intp_fine_to_coarse( cOmegaKine, fOmegaKine) ! interpolate all pdfs by average and scale accordingly do iDir = 1, QQ t_f_eq(iDir) = t_f_eq(iDir) * inv_nSourceElems t_f_neq(iDir) = t_f_neq(iDir) * inv_nSourceElems * fac f(iDir) = t_f_eq(iDir) + t_f_neq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo end do end subroutine fillMyGhostsFromFiner_avg2D_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> Fill coarse target ghost from fine source fluid by average interpolation. !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillMyGhostsFromFiner_avg2DIncomp_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f( layout%fStencil%QQ ), f_eq( layout%fStencil%QQ ) integer :: QQ real(kind=rk) :: inv_nSourceElems ! source elements' pdf real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation type(mus_fluid_type), pointer :: fluid real(kind=rk) :: cVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- ! for initialization do something else vel = 0._rk fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level + 1 targetLevel = level dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all coarse target elements do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element treeId targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner) ! Find out how many fine source elements we have for interpolation. nSourceELems = tLevelDesc%depFromFiner( iElem )%elem%nVals inv_nSourceElems = 1.0_rk / dble(nSourceElems) ! Now loop over all fine source elements for this target: t_f_eq = 0._rk t_f_neq = 0._rk do iSourceElem = 1, nSourceElems ! Get the source element position sourceElem = tLevelDesc%depFromFiner( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) do iDir = 1, QQ f_eq(idir) = layout%weight( idir ) & & * ( rho + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:) * vel(:)) * 0.5_rk * cs2inv ) ) t_f_eq(iDir) = t_f_eq(iDir) + f_eq(iDir) t_f_neq(iDir) = t_f_neq(iDir) + (f(iDir) - f_eq(iDir)) enddo end do ! iSourceElem ! get normalized kinematic viscosity on target (coarse) element cVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_f = 2 v^s_c ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(2_rk*cVisc) cOmegaKine = mus_calcOmegaFromVisc(cVisc) !evaluate scaling factor fac = getNonEqFac_intp_fine_to_coarse( cOmegaKine, fOmegaKine) ! interpolate all pdfs by average and scale accordingly do iDir = 1, QQ t_f_eq(iDir) = t_f_eq(iDir) * inv_nSourceElems t_f_neq(iDir) = t_f_neq(iDir) * inv_nSourceElems * fac f(iDir) = t_f_eq(iDir) + t_f_neq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo end do end subroutine fillMyGhostsFromFiner_avg2Dincomp_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> Interpolate auxiliary field from coarse source to fine target !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine_arbitraryVal]] in intp/[[mus_interpolate_header_module]].f90 !! in order to be callable via !! [[mus_interpolation_method_type:do_intpArbiVal]] function pointer. subroutine fillArbiFinerGhostsFromMe_weighAvg( method, tLevelDesc, level, & & stencil, sVal, tVal, nTargets,& & targetList, nScalars ) ! ------------------------------------------------------------------ ! class(mus_interpolation_method_type), intent(inout) :: method !> my refinement level integer, intent(in) :: level !> stencil header type(tem_stencilHeader_type), intent(in) :: stencil !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sVal(:) !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tVal(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> number of scalars to interpolate integer, intent(in) :: nScalars ! ------------------------------------------------------------------ ! ! -------------------------------------------------------------------------- integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: weight( stencil%QQ ) real(kind=rk) :: sArbi(nScalars, stencil%QQ) ! temp source ArbiField integer :: iVar ! -------------------------------------------------------------------------- sourceLevel = level targetLevel = level + 1 ! Treat all coarse target elements do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser) ! Get how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals ! Read the pre-calculated weights (like in Average intp.) weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) & & %weight(1:nSourceElems) ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element sourceElem = tLevelDesc%depFromCoarser( iElem ) & & %elem%val( iSourceElem ) ! Get souce auxilary variables sArbi(:, iSourceElem) = sVal( (sourceElem-1)*nScalars+1 & & : sourceElem*nScalars ) end do ! iSourceElem ! interpolate all auxiliary variables by average do iVar = 1, nScalars tVal((targetElem-1)*nScalars+iVar) & & = sum( weight(1:nSourceElems)*sArbi(iVar,1:nSourceElems) ) end do enddo end subroutine fillArbiFinerGhostsFromMe_weighAvg ! **************************************************************************** ! ! **************************************************************************** ! !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillFinerGhostsFromMe_weighAvg_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f_neq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: f_eq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: weight(layout%fStencil%QQ) real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation ! pdf to reconstruct from real(kind=rk) :: f( layout%fStencil%QQ ) integer :: QQ type(mus_fluid_type), pointer :: fluid real(kind=rk) :: fVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level targetLevel = level + 1 dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all fine target elements: do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element treeId targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser) ! Find out how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals ! Read the pre-calculated weights (like in Average intp.) weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) & & %weight(1:nSourceElems) ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element position sourceElem = tLevelDesc%depFromCoarser( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) vel(3) = sAuxField(elemOff + vel_pos(3)) do iDir = 1, QQ ! calculate equilibrium density f_eq(idir,isourceelem) = layout%weight( idir ) & & * rho & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) f_neq(iDir, iSourceElem) = f(iDir) - f_eq(iDir, iSourceElem) enddo end do ! iSourceElem ! interpolate all pdfs by weighted average do iDir = 1,QQ t_f_eq(iDir) = sum(weight(1:nSourceElems)*f_eq(iDir,1:nSourceElems)) t_f_neq(iDir) = sum(weight(1:nSourceElems)*f_neq(iDir,1:nSourceElems)) end do ! get normalized kinematic viscosity on target (fine) element fVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_c = 0.5 v^s_f ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(fVisc) cOmegaKine = mus_calcOmegaFromVisc(0.5_rk*fVisc) !evaluate scaling factor fac = getNonEqFac_intp_coarse_to_fine( cOmegaKine, fOmegaKine ) ! Rescale the non eq pdfs do iDir=1, QQ t_f_neq(iDir) = t_f_neq(iDir) * fac f(iDir) = t_f_neq(iDir) + t_f_eq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo end do ! indElem end subroutine fillFinerGhostsFromMe_weighAvg_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillFinerGhostsFromMe_weighAvgLES_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f_neq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: f_eq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: weight(layout%fStencil%QQ), sTurbVisc(layout%fStencil%QQ) real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation ! pdf to reconstruct from real(kind=rk) :: f( layout%fStencil%QQ ), tTurbVisc integer :: QQ type(mus_fluid_type), pointer :: fluid real(kind=rk) :: fVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level targetLevel = level + 1 dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all fine target elements: do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element treeId targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser) ! Find out how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals ! Read the pre-calculated weights (like in Average intp.) weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) & & %weight(1:nSourceElems) ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element position sourceElem = tLevelDesc%depFromCoarser( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) vel(3) = sAuxField(elemOff + vel_pos(3)) do iDir = 1, QQ ! calculate equilibrium density f_eq(idir,isourceelem) = layout%weight( idir ) & & * rho & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) f_neq(iDir, iSourceElem) = f(iDir) - f_eq(iDir, iSourceElem) enddo ! get turbulent viscosity sTurbVisc(iSourceElem) = fluid%turbulence%dataOnLvl(sourceLevel) & & %visc(sourceElem) end do ! iSourceElem ! interpolate all pdfs by weighted average tTurbVisc = sum(weight(1:nSourceElems)*sTurbVisc(1:nSourceElems)) do iDir = 1,QQ t_f_eq(iDir) = sum(weight(1:nSourceElems)*f_eq(iDir,1:nSourceElems)) t_f_neq(iDir) = sum(weight(1:nSourceElems)*f_neq(iDir,1:nSourceElems)) end do ! scale interpolated turbulent viscosity to target element fluid%turbulence%dataOnLvl(targetLevel)%visc(targetElem) & & = fluid%turbulence%fac_c2f*tTurbVisc ! get normalized kinematic viscosity on target (fine) element fVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_c = 0.5 v^s_f ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(fVisc + fluid%turbulence%fac_c2f*tTurbVisc) cOmegaKine = mus_calcOmegaFromVisc(0.5_rk*fVisc + tTurbVisc) !evaluate scaling factor fac = getNonEqFac_intp_coarse_to_fine( cOmegaKine, fOmegaKine ) ! Rescale the non eq pdfs do iDir=1, QQ t_f_neq(iDir) = t_f_neq(iDir) * fac f(iDir) = t_f_neq(iDir) + t_f_eq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo end do ! indElem end subroutine fillFinerGhostsFromMe_weighAvgLES_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillFinerGhostsFromMe_weighAvgIncomp_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f_neq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: f_eq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: weight(layout%fStencil%QQ) real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation ! pdf to reconstruct from real(kind=rk) :: f( layout%fStencil%QQ ) integer :: QQ type(mus_fluid_type), pointer :: fluid real(kind=rk) :: fVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level targetLevel = level + 1 dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all fine target elements: do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element treeId targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser) ! Find out how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals ! Read the pre-calculated weights (like in Average intp.) weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) & & %weight(1:nSourceElems) ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element position sourceElem = tLevelDesc%depFromCoarser( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) vel(3) = sAuxField(elemOff + vel_pos(3)) do iDir = 1, QQ f_eq(idir,isourceelem) = layout%weight( idir ) & & * ( rho + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:) * vel(:)) * 0.5_rk * cs2inv ) ) f_neq(iDir, iSourceElem) = f(iDir) - f_eq(iDir, iSourceElem) enddo end do ! iSourceElem ! interpolate all pdfs by weighted average do iDir = 1,QQ t_f_eq(iDir) = sum(weight(1:nSourceElems)*f_eq(iDir,1:nSourceElems)) t_f_neq(iDir) = sum(weight(1:nSourceElems)*f_neq(iDir,1:nSourceElems)) end do ! get normalized kinematic viscosity on target (fine) element fVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_c = 0.5 v^s_f ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(fVisc) cOmegaKine = mus_calcOmegaFromVisc(0.5_rk*fVisc) !evaluate scaling factor fac = getNonEqFac_intp_coarse_to_fine( cOmegaKine, fOmegaKine ) ! Rescale the non eq pdfs do iDir=1, QQ t_f_neq(iDir) = t_f_neq(iDir) * fac f(iDir) = t_f_neq(iDir) + t_f_eq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo end do ! indElem end subroutine fillFinerGhostsFromMe_weighAvgIncomp_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillFinerGhostsFromMe_weighAvgIncompLES_feq_fneq( method, & & fieldProp, tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, varSys, & & derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f_neq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: f_eq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: weight(layout%fStencil%QQ), sTurbVisc(layout%fStencil%QQ) real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation ! pdf to reconstruct from real(kind=rk) :: f( layout%fStencil%QQ ), tTurbVisc integer :: QQ type(mus_fluid_type), pointer :: fluid real(kind=rk) :: fVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level targetLevel = level + 1 dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all fine target elements: do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element treeId targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser) ! Find out how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals ! Read the pre-calculated weights (like in Average intp.) weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) & & %weight(1:nSourceElems) ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element position sourceElem = tLevelDesc%depFromCoarser( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) vel(3) = sAuxField(elemOff + vel_pos(3)) do iDir = 1, QQ f_eq(idir,isourceelem) = layout%weight( idir ) & & * ( rho + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:) * vel(:)) * 0.5_rk * cs2inv ) ) f_neq(iDir, iSourceElem) = f(iDir) - f_eq(iDir, iSourceElem) enddo ! get turbulent viscosity sTurbVisc(iSourceElem) = fluid%turbulence%dataOnLvl(sourceLevel) & & %visc(sourceElem) end do ! iSourceElem ! interpolate all pdfs by weighted average tTurbVisc = sum(weight(1:nSourceElems)*sTurbVisc(1:nSourceElems)) do iDir = 1,QQ t_f_eq(iDir) = sum(weight(1:nSourceElems)*f_eq(iDir,1:nSourceElems)) t_f_neq(iDir) = sum(weight(1:nSourceElems)*f_neq(iDir,1:nSourceElems)) end do ! scale interpolated turbulent viscosity to target element fluid%turbulence%dataOnLvl(targetLevel)%visc(targetElem) & & = fluid%turbulence%fac_c2f*tTurbVisc ! get normalized kinematic viscosity on target (fine) element fVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_c = 0.5 v^s_f ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(fVisc + fluid%turbulence%fac_c2f*tTurbVisc) cOmegaKine = mus_calcOmegaFromVisc(0.5_rk*fVisc + tTurbVisc) !evaluate scaling factor fac = getNonEqFac_intp_coarse_to_fine( cOmegaKine, fOmegaKine ) ! Rescale the non eq pdfs do iDir=1, QQ t_f_neq(iDir) = t_f_neq(iDir) * fac f(iDir) = t_f_neq(iDir) + t_f_eq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo end do ! indElem end subroutine fillFinerGhostsFromMe_weighAvgIncompLES_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillFinerGhostsFromMe_weighAvg2D_feq_fneq( method, fieldProp, & & tLevelDesc, level, sState, sNeigh, snSize, sAuxField, tState, & & tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f_neq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: f_eq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: weight(layout%fStencil%QQ) real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation ! pdf to reconstruct from real(kind=rk) :: f( layout%fStencil%QQ ) integer :: QQ type(mus_fluid_type), pointer :: fluid real(kind=rk) :: fVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- vel = 0._rk fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level targetLevel = level + 1 dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all fine target elements: do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element treeId targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser) ! Find out how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals ! Read the pre-calculated weights (like in Average intp.) weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) & & %weight(1:nSourceElems) ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element position sourceElem = tLevelDesc%depFromCoarser( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) do iDir = 1, QQ ! calculate equilibrium density f_eq(idir,isourceelem) = layout%weight( idir ) & & * rho & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) f_neq(iDir, iSourceElem) = f(iDir) - f_eq(iDir, iSourceElem) enddo end do ! iSourceElem ! interpolate all pdfs by weighted average do iDir = 1,QQ t_f_eq(iDir) = sum(weight(1:nSourceElems)*f_eq(iDir,1:nSourceElems)) t_f_neq(iDir) = sum(weight(1:nSourceElems)*f_neq(iDir,1:nSourceElems)) end do ! get normalized kinematic viscosity on target (fine) element fVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_c = 0.5 v^s_f ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(fVisc) cOmegaKine = mus_calcOmegaFromVisc(0.5_rk*fVisc) !evaluate scaling factor fac = getNonEqFac_intp_coarse_to_fine( cOmegaKine, fOmegaKine ) ! Rescale the non eq pdfs do iDir=1, QQ t_f_neq(iDir) = t_f_neq(iDir) * fac f(iDir) = t_f_neq(iDir) + t_f_eq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo end do ! indElem end subroutine fillFinerGhostsFromMe_weighAvg2D_feq_fneq ! **************************************************************************** ! ! **************************************************************************** ! !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner !! 1. Calculate Equilibrium and nonEquilibrium !! 2. Compute scaling factor !! 3. calculate target: Eq + Scale * nonEquilibrium !! !! This subroutine's interface must match the abstract interface definition !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer. subroutine fillFinerGhostsFromMe_weighAvg2DIncomp_feq_fneq( method, & & fieldProp, tLevelDesc, level, sState, sNeigh, snSize, sAuxField, & & tState, tNeigh, tnSize, layout, nTargets, targetList, physics, time, & & varSys, derVarPos ) ! -------------------------------------------------------------------- ! class(mus_interpolation_method_type), intent(inout) :: method !> Array of field properties (fluid or species) type(mus_field_prop_type), target, intent(in) :: fieldProp(:) !> level descriptor on target level type( tem_levelDesc_type ), intent(in) :: tLevelDesc !> my refinement level integer, intent(in) :: level !> State vector of SOURCE FLUID elements real(kind=rk), intent(in) :: sState(:) integer, intent(in) :: sNeigh(:) integer, intent(in) :: snSize !> State vector of TARGET GHOST elements real(kind=rk), intent(inout) :: tState(:) integer, intent(in) :: tNeigh(:) integer, intent(in) :: tnSize !> AuxField variable to read rho and vel from source elements real(kind=rk), intent(inout) :: sAuxField(:) !> the layout used type( mus_scheme_layout_type ), intent(in) :: layout !> List of target elements ( their position in depSource list ) integer, intent(in) :: nTargets integer, intent(in) :: targetList(nTargets) !> physics type to convert lattice to physics SI unit and vice versa !! @todo: This can be replaced by scale factor by level type( mus_physics_type ), intent(in) :: physics !> time required to compute viscosity on target element barycenter type(tem_time_type), intent(in) :: time !> scheme variable system type( tem_varSys_type ), intent(in) :: varSys !> position of all derive variable in varSys for all fields type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: sourceLevel ! level of source elements integer :: sourceElem ! treeId of current source element integer :: targetLevel ! level of target elements integer :: targetElem ! treeId of current source element integer :: iDir ! current direction (discrete velocity) for loop integer :: iElem ! current target element (for outer loop) integer :: indElem ! element counter for indirection list integer :: iSourceElem ! current source element (for inner loop) integer :: nSourceElems ! number of source elements for the current target real(kind=rk) :: f_neq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: f_eq( layout%fStencil%QQ, layout%fStencil%QQ ) real(kind=rk) :: weight(layout%fStencil%QQ) real(kind=rk) :: t_f_eq( layout%fStencil%QQ ) ! temp pdf calculation real(kind=rk) :: t_f_neq( layout%fStencil%QQ ) ! temp pdf calculation ! pdf to reconstruct from real(kind=rk) :: f( layout%fStencil%QQ ) integer :: QQ type(mus_fluid_type), pointer :: fluid real(kind=rk) :: fVisc, cOmegaKine, fOmegaKine, rho, vel(3), fac integer :: elemOff, dens_pos, vel_pos(3), nScalars ! -------------------------------------------------------------------------- vel = 0._rk fluid => fieldProp(1)%fluid nScalars = varSys%nScalars QQ = layout%fStencil%QQ sourceLevel = level targetLevel = level + 1 dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! Treat all fine target elements: do indElem = 1, nTargets iElem = targetList( indElem ) ! Read the target element treeId targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser) ! Find out how many fine source elements we have for interpolation. nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals ! Read the pre-calculated weights (like in Average intp.) weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) & & %weight(1:nSourceElems) ! Now loop over all fine source elements for this target: do iSourceElem = 1, nSourceElems ! Get the source element position sourceElem = tLevelDesc%depFromCoarser( iElem ) & & %elem%val( iSourceElem ) do iDir = 1, QQ ! this is post collision f(iDir) = sState( ( sourceelem-1)* nscalars+ idir+( 1-1)* qq) enddo ! element offset for auxField array elemOff = (sourceElem-1)*varSys%nAuxScalars ! local density rho = sAuxField(elemOff + dens_pos) ! local x-, y-, z-velocity vel(1) = sAuxField(elemOff + vel_pos(1)) vel(2) = sAuxField(elemOff + vel_pos(2)) do iDir = 1, QQ f_eq(idir,isourceelem) = layout%weight( idir ) & & * ( rho + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:) * vel(:)) * 0.5_rk * cs2inv ) ) f_neq(iDir, iSourceElem) = f(iDir) - f_eq(iDir, iSourceElem) enddo end do ! iSourceElem ! interpolate all pdfs by weighted average do iDir = 1,QQ t_f_eq(iDir) = sum(weight(1:nSourceElems)*f_eq(iDir,1:nSourceElems)) t_f_neq(iDir) = sum(weight(1:nSourceElems)*f_neq(iDir,1:nSourceElems)) end do ! get normalized kinematic viscosity on target (fine) element fVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem) ! relation between coarse and fine grid kinematic viscosity: ! v^s_c = 0.5 v^s_f ! calculate omega on source and target level fOmegaKine = mus_calcOmegaFromVisc(fVisc) cOmegaKine = mus_calcOmegaFromVisc(0.5_rk*fVisc) !evaluate scaling factor fac = getNonEqFac_intp_coarse_to_fine( cOmegaKine, fOmegaKine ) ! Rescale the non eq pdfs do iDir=1, QQ t_f_neq(iDir) = t_f_neq(iDir) * fac f(iDir) = t_f_neq(iDir) + t_f_eq(iDir) ! Now write the resulting pdf in the current direction to the target tState( ( targetelem-1)* nscalars+ idir+( 1-1)* qq) = f(iDir) enddo end do ! indElem end subroutine fillFinerGhostsFromMe_weighAvg2DIncomp_feq_fneq ! **************************************************************************** ! end module mus_interpolate_average_module ! **************************************************************************** !