! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011-2013 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2011, 2017, 2019 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2011-2016, 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2011 Jan Hueckelheim <j.hueckelheim@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) 2012 Sathish Krishnan P S <s.krishnan@grs-sim.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 2021-2022 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.de> ! Copyright (c) 2021 Tobias Horstmann <tobias.horstmann@dlr.de> ! Copyright (c) 2021 Jana Gericke <jana.gericke@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: Jiaxing Qi !! Routines and parameter definitions for the standard D2Q9 model !! ! 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. module mus_d2q9_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: cs2inv, div1_12, cs4inv, cs6inv, cs8inv, & & div1_4, div1_9, div1_6, div1_18, div1_36, & & div1_2, div4_9, div2_9, div1_3, div2_3, & & cs2, rho0 use tem_aux_module, only: tem_abort use tem_property_module, only: prp_fluid, prp_fineGhostClosestToFluid use tem_construction_module, only: tem_levelDesc_type ! 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_derVarPos_module, only: mus_derVarPos_type use mus_param_module, only: mus_param_type use mus_varSys_module, only: mus_varSys_data_type use mus_directions_module, only: qN0, q0N, q10, q01, qNN, qN1, q1N, q11, & & q__W, q__S, q__E, q__N, & & q_SW, q_NW, q_SE, q_NE use mus_gradData_module, only: mus_gradData_type use mus_derivedQuantities_module2, only: second_order_moments_2D, & & getHermitepolynomials use mus_scheme_type_module, only: mus_scheme_type implicit none private public :: bgk_improved_advRel_d2q9 public :: bgk_advRel_d2q9 public :: mrt_advRel_d2q9 public :: mrt_advRel_d2q9_incomp public :: bgk_advRel_d2q9_incomp public :: f_f_eq_regularized_2nd_ord_d2q9 public :: f_f_eq_regularized_4th_ord_d2q9 public :: bgk_Regularized_d2q9 public :: bgk_RecursiveRegularized_d2q9 public :: bgk_HybridRecursiveRegularized_d2q9 public :: bgk_ProjectedRecursiveRegularized_d2q9 public :: bgk_DualRelaxationTime_RR_d2q9 ! ============================================================================ ! D2Q9 flow model ! ============================================================================ !> Definition of the discrete velocity set integer, parameter :: QQ = 9 !< number of pdf directions ! General direction index integer, parameter :: q00 = 9 !< rest integer, parameter :: q__0 = 9 !< rest density is last contains ! **************************************************************************** ! !> Improved BGK model (with Galilean correction term) !! taken from Martin Geier cumulent paper 2015 !! Geier, M., Schönherr, M., Pasquali, A., & Krafczyk, M. (2015). !! The cumulant lattice Boltzmann equation in three dimensions : Theory and !! validation. Computers and Mathematics with Applications. !! !! 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_improved_advRel_d2q9( 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 real(kind=rk) :: f(-1:1,-1:1) real(kind=rk) :: u,v,u2,v2 real(kind=rk) :: rho, rho_omg real(kind=rk) :: inv_rho real(kind=rk) :: omega, cmpl_o, fac real(kind=rk) :: m20, m02 real(kind=rk) :: sumX1, sumXN, sumY1, sumYN real(kind=rk) :: Gx, Gy real(kind=rk) :: x0, x1, xn, y0, y1, yn integer :: dens_pos, vel_pos(2), elemOff ! -------------------------------------------------------------------- ! dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve f(-1, 0) = inState( neigh(( qn0-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f( 0,-1) = inState( neigh(( q0n-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f( 1, 0) = inState( neigh(( q10-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f( 0, 1) = inState( neigh(( q01-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f(-1,-1) = inState( neigh(( qnn-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f(-1, 1) = inState( neigh(( qn1-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f( 1,-1) = inState( neigh(( q1n-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f( 1, 1) = inState( neigh(( q11-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f( 0, 0) = inState( neigh(( q00-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) ! element offset for auxField array elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + dens_pos) ! local x-, y- and z-velocity u = auxField(elemOff + vel_pos(1)) v = auxField(elemOff + vel_pos(2)) ! u v square u2 = u*u v2 = v*v sumX1 = sum(f( 1,:)) sumXN = sum(f(-1,:)) sumY1 = sum(f(:, 1)) sumYN = sum(f(:,-1)) inv_rho = 1.0_rk / rho ! second moments, by equation A.7 and A.8 m20 = ( sumX1 + sumXN ) * inv_rho m02 = ( sumY1 + sumYN ) * inv_rho ! relaxation rate omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) cmpl_o = 1._rk - omega fac = 4.5_rk - 2.25_rk * omega ! calculate Galilean correction term, by equation A.13 and A.14 Gx = fac * u2 * ( m20 - div1_3 - u2 ) Gy = fac * v2 * ( m02 - div1_3 - v2 ) ! X Y components of eq ! by equation A.19 - A.21 X0 = -div2_3 + u2 + Gx X1 = - ( X0 + 1.0_rk + u ) * 0.5_rk XN = X1 + u Y0 = -div2_3 + v2 + Gy Y1 = - ( Y0 + 1.0_rk + v ) * 0.5_rk YN = Y1 + v rho_omg = rho * omega X0 = X0 * rho_omg X1 = X1 * rho_omg XN = XN * rho_omg outState( ( ielem-1)* qq+ q00+( 1-1)* qq) & & = cmpl_o * f(0,0) + X0 * Y0 outState( ( ielem-1)* qq+ qn0+( 1-1)* qq) & & = cmpl_o * f(-1,0) + XN * Y0 outState( ( ielem-1)* qq+ q10+( 1-1)* qq) & & = cmpl_o * f(1,0) + X1 * Y0 outState( ( ielem-1)* qq+ q0n+( 1-1)* qq) & & = cmpl_o * f(0,-1) + X0 * YN outState( ( ielem-1)* qq+ q01+( 1-1)* qq) & & = cmpl_o * f(0,1) + X0 * Y1 outState( ( ielem-1)* qq+ qnn+( 1-1)* qq) & & = cmpl_o * f(-1,-1) + XN * YN outState( ( ielem-1)* qq+ q11+( 1-1)* qq) & & = cmpl_o * f(1,1) + X1 * Y1 outState( ( ielem-1)* qq+ q1n+( 1-1)* qq) & & = cmpl_o * f(1,-1) + X1 * YN outState( ( ielem-1)* qq+ qn1+( 1-1)* qq) & & = cmpl_o * f(-1,1) + XN * Y1 end do nodeloop !$omp end do nowait end subroutine bgk_improved_advRel_d2q9 ! **************************************************************************** ! ! **************************************************************************** ! !> No comment yet! !! !! TODO add comment !! !! 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_d2q9( 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 ! temporary distribution variables real(kind=rk) :: f1, f2, f3, f4, f5, f6, f7, f8, f9 real(kind=rk) :: fEq1, fEq2, fEq3, fEq4, fEq5, fEq6, fEq7, fEq8, fEq9 real(kind=rk) :: u_x, u_y real(kind=rk) :: rho, usq, ucx real(kind=rk) :: omega, cmpl_o real(kind=rk) :: c0, c1, c2, c3 integer :: dens_pos, vel_pos(2), elemOff !type(mus_varSys_data_type), pointer :: fPtr !type( tem_levelDesc_type ), pointer :: levelDesc ! -------------------------------------------------------------------- ! dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) !call c_f_pointer( varSys%method%val( 1 )%method_data, fPtr ) !levelDesc => fPtr%solverdata%scheme%levelDesc(level) !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve f1 = inState( neigh(( 1-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f2 = inState( neigh(( 2-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f3 = inState( neigh(( 3-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f4 = inState( neigh(( 4-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f5 = inState( neigh(( 5-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f6 = inState( neigh(( 6-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f7 = inState( neigh(( 7-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f8 = inState( neigh(( 8-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f9 = inState( neigh(( 9-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) !! no collision for coarse ghosts and for fine ghost cells far from fluid cell !if ( .not. btest(levelDesc%property( iElem ), prp_fluid) ) then ! if ( .not. btest(levelDesc%property( iElem ), prp_fineGhostClosestToFluid) ) then ! outState( ( ielem-1)* 9+ 9+( 1-1)* 9) = f9 ! outState( ( ielem-1)* 9+ 1+( 1-1)* 9) = f1 ! outState( ( ielem-1)* 9+ 2+( 1-1)* 9) = f2 ! outState( ( ielem-1)* 9+ 3+( 1-1)* 9) = f3 ! outState( ( ielem-1)* 9+ 4+( 1-1)* 9) = f4 ! outState( ( ielem-1)* 9+ 5+( 1-1)* 9) = f5 ! outState( ( ielem-1)* 9+ 6+( 1-1)* 9) = f6 ! outState( ( ielem-1)* 9+ 7+( 1-1)* 9) = f7 ! outState( ( ielem-1)* 9+ 8+( 1-1)* 9) = f8 ! CYCLE ! ! fine ghost cells closest to fluid cells should also skip the second collision ! else if ( params%iNesting( level ) == 2 ) then ! outState( ( ielem-1)* 9+ 9+( 1-1)* 9) = f9 ! outState( ( ielem-1)* 9+ 1+( 1-1)* 9) = f1 ! outState( ( ielem-1)* 9+ 2+( 1-1)* 9) = f2 ! outState( ( ielem-1)* 9+ 3+( 1-1)* 9) = f3 ! outState( ( ielem-1)* 9+ 4+( 1-1)* 9) = f4 ! outState( ( ielem-1)* 9+ 5+( 1-1)* 9) = f5 ! outState( ( ielem-1)* 9+ 6+( 1-1)* 9) = f6 ! outState( ( ielem-1)* 9+ 7+( 1-1)* 9) = f7 ! outState( ( ielem-1)* 9+ 8+( 1-1)* 9) = f8 ! CYCLE ! endif !end if ! 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)) ! calculate fEq usq = u_x * u_x + u_y * u_y c1 = rho - rho * usq * div1_2 * cs2inv feq9 = div4_9 * c1 c0 = rho * cs2inv * cs2inv * div1_2 c2 = rho * cs2inv * u_x c3 = rho * cs2inv * u_y feq1 = div1_9 * ( c1 - c2 + u_x * u_x * c0 ) feq3 = feq1 + div2_9 * c2 feq2 = div1_9 * ( c1 - c3 + u_y * u_y * c0 ) feq4 = feq2 + div2_9 * c3 ucx = u_x + u_y feq5 = div1_36 * ( c1 - c2 - c3 + ucx * ucx * c0 ) feq8 = feq5 + div1_18 * ( c2 + c3 ) ucx = u_x - u_y feq6 = div1_36 * ( c1 - c2 + c3 + ucx * ucx * c0 ) feq7 = feq6 + div1_18 * ( c2 - c3 ) omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) cmpl_o = 1._rk - omega outState( ( ielem-1)* 9+ 9+( 1-1)* 9) & & = cmpl_o * f9 + omega * fEq9 outState( ( ielem-1)* 9+ 1+( 1-1)* 9) & & = cmpl_o * f1 + omega * fEq1 outState( ( ielem-1)* 9+ 2+( 1-1)* 9) & & = cmpl_o * f2 + omega * fEq2 outState( ( ielem-1)* 9+ 3+( 1-1)* 9) & & = cmpl_o * f3 + omega * fEq3 outState( ( ielem-1)* 9+ 4+( 1-1)* 9) & & = cmpl_o * f4 + omega * fEq4 outState( ( ielem-1)* 9+ 5+( 1-1)* 9) & & = cmpl_o * f5 + omega * fEq5 outState( ( ielem-1)* 9+ 6+( 1-1)* 9) & & = cmpl_o * f6 + omega * fEq6 outState( ( ielem-1)* 9+ 7+( 1-1)* 9) & & = cmpl_o * f7 + omega * fEq7 outState( ( ielem-1)* 9+ 8+( 1-1)* 9) & & = cmpl_o * f8 + omega * fEq8 end do nodeloop !$omp end do nowait end subroutine bgk_advRel_d2q9 ! **************************************************************************** ! ! **************************************************************************** ! !> Advection relaxation routine for the D2Q9 MRT model !! f( x+c, t+1 ) = f(x,t) - M^(-1)S( m - meq ) !! !! 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 mrt_advRel_d2q9( 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 real(kind=rk) :: f1, f2, f3, f4, f5, f6, f7, f8, f9 real(kind=rk) :: fneq1, fneq2, fneq3, fneq4, fneq5, fneq6, fneq7, fneq8, & & fneq9 real(kind=rk) :: mom2, mom3, mom5, mom7, mom8, mom9 real(kind=rk) :: meq2, meq3, meq8, meq9 real(kind=rk) :: mneq2, mneq3, mneq5, mneq7, mneq8, mneq9 real(kind=rk) :: s_mrt(9), omegaBulk real(kind=rk) :: rho ! local density real(kind=rk) :: ux, uy, u2, v2 integer :: dens_pos, vel_pos(2), elemOff ! -------------------------------------------------------------------- ! dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) ! overwrite omegaKine indices 8 and 9 in element loop omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level) s_mrt = fieldProp(1)%fluid%mrtPtr( omegaKine = 1.0_rk, & & omegaBulk = omegaBulk, & & QQ = 9 ) !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve ! First load all local values into temp array f1 = inState( neigh(( 1-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f2 = inState( neigh(( 2-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f3 = inState( neigh(( 3-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f4 = inState( neigh(( 4-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f5 = inState( neigh(( 5-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f6 = inState( neigh(( 6-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f7 = inState( neigh(( 7-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f8 = inState( neigh(( 8-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f9 = inState( neigh(( 9-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) ! element offset for auxField array elemOff = (iElem - 1) * varSys%nAuxScalars ! local density rho = auxField(elemOff + dens_pos) ! local x-, y- and z-velocity ux = auxField(elemOff + vel_pos(1)) uy = auxField(elemOff + vel_pos(2)) ! square of velocity u2 = ux * ux v2 = uy * uy ! meq1 = rho ! rho meq2 = -2._rk * rho + 3._rk * rho * (u2 + v2) ! e meq3 = rho - 3._rk * rho * (u2 + v2) ! eps ! meq4 = ux ! jx ! meq5 = -ux ! qx ! meq6 = uy ! jy ! meq7 = -uy ! qy meq8 = rho * (u2 - v2) ! pxx meq9 = rho * ux * uy ! pxy ! convert pdf into moment mom2 = - f1 - f2 - f3 - f4 + 2.0_rk * ( f5 + f6 + f7 + f8 ) - 4.0_rk * f9 mom3 = -2.0_rk * (f1 + f2 + f3 + f4) + f5 + f6 + f7 + f8 + 4.0_rk * f9 mom5 = 2.0_rk * ( f1 - f3 ) - f5 - f6 + f7 + f8 mom7 = 2.0_rk * ( f2 - f4 ) - f5 + f6 - f7 + f8 mom8 = f1 - f2 + f3 - f4 mom9 = f5 - f6 - f7 + f8 ! omega s_mrt(8) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) s_mrt(9) = s_mrt(8) ! compute neq moment ! mneq1 = s_mrt1 * ( mom1 - meq1 ) = 0 mneq2 = s_mrt(2) * ( mom2 - meq2 ) mneq3 = s_mrt(3) * ( mom3 - meq3 ) ! mneq4 = s_mrt4 * ( mom4 - meq4 ) = 0 mneq5 = s_mrt(5) * ( mom5 + rho*ux ) ! meq5 = -ux ! mneq6 = s_mrt6 * ( mom6 - meq6 ) = 0 mneq7 = s_mrt(7) * ( mom7 + rho*uy ) ! meq7 = -uy mneq8 = s_mrt(8) * ( mom8 - meq8 ) mneq9 = s_mrt(9) * ( mom9 - meq9 ) ! compute fNeq ! do iDir = 1, 9 ! fneq(iDir) = sum( MMIvD2Q9(iDir,:) * mneq(:) ) ! end do ! mneq1 = mneq4 = mneq6 = 0 fNeq1 = -div1_36 * mneq2 - div1_18 * mneq3 + div1_6 * mneq5 & & + div1_4 * mneq8 fNeq2 = -div1_36 * mneq2 - div1_18 * mneq3 + div1_6 * mneq7 & & - div1_4 * mneq8 fNeq3 = -div1_36 * mneq2 - div1_18 * mneq3 - div1_6 * mneq5 & & + div1_4 * mneq8 fNeq4 = -div1_36 * mneq2 - div1_18 * mneq3 - div1_6 * mneq7 & & - div1_4 * mneq8 fNeq5 = div1_18 * mneq2 + div1_36 * mneq3 - div1_12 * (mneq5 + mneq7) & & + div1_4 * mneq9 fNeq6 = div1_18 * mneq2 + div1_36 * mneq3 - div1_12 * (mneq5 - mneq7) & & - div1_4 * mneq9 fNeq7 = div1_18 * mneq2 + div1_36 * mneq3 + div1_12 * (mneq5 - mneq7) & & - div1_4 * mneq9 fNeq8 = div1_18 * mneq2 + div1_36 * mneq3 + div1_12 * (mneq5 + mneq7) & & + div1_4 * mneq9 fNeq9 = div1_9 * (-mneq2 + mneq3) outState( ( ielem-1)* 9+ 1+( 1-1)* 9) = f1 - fNeq1 outState( ( ielem-1)* 9+ 2+( 1-1)* 9) = f2 - fNeq2 outState( ( ielem-1)* 9+ 3+( 1-1)* 9) = f3 - fNeq3 outState( ( ielem-1)* 9+ 4+( 1-1)* 9) = f4 - fNeq4 outState( ( ielem-1)* 9+ 5+( 1-1)* 9) = f5 - fNeq5 outState( ( ielem-1)* 9+ 6+( 1-1)* 9) = f6 - fNeq6 outState( ( ielem-1)* 9+ 7+( 1-1)* 9) = f7 - fNeq7 outState( ( ielem-1)* 9+ 8+( 1-1)* 9) = f8 - fNeq8 outState( ( ielem-1)* 9+ 9+( 1-1)* 9) = f9 - fNeq9 enddo nodeloop !$omp end do nowait end subroutine mrt_advRel_d2q9 ! **************************************************************************** ! ! **************************************************************************** ! !> Advection relaxation routine for the D2Q9 MRT model !! f( x+c, t+1 ) = f(x,t) - M^(-1)S( m - meq ) !! !! 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 mrt_advRel_d2q9_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 real(kind=rk) :: f1, f2, f3, f4, f5, f6, f7, f8, f9 real(kind=rk) :: fneq1, fneq2, fneq3, fneq4, fneq5, fneq6, fneq7, fneq8, & & fneq9 real(kind=rk) :: mom2, mom3, mom5, mom7, mom8, mom9 real(kind=rk) :: meq2, meq3, meq8, meq9 real(kind=rk) :: mneq2, mneq3, mneq5, mneq7, mneq8, mneq9 real(kind=rk) :: s_mrt(9), omegaBulk real(kind=rk) :: rho ! local density real(kind=rk) :: ux, uy, u2, v2 integer :: dens_pos, vel_pos(2), elemOff ! -------------------------------------------------------------------- ! dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) ! overwrite omegaKine indices 8 and 9 in element loop omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level) s_mrt = fieldProp(1)%fluid%mrtPtr( omegaKine = 1.0_rk, & & omegaBulk = omegaBulk, & & QQ = 9 ) !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve ! First load all local values into temp array f1 = inState( neigh(( 1-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f2 = inState( neigh(( 2-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f3 = inState( neigh(( 3-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f4 = inState( neigh(( 4-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f5 = inState( neigh(( 5-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f6 = inState( neigh(( 6-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f7 = inState( neigh(( 7-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f8 = inState( neigh(( 8-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) f9 = inState( neigh(( 9-1)* nelems+ ielem)+( 1-1)* qq+ qq*0) ! element offset for auxField array elemOff = (iElem - 1) * varSys%nAuxScalars ! local density rho = auxField(elemOff + dens_pos) ! local x-, y- and z-velocity ux = auxField(elemOff + vel_pos(1)) uy = auxField(elemOff + vel_pos(2)) u2 = ux * ux v2 = uy * uy ! meq1 = rho ! rho meq2 = -2._rk * rho + 3._rk * rho0 * (u2 + v2) ! e meq3 = rho - 3._rk * rho0 * (u2 + v2) ! eps ! meq4 = ux ! jx ! meq5 = -ux ! qx ! meq6 = uy ! jy ! meq7 = -uy ! qy meq8 = rho0 * (u2 - v2) ! pxx meq9 = rho0 * ux * uy ! pxy ! convert pdf into moment mom2 = - f1 - f2 - f3 - f4 + 2.0_rk * (f5 + f6 + f7 + f8) - 4.0_rk * f9 mom3 = -2.0_rk * (f1 + f2 + f3 + f4) + f5 + f6 + f7 + f8 + 4.0_rk * f9 mom5 = 2.0_rk * (f1 - f3) - f5 - f6 + f7 + f8 mom7 = 2.0_rk * (f2 - f4) - f5 + f6 - f7 + f8 mom8 = f1 - f2 + f3 - f4 mom9 = f5 - f6 - f7 + f8 ! omega s_mrt(8) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) s_mrt(9) = s_mrt(8) ! compute neq moment ! mneq1 = s_mrt1 * (mom1 - meq1) = 0 mneq2 = s_mrt(2) * (mom2 - meq2) mneq3 = s_mrt(3) * (mom3 - meq3) ! mneq4 = s_mrt4 * (mom4 - meq4) = 0 mneq5 = s_mrt(5) * (mom5 + rho0 * ux) ! meq5 = -ux ! mneq6 = s_mrt6 * (mom6 - meq6) = 0 mneq7 = s_mrt(7) * (mom7 + rho0 * uy) ! meq7 = -uy mneq8 = s_mrt(8) * (mom8 - meq8) mneq9 = s_mrt(9) * (mom9 - meq9) ! compute fNeq ! do iDir = 1, 9 ! fneq(iDir) = sum( MMIvD2Q9(iDir,:) * mneq(:) ) ! end do ! mneq1 = mneq4 = mneq6 = 0 fNeq1 = -div1_36 * mneq2 - div1_18 * mneq3 + div1_6 * mneq5 & & + div1_4 * mneq8 fNeq2 = -div1_36 * mneq2 - div1_18 * mneq3 + div1_6 * mneq7 & & - div1_4 * mneq8 fNeq3 = -div1_36 * mneq2 - div1_18 * mneq3 - div1_6 * mneq5 & & + div1_4 * mneq8 fNeq4 = -div1_36 * mneq2 - div1_18 * mneq3 - div1_6 * mneq7 & & - div1_4 * mneq8 fNeq5 = div1_18 * mneq2 + div1_36 * mneq3 - div1_12 * (mneq5 + mneq7) & & + div1_4 * mneq9 fNeq6 = div1_18 * mneq2 + div1_36 * mneq3 - div1_12 * (mneq5 - mneq7) & & - div1_4 * mneq9 fNeq7 = div1_18 * mneq2 + div1_36 * mneq3 + div1_12 * (mneq5 - mneq7) & & - div1_4 * mneq9 fNeq8 = div1_18 * mneq2 + div1_36 * mneq3 + div1_12 * (mneq5 + mneq7) & & + div1_4 * mneq9 fNeq9 = div1_9 * (-mneq2 + mneq3) outState( ( ielem-1)* 9+ 1+( 1-1)* 9) = f1 - fNeq1 outState( ( ielem-1)* 9+ 2+( 1-1)* 9) = f2 - fNeq2 outState( ( ielem-1)* 9+ 3+( 1-1)* 9) = f3 - fNeq3 outState( ( ielem-1)* 9+ 4+( 1-1)* 9) = f4 - fNeq4 outState( ( ielem-1)* 9+ 5+( 1-1)* 9) = f5 - fNeq5 outState( ( ielem-1)* 9+ 6+( 1-1)* 9) = f6 - fNeq6 outState( ( ielem-1)* 9+ 7+( 1-1)* 9) = f7 - fNeq7 outState( ( ielem-1)* 9+ 8+( 1-1)* 9) = f8 - fNeq8 outState( ( ielem-1)* 9+ 9+( 1-1)* 9) = f9 - fNeq9 enddo nodeloop !$omp end do nowait end subroutine mrt_advRel_d2q9_incomp ! **************************************************************************** ! ! **************************************************************************** ! !> No comment yet! !! !! TODO add comment !! !! 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_d2q9_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 ! temporary distribution variables real(kind=rk) :: f1, f2, f3, f4, f5, f6, f7, f8, f9 real(kind=rk) :: fEq1, fEq2, fEq3, fEq4, fEq5, fEq6, fEq7, fEq8, fEq9 real(kind=rk) :: u_x, u_y real(kind=rk) :: rho, usq, ucx real(kind=rk) :: omega, cmpl_o real(kind=rk) :: c0, c1, c2, c3 integer :: dens_pos, vel_pos(2), elemOff ! -------------------------------------------------------------------- ! dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve f1 = inState( neigh(( 1-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f2 = inState( neigh(( 2-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f3 = inState( neigh(( 3-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f4 = inState( neigh(( 4-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f5 = inState( neigh(( 5-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f6 = inState( neigh(( 6-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f7 = inState( neigh(( 7-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f8 = inState( neigh(( 8-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) f9 = inState( neigh(( 9-1)* nelems+ ielem)+( 1-1)* 9+ 9*0) ! 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)) ! calculate fEq usq = u_x * u_x + u_y * u_y c1 = rho - rho0 * usq * div1_2 * cs2inv feq9 = div4_9 * c1 c0 = rho0 * cs2inv * cs2inv * div1_2 c2 = rho0 * cs2inv * u_x c3 = rho0 * cs2inv * u_y feq1 = div1_9 * ( c1 - c2 + u_x*u_x*c0 ) feq3 = feq1 + div2_9 * c2 feq2 = div1_9 * ( c1 - c3 + u_y*u_y*c0 ) feq4 = feq2 + div2_9 * c3 ucx = u_x+u_y feq5 = div1_36 * ( c1 - c2 - c3 + ucx*ucx*c0 ) feq8 = feq5 + div1_18 * ( c2 + c3 ) ucx = u_x-u_y feq6 = div1_36 * ( c1 - c2 + c3 + ucx*ucx*c0 ) feq7 = feq6 + div1_18 * ( c2 - c3 ) omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) cmpl_o = 1._rk - omega outState(( ielem-1)* 9+9+( 1-1)* 9) & & = cmpl_o * f9 + omega * fEq9 outState(( ielem-1)* 9+1+( 1-1)* 9) & & = cmpl_o * f1 + omega * fEq1 outState(( ielem-1)* 9+2+( 1-1)* 9) & & = cmpl_o * f2 + omega * fEq2 outState(( ielem-1)* 9+3+( 1-1)* 9) & & = cmpl_o * f3 + omega * fEq3 outState(( ielem-1)* 9+4+( 1-1)* 9) & & = cmpl_o * f4 + omega * fEq4 outState(( ielem-1)* 9+5+( 1-1)* 9) & & = cmpl_o * f5 + omega * fEq5 outState(( ielem-1)* 9+6+( 1-1)* 9) & & = cmpl_o * f6 + omega * fEq6 outState(( ielem-1)* 9+7+( 1-1)* 9) & & = cmpl_o * f7 + omega * fEq7 outState(( ielem-1)* 9+8+( 1-1)* 9) & & = cmpl_o * f8 + omega * fEq8 end do nodeloop !$omp end do nowait end subroutine bgk_advRel_d2q9_incomp ! **************************************************************************** ! ! ****************************************************************************** ! ! > Todo add description ! ! Based on: ! Latt, J. and Chopard, B. (2005) 'Lattice Boltzmann Method with regularized ! non-equilibrium distribution functions' pure subroutine f_f_eq_regularized_2nd_ord_d2q9 ( weight, rho, u_x, u_y, feq, & & f1, a12xx, a12yy, a12xy ) ! -------------------------------------------------------------------- ! !> Weights of the stencil real(kind=rk), intent(in) :: weight(QQ) !> Density, velocity components real(kind=rk), intent(in) :: rho real(kind=rk), intent(in) :: u_x real(kind=rk), intent(in) :: u_y !> Equilibrium pdf and full pdf real(kind=rk), intent(out) :: feq(QQ) real(kind=rk), intent(out) :: f1(QQ) !> Coefficients of f1: a12xx, a12yy, a12xy real(kind=rk), intent(in) :: a12xx, a12yy, a12xy ! -------------------------------------------------------------------- ! real(kind=rk) :: u_x_sqr, u_y_sqr, u_x_u_y, f00, f01, f02, f12 ! -------------------------------------------------------------------------- u_x_sqr = u_x**2 u_y_sqr = u_y**2 u_x_u_y = u_x * u_y f00 = 1.0_rk !iDir = 1 f01 = -cs2inv*u_x f02 = div1_6*cs4inv*(2._rk*u_x_sqr - u_y_sqr) feq(1) = weight(1) * rho * (f00 + f01 + f02) f12 = div1_6*cs4inv*(2._rk*a12xx - a12yy) f1(1) = weight(1) * f12 !iDir = 3 f01 = -f01 feq(3) = weight(3) * rho * (f00 + f01 + f02) f1(3) = weight(3) * f12 !iDir = 2 f01 = -cs2inv*u_y f02 = div1_6*cs4inv*(-u_x_sqr + 2._rk*u_y_sqr) feq(2) = weight(2) * rho * (f00 + f01 + f02) f12 = div1_6*cs4inv*(-a12xx + 2._rk*a12yy) f1(2) = weight(2) * f12 !iDir = 4 f01 = -f01 feq(4) = weight(4) * rho * (f00 + f01 + f02) f1(4) = weight(4) * f12 !iDir = 5 f01 = cs2inv*(-u_x - u_y) f02 = cs4inv*(div1_3*(u_x_sqr + u_y_sqr) + u_x_u_y) feq(5) = weight(5) * rho * (f00 + f01 + f02) f12 = cs4inv*(div1_3*(a12xx + a12yy) + a12xy) f1(5) = weight(5) * f12 !iDir = 8 f01 = -f01 feq(8) = weight(8) * rho * (f00 + f01 + f02) f1(8) = weight(8) * f12 !iDir = 6 f01 = cs2inv*(-u_x + u_y) f02 = f02 - 2._rk*cs4inv*u_x_u_y feq(6) = weight(6) * rho * (f00 + f01 + f02) f12 = f12 - 2._rk*cs4inv*a12xy f1(6) = weight(6) * f12 !iDir = 7 f01 = -f01 feq(7) = weight(7) * rho * (f00 + f01 + f02) f1(7) = weight(7) * f12 !iDir = 9 !f01 = 0._rk f02 = -div1_6*cs4inv*(u_x_sqr + u_y_sqr) feq(9) = weight(9) * rho * (f00 + f02) f12 = -div1_6*cs4inv*(a12xx + a12yy) f1(9) = weight(9) * f12 end subroutine f_f_eq_regularized_2nd_ord_d2q9 ! ****************************************************************************** ! ! ****************************************************************************** ! ! > Modified by GGSpinelli THorstmann ! based on Lattice Boltzmann Method with regularized non-equilibrium distribution ! functions, Jonas Latt and Bastien Chopard 2005 pure subroutine f_f_eq_regularized_4th_ord_d2q9 ( weight, rho, u_x, u_y, feq, & & f1, a12xx, a12yy, a12xy ) ! -------------------------------------------------------------------- ! !> weights of the stencil real(kind=rk), intent(in) :: weight(QQ) !> density, velocity components real(kind=rk), intent(in) :: rho real(kind=rk), intent(in) :: u_x real(kind=rk), intent(in) :: u_y !> equilibrium pdf and full pdf real(kind=rk), intent(out) :: feq(QQ) real(kind=rk), intent(out) :: f1(QQ) !> coefficients of f1: a12xx, a12yy, a12xy real(kind=rk), intent(in) :: a12xx, a12yy, a12xy ! -------------------------------------------------------------------- ! real(kind=rk) :: u_x_sqr_u_y, u_y_sqr_u_x, u_x_sqr_u_y_sqr, f03, f04 real(kind=rk) :: u_x_sqr, u_y_sqr, a13xyy, a13xxy, a14xxyy, f13, f14 ! --------------------------------------------------------------------------- call f_f_eq_regularized_2nd_ord_d2q9 ( weight, rho, u_x, u_y, feq, f1, & & a12xx, a12yy, a12xy ) u_x_sqr = u_x**2 u_y_sqr = u_y**2 u_x_sqr_u_y = u_x_sqr * u_y u_y_sqr_u_x = u_y_sqr * u_x u_x_sqr_u_y_sqr = u_x_sqr * u_y_sqr a13xyy = 2.0_rk * u_y * a12xy + u_x * a12yy a13xxy = 2.0_rk * u_x * a12xy + u_y * a12xx a14xxyy = u_y*a13xxy + u_x_sqr*a12yy + 2.0_rk*u_x*u_y*a12xy !iDir = 1 f03 = div1_6*cs6inv*(u_y_sqr_u_x) f04 = div1_18*cs8inv*(-u_x_sqr_u_y_sqr) feq(1) = feq(1) + weight(1) * rho * (f03 + f04) f13 = div1_6*cs6inv*(a13xyy) f14 = div1_18*cs8inv*(-a14xxyy) f1(1) = f1(1) + weight(1) * (f13 + f14) !iDir = 3 f03 = -f03 feq(3) = feq(3) + weight(3) * rho * (f03 + f04) f13 = -f13 f1(3) = f1(3) + weight(3) * (f13 + f14) !iDir = 2 f03 = div1_6*cs6inv*(u_x_sqr_u_y) feq(2) = feq(2) + weight(2) * rho * (f03 + f04) f13 = div1_6*cs6inv*(a13xxy) f1(2) = f1(2) + weight(2) * (f13 + f14) !iDir = 4 f03 = -f03 feq(4) = feq(4) + weight(4) * rho * (f03 + f04) f13 = -f13 f1(4) = f1(4) + weight(4) * (f13 + f14) !iDir = 5 f03 = div1_3*cs6inv*(-(u_x_sqr_u_y + u_y_sqr_u_x)) f04 = -2._rk*f04 feq(5) = feq(5) + weight(5) * rho * (f03 + f04) f13 = div1_3*cs6inv*(-(a13xxy + a13xyy)) f14 = -2._rk*f14 f1(5) = f1(5) + weight(5) * (f13 + f14) !iDir = 8 f03 = -f03 feq(8) = feq(8) + weight(8) * rho * (f03 + f04) f13 = -f13 f1(8) = f1(8) + weight(8) * (f13 + f14) !iDir = 6 f03 = div1_3*cs6inv*((u_x_sqr_u_y - u_y_sqr_u_x)) feq(6) = feq(6) + weight(6) * rho * (f03 + f04) f13 = div1_3*cs6inv*((a13xxy - a13xyy)) f1(6) = f1(6) + weight(6) * (f13 + f14) !iDir = 7 f03 = -f03 feq(7) = feq(7) + weight(7) * rho * (f03 + f04) f13 = -f13 f1(7) = f1(7) + weight(7) * (f13 + f14) !iDir = 9 !f03 = 0._rk f04 = div1_4*f04 feq(9) = feq(9) + weight(9) * rho * (f04) !f03 + !f13 = 0._rk f14 = div1_4*f14 f1(9) = f1(9) + weight(9) * (f14) !f13 + end subroutine f_f_eq_regularized_4th_ord_d2q9 ! ****************************************************************************** ! ! ****************************************************************************** ! ! > Modified by GGSpinelli THorstmann ! based on Lattice Boltzmann Method with regularized non-equilibrium distribution ! functions, Jonas Latt and Bastien Chopard 2005 subroutine bgk_Regularized_d2q9 ( fieldProp, inState, outState, auxField, & & neigh, nElems, nSolve, level, layout, & & params, varSys, derVarPos) !, gradData ) ! -------------------------------------------------------------------- ! !> 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(:) ! -------------------------------------------------------------------- ! ! indeces integer :: iElem, iDir ! temporary distribution variables real(kind=rk) :: f( QQ ), SOM(3), SOM_neq(3) real(kind=rk) :: rho, u_x, u_y, a12xx, a12xy, a12yy real(kind=rk) :: omega, cmpl_o, feq(QQ), f1(QQ) integer :: denspos, velpos(3), elemOff, nScalars ! --------------------------------------------------------------------------- !cdir nodep !ibm* novector !dir$ novector denspos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) velpos(1:2) = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) nScalars = varSys%nScalars !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve do iDir = 1, QQ f(iDir) = inState( neigh(( idir-1)* nelems+ ielem)+( 1-1)* qq+ nscalars*0) enddo ! element offset for auxField array elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + denspos) ! local x-, y- and z-velocity u_x = auxField(elemOff + velpos(1)) u_y = auxField(elemOff + velpos(2)) ! non equilibrium second-order moments ! SOM_neq = SOM - SOM_eq SOM = second_order_moments_2D( f=f, layout=layout ) SOM_neq(1) = SOM(1) - rho * (cs2 + (u_x * u_x)) SOM_neq(2) = SOM(2) - rho * (cs2 + (u_y * u_y)) SOM_neq(3) = SOM(3) - rho * u_x * u_y ! Hermitian coefficients omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) cmpl_o = 1.0_rk - omega a12xx = SOM_neq(1) a12yy = SOM_neq(2) a12xy = SOM_neq(3) call f_f_eq_regularized_2nd_ord_d2q9( layout%weight(:), rho, u_x, u_y, feq, f1, & & a12xx, a12yy, a12xy ) do iDir = 1, QQ outState( ( ielem-1)* nscalars+ idir+(1-1)* qq) = feq(iDir) & & + cmpl_o * f1(iDir) enddo enddo nodeloop !$omp end do nowait end subroutine bgk_Regularized_d2q9 ! ****************************************************************************** ! ! ****************************************************************************** ! ! > Modified by GGSpinelli THorstmann ! based on Lattice Boltzmann Method with regularized non-equilibrium distribution ! functions, Jonas Latt and Bastien Chopard 2005 subroutine bgk_RecursiveRegularized_d2q9 ( fieldProp, inState, outState, auxField, & & neigh, nElems, nSolve, level, layout, & & params, varSys, derVarPos) !, gradData ) ! -------------------------------------------------------------------- ! !> 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(:) ! -------------------------------------------------------------------- ! ! indeces integer :: iElem, iDir ! temporary distribution variables real(kind=rk) :: f( QQ ), SOM(3), SOM_neq(3) real(kind=rk) :: rho, u_x, u_y, a12xx, a12xy, a12yy real(kind=rk) :: omega, cmpl_o, feq(QQ), f1(QQ) integer :: denspos, velpos(3), elemOff, nScalars ! --------------------------------------------------------------------------- !cdir nodep !ibm* novector !dir$ novector denspos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) velpos(1:2) = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) nScalars = varSys%nScalars !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve do iDir = 1, QQ f(iDir) = inState( neigh(( idir-1)* nelems+ ielem)+( 1-1)* qq+ nscalars*0) enddo ! element offset for auxField array elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + denspos) ! local x-, y- and z-velocity u_x = auxField(elemOff + velpos(1)) u_y = auxField(elemOff + velpos(2)) ! non equilibrium second-order moments ! SOM_neq = SOM - SOM_eq SOM = second_order_moments_2D( f=f, layout=layout ) SOM_neq(1) = SOM(1) - rho * (cs2 + (u_x * u_x)) SOM_neq(2) = SOM(2) - rho * (cs2 + (u_y * u_y)) SOM_neq(3) = SOM(3) - rho * u_x * u_y ! Hermitian coefficients omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) cmpl_o = 1.0_rk - omega a12xx = SOM_neq(1) a12yy = SOM_neq(2) a12xy = SOM_neq(3) call f_f_eq_regularized_4th_ord_d2q9( layout%weight(:), rho, u_x, u_y, feq, f1, & & a12xx, a12yy, a12xy ) do iDir = 1, QQ outState( ( ielem-1)* nscalars+ idir+(1-1)* qq) = feq(iDir) & & + cmpl_o * f1(iDir) enddo enddo nodeloop !$omp end do nowait end subroutine bgk_RecursiveRegularized_d2q9 ! ****************************************************************************** ! ! ****************************************************************************** ! ! > Modified by GGSpinelli THorstmann ! based on High-order extension of the recursive regularized lattice ! Boltzmann method, PhD Thesis, COREIXAS 2018 ! Correction term from but doesnt work! ! Solid wall and open boundary conditionsin hybrid recursive regularized ! latticeBoltzmann method for compressible flows, Feng 2019, PoF subroutine bgk_ProjectedRecursiveRegularized_d2q9 ( fieldProp, inState, outState, auxField, & & neigh, nElems, nSolve, level, layout, & & params, varSys, derVarPos) !, gradData ) ! -------------------------------------------------------------------- ! !> 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(:) ! -------------------------------------------------------------------- ! !> gradient data type(mus_gradData_type), pointer :: gradData type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme ! indeces integer :: iElem, iDir ! temporary distribution variables real(kind=rk) :: f( QQ ), SR(3), gradU(2,2,1)!TODO::, tr_SR !real(kind=rk) :: gradRHOU3(2,1), Scorr real(kind=rk) :: rho, u_x, u_y, a12xx, a12xy, a12yy real(kind=rk) :: omega, cmpl_o, feq(QQ), f1(QQ), taup integer :: denspos, velpos(3), elemOff, nScalars ! --------------------------------------------------------------------------- !cdir nodep !ibm* novector !dir$ novector denspos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) velpos(1:2) = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) nScalars = varSys%nScalars ! access gradData ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( 1 )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme gradData => scheme%gradData(level) !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve do iDir = 1, QQ f(iDir) = inState( neigh(( idir-1)* nelems+ ielem)+( 1-1)* qq+ nscalars*0) enddo ! element offset for auxField array elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + denspos) ! local x-, y- and z-velocity u_x = auxField(elemOff + velpos(1)) u_y = auxField(elemOff + velpos(2)) ! Stress tensor components gradU(:,:,1:1) = scheme%Grad%U_ptr( & & auxField = auxField, & & gradData = gradData, & & velPos = velPos, & & nAuxScalars = varSys%nAuxScalars, & & nDims = 2, & & nSolve = 1, & & elemOffset = iElem -1 ) ! symmetric strain rate tensors ! transformed inro RHS of a1 FD equation ! the trace is needed only for energy conservation, which is not ! done in Musubi at the moment. !TODO::tr_SR = (gradU(1, 1, 1) + gradU(2, 2, 1)) ! * 2 / D = 1 SR(1) = 2._rk * gradU(1, 1, 1) !S_XX !TODO:: - tr_SR SR(2) = 2._rk * gradU(2, 2, 1) !S_YY !TODO:: - tr_SR SR(3) = gradU(1, 2, 1) + gradU(2, 1, 1) !S_XY ! Hermitian coefficients omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) cmpl_o = 1.0_rk - omega taup = rho * cs2 / omega a12xx = -taup * SR(1) a12yy = -taup * SR(2) a12xy = -taup * SR(3) call f_f_eq_regularized_4th_ord_d2q9( layout%weight(:), rho, u_x, u_y, feq, f1, & & a12xx, a12yy, a12xy ) do iDir = 1, QQ outState( ( ielem-1)* nscalars+ idir+(1-1)* qq) = feq(iDir) & & + cmpl_o * f1(iDir) enddo enddo nodeloop !$omp end do nowait end subroutine bgk_ProjectedRecursiveRegularized_d2q9 ! **************************************************************************** ! ! **************************************************************************** ! !> Hybrid recursive regularization relaxation routine for the BGK model. !! todo add paper !! todo add some more information if needed !! !! 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_HybridRecursiveRegularized_d2q9( 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 !> Vurrent 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(:) ! -------------------------------------------------------------------- ! !> Gradient data type(mus_gradData_type), pointer :: gradData type(mus_varSys_data_type), pointer :: fPtr !type( tem_levelDesc_type ), pointer :: levelDesc ! Self-describing variables, loop indices, etc. integer :: iElem, iDir ! Temporary distribution variables real(kind=rk) :: pdfTmp(QQ) !symmetric strain rate tensor (SR) and it's trace (tr) real(kind=rk) :: SR(3)!TODO::, tr_SR ! Velocity gradient real(kind=rk) :: gradU(2,2,1) ! Second-order moments (SOM) and it's equilibrium real(kind=rk) :: SOM(3), SOM_neq(3) !real(kind=rk) :: gradRHOU3(2,1), Scorr ! Density and velocity real(kind=rk) :: rho, u_x, u_y ! Hermitian coefficents real(kind=rk) :: a12xx, a12xy, a12yy ! Relaxation parameter omega and it's complementary real(kind=rk) :: omega, cmpl_o ! todo add comment, maybe split in multiple lines real(kind=rk) :: taup, feq(QQ), f1(QQ) ! todo add comment .... sigma and it's complementary real(kind=rk) :: sigma, cmpl_sgm ! todo add comment integer :: denspos, velpos(3), elemOff, nScalars ! -------------------------------------------------------------------------- !cdir nodep !ibm* novector !dir$ novector ! Get position of density and velocity in auxField to determine them later denspos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) velpos(1:2) = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) nScalars = varSys%nScalars ! todo can we delete the commented out part? !call getHermitepolynomials( 2, QQ, layout, 4) ! Access data via variable method data which is a state variable call c_f_pointer( varSys%method%val( 1 )%method_data, fPtr ) ! Set pointers for gradData via method data ! todo remove commented out. it's just for explanation. we can reduce it to ! one pointer !scheme => fPtr%solverData%scheme !gradData => scheme%gradData(level) gradData => fPtr%solverData%scheme%gradData(level) ! todo can we delete levelDesc as it is not used anymore? see commented ! block below. !levelDesc => fPtr%solverdata%scheme%levelDesc(level) ! Sigma value, ideally read from input ! todo what happens if not from input --> add information sigma = fieldProp(1)%fluid%HRR_sigma ! if sigma == 1 --> no HRR !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve ! Set pdf for all QQ-directions !> Generic fetching step: !! Streaming for pull !! Local copy for push do iDir = 1, QQ pdfTmp(iDir) = inState( neigh(( idir-1)* nelems+ ielem)+( 1-1)* qq+ nscalars*0) enddo ! todo can we delete the commented out part? !! no collision for coarse ghosts and for fine ghost cells far from fluid cell !if ( .not. btest(levelDesc%property( iElem ), prp_fluid) ) then ! if ( .not. btest(levelDesc%property( iElem ), prp_fineGhostClosestToFluid) ) then ! do iDir = 1, QQ ! outState( ( ielem-1)* 9+ idir+( 1-1)* 9) = f(iDir) ! enddo ! CYCLE ! ! fine ghost cells closest to fluid cells should also skip the second collision ! else if ( params%iNesting( level ) == 2 ) then ! do iDir = 1, QQ ! outState( ( ielem-1)* 9+ idir+( 1-1)* 9) = f(iDir) ! enddo ! CYCLE ! endif !end if ! Calculate element offset for auxField array elemOff = (iElem-1)*varSys%nAuxScalars ! Set local density rho = auxField(elemOff + denspos) ! Set local x- and y-velocity u_x = auxField(elemOff + velpos(1)) u_y = auxField(elemOff + velpos(2)) ! Calculate velocity gradient using central or forward difference, latter ! one if element has a boundary as neighbor gradU(:,:,1:1) = fPtr%solverData%scheme%Grad%U_ptr( & & auxField = auxField, & & gradData = gradData, & & velPos = velPos, & & nAuxScalars = varSys%nAuxScalars, & & nDims = 2, & & nSolve = 1, & & elemOffset = iElem-1 ) ! Calculate symmetric strain rate tensor (SR) and it's trace (tr) ! the trace is needed only for energy conservation, which is not ! done in Musubi at the moment !TODO::tr_SR = (gradU(1, 1, 1) + gradU(2, 2, 1)) ! * 2 / D = 1 SR(1) = 2._rk * gradU(1, 1, 1) !S_XX !TODO:: - tr_SR SR(2) = 2._rk * gradU(2, 2, 1) !S_YY !TODO:: - tr_SR SR(3) = gradU(1, 2, 1) + gradU(2, 1, 1) !S_XY ! Determine the non-equilibrium second-order moments via ! SOM_neq = SOM - SOM_eq ! First, get the second-order moments SOM = second_order_moments_2D( f=pdfTmp, layout=layout ) ! Second, subtract it's equilibrium: SOM_eq = rho*(cs²+u²) SOM_neq(1) = SOM(1) - rho * (cs2 + (u_x * u_x)) SOM_neq(2) = SOM(2) - rho * (cs2 + (u_y * u_y)) SOM_neq(3) = SOM(3) - rho * u_x * u_y ! Hermitian coefficients ! Relaxation parameter omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) ! Complmentary of relaxation parameter (cmpl_o) cmpl_o = 1.0_rk - omega ! Complmentary of sigma (cmpl_sigma) cmpl_sgm = 1.0_rk - sigma ! todo add comment taup = rho * cs2 / omega ! todo add comment a12xx = SOM_neq(1) * sigma + cmpl_sgm * (-taup * SR(1)) a12yy = SOM_neq(2) * sigma + cmpl_sgm * (-taup * SR(2)) a12xy = SOM_neq(3) * sigma + cmpl_sgm * (-taup * SR(3)) ! todo add info call f_f_eq_regularized_4th_ord_d2q9( layout%weight(:), rho, u_x, u_y, & & feq, f1, a12xx, a12yy, a12xy ) ! todo rename f1 ! Relaxation do iDir = 1, QQ outState( ( ielem-1)* nscalars+ idir+(1-1)* qq) = & & feq(iDir) + cmpl_o * f1(iDir) enddo enddo nodeloop !$omp end do nowait end subroutine bgk_HybridRecursiveRegularized_d2q9 ! ****************************************************************************** ! ! ****************************************************************************** ! subroutine bgk_DualRelaxationTime_RR_d2q9 ( 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(:) ! -------------------------------------------------------------------- ! ! indeces integer :: iElem, iDir ! temporary distribution variables real(kind=rk) :: f( QQ ), SOM(3), SOM_neq(3) real(kind=rk) :: rho, u_x, u_y, a12xx, a12xy, a12yy real(kind=rk) :: omega, tau, tauN, CoefTauNTau real(kind=rk) :: feq(QQ), f1(QQ), f_temp !real(kind=rk) :: sigma integer :: denspos, velpos(3), elemOff, nScalars ! --------------------------------------------------------------------------- !cdir nodep !ibm* novector !dir$ novector denspos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) velpos(1:2) = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:2) nScalars = varSys%nScalars tauN = fieldProp(1)%fluid%DRT_tauN !$omp do schedule(static) !NEC$ ivdep !DIR$ NOVECTOR nodeloop: do iElem = 1, nSolve do iDir = 1, QQ f(iDir) = inState( neigh(( idir-1)* nelems+ ielem)+( 1-1)* qq+ nscalars*0) enddo ! element offset for auxField array elemOff = (iElem-1)*varSys%nAuxScalars ! local density rho = auxField(elemOff + denspos) ! local x-, y- and z-velocity u_x = auxField(elemOff + velpos(1)) u_y = auxField(elemOff + velpos(2)) ! non equilibrium second-order moments ! SOM_neq = SOM - SOM_eq SOM = second_order_moments_2D( f=f, layout=layout ) SOM_neq(1) = SOM(1) - rho * (cs2 + (u_x * u_x)) SOM_neq(2) = SOM(2) - rho * (cs2 + (u_y * u_y)) SOM_neq(3) = SOM(3) - rho * u_x * u_y !Relaxation coefficients ! remains constant on uniform mesh. Does loop over elements include voxels off different size? omega = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem) tau = 1.0_rk / omega CoefTauNTau = ( tau - tauN ) / ( tau * tauN ) a12xx = SOM_neq(1) a12yy = SOM_neq(2) a12xy = SOM_neq(3) call f_f_eq_regularized_2nd_ord_d2q9( layout%weight(:), rho, u_x, u_y, feq, f1, & & a12xx, a12yy, a12xy ) do iDir = 1, QQ f_temp = f(iDir) - 1.0_rk/tauN * ( f(iDir) - feq(iDir) ) outState( ( ielem-1)* nscalars+ idir+(1-1)* qq) = f_temp & & + CoefTauNTau * f1(iDir) enddo enddo nodeloop !$omp end do nowait end subroutine bgk_DualRelaxationTime_RR_d2q9 ! ****************************************************************************** ! end module mus_d2q9_module ! **************************************************************************** !