! Copyright (c) 2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 2020 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2021-2022 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! **************************************************************************** ! !> author: Kannan Masilamani !! Non-equilibrium extrapolatioon boundary condition treatment routines for !! fluid simulation. The boundary condition for straight boundary with qVal=0 !! is based on: Guo, Z., & Shi, B. (2002). "Non-equilibrium extrapolation method !! for velocity and pressure boundary conditions in the lattice Boltzmann !! method." Chinese Physics. !! The boundary condition for curved boundary is based on !! Guo, Z., Shi, B., & Zheng, C. (2002). An extrapolation method for boundary !! conditions in lattice Boltzmann method. Physics of Fluids, 14(May), 10–14. !! !! Details on implementation can be found in [[tem_bc_module]] ! 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_bc_fluid_nonEqExpol_module ! include treelm modules use env_module, only: rk use tem_param_module, only: div1_3, div4_3, csInv, cs, cs2, cs2inv, & & rho0, rho0Inv use tem_time_module, only: tem_time_type use treelmesh_module, only: treelmesh_type use tem_varSys_module, only: tem_varSys_type use tem_construction_module, only: tem_levelDesc_type ! include musubi modules use mus_bc_header_module, only: boundary_type, & & glob_boundary_type use mus_scheme_layout_module, only: mus_scheme_layout_type use mus_field_prop_module, only: mus_field_prop_type use mus_derVarPos_module, only: mus_derVarPos_type use mus_physics_module, only: mus_physics_type use mus_mixture_module, only: mus_mixture_type use mus_derivedQuantities_module2, only: getVelocity, & & getVelocity_incomp implicit none private public :: velocity_nonEqExpol public :: velocity_nonEqExpol_incomp public :: velocity_nonEqExpol_curved, velocity_nonEqExpol_curved_incomp public :: pressure_nonEqExpol public :: pressure_nonEqExpol_incomp contains ! **************************************************************************** ! !> Element-wise Dirichlet velocity non-equilibrium boundary condition for !! straight boundary to update all directions. !! For straight wall, values are extrapolated along !! boundary normal instead of along the link and qVal =0.0 for straight wall. !! !! Notation: b (fictious boundary) w (physical bolundary or surface) !! f (local element) ff (first neighbor fluid element) !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'inlet', !! kind = 'velocity_nonEqExpol', !! velocity = 'inlet_vel, !! } !!} !!variable = { !! name = 'inlet_vel', !! ncomponents = 3, !! vartype = 'st_fun', !! st_fun = {0.06, 0.0, 0.0} !!} !!``` !! This is described in the paper: !! Non-equilibrium extrapolation method !! for velocity and pressure boundary conditions in the lattice Boltzmann !! method." Chinese Physics. !! !! More informations concerning the nonEqExpol can be found in the !! corresponding subroutine. !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine velocity_nonEqExpol( me, state, bcBuffer, globBC, levelDesc, & & tree, nSize, iLevel, sim_time, neigh, & & layout, fieldProp, varPos, nScalars, varSys, & & derVarPos, physics, iField, mixture ) ! -------------------------------------------------------------------- ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! -------------------------------------------------------------------- ! ! also determined in mus_bc_header_module integer :: iDir ! variables for fictious boundary element real(kind=rk) :: feq_b, vel_b(3) ! variables for surface (density and link-wise velocity) real(kind=rk) :: vel_w(globBC%nElems(iLevel)*3) ! variables for overnext fluid element real(kind=rk) :: feq_ff, rho_ff, rhoInv, vel_ff(3) real(kind=rk) :: inv_vel ! temporary local pdf values and the ones of overnext fluid real(kind=rk) :: pdf_ff(layout%fStencil%QQ) integer :: bcVel_pos, QQ, iElem, elemPos ! -------------------------------------------------------------------- ! QQ = layout%fstencil%QQ inv_vel = 1._rk / physics%fac(iLevel)%vel ! position of boundary velocity in varSys bcVel_pos = me%bc_states%velocity%varPos ! get velocity_phy on boundary call varSys%method%val(bcVel_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%velocity%pntIndex & & %indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = vel_w ) ! convert physical velocity into LB velocity vel_w = vel_w * inv_vel ! Update all directions of boundary elements do iElem = 1, globBC%nElems(iLevel) ! determine needed quantities of the overnext fluid neighbor element ! x_ff pdf_ff = me%neigh(iLevel)%neighBufferPre_nNext(1, (iElem-1) * QQ + 1 & & : (iElem-1) * QQ + QQ ) ! calculate density rho_ff = sum(pdf_ff) rhoInv = 1.0_rk / rho_ff ! calulate velocity vel_ff(1) = sum( pdf_ff * layout%fstencil%cxdirrk(1, :) ) * rhoinv vel_ff(2) = sum( pdf_ff * layout%fstencil%cxdirrk(2, :) ) * rhoinv vel_ff(3) = sum( pdf_ff * layout%fstencil%cxdirrk(3, :) ) * rhoinv ! determine needed quantities of fictious boundary (_b) ! compute velocity on boundary (eq.18) ! qVal = 0.0, boundary node is one lattice away from fluid node! vel_b = vel_w( (iElem-1)*3+1:iElem*3) ! position of this boundary element in state array elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) do iDir = 1, layout%fStencil%QQ ! derive feq from velocity and density for all direction ! calculate equilibrium density feq_ff = layout%weight( idir ) & & * rho_ff & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_ff(:)*vel_ff(:)) * 0.5_rk * cs2inv ) ) ! extrapolate density on boundary ! derive feq from velocity and density per direction ! calculate equilibrium density feq_b = layout%weight( idir ) & & * rho_ff & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_b(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_b(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_b(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_b(:)*vel_b(:)) * 0.5_rk * cs2inv ) ) ! Update pre-collision PDF state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_b + (pdf_ff(iDir) - fEq_ff) end do end do !iElem end subroutine velocity_nonEqExpol ! **************************************************************************** ! ! **************************************************************************** ! !> Linkwise Dirichlet velocity non-equilibrium boundary condition for !! curved using the subroutine "mus_set_nonEqExpol". !! For curved wall, values are extrapolated along !! element normal !! !! Notation: b (fictious boundary) w (physical bolundary or surface) !! f (local element) ff (overnext fluidneighbor) !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'inlet', !! kind = 'velocity_nonEqExpol', !! velocity = 'inlet_vel, !! curved = true !! } !!} !!variable = { !! name = 'inlet_vel', !! ncomponents = 3, !! vartype = 'st_fun', !! st_fun = {0.06, 0.0, 0.0} !!} !!``` !! This is described in the paper: !! Guo, Z.; Zheng, C. & Shi B. (2002). An extrapolation method for boundary !! conditions in lattice Boltzmann method. !! Physics of Fluids 14, 2007 (2002); https://doi.org/10.1063/1.1471914 !! !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine velocity_nonEqExpol_curved( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, & & sim_time, neigh, layout, fieldProp, & & varPos, nScalars, varSys, derVarPos, & & physics, iField, mixture ) ! -------------------------------------------------------------------- ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! -------------------------------------------------------------------- ! ! coefficients which are calculated in mus_bc_header_module real(kind=rk) :: c_w, c_f, c_ff, c_neq_f, c_neq_ff ! also determined in mus_bc_header_module integer :: iDir, posInBuffer, posInNeighBuf ! variables for fictious boundary element real(kind=rk) :: feq_b, fneq_b, rho_b, vel_b(3) ! variables for surface (density and link-wise velocity) real(kind=rk) :: vel_w(me%links(iLevel)%nVals*3) ! variables for fluid element real(kind=rk) :: feq_f, fneq_f, rho_f, vel_f(3) ! variables for overnext fluid element real(kind=rk) :: feq_ff, fneq_ff, rho_ff, vel_ff(3) real(kind=rk) :: inv_vel ! temporary local pdf values and the ones of overnext fluid real(kind=rk) :: pdfTmp(layout%fStencil%QQ) real(kind=rk) :: pdf_ff(layout%fStencil%QQ) integer :: bcVel_pos, iLink, QQ, offset ! relaxation parameter !real(kind=rk) :: omega ! -------------------------------------------------------------------- ! QQ = layout%fstencil%QQ inv_vel = 1._rk / physics%fac(iLevel)%vel ! position of boundary velocity in varSys bcVel_pos = me%bc_states%velocity%varPos ! get velocity_phy on boundary (surface) call varSys%method%val(bcVel_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%velocity%pntIndex & & %indexLvl(ilevel) & & %val(1:me%links(iLevel)%nVals), & & nVals = me%links(iLevel)%nVals, & & res = vel_w ) ! convert physical velocity into LB velocity vel_w = vel_w * inv_vel ! loop over links do iLink = 1, me%links(iLevel)%nVals ! load pre-calculated coefficients from subroutine c_w = me%nonEqExpol(iLevel)% c_w(iLink) c_f = me%nonEqExpol(iLevel)% c_f(iLink) c_ff = me%nonEqExpol(iLevel)% c_ff(iLink) c_neq_f = me%nonEqExpol(iLevel)% c_neq_f(iLink) c_neq_ff = me%nonEqExpol(iLevel)%c_neq_ff(iLink) ! set link-wise direction iDir = me%nonEqExpol(iLevel)%iDir(iLink) ! determine needed quantities of the current element (_f) ! calculate density posInBuffer = me%nonEqExpol(iLevel)%posInBuffer(iLink) pdfTmp(1:QQ) = bcBuffer( (posInBuffer - 1) * nScalars + varPos(1) & & : ( posInBuffer - 1) * nScalars + varPos(QQ) ) rho_f = sum(pdfTmp) ! calulate velocity (facilitate it by using a pure function) vel_f = getVelocity( subset = pdfTmp, & stencil = layout%fStencil, & varPos = varPos ) ! derive feq from velocity and density per direction ! calculate equilibrium density feq_f = layout%weight( idir ) & & * rho_f & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_f(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_f(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_f(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_f(:)*vel_f(:)) * 0.5_rk * cs2inv ) ) ! determine needed quantities of the overnext fluid neighbor element (_ff) ! calculate density posInNeighBuf = me%nonEqExpol(iLevel)%posInNeighBuf(iLink) pdf_ff = me%neigh(iLevel)%neighBufferPost( & & 1, & & (posInNeighBuf - 1) * QQ + 1:(posInNeighBuf - 1) * QQ + QQ ) rho_ff = sum(pdf_ff) ! calulate velocity (facilitate it by using a pure function) vel_ff = getVelocity( subset = pdf_ff, & & stencil = layout%fStencil, & varPos = varPos ) ! derive feq from velocity and density per direction ! calculate equilibrium density feq_ff = layout%weight( idir ) & & * rho_ff & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_ff(:)*vel_ff(:)) * 0.5_rk * cs2inv ) ) ! need offset for velocity offset = (iLink-1)*3 ! determine needed quantities of fictious boundary (_b) ! compute velocity on boundary (eq.18) vel_b = (c_w * vel_w(offset+1:offset+3)) + (c_f * vel_f) + (c_ff * vel_ff) ! extrapolate density on boundary !rho_b = div4_3*rho_f - div1_3*rho_ff rho_b = rho_f ! derive feq from velocity and density per direction ! calculate equilibrium density feq_b = layout%weight( idir ) & & * rho_b & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_b(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_b(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_b(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_b(:)*vel_b(:)) * 0.5_rk * cs2inv ) ) ! use pdf_b = fEq_b + fneq_b (eq.16) to determine ! non-equilibrium components (eq.19) fneq_f = pdfTmp(iDir) - feq_f fneq_ff = pdf_ff(iDir) - feq_ff fneq_b = c_nEq_f * fneq_f + c_nEq_ff * fneq_ff ! write into state ! fneq is computed from post-collision pdf so no need to multiply ! (1-omega) as in the literature state( me%links(iLevel)%val(iLink) ) = feq_b + fneq_b end do !iLink end subroutine velocity_nonEqExpol_curved ! **************************************************************************** ! ! **************************************************************************** ! !> Element-wise Dirichlet velocity non-equilibrium boundary condition for !! straight boundary to update all directions. !! For straight wall, values are extrapolated along !! boundary normal instead of along the link and qVal=0. !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine velocity_nonEqExpol_incomp( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, & & sim_time, neigh, layout, fieldProp, & & varPos, nScalars, varSys, derVarPos, & & physics, iField, mixture ) ! -------------------------------------------------------------------- ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! -------------------------------------------------------------------- ! ! also determined in mus_bc_header_module integer :: iDir ! variables for fictious boundary element real(kind=rk) :: feq_b, vel_b(3) ! variables for surface (density and link-wise velocity) real(kind=rk) :: vel_w(globBC%nElems(iLevel)*3) ! variables for overnext fluid element real(kind=rk) :: feq_ff, rho_ff, vel_ff(3) real(kind=rk) :: inv_vel ! temporary local pdf values real(kind=rk) :: pdf_ff(layout%fStencil%QQ) integer :: bcVel_pos, QQ, iElem, elemPos ! -------------------------------------------------------------------- ! QQ = layout%fstencil%QQ inv_vel = 1._rk / physics%fac(iLevel)%vel ! position of boundary velocity in varSys bcVel_pos = me%bc_states%velocity%varPos ! get velocity_phy on boundary call varSys%method%val(bcVel_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%velocity%pntIndex & & %indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = vel_w ) ! convert physical velocity into LB velocity vel_w = vel_w * inv_vel ! Update all directions of boundary elements do iElem = 1, globBC%nElems(iLevel) ! determine needed quantities of the overnext fluid neighbor element ! x_ff pdf_ff = me%neigh(iLevel)%neighBufferPre_nNext(1, (iElem-1) * QQ + 1 & & : (iElem-1) * QQ + QQ ) ! calculate density rho_ff = sum(pdf_ff) ! calulate velocity vel_ff(1) = sum( pdf_ff * layout%fstencil%cxdirrk(1, :) ) * rho0inv vel_ff(2) = sum( pdf_ff * layout%fstencil%cxdirrk(2, :) ) * rho0inv vel_ff(3) = sum( pdf_ff * layout%fstencil%cxdirrk(3, :) ) * rho0inv ! determine needed quantities of fictious boundary (_b) ! compute velocity on boundary (eq.18) ! qVal = 0.0, boundary node is one lattice away from fluid node! vel_b = vel_w( (iElem-1)*3+1:iElem*3) ! position of this boundary element in state array elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) do iDir = 1, layout%fStencil%QQ ! derive feq from velocity and density for all direction feq_ff = layout%weight( idir ) & & * ( rho_ff + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_ff(:) * vel_ff(:)) * 0.5_rk * cs2inv ) ) ! extrapolate density on boundary ! derive feq from velocity and density per direction feq_b = layout%weight( idir ) & & * ( rho_ff + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_b(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_b(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_b(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_b(:) * vel_b(:)) * 0.5_rk * cs2inv ) ) ! Update pre-collision PDF state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_b + (pdf_ff(iDir) - fEq_ff) end do end do !iElem end subroutine velocity_nonEqExpol_incomp ! **************************************************************************** ! ! **************************************************************************** ! !> Linkwise Dirichlet velocity non-equilibrium boundary condition for !! curved wal using the subroutine "mus_set_nonEqExpol" for incompressible !! model. !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine velocity_nonEqExpol_curved_incomp( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, sim_time, neigh, layout, & & fieldProp, varPos, nScalars, varSys, derVarPos, physics, & & iField, mixture ) ! -------------------------------------------------------------------- ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! -------------------------------------------------------------------- ! ! coefficients which are calculated in mus_bc_header_module real(kind=rk) :: c_w, c_f, c_ff, c_neq_f, c_neq_ff ! also determined in mus_bc_header_module integer :: iDir, posInBuffer, posInNeighBuf ! variables for fictious boundary element real(kind=rk) :: feq_b, fneq_b, rho_b, vel_b(3) ! variables for surface (density and link-wise velocity) real(kind=rk) :: vel_w(me%links(iLevel)%nVals*3) ! variables for fluid element real(kind=rk) :: feq_f, fneq_f, rho_f, vel_f(3) ! variables for overnext fluid element real(kind=rk) :: feq_ff, fneq_ff, rho_ff, vel_ff(3) real(kind=rk) :: inv_vel ! temporary local pdf values and the ones of overnext fluid real(kind=rk) :: pdfTmp(layout%fStencil%QQ) real(kind=rk) :: pdf_ff(layout%fStencil%QQ) integer :: bcVel_pos, iLink, QQ, offset ! -------------------------------------------------------------------- ! QQ = layout%fstencil%QQ inv_vel = 1._rk / physics%fac(iLevel)%vel ! position of boundary velocity in varSys bcVel_pos = me%bc_states%velocity%varPos ! get velocity_phy on boundary (surface) call varSys%method%val(bcVel_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%velocity%pntIndex & & %indexLvl(ilevel) & & %val(1:me%links(iLevel)%nVals), & & nVals = me%links(iLevel)%nVals, & & res = vel_w ) ! convert physical velocity into LB velocity vel_w = vel_w * inv_vel ! loop over links do iLink = 1, me%links(iLevel)%nVals ! load pre-calculated coefficients from subroutine c_w = me%nonEqExpol(iLevel)% c_w(iLink) c_f = me%nonEqExpol(iLevel)% c_f(iLink) c_ff = me%nonEqExpol(iLevel)% c_ff(iLink) c_neq_f = me%nonEqExpol(iLevel)% c_neq_f(iLink) c_neq_ff = me%nonEqExpol(iLevel)%c_neq_ff(iLink) ! set link-wise direction iDir = me%nonEqExpol(iLevel)%iDir(iLink) ! determine needed quantities of the current element (_f) ! calculate density posInBuffer = me%nonEqExpol(iLevel)%posInBuffer(iLink) pdfTmp(1:QQ) = bcBuffer( (posInBuffer - 1) * nScalars + varPos(1) & & : (posInBuffer - 1) * nScalars + varPos(QQ) ) ! calulcate density rho_f = sum(pdfTmp) ! calulate velocity (facilitate it by using a pure function) vel_f = getVelocity_incomp( subset = pdfTmp, & stencil = layout%fStencil, & varPos = varPos ) ! derive feq from velocity and density per direction (using pure funct.) feq_f = layout%weight( idir ) & & * ( rho_f + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_f(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_f(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_f(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_f(:) * vel_f(:)) * 0.5_rk * cs2inv ) ) ! determine needed quantities of the overnext fluid neighbor element (_ff) ! calculate density posInNeighBuf = me%nonEqExpol(iLevel)%posInNeighBuf(iLink) pdf_ff = me%neigh(iLevel)%neighBufferPost( & & 1, & & (posInNeighBuf - 1) * QQ + 1:(posInNeighBuf - 1) * QQ + QQ ) rho_ff = sum(pdf_ff) ! calulate velocity (facilitate it by using a pure function) vel_ff = getVelocity_incomp( subset = pdf_ff, & & stencil = layout%fStencil, & varPos = varPos ) ! derive feq from velocity and density per direction (using pure funct.) feq_ff = layout%weight( idir ) & & * ( rho_ff + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_ff(:) * vel_ff(:)) * 0.5_rk * cs2inv ) ) ! need offset for velocity offset = (iLink-1)*3 ! determine needed quantities of fictious boundary (_b) ! compute velocity on boundary (eq.18) vel_b = (c_w * vel_w(offset+1:offset+3)) + (c_f * vel_f) + (c_ff * vel_ff) ! extrapolate density on boundary !rho_b = div4_3*rho_f - div1_3*rho_ff rho_b = rho_f ! derive feq from velocity and density per direction (using pure funct.) feq_b = layout%weight( idir ) & & * ( rho_b + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_b(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_b(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_b(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_b(:) * vel_b(:)) * 0.5_rk * cs2inv ) ) ! use pdf_b = fEq_b + fneq_b (eq.16) to determine ! non-equilibrium components (eq.19) fneq_f = pdfTmp(iDir) - feq_f fneq_ff = pdf_ff(iDir) - feq_ff fneq_b = (c_nEq_f * fneq_f) + (c_nEq_ff * fneq_ff) ! write into state ! KM: Here dividing nonequilibrium by 1-omega does not work state( me%links(iLevel)%val(iLink) ) = feq_b + fneq_b end do !iLink end subroutine velocity_nonEqExpol_curved_incomp ! **************************************************************************** ! ! **************************************************************************** ! !> Element-wise Dirichlet pressure non-equilibrium boundary condition for !! straight boundary and updates all directions. !! For straight wall, values are extrapolated along !! boundary normal instead of along the link and qVal=0. !! !! Notation: b (fictious boundary) w (physical bolundary or surface) !! f (local element) ff (overnext fluidneighbor) !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'pressure_nonEqExpol', !! pressure = 'press, !! curved = true !! } !!} !!variable = { !! name = 'press', !! ncomponents = 1, !! vartype = 'st_fun', !! st_fun = 0.0 !!} !!``` !! It is based on the paper: !! Guo, Z.; Zheng, C. & Shi B. (2002). An extrapolation method for boundary !! conditions in lattice Boltzmann method. !! Physics of Fluids 14, 2007 (2002); https://doi.org/10.1063/1.1471914 !! !! More informations concerning the nonEqExpol can be found in the !! corresponding subroutine. !! !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine pressure_nonEqExpol( me, state, bcBuffer, globBC, levelDesc, & & tree, nSize, iLevel, sim_time, neigh, & & layout, fieldProp, varPos, nScalars, varSys, & & derVarPos, physics, iField, mixture ) ! -------------------------------------------------------------------- ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! -------------------------------------------------------------------- ! ! also determined in mus_bc_header_module integer :: iDir ! variables for fictious boundary element real(kind=rk) :: feq_b, rho_b ! variables for surface (density and link-wise velocity) real(kind=rk) :: rho_w(globBC%nElems(iLevel)) ! variables for overnext fluid element real(kind=rk) :: feq_ff, rho_ff, rhoInv, vel_ff(3) ! temporary local pdf values and the ones of overnext fluid real(kind=rk) :: pdf_ff(layout%fStencil%QQ) integer :: bcPress_pos, QQ, iElem, elemPos ! relaxation parameter ! real(kind=rk) :: omega real(kind=rk) :: inv_rho_phy ! -------------------------------------------------------------------- ! QQ = layout%fstencil%QQ inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv ! position of boundary velocity in varSys bcPress_pos = me%bc_states%pressure%varPos ! get pressure_phy on boundary (surface) and store it to rho_w call varSys%method%val(bcPress_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%pressure%pntIndex & & %indexLvl(ilevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = rho_w ) ! convert physical physical pressure into LB density rho_w = rho_w * inv_rho_phy ! Update all directions of boundary elements do iElem = 1, globBC%nElems(iLevel) ! determine needed quantities of the overnext fluid neighbor element ! x_ff pdf_ff = me%neigh(iLevel)%neighBufferPre_nNext(1, (iElem-1) * QQ + 1 & & : (iElem-1) * QQ + QQ ) ! calculate density rho_ff = sum(pdf_ff) rhoInv = 1.0_rk / rho_ff ! calulate velocity vel_ff(1) = sum( pdf_ff * layout%fstencil%cxdirrk(1, :) ) * rhoinv vel_ff(2) = sum( pdf_ff * layout%fstencil%cxdirrk(2, :) ) * rhoinv vel_ff(3) = sum( pdf_ff * layout%fstencil%cxdirrk(3, :) ) * rhoinv ! density on boundary, qVal=0: boundary and fluid node overlaps rho_b = rho_w(iElem) ! position of this boundary element in state array elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) do iDir = 1, layout%fStencil%QQ ! derive feq from velocity and density for all direction ! calculate equilibrium density feq_ff = layout%weight( idir ) & & * rho_ff & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_ff(:)*vel_ff(:)) * 0.5_rk * cs2inv ) ) ! extrapolate velocity on boundary ! derive feq from velocity and density per direction ! calculate equilibrium density feq_b = layout%weight( idir ) & & * rho_b & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_ff(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_ff(:)*vel_ff(:)) * 0.5_rk * cs2inv ) ) ! Update pre-collision PDF state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_b + (pdf_ff(iDir) - fEq_ff) end do end do !iElem end subroutine pressure_nonEqExpol ! **************************************************************************** ! ! **************************************************************************** ! !> Element-wise Dirichlet pressure non-equilibrium boundary condition for !! straight boundary and update all directions. !! For straight wall, values are extrapolated along !! boundary normal instead of along the link and qVal=0. !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine pressure_nonEqExpol_incomp( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, & & sim_time, neigh, layout, fieldProp, & & varPos, nScalars, varSys, derVarPos, & & physics, iField, mixture ) ! -------------------------------------------------------------------- ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! -------------------------------------------------------------------- ! ! also determined in mus_bc_header_module integer :: iDir ! variables for fictious boundary element real(kind=rk) :: feq_b, rho_b ! variables for surface (density and link-wise velocity) real(kind=rk) :: rho_w(globBC%nElems(iLevel)) ! variables for overnext fluid element real(kind=rk) :: feq_ff, rho_ff, vel_ff(3) ! temporary local pdf values and the ones of overnext fluid real(kind=rk) :: pdf_ff(layout%fStencil%QQ) integer :: bcPress_pos, QQ, iElem, elemPos ! relaxation parameter real(kind=rk) :: inv_rho_phy ! -------------------------------------------------------------------- ! QQ = layout%fstencil%QQ inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv ! position of boundary velocity in varSys bcPress_pos = me%bc_states%pressure%varPos ! get pressure_phy on boundary (surface) and store it to rho_w call varSys%method%val(bcPress_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%pressure%pntIndex & & %indexLvl(ilevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = rho_w ) ! convert physical physical pressure into LB density rho_w = rho_w * inv_rho_phy ! Update all directions of boundary elements do iElem = 1, globBC%nElems(iLevel) ! determine needed quantities of the overnext fluid neighbor element ! x_ff pdf_ff = me%neigh(iLevel)%neighBufferPre_nNext(1, (iElem-1) * QQ + 1 & & : (iElem-1) * QQ + QQ ) ! calculate density rho_ff = sum(pdf_ff) ! calulate velocity vel_ff(1) = sum( pdf_ff * layout%fstencil%cxdirrk(1, :) ) * rho0inv vel_ff(2) = sum( pdf_ff * layout%fstencil%cxdirrk(2, :) ) * rho0inv vel_ff(3) = sum( pdf_ff * layout%fstencil%cxdirrk(3, :) ) * rho0inv ! density on boundary, qVal=0: boundary and fluid node overlaps rho_b = rho_w(iElem) ! position of this boundary element in state array elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) do iDir = 1, layout%fStencil%QQ ! derive feq from velocity and density for all direction feq_ff = layout%weight( idir ) & & * ( rho_ff + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_ff(:) * vel_ff(:)) * 0.5_rk * cs2inv ) ) ! extrapolate velocity on boundary ! derive feq from velocity and density per direction feq_b = layout%weight( idir ) & & * ( rho_b + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_ff(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_ff(:) * vel_ff(:)) * 0.5_rk * cs2inv ) ) ! Update pre-collision PDF state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_b + (pdf_ff(iDir) - fEq_ff) end do end do !iElem end subroutine pressure_nonEqExpol_incomp ! **************************************************************************** ! end module mus_bc_fluid_nonEqExpol_module