! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011-2014 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2011, 2013 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2011-2016, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2011 Jan Hueckelheim <j.hueckelheim@grs-sim.de> ! Copyright (c) 2012-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2012, 2014 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2018 Raphael Haupt <raphael.haupt@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) 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) 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. ! ****************************************************************************** ! !> This module provides the definition and methods for !! BGK advection relaxation scheme. module mus_bgk_module use iso_c_binding, only: c_f_pointer ! include treelm modules use env_module, only: rk use tem_varSys_module, only: tem_varSys_type !, tem_varSys_op_type use tem_param_module, only: div1_3, div1_36, div1_8, div3_4h, div1_4,& & div3_8, div9_16, div3_16, cs2inv, cs4inv,& & rho0 ! use tem_spacetime_fun_module, only: tem_st_fun_listElem_type, & ! & tem_spacetime_for ! use tem_logging_module, only: logUnit ! use tem_dyn_array_module, only: PositionOfVal ! use tem_aux_module, only: tem_abort ! include musubi modules use mus_field_prop_module, only: mus_field_prop_type use mus_scheme_layout_module, only: mus_scheme_layout_type use mus_scheme_type_module, only: mus_scheme_type use mus_derVarPos_module, only: mus_derVarPos_type use mus_param_module, only: mus_param_type use mus_varSys_module, only: mus_varSys_data_type implicit none private public :: bgk_advRel_generic public :: bgk_advRel_generic_incomp public :: bgk_advRel_flekkoy public :: bgk_advRel_flekkoy_noFluid contains ! **************************************************************************** ! !> Advection relaxation routine for the flekkoy diffusion model. !! !! This subroutine interface must match the abstract interface definition !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable !! via [[mus_scheme_type:compute]] function pointer. subroutine bgk_advRel_flekkoy( fieldProp, inState, outState, auxField, & & neigh, nElems, nSolve, level, layout, & & params, varSys, derVarPos ) ! -------------------------------------------------------------------- ! !> Array of field properties (fluid or species) type(mus_field_prop_type), intent(in) :: fieldProp(:) !> variable system definition type(tem_varSys_type), intent(in) :: varSys !> current layout type(mus_scheme_layout_type), intent(in) :: layout !> number of elements in state Array integer, intent(in) :: nElems !> input pdf vector real(kind=rk), intent(in) :: inState(nElems * varSys%nScalars) !> output pdf vector real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars) !> Auxiliary field computed from pre-collision state !! Is updated with correct velocity field for multicomponent models real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars) !> connectivity vector integer, intent(in) :: neigh(nElems * layout%fStencil%QQ) !> number of elements solved in kernel integer, intent(in) :: nSolve !> current level integer,intent(in) :: level !> global parameters type(mus_param_type),intent(in) :: params !> position of derived quantities in varsys for all fields type( mus_derVarPos_type ), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: iElem, iDir type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: pdfTmp( layout%fStencil%QQ ) ! temporary local pdf values real(kind=rk) :: rho, feq real(kind=rk) :: d_omega real(kind=rk) :: transVel( nSolve*3 ) ! velocity from the transport field real(kind=rk) :: uc ! u_i,fluid * c_i integer :: vel_varPos ! position of transport velocity variable in varSys real(kind=rk) :: inv_vel, u_fluid(3) ! -------------------------------------------------------------------------- ! access scheme via 1st variable method data which is a state variable call C_F_POINTER( varSys%method%val(derVarPos(1)%pdf)%method_Data, fPtr ) scheme => fPtr%solverData%scheme ! passive scalar has only one transport Variable vel_varPos = scheme%transVar%method(1)%data_varPos ! Get velocity field call varSys%method%val(vel_varPos)%get_valOfIndex( & & varSys = varSys, & & time = params%general%simControl%now, & & iLevel = level, & & idx = scheme%transVar%method(1) & & %pntIndex%indexLvl(level) & & %val(1:nSolve), & & nVals = nSolve, & & res = transVel ) ! convert physical velocity into LB velocity inv_vel = 1.0_rk / params%physics%fac( level )%vel transVel = transVel * inv_vel ! initialize and define some field wise parameters d_omega = 2._rk / ( 1._rk + 6._rk & & * fieldProp(1)%species%diff_coeff( 1 )) nodeloop: do iElem = 1, nSolve ! x-, y- and z-velocity from transport field u_fluid = transVel( (iElem-1)*3+1 : iElem*3 ) do iDir = 1, layout%fStencil%QQ pdfTmp( iDir ) = & & instate( neigh (( idir-1)* nelems+ ielem)+( 1-1)* layout%fstencil%qq+ varsys%nscalars*0 ) end do rho = sum( pdfTmp ) do iDir = 1, layout%fStencil%QQ ! compute c_i * u uc = dble( layout%fStencil%cxDir(1, iDir)) * u_fluid(1) + & & dble( layout%fStencil%cxDir(2, iDir)) * u_fluid(2) + & & dble( layout%fStencil%cxDir(3, iDir)) * u_fluid(3) ! compute the equilibrium (fi_eq = weight_i * rho * ( 1+c_i*u / cs^2)) feq = rho * layout%weight( iDir ) * (1._rk+3._rk*uc) outstate( & & ( ielem-1)* varsys%nscalars+ idir+( 1-1)* layout%fstencil%qq ) = & & pdfTmp( iDir ) + d_omega * ( feq - pdfTmp( idir )) end do end do nodeloop end subroutine bgk_advRel_flekkoy ! ****************************************************************************** ! ! ****************************************************************************** ! !> Advection relaxation routine for the flekkoy diffusion model without prior !! definition of a dependent fluid scheme. Velocity is read from the depend !! table. !! !! This subroutine interface must match the abstract interface definition !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable !! via [[mus_scheme_type:compute]] function pointer. subroutine bgk_advRel_flekkoy_noFluid( fieldProp, inState, outState, & & auxField, neigh, nElems, nSolve, & & level, layout, params, varSys, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Array of field properties (fluid or species) type(mus_field_prop_type), intent(in) :: fieldProp(:) !> variable system definition type(tem_varSys_type), intent(in) :: varSys !> current layout type(mus_scheme_layout_type), intent(in) :: layout !> number of elements in state Array integer, intent(in) :: nElems !> input pdf vector real(kind=rk), intent(in) :: inState(nElems * varSys%nScalars) !> output pdf vector real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars) !> Auxiliary field computed from pre-collision state !! Is updated with correct velocity field for multicomponent models real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars) !> connectivity vector integer, intent(in) :: neigh(nElems * layout%fStencil%QQ) !> number of elements solved in kernel integer, intent(in) :: nSolve !> current level integer,intent(in) :: level !> global parameters type(mus_param_type),intent(in) :: params !> position of derived quantities in varsys for all fields type( mus_derVarPos_type ), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: iElem integer :: iDir integer :: QQ, nScalars real(kind=rk) :: pdfTmp( layout%fStencil%QQ ) ! temporary local pdf values real(kind=rk) :: rho real(kind=rk) :: feq real(kind=rk) :: d_omega ! --------------------------------------------------------------------------- QQ = layout%fStencil%QQ nScalars = varSys%nScalars ! initialize and define some field wise parameters d_omega = 2._rk / ( 1._rk + 6._rk * fieldProp(1)%species%diff_coeff(1) ) nodeloop: do iElem = 1, nSolve rho = 0._rk do iDir = 1, QQ ! store the pdfs for this element and compute the density locally pdfTmp( iDir ) = & & instate( neigh (( idir-1)* nelems+ ielem)+( 1-1)* qq+ nscalars*0) rho = rho + pdfTmp( iDir ) end do directionloop: do iDir = 1, QQ ! compute the equilibrium ! fi_eq = weight_i * rho * ( 1+c_i*u / cs^2)) ! = weight_i * rho feq = rho * layout%weight( iDir ) outstate( ( ielem-1)* nscalars+ idir+( 1-1)* qq ) = & & pdfTmp( iDir ) + d_omega * ( feq - pdfTmp( idir )) end do directionloop end do nodeloop end subroutine bgk_advRel_flekkoy_noFluid ! ****************************************************************************** ! ! ****************************************************************************** ! !> Advection relaxation routine for the !! BGK model with an explicit calculation of all equilibrium !! quantities. Slow and simple. This routine should only be !! used for testing purposes !! !! This subroutine interface must match the abstract interface definition !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable !! via [[mus_scheme_type:compute]] function pointer. subroutine bgk_advRel_generic( fieldProp, inState, outState, auxField, & & neigh, nElems, nSolve, level, layout, & & params, varSys, derVarPos ) ! -------------------------------------------------------------------- ! !> Array of field properties (fluid or species) type(mus_field_prop_type), intent(in) :: fieldProp(:) !> variable system definition type(tem_varSys_type), intent(in) :: varSys !> current layout type(mus_scheme_layout_type), intent(in) :: layout !> number of elements in state Array integer, intent(in) :: nElems !> input pdf vector real(kind=rk), intent(in) :: inState(nElems * varSys%nScalars) !> output pdf vector real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars) !> Auxiliary field computed from pre-collision state !! Is updated with correct velocity field for multicomponent models real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars) !> connectivity vector integer, intent(in) :: neigh(nElems * layout%fStencil%QQ) !> number of elements solved in kernel integer, intent(in) :: nSolve !> current level integer,intent(in) :: level !> global parameters type(mus_param_type),intent(in) :: params !> position of derived quantities in varsys for all fields type( mus_derVarPos_type ), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: iElem, iDir ! voxel element counter integer :: QQ, nScalars ! temporary distribution variables real(kind=rk) pdfTmp(layout%fStencil%QQ) real(kind=rk) ux(3) ! local velocity real(kind=rk) rho ! local density real(kind=rk) usq ! square velocity ! derived constants ! equilibrium calculation variables real(kind=rk) ucx real(kind=rk) eqState(layout%fStencil%QQ) real(kind=rk) omega integer :: dens_pos, vel_pos(3), elemOff ! --------------------------------------------------------------------------- dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) QQ = layout%fStencil%QQ ! nElems = size(neigh)/QQ nScalars = varSys%nScalars nodeloop: do iElem=1,nSolve !> Generic fetching step: !! Streaming for pull !! Local copy for push ux = 0._rk rho = 0._rk do iDir = 1, QQ pdfTmp( iDir ) = inState( neigh(( idir-1)* nelems+ ielem)+( 1-1)* qq+ nscalars*0) end do ! element offset for auxField array elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + dens_pos) ! local x-, y- and z-velocity ux(1) = auxField(elemOff + vel_pos(1)) ux(2) = auxField(elemOff + vel_pos(2)) ux(3) = auxField(elemOff + vel_pos(3)) ! square velocity and derived constants usq = ux(1)*ux(1) + ux(2)*ux(2) + ux(3)*ux(3) !> relaxation parameter omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) do iDir = 1, QQ !> Pre-calculate velocitiy terms ucx = layout%fStencil%cxDirRK( 1, iDir )*ux(1) & & + layout%fStencil%cxDirRK( 2, iDir )*ux(2) & & + layout%fStencil%cxDirRK( 3, iDir )*ux(3) !> Calculate equilibrium distribution functions fEq eqState(iDir) = layout%weight(iDir)*rho*( 1.0_rk & & + ucx*cs2inv & & + ucx*ucx*cs2inv*cs2inv*0.5_rk & & - usq*0.5_rk*cs2inv ) !> Relaxation outState( ( ielem-1)* nscalars+ idir+( 1-1)* qq) = & & pdfTmp( iDir ) - omega*( pdfTmp(iDir ) - eqState( iDir )) end do !< iDir end do nodeloop end subroutine bgk_advRel_generic ! ****************************************************************************** ! ! ****************************************************************************** ! !> Advection relaxation routine for the !! BGK model with an explicit calculation of all equilibrium !! quantities. Slow and simple. This routine should only be !! used for testing purposes !! !! This subroutine interface must match the abstract interface definition !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable !! via [[mus_scheme_type:compute]] function pointer. subroutine bgk_advRel_generic_incomp( fieldProp, inState, outState, & & auxField, neigh, nElems, nSolve, & & level, layout, params, varSys, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Array of field properties (fluid or species) type(mus_field_prop_type), intent(in) :: fieldProp(:) !> variable system definition type(tem_varSys_type), intent(in) :: varSys !> current layout type(mus_scheme_layout_type), intent(in) :: layout !> number of elements in state Array integer, intent(in) :: nElems !> input pdf vector real(kind=rk), intent(in) :: inState(nElems * varSys%nScalars) !> output pdf vector real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars) !> Auxiliary field computed from pre-collision state !! Is updated with correct velocity field for multicomponent models real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars) !> connectivity vector integer, intent(in) :: neigh(nElems * layout%fStencil%QQ) !> number of elements solved in kernel integer, intent(in) :: nSolve !> current level integer,intent(in) :: level !> global parameters type(mus_param_type),intent(in) :: params !> position of derived quantities in varsys for all fields type( mus_derVarPos_type ), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! integer :: iElem, iDir ! voxel element counter integer :: QQ, nScalars ! temporary distribution variables real(kind=rk) pdfTmp(layout%fStencil%QQ) real(kind=rk) ux(3) ! local velocity real(kind=rk) rho ! local density real(kind=rk) usq ! square velocity ! derived constants ! equilibrium calculation variables real(kind=rk) ucx real(kind=rk) eqState(layout%fStencil%QQ) real(kind=rk) omega integer :: dens_pos, vel_pos(3), elemOff ! --------------------------------------------------------------------------- dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) QQ = layout%fStencil%QQ ! nElems = size(neigh)/QQ nScalars = varSys%nScalars nodeloop: do iElem = 1, nSolve !> Generic fetching step: !! Streaming for pull !! Local copy for push ux = 0._rk rho = 0._rk do iDir = 1, QQ pdfTmp( iDir ) = inState( neigh((idir-1)* nelems+ ielem)+( 1-1)* qq+ nscalars*0) end do ! element offset for auxField array elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + dens_pos) ! local x-, y- and z-velocity ux(1) = auxField(elemOff + vel_pos(1)) ux(2) = auxField(elemOff + vel_pos(2)) ux(3) = auxField(elemOff + vel_pos(3)) ! square velocity and derived constants usq = ux(1)*ux(1) + ux(2)*ux(2) + ux(3)*ux(3) !> relaxation parameter omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) do iDir = 1,QQ !> Pre-calculate velocitiy terms ucx = layout%fStencil%cxDir( 1, iDir )*ux(1) & & + layout%fStencil%cxDir( 2, iDir )*ux(2) & & + layout%fStencil%cxDir( 3, iDir )*ux(3) !> Calculate equilibrium distribution functions fEq eqState(iDir) = layout%weight(iDir)*( rho + rho0*( & & + ucx*cs2inv & & + ucx*ucx*cs2inv*cs2inv*0.5_rk & & - usq*0.5_rk*cs2inv )) !> Relaxation outState( (ielem-1)*nscalars+ idir+(1-1)*qq) = & & pdfTmp( iDir ) - omega*( pdfTmp(iDir) - eqState(iDir) ) end do ! iDir end do nodeloop end subroutine bgk_advRel_generic_incomp ! ****************************************************************************** ! end module mus_bgk_module ! ****************************************************************************** !