! Copyright (c) 2019-2021 Kannan Masilamani <kannan.masilamani@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. ! ****************************************************************************** ! !> This module contains function to compute eddy viscosity for !! Smagorinsky les turbulence model. !! author: Kannan Masilamani ! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de> ! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2014-2015, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. module mus_Smagorinsky_module ! treelm modules use env_module, only: rk use tem_param_module, only: div1_2, div1_3, div2_3, & & cs2, rho0, rho0Inv use tem_compileconf_module, only: vlen ! musubi modules use mus_turbulence_module, only: mus_turbulence_config_type use mus_gradData_module, only: mus_gradData_type, mus_Grad_type use mus_scheme_layout_module, only: mus_scheme_layout_type use mus_derivedQuantities_module2, only: secondMom, getEquilibrium, & & getEquilibriumIncomp implicit none private public :: mus_turbVisc_Smagorinsky_fromGradU3D public :: mus_turbVisc_Smagorinsky_fromGradU2D public :: mus_turbVisc_Smagorinsky_fromGradU3D_incomp public :: mus_turbVisc_Smagorinsky_fromGradU2D_incomp public :: mus_turbVisc_Smagorinsky_fromPreColPDF public :: mus_turbVisc_Smagorinsky_fromPreColPDF_incomp contains ! ************************************************************************** ! !> Calculate eddy viscosity with smagorinsky model for compressible model !! using gradient of velocity !! The formula is taken from https://caefn.com/openfoam/smagorinsky-sgs-model !! nu_t = C_k delta sqrt(k_sgs) !! k_sgs = ((-b+sqrt(b^2+4ac))/2a)^2 !! a = C_e/delta, b=2/3 tr(dev(Strain)), c = 2 C_k delta (dev(Strain):Strain) subroutine mus_turbVisc_Smagorinsky_fromGradU3D(turbVisc, turbConfig, & & gradData, auxField, velPos, nSolve, nAuxScalars, dxL, dtL, Grad) ! -------------------------------------------------------------------------- !> output: turbulent viscosity real(kind=rk), intent(out) :: turbVisc(:) !> Contains turbulenct coefficients type(mus_turbulence_config_type), intent(in) :: turbConfig !> gradient data type(mus_gradData_type), intent(in) :: gradData !> Auxiliary field variable array real(kind=rk), intent(in) :: auxField(:) !> position of velocity components in auxField integer, intent(in) :: velPos(3) !> Number of element to solve in this level integer, intent(in) :: nSolve !> number of scalars in auxField array integer, intent(in) :: nAuxScalars !> current level lattice element size real(kind=rk), intent(in) :: dxL !> current level lattice time step size real(kind=rk), intent(in) :: dtL !> Object that contains pointers to calculate gradients type(mus_Grad_type), intent(in) :: Grad ! -------------------------------------------------------------------------- integer :: iElem real(kind=rk) :: SR(6), devSR(3), oneThird_trSR, devSR_SR, tr_devSR real(kind=rk) :: sqrt_k_sgs, visc_coeff real(kind=rk) :: a, b, c !> gradient of velocity real(kind=rk) :: gradU(3,3,vlen) integer :: nChunks, iChunks, nChunkElems, low_bound, elemPos ! -------------------------------------------------------------------------- ! viscosity coeff visc_coeff = turbConfig%coeff%C_k * dxL / dtL ! gradU = getGradU(iElem, auxField, gradData, velPos, nAuxScalars, 3) nChunks = ceiling(real(nSolve, kind=rk) / real(vlen, kind=rk)) do iChunks = 1, nChunks ! calculate the end number of iElem loop nChunkElems = min(vlen, nSolve - ((iChunks - 1) * vlen)) low_bound = (iChunks - 1) * vlen gradU(:,:,1:nChunkElems) = Grad%U_ptr( auxField = auxField, & & gradData = gradData, & & velPos = velPos, & & nAuxScalars = nAuxScalars, & & nDims = 3, & & nSolve = nChunkElems, & & elemOffset = low_bound ) do iElem = 1, nChunkElems ! symmetric strain rate tensors SR(1) = gradU(1,1,iElem) !S_XX SR(2) = gradU(2,2,iElem) !S_YY SR(3) = gradU(3,3,iElem) !S_ZZ SR(4) = (gradU(1,2,iElem)+gradU(2,1,iElem))*0.5_rk !S_XY SR(5) = (gradU(2,3,iElem)+gradU(3,2,iElem))*0.5_rk !S_YZ SR(6) = (gradU(1,3,iElem)+gradU(3,1,iElem))*0.5_rk !S_XZ ! onethird of trace of strain rate oneThird_trSR = (SR(1) + SR(2) + SR(3))*div1_3 ! deviatoric strain rate differs from strain rate only in diagonal terms devSR(1) = SR(1) - oneThird_trSR devSR(2) = SR(2) - oneThird_trSR devSR(3) = SR(3) - oneThird_trSR ! inner product of devSR:SR devSR_SR = SR(1)*devSR(1) + SR(2)*devSR(2) + SR(3)*devSR(3) & & + 2.0_rk*( SR(4)**2 + SR(5)**2 + SR(6)**2) ! trace of deviatoric strain rate tr_devSR = devSR(1) + devSR(2) + devSR(3) ! parameters to compute subgrid scale kinetic energy a = turbConfig%coeff%C_e / dxL b = tr_devSR * div2_3 c = 2.0_rk * turbConfig%coeff%C_k * dxL * devSR_SR ! subgrid scale kinetic energy sqrt_k_sgs = (-b + sqrt(b**2 + 4.0_rk*a*c))/(2.0_rk*a) elemPos = low_bound + iElem ! subgrid scale turbulent viscosity normalized to current level turbVisc(elemPos) = visc_coeff * sqrt_k_sgs end do!iELem end do!iChunks end subroutine mus_turbVisc_Smagorinsky_fromGradU3D ! ************************************************************************** ! ! ************************************************************************** ! !> Calculate eddy viscosity with smagorinsky model for compressible model !! using gradient of velocity for 2D layout subroutine mus_turbVisc_Smagorinsky_fromGradU2D(turbVisc, turbConfig, & & gradData, auxField, velPos, nSolve, nAuxScalars, dxL, dtL, Grad) ! -------------------------------------------------------------------------- !> output: turbulent viscosity real(kind=rk), intent(out) :: turbVisc(:) !> Contains turbulenct coefficients type(mus_turbulence_config_type), intent(in) :: turbConfig !> gradient data type(mus_gradData_type), intent(in) :: gradData !> Auxiliary field variable array real(kind=rk), intent(in) :: auxField(:) !> position of velocity components in auxField integer, intent(in) :: velPos(3) !> Number of element to solve in this level integer, intent(in) :: nSolve !> number of scalars in auxField array integer, intent(in) :: nAuxScalars !> current level lattice element size real(kind=rk), intent(in) :: dxL !> current level lattice time step size real(kind=rk), intent(in) :: dtL !> Object that contains pointers to calculate gradients type(mus_Grad_type), intent(in) :: Grad ! -------------------------------------------------------------------------- integer :: iElem real(kind=rk) :: SR(3), devSR(2), oneHalf_trSR, devSR_SR, tr_devSR real(kind=rk) :: sqrt_k_sgs, visc_coeff real(kind=rk) :: a, b, c !> gradient of velocity real(kind=rk) :: gradU(2,2,vlen) integer :: nChunks, iChunks, nChunkElems, low_bound, elemPos ! -------------------------------------------------------------------------- ! viscosity coeff visc_coeff = turbConfig%coeff%C_k * dxL / dtL nChunks = ceiling(real(nSolve, kind=rk) / real(vlen, kind=rk)) do iChunks = 1, nChunks ! calculate the end number of iElem loop nChunkElems = min(vlen, nSolve - ((iChunks - 1) * vlen)) low_bound = (iChunks - 1) * vlen gradU(:,:,1:nChunkElems) = Grad%U_ptr( auxField = auxField, & & gradData = gradData, & & velPos = velPos, & & nAuxScalars = nAuxScalars, & & nDims = 2, & & nSolve = nChunkElems, & & elemOffset = low_bound ) do iElem = 1, nChunkElems ! symmetric strain rate tensors SR(1) = gradU(1,1,iElem) !S_XX SR(2) = gradU(2,2,iElem) !S_YY SR(3) = (gradU(1,2,iElem)+gradU(2,1,iElem))*0.5_rk !S_XY ! one half of trace of strain rate oneHalf_trSR = (SR(1) + SR(2))*div1_2 ! deviatoric strain rate differs from strain rate only in diagonal terms devSR(1) = SR(1) - oneHalf_trSR devSR(2) = SR(2) - oneHalf_trSR ! inner product of devSR:SR devSR_SR = SR(1)*devSR(1) + SR(2)*devSR(2) & & + 2.0_rk*( SR(3)**2 ) ! trace of deviatoric strain rate tr_devSR = devSR(1) + devSR(2) ! parameters to compute subgrid scale kinetic energy a = turbConfig%coeff%C_e / dxL b = tr_devSR c = 2.0_rk * turbConfig%coeff%C_k * dxL * devSR_SR ! subgrid scale kinetic energy sqrt_k_sgs = (-b + sqrt(b**2 + 4.0_rk*a*c))/(2.0_rk*a) elemPos = low_bound + iElem ! subgrid scale turbulent viscosity turbVisc(elemPos) = visc_coeff * sqrt_k_sgs end do end do end subroutine mus_turbVisc_Smagorinsky_fromGradU2D ! ************************************************************************** ! ! ************************************************************************** ! !> Calculate eddy viscosity with smagorinsky model for incompressible model !! using gradient of velocity subroutine mus_turbVisc_Smagorinsky_fromGradU3D_incomp(turbVisc, turbConfig, & & gradData, auxField, velPos, nSolve, nAuxScalars, dxL, dtL, Grad) ! -------------------------------------------------------------------------- !> output: turbulent viscosity real(kind=rk), intent(out) :: turbVisc(:) !> Contains turbulenct coefficients type(mus_turbulence_config_type), intent(in) :: turbConfig !> gradient data type(mus_gradData_type), intent(in) :: gradData !> Auxiliary field variable array real(kind=rk), intent(in) :: auxField(:) !> position of velocity components in auxField integer, intent(in) :: velPos(3) !> Number of element to solve in this level integer, intent(in) :: nSolve !> number of scalars in auxField array integer, intent(in) :: nAuxScalars !> current level lattice element size real(kind=rk), intent(in) :: dxL !> current level lattice time step size real(kind=rk), intent(in) :: dtL !> Object that contains pointers to calculate gradients type(mus_Grad_type), intent(in) :: Grad ! -------------------------------------------------------------------------- integer :: iElem real(kind=rk) :: SR(6), SR_mag real(kind=rk) :: k_sgs, visc_coeff !> gradient of velocity real(kind=rk) :: gradU(3,3,vlen) integer :: nChunks, iChunks, nChunkElems, low_bound, elemPos ! -------------------------------------------------------------------------- ! viscosity coeff visc_coeff = turbConfig%coeff%C_k * dxL / dtL nChunks = ceiling(real(nSolve, kind=rk) / real(vlen, kind=rk)) do iChunks = 1, nChunks ! calculate the end number of iElem loop nChunkElems = min(vlen, nSolve - ((iChunks - 1) * vlen)) low_bound = (iChunks - 1) * vlen gradU(:,:,1:nChunkElems) = Grad%U_ptr( auxField = auxField, & & gradData = gradData, & & velPos = velPos, & & nAuxScalars = nAuxScalars, & & nDims = 3, & & nSolve = nChunkElems, & & elemOffset = low_bound ) do iElem = 1, nChunkElems ! symmetric strain rate tensors SR(1) = gradU(1,1,iElem) !S_XX SR(2) = gradU(2,2,iElem) !S_YY SR(3) = gradU(3,3,iElem) !S_ZZ SR(4) = (gradU(1,2,iElem)+gradU(2,1,iElem))*0.5_rk !S_XY SR(5) = (gradU(2,3,iElem)+gradU(3,2,iElem))*0.5_rk !S_YZ SR(6) = (gradU(1,3,iElem)+gradU(3,1,iElem))*0.5_rk !S_XZ ! magnitude of strain rate tensor, sqrt(2 S:S) SR_mag = sqrt(2.0_rk*( SR(1)**2 + SR(2)**2 + SR(3)**2 & & + 2.0_rk*(SR(4)**2 + SR(5)**2 + SR(6)**2))) ! subgrid scale kinetic energy of incompressible model k_sgs = turbConfig%coeff%C_k * dxL**2 * SR_mag**2 / turbConfig%coeff%C_e elemPos = low_bound + iElem ! subgrid scale turbulent viscosity turbVisc(elemPos) = visc_coeff * sqrt(k_sgs) end do end do end subroutine mus_turbVisc_Smagorinsky_fromGradU3D_incomp ! ************************************************************************** ! ! ************************************************************************** ! !> Calculate eddy viscosity with smagorinsky model for incompressible model !! using gradient of velocity for 2D layout subroutine mus_turbVisc_Smagorinsky_fromGradU2D_incomp(turbVisc, turbConfig, & & gradData, auxField, velPos, nSolve, nAuxScalars, dxL, dtL, Grad) ! -------------------------------------------------------------------------- !> output: turbulent viscosity real(kind=rk), intent(out) :: turbVisc(:) !> Contains turbulenct coefficients type(mus_turbulence_config_type), intent(in) :: turbConfig !> gradient data type(mus_gradData_type), intent(in) :: gradData !> Auxiliary field variable array real(kind=rk), intent(in) :: auxField(:) !> position of velocity components in auxField integer, intent(in) :: velPos(3) !> Number of element to solve in this level integer, intent(in) :: nSolve !> number of scalars in auxField array integer, intent(in) :: nAuxScalars !> current level lattice element size real(kind=rk), intent(in) :: dxL !> current level lattice time step size real(kind=rk), intent(in) :: dtL !> Object that contains pointers to calculate gradients type(mus_Grad_type), intent(in) :: Grad ! -------------------------------------------------------------------------- integer :: iElem real(kind=rk) :: SR(3), SR_mag real(kind=rk) :: k_sgs, visc_coeff !> gradient of velocity real(kind=rk) :: gradU(2,2,vlen) integer :: nChunks, iChunks, nChunkElems, low_bound, elemPos ! -------------------------------------------------------------------------- ! viscosity coeff visc_coeff = turbConfig%coeff%C_k * dxL / dtL nChunks = ceiling(real(nSolve, kind=rk) / real(vlen, kind=rk)) do iChunks = 1, nChunks ! calculate the end number of iElem loop nChunkElems = min(vlen, nSolve - ((iChunks - 1) * vlen)) low_bound = (iChunks - 1) * vlen gradU(:,:,1:nChunkElems) = Grad%U_ptr( auxField = auxField, & & gradData = gradData, & & velPos = velPos, & & nAuxScalars = nAuxScalars, & & nDims = 2, & & nSolve = nChunkElems, & & elemOffset = low_bound ) do iElem = 1, nChunkElems ! symmetric strain rate tensors SR(1) = gradU(1,1,iElem) !S_XX SR(2) = gradU(2,2,iElem) !S_YY SR(3) = (gradU(1,2,iElem)+gradU(2,1,iElem))*0.5_rk !S_XY ! magnitude of strain rate tensor, sqrt(2 S:S) SR_mag = sqrt(2.0_rk*( SR(1)**2 + SR(2)**2 + 2.0_rk*(SR(3)**2) )) ! subgrid scale kinetic energy of incompressible model k_sgs = turbConfig%coeff%C_k * dxL**2 * SR_mag**2 / turbConfig%coeff%C_e elemPos = low_bound + iElem ! subgrid scale turbulent viscosity turbVisc(elemPos) = visc_coeff * sqrt(k_sgs) end do end do end subroutine mus_turbVisc_Smagorinsky_fromGradU2D_incomp ! ************************************************************************** ! ! ************************************************************************** ! !> Calculate eddy viscosity with smagorinsky model for compressible model !! using pre-collision PDF. !! Schneider, A. (2015). A Consistent Large Eddy Approach for Lattice !! Boltzmann Methods and its Application to Complex Flows. !! Technical University Kaiserslautern. subroutine mus_turbVisc_Smagorinsky_fromPreColPDF(turbVisc, turbConfig, & & state, neigh, auxField, densPos, velPos, nSize, nSolve, nScalars, & & nAuxScalars, layout, dxL, dtL, viscKine) ! -------------------------------------------------------------------------- !> output: turbulent viscosity real(kind=rk), intent(out) :: turbVisc(:) !> Contains turbulenct coefficients type(mus_turbulence_config_type), intent(in) :: turbConfig !> state array real(kind=rk), intent(in) :: state(:) !> neigh array to obtain precollision pdf integer, intent(in) :: neigh(:) !> Auxiliary field variable array real(kind=rk), intent(in) :: auxField(:) !> position of density in auxField integer, intent(in) :: densPos !> position of velocity components in auxField integer, intent(in) :: velPos(3) !> number of elements in state array integer, intent(in) :: nSize !> Number of element to solve in this level integer, intent(in) :: nSolve !> number of scalars in state array integer, intent(in) :: nScalars !> number of scalars in auxField array integer, intent(in) :: nAuxScalars !> scheme layout type(mus_scheme_layout_type), intent(in) :: layout !> current level lattice element size real(kind=rk), intent(in) :: dxL !> current level lattice time step size real(kind=rk), intent(in) :: dtL !> Background kinematic viscosity in lattice divided by dtL real(kind=rk), intent(in) :: viscKine(:) ! -------------------------------------------------------------------------- integer :: iElem, iDir, QQ, elemOff real(kind=rk) :: visc_coeff, rho, inv_rho, vel(3) !> precollision PDF real(kind=rk) :: f_preCol(layout%fStencil%QQ) real(kind=rk) :: fEq(layout%fStencil%QQ), nEq(layout%fStencil%QQ) real(kind=rk) :: nEqTens(6), nEqTensMag, viscKineTerm, nEqTensTerm ! -------------------------------------------------------------------------- QQ = layout%fStencil%QQ ! viscosity coeff visc_coeff = (turbConfig%coeff%C_s * dxL )**2 do iElem = 1, nSolve ! Get pre-collisiton PDF do iDir = 1, QQ f_preCol(iDir) = state ( & & neigh((idir-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0) end do elemOff = (iElem-1)*nAuxScalars ! density rho = auxField( elemOff + densPos) inv_rho = 1.0_rk/rho ! velocity vel(1) = auxField( elemOff + velPos(1) ) vel(2) = auxField( elemOff + velPos(2) ) vel(3) = auxField( elemOff + velPos(3) ) ! Calculate the equilibrium distribution function fEq(:) = getEquilibrium( rho, vel, layout) ! Calculate the non-equilibrium part nEq(:) = f_preCol(:) - fEq(:) ! Now calculate the symmetric deviatoric second-order tensor of ! nonEquilibrium part ! the static part cs2 I is usually neglected for weakly compressible flows ! however, in current implementation it is considered nEqTens(1) = sum( (layout%fStencil%cxcx(1,:) - cs2) * nEq) !XX nEqTens(2) = sum( (layout%fStencil%cxcx(2,:) - cs2) * nEq) !YY nEqTens(3) = sum( (layout%fStencil%cxcx(3,:) - cs2) * nEq) !ZZ nEqTens(4) = sum( (layout%fStencil%cxcx(4,:) ) * nEq) !XY nEqTens(5) = sum( (layout%fStencil%cxcx(5,:) ) * nEq) !YZ nEqTens(6) = sum( (layout%fStencil%cxcx(6,:) ) * nEq) !XZ ! magnitude of second-order tensor nEqTensMag = sqrt(nEqTens(1)**2 + nEqTens(2)**2 + nEqTens(3)**2 & & + 2.0_rk*(nEqTens(4)**2 + nEqTens(5)**2 + nEqTens(6)**2) ) ! turbulent viscosity !nu_t = (sqrt((v+cs2 dt/2)^2+2(Cs dx)^2|nEQ|/rho) - (v+cs2 dt/2))/2 ! viscKine is scaled to current level dt so multiply viscKine with dt ! to get unit consistent viscKineTerm = (viscKine(iElem) + div1_2 * cs2) * dtL nEqTensTerm = 2.0_rk * visc_coeff * nEqTensMag * inv_rho turbVisc(iElem) = ( sqrt(viscKineTerm**2 + nEqTensTerm) & & - viscKineTerm ) * div1_2 / dtL end do end subroutine mus_turbVisc_Smagorinsky_fromPreColPDF ! ************************************************************************** ! ! ************************************************************************** ! !> Calculate eddy viscosity with smagorinsky model for incompressible model !! using pre-collision PDF subroutine mus_turbVisc_Smagorinsky_fromPreColPDF_incomp(turbVisc, & & turbConfig, state, neigh, auxField, densPos, velPos, nSize, nSolve, & & nScalars, nAuxScalars, layout, dxL, dtL, viscKine) ! -------------------------------------------------------------------------- !> output: turbulent viscosity real(kind=rk), intent(out) :: turbVisc(:) !> Contains turbulenct coefficients type(mus_turbulence_config_type), intent(in) :: turbConfig !> state array real(kind=rk), intent(in) :: state(:) !> neigh array to obtain precollision pdf integer, intent(in) :: neigh(:) !> Auxiliary field variable array real(kind=rk), intent(in) :: auxField(:) !> position of density in auxField integer, intent(in) :: densPos !> position of velocity components in auxField integer, intent(in) :: velPos(3) !> number of elements in state array integer, intent(in) :: nSize !> Number of element to solve in this level integer, intent(in) :: nSolve !> number of scalars in state array integer, intent(in) :: nScalars !> number of scalars in auxField array integer, intent(in) :: nAuxScalars !> scheme layout type(mus_scheme_layout_type), intent(in) :: layout !> current level lattice element size real(kind=rk), intent(in) :: dxL !> current level lattice time step size real(kind=rk), intent(in) :: dtL !> Background kinematic viscosity in lattice divided by dtL real(kind=rk), intent(in) :: viscKine(:) ! -------------------------------------------------------------------------- integer :: iElem, iDir, QQ, elemOff real(kind=rk) :: visc_coeff, rho, vel(3) !> precollision PDF real(kind=rk) :: f_preCol(layout%fStencil%QQ) real(kind=rk) :: fEq(layout%fStencil%QQ), nEq(layout%fStencil%QQ) real(kind=rk) :: nEqTens(6), nEqTensMag, viscKineTerm, nEqTensTerm ! -------------------------------------------------------------------------- QQ = layout%fStencil%QQ ! viscosity coeff visc_coeff = (turbConfig%coeff%C_s * dxL )**2 do iElem = 1, nSolve ! Get pre-collisiton PDF do iDir = 1, QQ f_preCol(iDir) = state ( & & neigh((idir-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0) end do elemOff = (iElem-1)*nAuxScalars ! density rho = auxField( elemOff + densPos) ! velocity vel(1) = auxField( elemOff + velPos(1) ) vel(2) = auxField( elemOff + velPos(2) ) vel(3) = auxField( elemOff + velPos(3) ) ! Calculate the equilibrium distribution function fEq(:) = getEquilibriumIncomp( rho, vel, layout, rho0) ! Calculate the non-equilibrium part nEq(:) = f_preCol(:) - fEq(:) ! Now calculate the symmetric second-order tensor of nonEquilibrium part nEqTens = secondMom(layout%fStencil%cxcx, nEq, layout%fStencil%QQ) ! magnitude of second-order tensor nEqTensMag = sqrt(nEqTens(1)**2 + nEqTens(2)**2 + nEqTens(3)**2 & & + 2.0_rk*(nEqTens(4)**2 + nEqTens(5)**2 + nEqTens(6)**2) ) ! turbulent viscosity !nu_t = (sqrt((v+cs2 dt/2)^2+2(Cs dx)^2|nEQ|/rho) - (v+cs2 dt/2))/2 ! viscKine is scaled to current level dt so multiply viscKine with dt ! to get unit consitent viscKineTerm = (viscKine(iElem) + div1_2 * cs2) * dtL ! use reference density for incompressible model nEqTensTerm = 2.0_rk * visc_coeff * nEqTensMag * rho0Inv turbVisc(iElem) = ( sqrt(viscKineTerm**2 + nEqTensTerm) & & - viscKineTerm ) * div1_2 / dtL end do end subroutine mus_turbVisc_Smagorinsky_fromPreColPDF_incomp ! ************************************************************************** ! end module mus_Smagorinsky_module