! Copyright (c) 2017, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2019 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 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. ! 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. ! **************************************************************************** ! !> author: Kannan Masilamani !! This module provides the definition and methods for !! MRT advection relaxation scheme for D3Q27 stencil. !! The weighted MRT is based on the following paper !! Abbas Fakhari, Diogo Bolster, Li-Shi Luo !! "A weighted multiple-relaxation-time lattice Boltzmann method for multiphase !! flows and its application to partial coalescence cascades" !! Journal of Computational Physics, 2017 !! !! Density (rho) and velocity (ux, uy, uz) are conserved during collision. !! i.e. m(1) = meq(1) --> mneq(1) = 0 !! m(2) = meq(2) --> mneq(2) = 0 !! m(3) = meq(3) --> mneq(3) = 0 !! m(4) = meq(4) --> mneq(4) = 0 !! !! Collision parameters are chosen as !! s(1:4) = 0 !! s(5:9) = omega !! s(10) = bulk_omega !! s(11:27) = 1.0 module mus_mrt_d3q27_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: rho0, cs2inv, cs4inv, t2cs2inv, t2cs4inv ! 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_param_module, only: mus_param_type use mus_varSys_module, only: mus_varSys_data_type use mus_derVarPos_module, only: mus_derVarPos_type use mus_mrtInit_module, only: WMMtrD3Q27, WMMIvD3Q27 implicit none private public :: weighted_mrt_advRel_d3q27_generic public :: weighted_mrt_advRel_d3q27_incomp_generic public :: weighted_mrt_advRel_d3q27 public :: weighted_mrt_advRel_d3q27_incomp !============================================================================= ! D3Q27 flow model !============================================================================= !> Definition of the discrete velocity set integer,parameter :: QQ = 27 !< number of pdf directions integer,parameter :: q000 = 27 !< rest density is last contains ! **************************************************************************** ! !> Unoptimized explicit implementation !! !! 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 weighted_mrt_advRel_d3q27_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 real(kind=rk) :: pdfTmp( QQ ) ! temporary local pdf values real(kind=rk) :: rho ! local density real(kind=rk) :: u_x ! local x-velocity real(kind=rk) :: u_y ! local y-velocity real(kind=rk) :: u_z ! local z-velocity real(kind=rk) :: usq, ucx real(kind=rk) :: fEq( QQ ) !< equilibrium distribution ! MRT Variables real(kind=rk) :: s_mrt( QQ ) real(kind=rk) :: mneq( QQ ) real(kind=rk) :: fneq( QQ ) real(kind=rk) :: omegaBulk 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) omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level) ! collision matrix s_mrt = fieldProp(1)%fluid%mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ) !NEC$ ivdep !cdir nodep !ibm* novector !dir$ novector nodeloop: do iElem = 1,nSolve !> First load all local values into temp array do iDir = 1, QQ pdfTmp( iDir ) = inState( neigh (( idir-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) end do elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + dens_pos) ! local x-, y- and z-velocity u_x = auxField(elemOff + vel_pos(1)) u_y = auxField(elemOff + vel_pos(2)) u_z = auxField(elemOff + vel_pos(3)) usq = u_x*u_x + u_y*u_y + u_z*u_z ! Calculate the equilibrium distribution function do iDir = 1, QQ ! velocity times lattice unit velocity ucx = layout%fStencil%cxDirRK(1,iDir) * u_x & & + layout%fStencil%cxDirRK(2,iDir) * u_y & & + layout%fStencil%cxDirRK(3,iDir) * u_z ! calculate equilibrium density fEq( iDir ) = layout%weight( iDir ) * rho * ( 1.d0 + ucx*cs2inv & & + ucx*ucx*t2cs4inv - usq*t2cs2inv ) enddo ! Calculate the non-equilibrium part fnEq(:) = pdfTmp(:) - fEq(:) ! convert non-equilibrium PDF into moments do iDir = 1, QQ mneq(iDir) = sum( fnEq(:) * WMMtrD3Q27(iDir,:) ) end do ! viscosity omega s_mrt( 5:9 ) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) ! multiply it with collision matrix mneq(:) = s_mrt(:) * mneq(:) ! compute fNeq do iDir = 1, QQ fneq(iDir) = sum( WMMIvD3Q27(iDir,:) * mneq(:) ) outState((ielem-1)*qq+idir+(1-1)*qq) = pdfTmp(iDir) - fneq(iDir) end do enddo nodeloop end subroutine weighted_mrt_advRel_d3q27_generic ! **************************************************************************** ! ! **************************************************************************** ! !> Unoptimized explicit implementation !! !! 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 weighted_mrt_advRel_d3q27_incomp_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 real(kind=rk) :: pdfTmp( QQ ) ! temporary local pdf values real(kind=rk) :: rho ! local density real(kind=rk) :: u_x ! local x-velocity real(kind=rk) :: u_y ! local y-velocity real(kind=rk) :: u_z ! local z-velocity real(kind=rk) :: usq, ucx real(kind=rk) :: fEq( QQ ) !< equilibrium distribution ! MRT Variables real(kind=rk) :: s_mrt( QQ ) real(kind=rk) :: mneq( QQ ) real(kind=rk) :: fneq( QQ ) real(kind=rk) :: omegaBulk 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) omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level) ! collision matrix s_mrt = fieldProp(1)%fluid%mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ) !NEC$ ivdep !cdir nodep !ibm* novector !dir$ novector nodeloop: do iElem = 1,nSolve !> First load all local values into temp array do iDir = 1, QQ pdfTmp( iDir ) = inState( neigh (( idir-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) end do elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + dens_pos) ! local x-, y- and z-velocity u_x = auxField(elemOff + vel_pos(1)) u_y = auxField(elemOff + vel_pos(2)) u_z = auxField(elemOff + vel_pos(3)) usq = u_x*u_x + u_y*u_y + u_z*u_z ! Calculate the equilibrium distribution function do iDir = 1, QQ ! velocity times lattice unit velocity ucx = layout%fStencil%cxDirRK(1,iDir) * u_x & & + layout%fStencil%cxDirRK(2,iDir) * u_y & & + layout%fStencil%cxDirRK(3,iDir) * u_z ! calculate equilibrium density fEq( iDir ) = layout%weight( iDir ) * ( rho + rho0 * ( ucx*cs2inv & & + ucx*ucx*t2cs4inv - usq*t2cs2inv ) ) enddo ! Calculate the non-equilibrium part fnEq(:) = pdfTmp(:) - fEq(:) ! convert non-equilibrium PDF into moments do iDir = 1, QQ mneq(iDir) = sum( fnEq(:) * WMMtrD3Q27(iDir,:) ) end do ! viscosity omega s_mrt( 5:9 ) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) ! multiply it with collision matrix mneq(:) = s_mrt(:) * mneq(:) ! compute fNeq do iDir = 1, QQ fneq(iDir) = sum( WMMIvD3Q27(iDir,:) * mneq(:) ) outState((ielem-1)*qq+idir+(1-1)*qq) = pdfTmp(iDir) - fneq(iDir) end do enddo nodeloop end subroutine weighted_mrt_advRel_d3q27_incomp_generic ! **************************************************************************** ! ! **************************************************************************** ! !> Semi-optimized explicit implementation !! !! 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 weighted_mrt_advRel_d3q27( 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 integer :: nScalars real(kind=rk) :: f(QQ) real(kind=rk) :: rho ! local density real(kind=rk) :: u_x ! local x-velocity real(kind=rk) :: u_y ! local y-velocity real(kind=rk) :: u_z ! local z-velocity real(kind=rk) :: sum_19_22, sum_23_26, sum_19_26, sum_21_22_23_24 real(kind=rk) :: sum_20_21_24_25, sum_20_22_23_25, sum_7_10, sum_11_14 real(kind=rk) :: sum_15_18, sum_11_18, sum_19_20_23_24, sum_19_21_23_25 ! MRT Variables real(kind=rk) :: s_mrt( QQ ), meq( QQ ) real(kind=rk) :: mneq( QQ ), mom( QQ ) integer :: dens_pos, vel_pos(3), elemOff real(kind=rk) :: omegaBulk ! --------------------------------------------------------------------------- dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) nScalars = varSys%nScalars omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level) ! overwrite omegaKine term in the element loop s_mrt = fieldProp(1)%fluid%mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ) !NEC$ ivdep !cdir nodep !ibm* novector !dir$ novector nodeloop: do iElem = 1,nSolve !> First load all local values into temp array do iDir = 1, QQ f( iDir ) = inState( neigh (( idir-1)* nelems+ ielem)+( 1-1)* qq+ qq*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 u_x = auxField(elemOff + vel_pos(1)) u_y = auxField(elemOff + vel_pos(2)) u_z = auxField(elemOff + vel_pos(3)) ! Zero moment mom( 1) = rho ! First moments mom( 2) = rho*u_x mom( 3) = rho*u_y mom( 4) = rho*u_z sum_19_22 = f(19) + f(20) + f(21) + f(22) sum_23_26 = f(23) + f(24) + f(25) + f(26) sum_19_26 = sum_19_22 + sum_23_26 sum_21_22_23_24 = f(21) + f(22) + f(23) + f(24) mom( 5) = f(15) - f(16) - f(17) + f(18) + sum_19_26 - 2._rk*sum_21_22_23_24 sum_20_21_24_25 = f(20) + f(21) + f(24) + f(25) mom( 6) = f(7) - f(8) - f(9) + f(10) + sum_19_26 - 2._rk*sum_20_21_24_25 sum_20_22_23_25 = f(20) + f(22) + f(23) + f(25) mom( 7) = f(11) - f(12) - f(13) + f(14) + sum_19_26 - 2._rk*sum_20_22_23_25 sum_7_10 = f(7) + f(8) + f(9) + f(10) sum_11_14 = f(11) + f(12) + f(13) + f(14) sum_15_18 = f(15) + f(16) + f(17) + f(18) sum_11_18 = sum_11_14 + sum_15_18 mom( 8) = 2.0_rk*( f(1) + f(4) - sum_7_10 ) - f(2) - f(3) - f(5) - f(6) & & + sum_11_18 mom( 9) = f(2) - f(3) + f(5) - f(6) - sum_11_14 + sum_15_18 mom(10) = sum_7_10 + sum_11_18 + 2._rk*( sum_19_26 ) - f(27) mom(11) = 2.0_rk*( f(1) - f(4) ) - f(11) + f(12) - f(13) + f(14) - f(15) & & - f(16) + f(17) + f(18) + 4._rk*( -sum_19_22 + sum_23_26 ) sum_19_20_23_24 = f(19) + f(20) + f(23) + f(24) mom(12) = 2.0_rk*( f(2) - f(5) ) - f(7) - f(8) + f(9) + f(10) - f(15) + f(16) & & - f(17) + f(18) + 4._rk*( sum_19_26 - 2._rk*sum_19_20_23_24 ) sum_19_21_23_25 = f(19) + f(21) + f(23) + f(25) mom(13) = 2.0_rk*( f(3) - f(6) ) - f(7) + f(8) - f(9) + f(10) - f(11) - f(12) & & + f(13) + f(14) + 4._rk*( sum_19_26 - 2._rk*sum_19_21_23_25 ) mom(14) = f(11) - f(12) + f(13) - f(14) - f(15) - f(16) + f(17) + f(18) mom(15) = -f(7) - f(8) + f(9) + f(10) + f(15) - f(16) + f(17) - f(18) mom(16) = f(7) - f(8) + f(9) - f(10) - f(11) - f(12) + f(13) + f(14) mom(17) = -f(19) + f(20) + f(21) - f(22) + f(23) - f(24) - f(25) + f(26) mom(18) = -f(1) - f(2) - f(3) - f(4) - f(5) - f(6) & & + 4._rk*( sum_19_26 ) + f(27) mom(19) = 2.0_rk*( -f(1) - f(4) ) + f(2) + f(3) + f(5) + f(6) - 4._rk*sum_7_10 & & + 2._rk*( sum_11_18 ) mom(20) = -f(2) + f(3) - f(5) + f(6) + 2._rk*( -sum_11_14 + sum_15_18 ) mom(21) = -f(15) + f(16) + f(17) - f(18) + 2.0_rk*( sum_19_26 - 2._rk*sum_21_22_23_24 ) mom(22) = -f(7) + f(8) + f(9) - f(10) + 2.0_rk*( sum_19_26 - 2._rk*sum_20_21_24_25 ) mom(23) = -f(11) + f(12) + f(13) - f(14) + 2.0_rk*( sum_19_26 - 2._rk*sum_20_22_23_25 ) ! third order pseudo-vector mom(24) = -f(1) + f(4) + 2._rk*( f(11) - f(12) + f(13) - f(14) + f(15) + f(16) - f(17) & & -f(18) ) + 4._rk*( -sum_19_22 + sum_23_26 ) mom(25) = -f(2) + f(5) + 2._rk*( f(7) + f(8) - f(9) - f(10) + f(15) - f(16) + f(17) & & -f(18) ) + 4._rk*( sum_19_26 - 2._rk*sum_19_20_23_24 ) mom(26) = -f(3) + f(6) + 2._rk*( f(7) - f(8) + f(9) - f(10) + f(11) + f(12) - f(13) & & -f(14) ) + 4._rk*( sum_19_26 - 2._rk*sum_19_21_23_25 ) mom(27) = 2._rk*( f(1) + f(2) + f(3) + f(4) + f(5) + f(6) ) + 4._rk*( -sum_7_10 & & -sum_11_18 ) + 8._rk*( sum_19_26 ) - f(27) ! ------------------------------------------------------------------------- ! Equilibrium moments meq(1:QQ) = 0.0_rk meq( 1) = rho meq( 2) = rho*u_x meq( 3) = rho*u_y meq( 4) = rho*u_z meq( 5) = meq(2) * u_y meq( 6) = meq(3) * u_z meq( 7) = meq(4) * u_x meq( 8) = rho * (2.0_rk*u_x*u_x - u_y*u_y - u_z*u_z) meq( 9) = rho * (u_y*u_y - u_z*u_z) meq(10) = rho * (u_x*u_x + u_y*u_y + u_z*u_z) ! update kinematic omega s_mrt(5:9) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) ! compute neq moment mneq(:) = s_mrt(:) * ( mom(:) - meq(:) ) ! compute fNeq do iDir = 1, QQ f(iDir) = f(iDir) - sum( WMMIvD3Q27(iDir,:) * mneq(:) ) end do do iDir = 1, QQ outState(( ielem-1)* nscalars+ idir+( 1-1)* qq) = & & f( iDir ) end do enddo nodeloop end subroutine weighted_mrt_advRel_d3q27 ! **************************************************************************** ! ! **************************************************************************** ! !> Semi-optimized explicit implementation !! !! 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 weighted_mrt_advRel_d3q27_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 integer :: nScalars real(kind=rk) :: f(QQ) real(kind=rk) :: rho ! local density real(kind=rk) :: u_x ! local x-velocity real(kind=rk) :: u_y ! local y-velocity real(kind=rk) :: u_z ! local z-velocity real(kind=rk) :: sum_19_22, sum_23_26, sum_19_26, sum_21_22_23_24 real(kind=rk) :: sum_20_21_24_25, sum_20_22_23_25, sum_7_10, sum_11_14 real(kind=rk) :: sum_15_18, sum_11_18, sum_19_20_23_24, sum_19_21_23_25 ! MRT Variables real(kind=rk) :: s_mrt( QQ ), meq( QQ ) real(kind=rk) :: mneq( QQ ), mom( QQ ) integer :: dens_pos, vel_pos(3), elemOff real(kind=rk) :: omegaBulk ! --------------------------------------------------------------------------- dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) nScalars = varSys%nScalars omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level) ! overwrite omegaKine term in the element loop s_mrt = fieldProp(1)%fluid%mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ) !NEC$ ivdep !cdir nodep !ibm* novector !dir$ novector nodeloop: do iElem = 1,nSolve !> First load all local values into temp array do iDir = 1, QQ f( iDir ) = inState( neigh (( idir-1)* nelems+ ielem)+( 1-1)* qq+ qq*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 u_x = auxField(elemOff + vel_pos(1)) u_y = auxField(elemOff + vel_pos(2)) u_z = auxField(elemOff + vel_pos(3)) ! Zero moment mom( 1) = rho ! First moments mom( 2) = rho*u_x mom( 3) = rho*u_y mom( 4) = rho*u_z sum_19_22 = f(19) + f(20) + f(21) + f(22) sum_23_26 = f(23) + f(24) + f(25) + f(26) sum_19_26 = sum_19_22 + sum_23_26 sum_21_22_23_24 = f(21) + f(22) + f(23) + f(24) mom( 5) = f(15) - f(16) - f(17) + f(18) + sum_19_26 - 2._rk*sum_21_22_23_24 sum_20_21_24_25 = f(20) + f(21) + f(24) + f(25) mom( 6) = f(7) - f(8) - f(9) + f(10) + sum_19_26 - 2._rk*sum_20_21_24_25 sum_20_22_23_25 = f(20) + f(22) + f(23) + f(25) mom( 7) = f(11) - f(12) - f(13) + f(14) + sum_19_26 - 2._rk*sum_20_22_23_25 sum_7_10 = f(7) + f(8) + f(9) + f(10) sum_11_14 = f(11) + f(12) + f(13) + f(14) sum_15_18 = f(15) + f(16) + f(17) + f(18) sum_11_18 = sum_11_14 + sum_15_18 mom( 8) = 2.0_rk*( f(1) + f(4) - sum_7_10 ) - f(2) - f(3) - f(5) - f(6) & & + sum_11_18 mom( 9) = f(2) - f(3) + f(5) - f(6) - sum_11_14 + sum_15_18 mom(10) = sum_7_10 + sum_11_18 + 2._rk*( sum_19_26 ) - f(27) mom(11) = 2.0_rk*( f(1) - f(4) ) - f(11) + f(12) - f(13) + f(14) - f(15) & & - f(16) + f(17) + f(18) + 4._rk*( -sum_19_22 + sum_23_26 ) sum_19_20_23_24 = f(19) + f(20) + f(23) + f(24) mom(12) = 2.0_rk*( f(2) - f(5) ) - f(7) - f(8) + f(9) + f(10) - f(15) + f(16) & & - f(17) + f(18) + 4._rk*( sum_19_26 - 2._rk*sum_19_20_23_24 ) sum_19_21_23_25 = f(19) + f(21) + f(23) + f(25) mom(13) = 2.0_rk*( f(3) - f(6) ) - f(7) + f(8) - f(9) + f(10) - f(11) - f(12) & & + f(13) + f(14) + 4._rk*( sum_19_26 - 2._rk*sum_19_21_23_25 ) mom(14) = f(11) - f(12) + f(13) - f(14) - f(15) - f(16) + f(17) + f(18) mom(15) = -f(7) - f(8) + f(9) + f(10) + f(15) - f(16) + f(17) - f(18) mom(16) = f(7) - f(8) + f(9) - f(10) - f(11) - f(12) + f(13) + f(14) mom(17) = -f(19) + f(20) + f(21) - f(22) + f(23) - f(24) - f(25) + f(26) mom(18) = -f(1) - f(2) - f(3) - f(4) - f(5) - f(6) & & + 4._rk*( sum_19_26 ) + f(27) mom(19) = 2.0_rk*( -f(1) - f(4) ) + f(2) + f(3) + f(5) + f(6) - 4._rk*sum_7_10 & & + 2._rk*( sum_11_18 ) mom(20) = -f(2) + f(3) - f(5) + f(6) + 2._rk*( -sum_11_14 + sum_15_18 ) mom(21) = -f(15) + f(16) + f(17) - f(18) + 2.0_rk*( sum_19_26 - 2._rk*sum_21_22_23_24 ) mom(22) = -f(7) + f(8) + f(9) - f(10) + 2.0_rk*( sum_19_26 - 2._rk*sum_20_21_24_25 ) mom(23) = -f(11) + f(12) + f(13) - f(14) + 2.0_rk*( sum_19_26 - 2._rk*sum_20_22_23_25 ) ! third order pseudo-vector mom(24) = -f(1) + f(4) + 2._rk*( f(11) - f(12) + f(13) - f(14) + f(15) + f(16) - f(17) & & -f(18) ) + 4._rk*( -sum_19_22 + sum_23_26 ) mom(25) = -f(2) + f(5) + 2._rk*( f(7) + f(8) - f(9) - f(10) + f(15) - f(16) + f(17) & & -f(18) ) + 4._rk*( sum_19_26 - 2._rk*sum_19_20_23_24 ) mom(26) = -f(3) + f(6) + 2._rk*( f(7) - f(8) + f(9) - f(10) + f(11) + f(12) - f(13) & & -f(14) ) + 4._rk*( sum_19_26 - 2._rk*sum_19_21_23_25 ) mom(27) = 2._rk*( f(1) + f(2) + f(3) + f(4) + f(5) + f(6) ) + 4._rk*( -sum_7_10 & & -sum_11_18 ) + 8._rk*( sum_19_26 ) - f(27) ! ------------------------------------------------------------------------- ! Equilibrium moments meq(1:QQ) = 0.0_rk meq( 1) = rho meq( 2) = rho0*u_x meq( 3) = rho0*u_y meq( 4) = rho0*u_z meq( 5) = meq(2) * u_y meq( 6) = meq(3) * u_z meq( 7) = meq(4) * u_x meq( 8) = rho0 * (2.0_rk*u_x*u_x - u_y*u_y - u_z*u_z) meq( 9) = rho0 * (u_y*u_y - u_z*u_z) meq(10) = rho0 * (u_x*u_x + u_y*u_y + u_z*u_z) ! update kinematic omega s_mrt(5:9) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) ! compute neq moment mneq(:) = s_mrt(:) * ( mom(:) - meq(:) ) ! compute fNeq do iDir = 1, QQ f(iDir) = f(iDir) - sum( WMMIvD3Q27(iDir,:) * mneq(:) ) end do do iDir = 1, QQ outState(( ielem-1)* nscalars+ idir+( 1-1)* qq) = & & f( iDir ) end do enddo nodeloop end subroutine weighted_mrt_advRel_d3q27_incomp ! **************************************************************************** ! end module mus_mrt_d3q27_module