! Copyright (c) 2012-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2012-2014 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2017 Raphael Haupt <raphael.haupt@uni-siegen.de> ! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de> ! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ****************************************************************************** ! !> Boundary condition treatment routines for multispecies simulation !! module mus_bc_species_module use iso_c_binding, only: c_f_pointer ! include treelm modules use env_module, only: rk, labelLen use tem_param_module, only: div1_3, div1_36, div1_8, div3_4h, & & div1_4, div3_8, div3_2, div9_16, div3_16,& & cs2inv, div1_6, cs4inv, & & t2cs2inv, t2cs4inv, div1_18, cs2 use tem_time_module, only: tem_time_type use treelmesh_module, only: treelmesh_type use tem_varSys_module, only: tem_varSys_type use tem_math_module, only: invert_matrix use tem_spatial_module, only: tem_spatial_for use tem_spacetime_fun_module, only: tem_spacetime_for use tem_isNaN_module, only: tem_isNaN use tem_debug_module, only: dbgUnit use tem_property_module, only: prp_solid use tem_construction_module, only: tem_levelDesc_type use tem_stencil_module, only: tem_stencilHeader_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_param_module, only: mus_param_type use mus_physics_module, only: mus_physics_type use mus_mixture_module, only: mus_mixture_type use mus_species_module, only: mus_species_type use mus_eNRTL_module, only: mus_calc_thermFactor, & & mus_calc_MS_DiffMatrix use mus_varSys_module, only: mus_varSys_data_type use mus_derQuanMSLiquid_module, only: momentumFromMacroLSE implicit none private public :: spc_inlet_eq, spc_inlet, spc_vel_bb public :: spc_outlet_zero_prsgrd public :: spc_outlet_eq, spc_outlet_vel, spc_outlet_expol public :: spc_moleFrac, spc_moleDiff_Flux public :: spc_moleFlux, spc_moleFlux_eq public :: spc_moleDens_eq public :: spc_blackbox_mem_ion, spc_blackbox_mem_solvent public :: spc_moleFrac_wtdf, spc_moleFrac_eq public :: spc_outflow, spc_inflow public :: spc_solvent_outflow, spc_solvent_inflow public :: spc_velocity_noneq_expol, spc_mole_fraction_noneq_expol ! moments BC public :: spc_moments_moleFrac public :: spc_moments_moleFlux public :: spc_moments_wall public :: spc_moments_vel contains ! ****************************************************************************** ! !> author: Kannan Masilamani !! Outlet boundary conditions with zero pressure gradient. !! !! These boundary conditions use the neighbor cells density and velocity in !! the equilibrium function. It is not necessary to specify density at !! boundary in the lua configuration file !! \( f = f^{eq}(\rho_{b-1},u_{b-1}) \) !! !! 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 spc_outlet_zero_prsgrd( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: rho, ux(3), usq, m_ratio integer :: iElem, iDir, iField_2, pos, QQ real(kind=rk) :: fTmp( layout%fStencil%QQ * varSys%nStateVars ) integer :: elemPos, posInBuffer ! --------------------------------------------------------------------------- QQ = layout%fStencil%QQ m_ratio = fieldProp%species%molWeigRatio do iElem = 1, globBC%nElems( iLevel ) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem) do iField_2 = 1, varSys%nStateVars do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iField_2)%state_varPos(iDir) fTmp(pos) = bcBuffer( & & ( posinbuffer-1)* nscalars+ pos ) end do end do !Calculate local density and velocity moments rho = 0._rk do iDir = 1,layout%fStencil%QQ rho = rho + fTmp(varPos(iDir)) end do call derVarPos%velFromState( state = fTmp, & & iField = iField, & & nElems = 1, & & varSys = varSys, & & layout = layout, & & res = ux ) usq = ux(1)*ux(1) + ux(2)*ux(2) + ux(3)*ux(3) elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) do iDir = 1,layout%fStencil%QQN ! Write the values if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then ! Only changed from save to fetch and added currT state( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) = & & layout%weight( iDir )*rho*(2.0_rk*m_ratio+9._rk*( ( & & layout%fStencil%cxDir( 1, layout%fStencil & & %cxDirInv( iDir ))*ux(1) & & + layout%fStencil%cxDir( 2, layout%fStencil & & %cxDirInv( iDir ))*ux(2) & & + layout%fStencil%cxDir( 3, layout%fStencil & & %cxDirInv( iDir ))*ux(3) & & )**2 - div1_3*usq)) & & - state( & & neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) end if end do end do end subroutine spc_outlet_zero_prsgrd ! ****************************************************************************** ! ! ****************************************************************************** ! !> Inlet boundary condition for defined species velocity and mole fraction !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_inlet', !! velocity = 'inlet_vel', !! mole_fraction = 'inlet_mole' !! } !!} !!variable = { !! { !! label = 'inlet_vel', !! ncomponents = 3, !! vartype = 'st_fun', !! st_fun = {0.1,0.0,0.0} !! }, !! { !! label = 'inlet_mole', !! ncomponents = 1, !! vartype = 'st_fun', !! st_fun = 0.1 !! } !!} !!``` !! !! 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 spc_inlet( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmp( layout%fStencil%QQ ) real(kind=rk) :: mass_dens real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel real(kind=rk) :: moleFrac(globBC%nElems(iLevel)) real(kind=rk) :: massFlux(3) integer :: iElem, iDir, QQ, posInBuffer, bcVel_pos, offset, elemPos integer :: bcMoleFrac_pos ! ------------------------------------------------------------------------ QQ = layout%fStencil%QQ inv_vel = 1.0_rk / physics%fac( iLevel )%vel ! position of boundary velocity in varSys bcVel_pos = me%bc_states%velocity%varPos ! Get velocity 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_b ) ! convert physical velocity into LB velocity vel_b = vel_b * inv_vel bcMoleFrac_pos = me%bc_states%moleFrac%varPos ! mole fraction call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFrac & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleFrac(:) ) do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*3 posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) fTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : & & (posInBuffer-1)*nScalars+varPos(QQ) ) ! local density !mass_dens = sum(fTmp) mass_dens = mixture%moleDens0LB * moleFrac(iElem) & & * fieldProp%species%molWeight massFlux = mass_dens * vel_b(offset+1:offset+3) elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1, layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then ! Depending on PUSH or pull, use + or - for cxDir, because directions ! are inverted state( neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fTmp(layout%fStencil%cxDirInv( iDir )) & & + layout%weight( iDir )*2._rk*cs2inv & & * ( layout%fStencil%cxDirRK( 1, iDir )*massFlux(1) & & + layout%fStencil%cxDirRK( 2, iDir )*massFlux(2) & & + layout%fStencil%cxDirRK( 3, iDir )*massFlux(3) ) end if end do end do end subroutine spc_inlet ! ****************************************************************************** ! ! ****************************************************************************** ! !> Inlet species velocity equilibrium boundary with specified !! mixture averaged mass velocity and its molefraction !! mixture kinematic pressure is extrapolated here. !! Density and velocity of all fields are used to compute equilibrium !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_inlet_eq', !! velocity = 'inlet_vel', !! mole_fraction = 'inlet_mole' !! } !!} !!variable = { !! { !! label = 'inlet_vel', !! ncomponents = 3, !! vartype = 'st_fun', !! st_fun = {0.1,0.0,0.0} !! }, !! { !! label = 'inlet_mole', !! ncomponents = 1, !! vartype = 'st_fun', !! st_fun = 0.1 !! } !!} !!``` !! !! 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 spc_inlet_eq( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) & & * varSys%nStateVars ) ! fEq uses AOS layout real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) ) real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars ) !real(kind=rk) :: uxB(3*globBC%nElems(iLevel)) !< Velocity on boundary element !> Velocity on boundary element real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3) real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars) real(kind=rk) :: spc_vel(globBC%nElems(iLevel)*3), inv_vel real(kind=rk) :: moleFrac(globBC%nElems(iLevel)) integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ integer :: offset, bcVel_pos, bcMoleFrac_pos, elemPos, posInBuffer ! ------------------------------------------------------------------------ QQ = layout%fStencil%QQ nFields = varSys%nStateVars inv_vel = 1.0_rk / physics%fac( iLevel )%vel ! position of boundary velocity in varSys bcVel_pos = me%bc_states%velocity%varPos ! Get velocity 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 = spc_vel ) ! convert physical velocity into LB velocity spc_vel = spc_vel * inv_vel ! position of molefraction spacetime variable in varSys bcMoleFrac_pos = me%bc_states%moleFrac%varPos ! mole fraction call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFrac & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleFrac ) ! copy state do iElem = 1, globBC%nElems(iLevel) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem) do iFieldLoc = 1, nFields do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) fTmp( ( ielem-1)* nscalars+ pos ) = & & bcBuffer( & & ( posinbuffer-1)* nscalars+ pos) end do end do !iField end do !iElem ! Derive all species velocities at once since ! it is inefficient to derive each species velocity call derVarPos%velocitiesFromState( state = fTmp, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = uxB ) ! get density and velocity. ! if current field, use velocity and density defined in lua file. ! for other fields, derive density and velocity from state mass_dens = 0.0_rk do iFieldLoc = 1, nFields if (iFieldLoc == iField) then ! store velocity in input_loc array do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc velocity(:,offset) = spc_vel((iElem-1)*3+1 : iElem*3) ! compute current species mass density from specified molefraction at ! boundary ! rho = n_t * chi_i * m_i + massFrac_i * rho0 * KinePress/(cs2*phi_i) ! massFrac_i = rho_i/rho ! rho = n_t * chi_i * m_i + rho_i * KinePress/(cs2*phi_i) ! 1st term is zero order density ! 2nd term is second order density with kinematic mixture pressure, ! p = cs2*(sum(phi_k*rho_k) - min_a(m_a)*n0)/rho0 !KM: using inital mixture number density rho0/mixtureMOlWeight !Using local tot_NuMdens increases density over time mass_dens( offset ) = mixture%moleDens0LB * moleFrac(iElem) & & * fieldProp%species%molWeight end do else do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc velocity(:,offset) = & & uxB((iElem-1)*nFields*3 + (iFieldLoc-1)*3 + 1 : & & (iElem-1)*nFields*3 + iFieldLoc*3 ) ! species density do iDir = 1, layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) mass_dens(offset) = mass_dens(offset) & & + fTmp(( ielem-1)* nscalars+ pos ) end do !iDir end do !iElem end if end do !iField ! compute equilibrium call derVarPos%equilFromMacro( density = mass_dens, & & velocity = velocity, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = fEq ) do iElem = 1, globBC%nElems(iLevel) elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1, layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then state(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq( iDir+(iElem-1)*QQ ) end if end do end do end subroutine spc_inlet_eq ! ****************************************************************************** ! ! ****************************************************************************** ! !> Inlet species velocity bounce back boundary with specified !! mixture averaged mass velocity and its molefraction !! mixture kinematic pressure is extrapolated here. !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_vel_bb', !! velocity = 'inlet_vel', !! } !!} !!variable = { !! name = 'inlet_vel', !! ncomponents = 3, !! vartype = 'st_fun', !! st_fun = {0.06, 0.0, 0.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 spc_vel_bb( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmp( layout%fStencil%QQ ) real(kind=rk) :: mass_dens real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel integer :: iElem, iDir, QQ, posInBuffer, bcVel_pos, offset, elemPos ! ------------------------------------------------------------------------ QQ = layout%fStencil%QQ inv_vel = 1.0_rk / physics%fac( iLevel )%vel ! position of boundary velocity in varSys bcVel_pos = me%bc_states%velocity%varPos ! Get velocity 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_b ) ! convert physical velocity into LB velocity vel_b = vel_b * inv_vel do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*3 posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) fTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : & & (posInBuffer-1)*nScalars+varPos(QQ) ) ! local density mass_dens = sum(fTmp) elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1,layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then ! Depending on PUSH or pull, use + or - for cxDir, because directions ! are inverted state( neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fTmp(layout%fStencil%cxDirInv( iDir )) & & + layout%weight( iDir )*2._rk*cs2inv*mass_dens & & * ( layout%fStencil%cxDirRK( 1, iDir )*vel_b(offset+1) & & + layout%fStencil%cxDirRK( 2, iDir )*vel_b(offset+2) & & + layout%fStencil%cxDirRK( 3, iDir )*vel_b(offset+3) ) end if end do end do end subroutine spc_vel_bb ! ****************************************************************************** ! ! ************************************************************************ ! !> species Outlet Pressure extrapolation boundary. NOT VERIFIED !! !! This is taken from the paper: !! M. Junk and Z. Yang, !! Asymptotic Analysis of Lattice Boltzmann Outflow Treatments !! Commun. Comput. Phys., pp. 1–11, 2011. !! !! * Pressure as defined in the configuration file !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_outlet_expol', !! pressure = 'p1' !! } !!} !!variable = { !! name = 'p1', !! ncomponents = 1, !! vartype = 'st_fun', !! st_fun = 1.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 spc_outlet_expol( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmpAll( layout%fStencil%QQ*globBC%nElems(iLevel) & & * varSys%nStateVars ) real(kind=rk) :: fTmp_1 real(kind=rk) :: fTmp_2 ! fEq uses AOS layout real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) ) integer :: iDir !< Direction counter integer :: iElem !< Element counter integer :: pos, ifieldloc, nFields, neighPos, QQ, elemPos ! ------------------------------------------------------------------------ QQ = layout%fStencil%QQ nFields = varSys%nStateVars ! Run over all elements in list for this boundary do iElem = 1, globBC%nElems(iLevel) neighPos = me%neigh(iLevel)%posInState(1,iElem) do iFieldLoc = 1, nFields do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) !neighbor element fTmpAll( ( ielem-1)* nscalars+ pos ) & & = state((neighpos-1)*nscalars+ idir+(ifieldloc-1)*qq) !local element !fTmpAll( ( ielem-1)* nscalars+ pos ) = bcBuffer( & ! & ( globbc%elemlvl( ilevel )%posinbcelembuf%val( ielem )-1)* nscalars+ pos ) end do end do !iField end do !iElem ! Calculate the equilibrium distribution call derVarPos%equilFromState( state = fTmpAll, & & iField = iField, & & nElems = globBC%nElems( iLevel ), & & varSys = varSys, & & layout = layout, & & res = fEq ) ! Run over all elements in list for this boundary do iElem = 1, globBC%nElems(iLevel) elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) ! Treat all directions, but actually we only want to treat ! non-normal directions ! Equation (3.2b) do iDir = 1, layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then ! KM: pre-collision of current time step must be used for extrapolation fTmp_1 = me%neigh( iLevel )%neighBufferPre_nNext( 1, iDir+(iElem-1)*QQ) fTmp_2 = me%neigh( iLevel )%neighBufferPre_nNext( 2, iDir+(iElem-1)*QQ) ! update the incoming velocity direction state(neigh (( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = 1.5_rk*fTmp_1 - 0.5_rk*fTmp_2 !& = feq( iDir+(iElem-1)*QQ ) end if end do ! then overwrite the normal direction with special treatment ! Equation (3.2a) iDir = globBC%elemLvl( iLevel )%normalInd%val( iElem ) state( neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) = & & feq( ( ielem-1)* qq+idir ) !& !& +(fTmpAll(( ielem-1)* nscalars+varpos(layout%fstencil%cxdirinv(idir)))& !& -feq( ( ielem-1)* layout%fstencil%qq+layout%fstencil%cxdirinv(idir) ) ) !& fTmpAll(( ielem-1)* nscalars+varpos(idir)) end do !iElem end subroutine spc_outlet_expol ! ****************************************************************************** ! ! ****************************************************************************** ! !> Outlet species velocity equilibrium boundary with specified !! mixture averaged mass velocity. !! molefraction is extrapolated here. !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_outlet_vel', !! velocity = 'inlet_vel', !! } !!} !!variable = { !! { !! label = 'inlet_vel', !! ncomponents = 3, !! vartype = 'st_fun', !! st_fun = {0.1,0.0,0.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 spc_outlet_vel( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) & & * varSys%nStateVars ) ! fEq uses AOS layout real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) ) real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars ) !> Velocity on boundary element real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3) real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars) real(kind=rk) :: spc_vel(globBC%nElems(iLevel)*3), inv_vel integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ integer :: offset, bcVel_pos, elemPos, posInBuffer ! ------------------------------------------------------------------------ QQ = layout%fStencil%QQ nFields = varSys%nStateVars inv_vel = 1.0_rk / physics%fac( iLevel )%vel ! position of boundary velocity in varSys bcVel_pos = me%bc_states%velocity%varPos ! Get velocity 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 = spc_vel ) ! convert physical velocity into LB velocity spc_vel = spc_vel * inv_vel ! copy state do iElem = 1, globBC%nElems(iLevel) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) do iFieldLoc = 1, nFields do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) fTmp( ( ielem-1)* nscalars+ pos ) & & = bcBuffer( & & ( posinbuffer-1)* nscalars+ pos ) end do end do !iField end do !iElem ! Derive all species velocities at once since ! it is inefficient to derive each species velocity call derVarPos%velocitiesFromState( state = fTmp, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = uxB ) ! use same velocity for all species and extrapolate density do iFieldLoc = 1, nFields do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc velocity(:,offset) = spc_vel((iElem-1)*3+1 : iElem*3) end do end do !iField ! species density mass_dens = 0.0_rk do iElem = 1, globBC%nElems(iLevel) do iFieldLoc = 1, nFields offset = (iElem-1)*nFields + iFieldLoc do iDir = 1, layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) mass_dens(offset) = mass_dens(offset) & & + fTmp(( ielem-1)* nscalars+ pos ) end do !iDir end do !iField end do !iElem ! compute equilibrium call derVarPos%equilFromMacro( density = mass_dens, & & velocity = velocity, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = fEq ) do iElem = 1, globBC%nElems(iLevel) elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1,layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then state(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq( iDir+(iElem-1)*QQ ) end if end do end do end subroutine spc_outlet_vel ! ****************************************************************************** ! ! ****************************************************************************** ! !> Outlet mixture pressure species equilibrium boundary !! kinematic pressure is computed from pressure !! species density and velocity are extrapolated !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_outlet_eq', !! pressure = 'press0' !! } !!} !!variable = { !! { !! label = 'press0', !! ncomponents = 1, !! vartype = 'st_fun', !! st_fun = 1.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 spc_outlet_eq( 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 ! -------------------------------------------------------------------- ! ! fEq uses AOS layout real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) ) !store pdf of all fields for derive function pointer input real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) & & * varSys%nStateVars ) integer :: iDir ! Direction counter integer :: iElem ! Element counter integer :: nFields, iFieldLoc, pos, offset, QQ!, bcPress_pos ! real(kind=rk) :: press(globBC%nElems(iLevel)) integer :: elemPos, posInBuffer, neighPos ! --------------------------------------------------------------------------- QQ = layout%fStencil%QQ nFields = varSys%nStateVars ! @todo KM: use pressure term to set species density at boundary ! get density value from the definition in the lua file ! position of boundary pressure in varSys !KM! bcPress_pos = me%bc_states%pressure%varPos !KM! ! get pressure variable from spacetime function !KM! call varSys%method%val(bcPress_pos)%get_valOfIndex( & !KM! & varSys = varSys, & !KM! & time = sim_time, & !KM! & iLevel = iLevel, & !KM! & idx = me%bc_states%pressure & !KM! & %pntIndex%indexLvl(iLevel) & !KM! & %val(1:globBC%nElems(iLevel)), & !KM! & nVals = globBC%nElems(iLevel), & !KM! & res = press ) !KM! !KM! ! If physical quantities are given, transform to lattice units by division !KM! ! with the conversion factor !KM! ! kinematic pressure = pressure/rho0 !KM! press = (press/mixture%rho0)/(physics%fac( iLevel )%press/physics%rho0) ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) neighPos = me%neigh(iLevel)%posInState(1,iElem) do iFieldLoc = 1, nFields offset = (iElem-1)*nFields + iFieldLoc do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) !fTmp( ( ielem-1)* nscalars+ pos ) & ! & = bcBuffer( & ! & ( posinbuffer-1)* nscalars+ pos ) fTmp( ( ielem-1)* nscalars+ pos ) & & = state((neighpos-1)*nscalars+ idir+(ifieldloc-1)*qq) end do end do end do !iElem ! Extrapolate equilibrium from local ! Calculate the equilibrium distribution call derVarPos%equilFromState( state = fTmp, & & iField = iField, & & nElems = globBC%nElems( iLevel ), & & varSys = varSys, & & layout = layout, & & res = fEq ) ! Run over all elements in list for this boundary do iElem = 1, globBC%nElems( iLevel ) elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) do iDir = 1,layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then ! Now assign the values state(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0 ) & & = fEq( iDir+(iElem-1)*QQ ) end if end do end do end subroutine spc_outlet_eq ! ****************************************************************************** ! ! ****************************************************************************** ! !> Mole fraction boundary condition !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_moleFrac', !! moleFraction = 0.0 !! } !!``` !! Post collision pdf of incoming link is updated with equilibrium functions !! which is similar to initial condition. !! \f$ \bar{f^c_k} = f^{eq}_k(\rho_k, u_k) + \frac{\lambda}{2} !! (f^{eq}_k(\rho_k, u_k) + f^{eq}_k(\rho_k, u^{eq}_k})) \f$ !! Here, !! $u_k$ - velocity of species k from original pdf computed from LSE !! $u^{eq}_k$ - equilibrium velocity of species k !! !! !! 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 spc_moleFrac( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: moleFrac(globBC%nElems(iLevel)) real(kind=rk) :: rho, press integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ integer :: bcMoleFrac_pos, elemPos, posInBuffer real(kind=rk) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars ) real(kind=rk) :: fEq, fEqStar, ucx, ucxStar, ucxQuad, ucxQuadStar real(kind=rk) :: usq, usqStar real(kind=rk) :: velAvg(3), velQuad(3), velQuadStar(3), eqVel(3) real(kind=rk) :: vel(3,varSys%nStateVars) real(kind=rk) :: mom(3,varSys%nStateVars) real(kind=rk) :: mass_dens(varSys%nStateVars), num_dens(varSys%nStateVars) real(kind=rk) :: totMassDens real(kind=rk) :: moleFrac_loc(varSys%nStateVars) real(kind=rk) :: molWeightInv(varSys%nStateVars), phi(varSys%nStateVars) real(kind=rk) :: resi_coeff( varSys%nStateVars, varSys%nStateVars ) real(kind=rk) :: weight0_inv, paramBInv type(mus_varSys_data_type), pointer :: fPtr !--------------------------------------------------------------------------- call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr ) ! write(*,*) 'boundary label ', trim(me%label) QQ = layout%fStencil%QQ nFields = varSys%nStateVars do iFieldLoc = 1, nFields ! species properties ! molecular weight inverse molWeightInv(iFieldLoc) = fPtr%solverData%scheme%field(iFieldLoc) & & %fieldProp%species%molWeightInv ! molecular weight ratio phi(iFieldLoc) = fPtr%solverData%scheme & & %field(iFieldLoc)%fieldProp%species%molWeigRatio ! resistivity coefficients resi_coeff(iFieldLoc, :) = fPtr%solverData%scheme & & %field(iFieldLoc)%fieldProp & & %species%resi_coeff(:) end do weight0_inv = 1.0_rk/layout%weight(layout%fStencil%restPosition) paramBInv = 1.0_rk / mixture%paramB ! position of molefraction spacetime variable in varSys bcMoleFrac_pos = me%bc_states%moleFrac%varPos ! mole fraction call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFrac & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleFrac ) ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem) ! local state vector to compute density and velocity of other species do iFieldLoc = 1, varSys%nStateVars do iDir = 1, QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) fTmp_all( pos ) = bcBuffer( & & ( posinbuffer-1)* nscalars+ pos) end do end do !KM: using inital mixture number density rho0/mixtureMOlWeight rho = mixture%moleDens0LB * moleFrac(iElem) & & * fieldProp%species%molWeight press = cs2 * rho * fieldProp%species%molWeigRatio !local density and momentum mass_dens = 0.0_rk vel = 0.0_rk do iFieldLoc = 1, nFields mass_dens(iFieldLoc) = & & sum( fTmp_all( varSys%method%val(iFieldLoc)%state_varPos(:) ) ) num_dens(iFieldLoc) = mass_dens(iFieldLoc) * molWeightInv(iFieldLoc) !velocity call derVarPos%momFromState( state = fTmp_all, & & iField = iFieldLoc, & & nElems = 1, & & varSys = varSys, & & layout = layout, & & res = mom(:, iFieldLoc) ) end do mass_dens(iField) = rho num_dens(iField) = rho / fieldProp%species%molWeight !mass flux do iFieldLoc = 1, nFields vel(:, iFieldLoc) = mom(:, iFieldLoc) / mass_dens(iFieldLoc) end do ! total mass density totMassDens = sum(mass_dens) ! mole fraction moleFrac_loc = num_dens/sum(num_dens) ! equilibrium velocity eqVel = vel(:, iField) do iFieldLoc = 1, nFields eqVel(:) = eqVel(:) + resi_coeff(iField, iFieldLoc) & & * phi(iField) * moleFrac_loc(iFieldLoc) & & * ( vel(:, iFieldLoc) - vel( :, iField ) ) & & / mixture%paramB end do ! mass averaged mixture velocity velAvg(1) = sum(mom(1,:)) / totMassDens velAvg(2) = sum(mom(2,:)) / totMassDens velAvg(3) = sum(mom(3,:)) / totMassDens velQuadStar(:) = mixture%theta_eq*velAvg(:) & & + (1.0_rk-mixture%theta_eq) * eqVel(:) velQuad(:) = mixture%theta_eq*velAvg(:) & & + (1.0_rk-mixture%theta_eq)*vel(:, iField) usqStar = dot_product(velQuadStar, velQuadStar)*t2cs2inv usq = dot_product(velQuad, velQuad)*t2cs2inv elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir =1,layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then ucx = dot_product(layout%fStencil%cxDir(:, iDir), & & vel(:,iField)) ucxQuad = dot_product(layout%fStencil%cxDir(:, iDir), & & velQuad) ucxStar = dot_product(layout%fStencil%cxDir(:, iDir), & & eqVel) ucxQuadStar = dot_product(layout%fStencil%cxDir(:, iDir), & & velQuadStar) ! eqVel is actually is rho_i*eqVel so ucxStar is not multiplied ! with rho in below equation fEqStar = layout%weight(iDir) * rho & & * ( phi(iField) + ucxStar * cs2inv & & + ucxQuadStar * ucxQuadStar * t2cs4inv - usqStar ) fEq = layout%weight(iDir) * rho & & * ( phi(iField) + ucx * cs2inv & & + ucxQuad * ucxQuad * t2cs4inv - usq ) if ( iDir == layout%fStencil%restPosition ) then ! equilibrium at rest fEq = layout%weight( iDir ) * rho * ( & & ( weight0_inv + (1.0_rk-weight0_inv) * phi(iField) ) - usq ) fEqStar = layout%weight( iDir ) * rho * ( & & ( weight0_inv + (1.0_rk-weight0_inv) * phi(iField) ) - usqStar ) end if ! set equilibrium ! setting transformed pdf in similar to initial condition i.e ! \bar{f^c} = f^{eq}(\rho_spc, u_spc) + \frac{\lambda}{2} ! (f^{eq}(\rho_spc, u_spc) + f^{eq}(\rho_spc, u^{eq}_spc})) state(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq + mixture%omega_diff*0.5_rk*( fEq - fEqStar ) end if end do! iDir end do end subroutine spc_moleFrac ! ****************************************************************************** ! ! ****************************************************************************** ! !> Mole fraction boundary condition !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_molefrac_eq', !! moleFraction = 0.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 spc_moleFrac_eq( 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 ! -------------------------------------------------------------------- ! ! fEq uses AOS layout real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) ) !store pdf of all fields for derive function pointer input real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) & & * varSys%nStateVars ) real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars ) real(kind=rk) :: moleFrac(globBC%nElems(iLevel)) integer :: iElem, iDir, iFieldLoc, nFields, pos !> Velocity on boundary element real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3) real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars) integer :: offset, QQ, bcMoleFrac_pos, elemPos, posInBuffer ! --------------------------------------------------------------------------- QQ = layout%fStencil%QQ nFields = varSys%nStateVars ! position of molefraction spacetime variable in varSys bcMoleFrac_pos = me%bc_states%moleFrac%varPos ! mole fraction call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFrac & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleFrac ) ! copy state do iElem = 1, globBC%nElems(iLevel) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) do iFieldLoc = 1, nFields do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) fTmp( ( ielem-1)* nscalars+ pos ) & & = bcBuffer( & & ( posinbuffer-1)* nscalars+ pos ) end do end do !iField end do !iElem ! get density and velocity. ! if current field, use velocity and density defined in lua file. ! for other fields, derive density and velocity from state mass_dens = 0.0_rk do iFieldLoc = 1, nFields if (iFieldLoc == iField) then ! store velocity in input_loc array do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc ! compute current species mass density from specified molefraction at ! boundary ! rho = n_t * chi_i * m_i + massFrac_i * rho0 * KinePress/(cs2*phi_i) ! massFrac_i = rho_i/rho ! rho = n_t * chi_i * m_i + rho_i * KinePress/(cs2*phi_i) ! 1st term is zero order density ! 2nd term is second order density with kinematic mixture pressure, ! p = cs2*(sum(phi_k*rho_k) - min_a(m_a)*n0)/rho0 !KM: using inital mixture number density rho0/mixtureMOlWeight !Using local tot_NuMdens increases density over time !kinePress = cs2 * ( dot_product(phi, mass_dens) & ! & - minval(molWeight)*mixture%moleDens0LB )/mixture%rho0LB mass_dens( offset ) = mixture%moleDens0LB * moleFrac(iElem) & & * fieldProp%species%molWeight end do else do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc ! species density do iDir = 1, layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) mass_dens(offset) = mass_dens(offset) & & + fTmp(( ielem-1)* nscalars+ pos ) end do !iDir end do !iElem end if end do !iField ! Derive all species velocities at once since ! it is inefficient to derive each species velocity call derVarPos%velocitiesFromState( state = fTmp, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = uxB ) do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields do iFieldLoc = 1, nFields velocity(:, offset+iFieldLoc) = & & uxB(offset*3 + (iFieldLoc-1)*3 + 1 : & & offset*3 + iFieldLoc*3 ) end do !iField end do !iElem fEq = 0.0_rk ! Calculate the equilibrium distribution call derVarPos%equilFromMacro( density = mass_dens, & & velocity = velocity, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = fEq ) do iElem = 1, globBC%nElems(iLevel) elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1,layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then state(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq( iDir+(iElem-1)*QQ ) end if end do end do end subroutine spc_moleFrac_eq ! ****************************************************************************** ! ! ****************************************************************************** ! !> Mole density boundary condition !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_moledens_eq', !! moleFraction = 0.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 spc_moleDens_eq( 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 ! -------------------------------------------------------------------- ! ! --------------------------------------------------------------------------- ! fEq uses AOS layout real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) ) !store pdf of all fields for derive function pointer input real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) & & * varSys%nStateVars ) real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars ) real(kind=rk) :: moleDens(globBC%nElems(iLevel)) integer :: iElem, iDir, iFieldLoc, nFields, pos !> Velocity on boundary element real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3) real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars) integer :: offset, QQ, bcMoleDens_pos, elemPos, posInBuffer ! --------------------------------------------------------------------------- QQ = layout%fStencil%QQ nFields = varSys%nStateVars ! position of molefraction spacetime variable in varSys bcMoleDens_pos = me%bc_states%moleDens%varPos ! mole fraction call varSys%method%val(bcMoleDens_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFrac & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleDens ) ! Convert to lattice Unit moleDens = moleDens / physics%moleDens0 ! copy state do iElem = 1, globBC%nElems(iLevel) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) do iFieldLoc = 1, nFields do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) fTmp( ( ielem-1)* nscalars+ pos ) & & = bcBuffer( & & ( posinbuffer-1)* nscalars+ pos ) end do end do !iField end do !iElem ! get density and velocity. ! if current field, use velocity and density defined in lua file. ! for other fields, derive density and velocity from state mass_dens = 0.0_rk do iFieldLoc = 1, nFields if (iFieldLoc == iField) then ! store velocity in input_loc array do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc ! compute current species mass density from specified molefraction at ! boundary ! rho = n_t * chi_i * m_i + massFrac_i * rho0 * KinePress/(cs2*phi_i) ! massFrac_i = rho_i/rho ! rho = n_t * chi_i * m_i + rho_i * KinePress/(cs2*phi_i) ! 1st term is zero order density ! 2nd term is second order density with kinematic mixture pressure, ! p = cs2*(sum(phi_k*rho_k) - min_a(m_a)*n0)/rho0 !KM: using inital mixture number density rho0/mixtureMOlWeight !Using local tot_NuMdens increases density over time !kinePress = cs2 * ( dot_product(phi, mass_dens) & ! & - minval(molWeight)*mixture%moleDens0LB )/mixture%rho0LB mass_dens( offset ) = moleDens(iElem) * fieldProp%species%molWeight end do else do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc ! species density do iDir = 1, layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) mass_dens(offset) = mass_dens(offset) & & + fTmp(( ielem-1)* nscalars+ pos ) end do !iDir end do !iElem end if end do !iField ! Derive all species velocities at once since ! it is inefficient to derive each species velocity call derVarPos%velocitiesFromState( state = fTmp, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = uxB ) do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields do iFieldLoc = 1, nFields velocity(:, offset+iFieldLoc) = & & uxB(offset*3 + (iFieldLoc-1)*3 + 1 : & & offset*3 + iFieldLoc*3 ) end do !iField end do !iElem fEq = 0.0_rk ! Calculate the equilibrium distribution call derVarPos%equilFromMacro( density = mass_dens, & & velocity = velocity, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = fEq ) do iElem = 1, globBC%nElems(iLevel) elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1,layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then state(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq( iDir+(iElem-1)*QQ ) end if end do end do end subroutine spc_moleDens_eq ! ****************************************************************************** ! ! ****************************************************************************** ! !> Mole fraction boundary condition with thermodynamic factor !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_moleFrac_wtdf', !! moleFraction = 0.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 spc_moleFrac_wtdf( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fEq !store pdf of all fields for derive function pointer input real(kind=rk) :: fTmp( varSys%nScalars ) real(kind=rk) :: mass_dens( varSys%nStateVars ) real(kind=rk) :: num_dens( varSys%nStateVars ) real(kind=rk) :: tot_numDens, tot_massDens real(kind=rk) :: moleFrac(globBC%nElems(iLevel)) integer :: iElem, iDir, iFieldLoc, iField_2, iField_3, nFields, pos, QQ real(kind=rk), dimension(varSys%nStateVars) :: molWeight, phi real(kind=rk), dimension(varSys%nStateVars) :: moleFrac_loc real(kind=rk), dimension(3, varSys%nStateVars) :: first_moments, velocity real(kind=rk), & & dimension(varSys%nStateVars, varSys%nStateVars) :: matA, invA, & & resi_coeff, thermodynamic_fac, inv_thermodyn_fac, diff_coeff real(kind=rk) :: temp, press, phy_moleDens_fac real(kind=rk) :: omega, paramB, theta_eq real(kind=rk) :: eqVel(3), velAvg(3), velQuad(3), usqr, ucx, ucxQuadTerm type(mus_varSys_data_type), pointer :: fPtr integer :: bcMoleFrac_pos, posInBuffer, posInState ! ------------------------------------------------------------------------ call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr ) QQ = layout%fStencil%QQ nFields = varSys%nStateVars do iFieldLoc = 1, nFields ! species properties ! molecular weight inverse molWeight(iFieldLoc) = fPtr%solverData%scheme%field(iFieldLoc) & & %fieldProp%species%molWeight ! molecular weight ratio phi(iFieldLoc) = fPtr%solverData%scheme%field(iFieldLoc) & & %fieldProp%species%molWeigRatio ! resistivity coefficients resi_coeff(iFieldLoc, :) = fPtr%solverData%scheme%field(iFieldLoc) & & %fieldProp%species%resi_coeff(:) end do !KM @todo check moleDens for multilevel phy_moleDens_fac = physics%moleDens0 !fixed parameter paramB = mixture%paramB !equilibrium theta theta_eq = mixture%theta_eq omega = mixture%relaxLvl(iLevel)%omega_diff ! temperature temp = mixture%temp0 ! atmospheric pressure press = mixture%atm_press !write(*,*) 'iField ', iField !write(*,*) 'resi_coeff ', resi_coeff !write(*,*) 'molWeight moleFrac BC', molWeight !write(*,*) 'phi ', phi !stop ! position of molefraction spacetime variable in varSys bcMoleFrac_pos = me%bc_states%moleFrac%varPos ! mole fraction call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFrac & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleFrac ) ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) !write(*,*) 'IElem ', iElem mass_dens = 0.0_rk first_moments = 0.0_rk posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) do iFieldLoc = 1, varSys%nStateVars do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) fTmp( pos ) = bcBuffer( & & ( posinbuffer-1)* nscalars+ pos ) mass_dens(iFieldLoc) = mass_dens(iFieldLoc) + fTmp( pos ) !field momentum (rho*u) first_moments( 1, iFieldLoc ) = first_moments( 1, iFieldLoc ) & & + fTmp( pos ) * layout%fStencil%cxDir( 1, iDir ) first_moments( 2, iFieldLoc ) = first_moments( 2, iFieldLoc ) & & + fTmp( pos ) * layout%fStencil%cxDir( 2, iDir ) first_moments( 3, iFieldLoc ) = first_moments( 3, iFieldLoc ) & & + fTmp( pos ) * layout%fStencil%cxDir( 3, iDir ) end do num_dens(iFieldLoc) = mass_dens(iFieldLoc)/molWeight(iFieldLoc) end do ! update num_dens, massDens and moleFrac of current species with specified ! molefraction num_dens(iField) = moleFrac(iElem) * mixture%moleDens0LB mass_dens(iField) = num_dens(iField) * molWeight(iField) !total number density tot_NumDens = sum(num_dens) !total mass density tot_massDens = sum(mass_dens) !mole fraction moleFrac_loc(:) = num_dens(:)/tot_NumDens !write(*,*) 'num_dens ', num_dens !write(*,*) 'num_dens*phy_moleDens_fac ', num_dens*phy_moleDens_fac ! MS-Diff coeff matrix from C++ code call mus_calc_MS_DiffMatrix( nFields, temp, press, & & num_dens*phy_moleDens_fac, diff_coeff ) !write(*,*) 'diff_coeff ', diff_coeff ! Convert to lattice unit resi_coeff = physics%fac(iLevel)%diffusivity/diff_coeff !write(*,*) 'resi_coeff ', resi_coeff call mus_calc_thermFactor( nFields, temp, press, moleFrac_loc, & & thermodynamic_fac ) !write(*,*) 'thermodyn_fac ', thermodynamic_fac inv_thermodyn_fac = invert_matrix( thermodynamic_fac ) !write(*,*) 'inv_thermodyn_fac ', inv_thermodyn_fac matA = 0.0_rk !build up matrix to solver LSE for actual velocity do iFieldLoc = 1, nFields !set diagonal part matA(iFieldLoc, iFieldLoc) = 1.0_rk do iField_2 = 1, nFields do iField_3 = 1, nFields matA(iFieldLoc, iField_2) = matA(iFieldLoc, iField_2) & & + omega * 0.5_rk & & * inv_thermodyn_fac(iFieldLoc, iField_2) & & * resi_coeff(iField_2, iField_3) & & * phi(iField_2) * moleFrac_loc(iField_3) & & / paramB end do end do !set non-diagonal part do iField_2 = 1, nFields do iField_3 = 1, nFields matA(iFieldLoc, iField_3) = matA(iFieldLoc, iField_3) & & - omega * 0.5_rk & & * inv_thermodyn_fac(iFieldLoc, iField_2) & & * resi_coeff(iField_2, iField_3) & & * phi(iField_3) * moleFrac_loc(iField_2) & & / paramB end do end do end do !write(*,*) 'matA ', matA ! invert matrix invA = invert_matrix( matA ) !write(*,*) 'invA ', invA !actual velocity of all species velocity(1, :) = matmul( invA, first_moments(1,:) ) / mass_dens(:) velocity(2, :) = matmul( invA, first_moments(2,:) ) / mass_dens(:) velocity(3, :) = matmul( invA, first_moments(3,:) ) / mass_dens(:) !write(*,*) 'velocity ', velocity ! equilibrium velocity with thermodynamic factor eqVel( : ) = mass_dens(iField)*velocity( :, iField ) do iField_2 = 1, nFields do iField_3 = 1, nFields eqVel( : ) = eqVel( : ) & & + inv_thermodyn_fac(iField, iField_2) & & * mass_dens(iField_2) & & * resi_coeff( iField_2, iField_3 ) * phi(iField_2)& & * moleFrac_loc(iField_3) & & * (velocity(:, iField_3) - velocity(:,iField_2)) & & / paramB end do end do !write(*,*) 'eqVel ', eqVel !compute mass averaged mixture velocity velAvg(1) = dot_product( mass_dens, velocity(1,:) )/tot_massDens velAvg(2) = dot_product( mass_dens, velocity(2,:) )/tot_massDens velAvg(3) = dot_product( mass_dens, velocity(3,:) )/tot_massDens velQuad = theta_eq*velAvg & & + (1.0_rk - theta_eq) * eqVel(:) / mass_dens(iField) usqr = dot_product( velQuad, velQuad ) * t2cs2inv posInState = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1,layout%fStencil%QQN ! Write the values ! The bitmask points into the incoming direction into the flow domain, ! which actually we want to update ! * For PUSH, we write to the incoming position, ! as the kernel reads it from there without propagation. ! * For PULL, we need to write to the inverse direction, as the kernel ! performs a bounce back before reading it. ! However, this bounced back direction actually comes from the ! non-existent boundary element and would point into the incoming ! direction, so the value has to be treated and set as if it points ! really into the incoming direction. if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then ! Depending on PUSH or pull, use + or - for cxDir, because directions ! are inverted ucx = dot_product( & & dble(layout%fStencil%cxDir(:, iDir)), eqVel(:) ) ucxQuadTerm = dot_product( & & dble(layout%fStencil%cxDir(:, iDir)), velQuad ) ! Equilibrium feq = layout%weight(iDir) * ( mass_dens(iField) * ( phi(iField) & & + ucxQuadTerm * ucxQuadTerm * t2cs4inv - usqr ) + ucx * cs2inv ) if ( iDir == layout%fStencil%restPosition ) then ! equilibrium at rest select case( trim(layout%fStencil%label) ) case('d2q9') feq = layout%weight( iDir ) * mass_dens(iField) * ( & & ( 9._rk - 5._rk * phi(iField) )/4._rk - usqr ) case('d3q19') feq = layout%weight( iDir ) * mass_dens(iField) * ( & & ( 3._rk - 2._rk * phi(iField) ) - usqr ) end select end if ! set equilibrium state( & & neigh (( idir-1)* nsize+ posinstate)+( ifield-1)* qq+ nscalars*0)& & = fEq !write(*,*) 'fEq ', fEq end if end do end do end subroutine spc_moleFrac_wtdf ! ****************************************************************************** ! ! ****************************************************************************** ! !> molar flux equilibrium boundary condition !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_moleflux', !! moleflux = 0.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 spc_moleFlux_eq( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) ) !store pdf of all fields for derive function pointer input real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) & & * varSys%nStateVars ) real(kind=rk) :: moleFlux(globBC%nElems(iLevel)*3), inv_flux !> Velocity on boundary element real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3) real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars ) real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars) real(kind=rk) :: num_dens integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ integer :: offset, posInBuffer, bcMoleFlux_pos, elemPos ! ------------------------------------------------------------------------ QQ = layout%fStencil%QQ nFields = varSys%nStateVars inv_flux = 1.0_rk / physics%fac( iLevel )%moleFlux ! position of boundary moleflux in varSys bcMoleFlux_pos = me%bc_states%moleFlux%varPos ! Get moleflux call varSys%method%val(bcMoleFlux_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFlux & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleFlux ) ! If physical quantities are given, transform to lattice units by division ! with the conversion factor moleFlux = moleFlux * inv_flux ! copy state do iElem = 1, globBC%nElems(iLevel) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) do iFieldLoc = 1, nFields do iDir = 1,layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) fTmp( ( ielem-1)* nscalars+ pos ) = & & bcBuffer( & & ( posinbuffer-1)* nscalars+ pos ) end do end do !iField end do !iElem ! get density and velocity. ! if current field, use velocity and density defined in lua file. ! for other fields, derive density and velocity from state ! species density mass_dens = 0.0_rk do iElem = 1, globBC%nElems(iLevel) do iFieldLoc = 1, nFields offset = (iElem-1)*nFields + iFieldLoc do iDir = 1, layout%fStencil%QQ pos = varSys%method%val(iFieldLoc)%state_varPos(iDir) mass_dens(offset) = mass_dens(offset) & & + fTmp(( ielem-1)* nscalars+ pos ) end do !iDir end do !iField end do !iElem ! Derive all species velocities at once since ! it is inefficient to derive each species velocity call derVarPos%velocitiesFromState( state = fTmp, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = uxB ) do iFieldLoc = 1, nFields if (iFieldLoc == iField) then ! store velocity in input_loc array do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc num_dens = mass_dens(offset) * fieldProp%species%molWeightInv velocity(:,offset) = moleFlux((iElem-1)*3+1:iElem*3) / num_dens end do else do iElem = 1, globBC%nElems(iLevel) offset = (iElem-1)*nFields + iFieldLoc velocity(:,offset) = & & uxB((iElem-1)*nFields*3 + (iFieldLoc-1)*3 + 1 : & & (iElem-1)*nFields*3 + iFieldLoc*3 ) end do !iElem end if end do !iField ! Calculate the equilibrium distribution call derVarPos%equilFromMacro( density = mass_dens, & & velocity = velocity, & & iField = iField, & & nElems = globBC%nElems(iLevel), & & varSys = varSys, & & layout = layout, & & res = fEq ) do iElem = 1, globBC%nElems(iLevel) if( .not. btest( levelDesc%property( & & globBC%elemLvl(iLevel)%elem%val(iElem)), prp_solid))then elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1, layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then state(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq( iDir+(iElem-1)*QQ ) end if end do end if end do end subroutine spc_moleFlux_eq ! ****************************************************************************** ! ! ****************************************************************************** ! !> molar flux boundary condition like velocity bounce back bc type !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_moleflux', !! mole_flux = 'mole_flux' !! } !!} !!variable = { !! name = 'mole_flux', !! ncomponents = 3, !! vartype = 'st_fun', !! st_fun = {0.06, 0.0, 0.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 spc_moleFlux( 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 ! -------------------------------------------------------------------- ! ! ------------------------------------------------------------------------ real(kind=rk) :: fTmp( layout%fStencil%QQ ) real(kind=rk) :: molWeight real(kind=rk) :: massFlux(3) real(kind=rk) :: moleFlux(globBC%nElems(iLevel)*3), inv_flux integer :: iElem, iDir, QQ, posInBuffer, offset, bcMoleFlux_pos integer :: elemPos ! ------------------------------------------------------------------------ QQ = layout%fStencil%QQ molWeight = fieldProp%species%molWeight inv_flux = 1.0_rk / physics%fac( iLevel )%moleFlux ! position of boundary velocity in varSys bcMoleFlux_pos = me%bc_states%moleFlux%varPos ! Get moleflux call varSys%method%val(bcMoleFlux_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFlux & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleFlux ) ! If physical quantities are given, transform to lattice units by division ! with the conversion factor moleFlux = moleFlux * inv_flux !moleFlux = moleFlux / physics%fac(iLevel)%flux ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) do iDir = 1,layout%fStencil%QQ fTmp( iDir ) = bcBuffer( & & ( posinbuffer-1)* nscalars+ varpos(idir) ) end do !caulate mass flux from moleFlux ! massFlux = moleflux * molWeight offset = (iElem-1)*3 massFlux(1) = moleFlux(offset+1) * molWeight massFlux(2) = moleFlux(offset+2) * molWeight massFlux(3) = moleFlux(offset+3) * molWeight !write(dbgUnit(1),*) 'massFlux ', massFlux !write(*,*) iElem, 'mass_flux ', massFlux elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) do iDir = 1,layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then state(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fTmp(layout%fStencil%cxDirInv( iDir )) & & + layout%weight( iDir )*2._rk*cs2inv & & * ( layout%fStencil%cxDir( 1, iDir )*massFlux(1) & & + layout%fStencil%cxDir( 2, iDir )*massFlux(2) & & + layout%fStencil%cxDir( 3, iDir )*massFlux(3) ) end if end do end do end subroutine spc_moleFlux ! ****************************************************************************** ! ! **************************************************************************** ! !> Inflow boundary condition for solvent based on non-Equilbrium !! extrapolation method. Similar to spc_velocity_noneq_expol except !! the mass density for solvent is enforced such that total moleDens0 is !! maintained. !! Default qVal=1.0. !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'inlet', !! kind = 'spc_solvent_inflow', !! velocity = {0.1,0.0,0.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 spc_solvent_inflow( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmpAll_1f(nScalars), fTmpAll_2f(nScalars) real(kind=rk) :: fEq_b, fEq_1f real(kind=rk) :: velAvg_1f(3), vel_b(3), eqVel_1f(3) real(kind=rk) :: usq_1f, ucx_1f, usq_b, ucx_b, ucxQuad_1f real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_1f, mass_dens_2f real(kind=rk), dimension(varSys%nStateVars) :: num_dens_1f, num_dens_2f real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars) real(kind=rk) :: mass_dens_b, num_dens_b, tot_NumDens_b real(kind=rk) :: omega_fac, paramBInv integer :: iElem, iDir, QQ, iFld, vPos integer :: bcVel_pos integer :: elemPos, elemOff, posInBuffer, neighPos real(kind=rk) :: vel_w(3*globBC%nElems(iLevel)), inv_vel type(mus_varSys_data_type), pointer :: fPtr integer :: nFields ! ------------------------------------------------------------------------- nFields = varSys%nStateVars ! 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:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = vel_w ) ! convert physical velocity into LB velocity inv_vel = 1.0_rk / physics%fac( iLevel )%vel vel_w = vel_w * inv_vel call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr ) ! resitivity coefficient of all species do iFld = 1, nFields !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:) resi_coeff(iFld, :) = fPtr%solverData%scheme & & %field(iFld)%fieldProp & & %species%resi_coeff(:) end do ! omega factor of LSE omega_fac = mixture%omega_diff * 0.5_rk paramBInv = 1.0_rk / mixture%paramB associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val, & & species => fPtr%solverData%scheme%field(:)%fieldProp%species, & & neighBufferPre_nNext => me%neigh(iLevel)%neighBufferPre_nNext, & & stencil => layout%fStencil ) QQ = layout%fStencil%QQ ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) ! Position of current boundary element in bcBuffer posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem ) elemOff = (posInBuffer-1)*nScalars ! Position of 2nd fluid neighbor in state array neighPos = me%neigh(iLevel)%posInState(1, iElem) mass_dens_2f = 0.0_rk do iFld = 1, nFields do iDir = 1, QQ ! Position of current field variable current dir in state array vPos = varSys%method%val(iFld)%state_varPos(iDir) ! State array of 1st fluid element fTmpAll_1f(vPos) = bcBuffer( elemOff + vPos ) ! Pre-Collision values of second fluid element fTmpAll_2f(vPos) = state( & & neigh((idir-1)* nsize+ neighpos)+( ifld-1)* qq+ nscalars*0 ) ! mass density of second neighbor mass_dens_2f(iFld) = mass_dens_2f(iFld) + fTmpAll_2f(vPos) end do num_dens_2f(iFld) = mass_dens_2f(iFld) * species(iFld)%molWeightInv end do ! calculate spc density and velAvg from PDF of first fluid call calcDensAndVelsFromPDF( nFields = nFields, & & iField = iField, & & mass_dens = mass_dens_1f, & & num_dens = num_dens_1f, & & velAvg = velAvg_1f, & & eqVel = eqVel_1f, & & varSys = varSys, & & pdf = fTmpAll_1f, & & stencil = layout%fStencil, & & species = species, & & resi_coeff = resi_coeff, & & omega_fac = omega_fac, & & paramBInv = paramBInv ) ! fluid node velocity square term of equilibrium function usq_1f = dot_product( velAvg_1f(1:3), velAvg_1f(1:3) ) * t2cs2inv ! Second-order extrapolation of number density of current species num_dens_b = ( 4.0_rk * num_dens_1f(iField) - num_dens_2f(iField) ) & & / 3.0_rk ! Second-order extrapolation of numerical total number density tot_NumDens_b = ( 4.0_rk * sum(num_dens_1f) - sum(num_dens_2f) ) & & / 3.0_rk ! Compute number density of current species at boundary element num_dens_b = mixture%moleDens0LB - ( tot_NumDens_b - num_dens_b ) ! Mass density of current species at boundary element mass_dens_b = num_dens_b * species(iField)%molWeight !mass_dens_b = mixture%rho0LB - ( sum(mass_dens_1f) - mass_dens_1f(iField) ) ! velocity at boundary node for qVal = 1.0 vel_b = vel_w( (iElem-1)*3 + 1: iElem*3 ) ! boundary node velocity square term of equilibrium function usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv ! Position of this boundary element in state array elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) do iDir = 1, layout%fStencil%QQN if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then ! Compute equilibrium on fluid node ucx_1f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_1f ) ucxQuad_1f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_1f ) fEq_1f = layout%weight(iDir) * mass_dens_1f(iField) & & * ( species(iField)%molWeigRatio + ucx_1f * cs2inv & & + ucxQuad_1f * ucxQuad_1f * t2cs4inv - usq_1f ) ! Compute equilibrium on boundary node ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b ) fEq_b = layout%weight(iDir) * mass_dens_b & & * ( species(iField)%molWeigRatio + ucx_b * cs2inv & & + ucx_b * ucx_b * t2cs4inv - usq_b ) state( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_b + ( fTmpAll_1f(varPos(iDir)) - fEq_1f ) end if end do end do end associate end subroutine spc_solvent_inflow ! **************************************************************************** ! ! **************************************************************************** ! !> Velocity boundary condition based on non-Equilbrium extrapolation method. !! Default qVal=1.0. !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'inlet', !! kind = 'spc_velocity_noneq_expol', !! velocity = {0.1,0.0,0.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 spc_velocity_noneq_expol( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmpAll_1f(nScalars), fTmpAll_2f(nScalars) real(kind=rk) :: fEq_b, fEq_1f real(kind=rk) :: velAvg_1f(3), vel_b(3), eqVel_1f(3) real(kind=rk) :: usq_1f, ucx_1f, usq_b, ucx_b real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_1f, mass_dens_2f, & & num_dens_1f real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars) real(kind=rk) :: mass_dens_b real(kind=rk) :: omega_fac, paramBInv integer :: iElem, iDir, QQ, iFld, vPos integer :: bcVel_pos integer :: elemPos, elemOff, posInBuffer, neighPos real(kind=rk) :: vel_w(3*globBC%nElems(iLevel)), inv_vel type(mus_varSys_data_type), pointer :: fPtr integer :: nFields ! ------------------------------------------------------------------------- nFields = varSys%nStateVars ! 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:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = vel_w ) ! convert physical velocity into LB velocity inv_vel = 1.0_rk / physics%fac( iLevel )%vel vel_w = vel_w * inv_vel call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr ) ! resitivity coefficient of all species do iFld = 1, nFields !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:) resi_coeff(iFld, :) = fPtr%solverData%scheme & & %field(iFld)%fieldProp & & %species%resi_coeff(:) end do ! omega factor of LSE omega_fac = mixture%omega_diff * 0.5_rk paramBInv = 1.0_rk / mixture%paramB associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val, & & species => fPtr%solverData%scheme%field(:)%fieldProp%species, & & neighBufferPre_nNext => me%neigh(iLevel)%neighBufferPre_nNext, & & stencil => layout%fStencil ) QQ = layout%fStencil%QQ ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) ! Position of current boundary element in bcBuffer posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem ) elemOff = (posInBuffer-1)*nScalars ! Position of 2nd fluid neighbor in state array neighPos = me%neigh(iLevel)%posInState(1, iElem) mass_dens_2f = 0.0_rk do iFld = 1, nFields do iDir = 1, QQ ! Position of current field variable current dir in state array vPos = varSys%method%val(iFld)%state_varPos(iDir) ! State array of 1st fluid element fTmpAll_1f(vPos) = bcBuffer( elemOff + vPos ) ! Pre-Collision values of second fluid element fTmpAll_2f(vPos) = state( & & neigh((idir-1)* nsize+ neighpos)+( ifld-1)* qq+ nscalars*0 ) ! mass density of second neighbor mass_dens_2f(iFld) = mass_dens_2f(iFld) + fTmpAll_2f(vPos) end do end do ! calculate spc density and velAvg from PDF of first fluid call calcDensAndVelsFromPDF( nFields = nFields, & & iField = iField, & & mass_dens = mass_dens_1f, & & num_dens = num_dens_1f, & & velAvg = velAvg_1f, & & eqVel = eqVel_1f, & & varSys = varSys, & & pdf = fTmpAll_1f, & & stencil = layout%fStencil, & & species = species, & & resi_coeff = resi_coeff, & & omega_fac = omega_fac, & & paramBInv = paramBInv ) ! fluid node velocity square term of equilibrium function usq_1f = dot_product( velAvg_1f(1:3), velAvg_1f(1:3) ) * t2cs2inv ! density of current species at boundary node extrapolated from ! fluid elements mass_dens_b = ( 4.0_rk * mass_dens_1f(iField) - mass_dens_2f(iField) ) & & / 3.0_rk ! velocity at boundary node for qVal = 1.0 vel_b = vel_w( (iElem-1)*3 + 1: iElem*3 ) ! boundary node velocity square term of equilibrium function usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv ! Position of this boundary element in state array elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) do iDir = 1, layout%fStencil%QQN if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then ! Compute equilibrium on fluid node ucx_1f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_1f ) fEq_1f = layout%weight(iDir) * mass_dens_1f(iField) & & * ( species(iField)%molWeigRatio + ucx_1f * cs2inv & & + ucx_1f * ucx_1f * t2cs4inv - usq_1f ) ! Compute equilibrium on boundary node ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b ) fEq_b = layout%weight(iDir) * mass_dens_b & & * ( species(iField)%molWeigRatio + ucx_b * cs2inv & & + ucx_b * ucx_b * t2cs4inv - usq_b ) state( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_b + ( fTmpAll_1f(varPos(iDir)) - fEq_1f ) end if end do end do end associate end subroutine spc_velocity_noneq_expol ! **************************************************************************** ! ! **************************************************************************** ! !> Mole fraction boundary condition for nonequilibrium extrapolation based. !! Default qVal=0.0. !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'spc_mole_fraction_noneq_expol', !! kind = 'mole_fraction', !! mole_fraction = 0.1 !! } !!} !!``` !! !! 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 spc_mole_fraction_noneq_expol( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmpAll_2f(nScalars) real(kind=rk) :: fEq_b, fEq_2f real(kind=rk) :: velAvg_2f(3), vel_b(3), eqVel_2f(3) real(kind=rk) :: usq_b, ucx_b, usq_2f, ucx_2f real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_2f real(kind=rk), dimension(varSys%nStateVars) :: num_dens_2f real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars) real(kind=rk) :: omega_fac, paramBInv integer :: iElem, iDir, QQ, iFld, vPos integer :: bcMoleFrac_pos integer :: elemPos, elemOff, neighPos, posInBuffer real(kind=rk) :: moleFrac_w(globBC%nElems(iLevel)), mass_dens_b type(mus_varSys_data_type), pointer :: fPtr integer :: nFields !> Inverse of lattice weight ar restPosition real(kind=rk) :: weight0Inv ! ------------------------------------------------------------------------- weight0Inv = 1.0_rk / layout%weight(layout%fStencil%restPosition) nFields = varSys%nStateVars ! position of molefraction spacetime variable in varSys bcMoleFrac_pos = me%bc_states%moleFrac%varPos ! mole fraction call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( & & varSys = varSys, & & time = sim_time, & & iLevel = iLevel, & & idx = me%bc_states%moleFrac & & %pntIndex%indexLvl(iLevel) & & %val(1:globBC%nElems(iLevel)), & & nVals = globBC%nElems(iLevel), & & res = moleFrac_w ) call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr ) ! resitivity coefficient of all species do iFld = 1, nFields !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:) resi_coeff(iFld, :) = fPtr%solverData%scheme & & %field(iFld)%fieldProp & & %species%resi_coeff(:) end do ! omega factor of LSE omega_fac = mixture%omega_diff * 0.5_rk paramBInv = 1.0_rk / mixture%paramB associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val, & & species => fPtr%solverData%scheme%field(:)%fieldProp%species, & & stencil => layout%fStencil ) QQ = layout%fStencil%QQ ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) ! Position of current boundary element in bcBuffer posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem ) elemOff = (posInBuffer-1)*nScalars ! Position of neighbor element in state array neighPos = me%neigh(iLevel)%posInState(1, iElem) do iFld = 1, nFields do iDir = 1, QQ ! Position of current field variable current dir in state array vPos = varSys%method%val(iFld)%state_varPos(iDir) ! Pre-Collision values of neighbor fluid element fTmpAll_2f(vPos) = state( & & neigh((idir-1)* nsize+ neighpos)+( ifld-1)* qq+ nscalars*0 ) ! State array of 1st fluid element !fTmpAll_1f(vPos) = bcBuffer( elemOff + vPos ) end do end do ! calculate spc density and velAvg from PDF of second fluid element call calcDensAndVelsFromPDF( nFields = nFields, & & iField = iField, & & mass_dens = mass_dens_2f, & & num_dens = num_dens_2f, & & velAvg = velAvg_2f, & & eqVel = eqVel_2f, & & varSys = varSys, & & pdf = fTmpAll_2f, & & stencil = layout%fStencil, & & species = species, & & resi_coeff = resi_coeff, & & omega_fac = omega_fac, & & paramBInv = paramBInv ) ! fluid node velocity square term of equilibrium function usq_2f = dot_product( velAvg_2f(1:3), velAvg_2f(1:3) ) * t2cs2inv ! density of current species at boundary node mass_dens_b = mixture%moleDens0LB * moleFrac_w(iElem) & & * species(iField)%molWeight ! velocity at boundary node extrapolated from 2nd fluid element vel_b = velAvg_2f ! boundary node velocity square term of equilibrium function usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv ! Position of this boundary element in state array elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) do iDir = 1, layout%fStencil%QQN if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then ! Compute equilibrium on fluid node ucx_2f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_2f ) fEq_2f = layout%weight(iDir) * mass_dens_2f(iField) & & * ( species(iField)%molWeigRatio + ucx_2f * cs2inv & & + ucx_2f * ucx_2f * t2cs4inv - usq_2f ) ! Compute equilibrium on boundary node ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b ) fEq_b = layout%weight(iDir) * mass_dens_b & & * ( species(iField)%molWeigRatio + ucx_b * cs2inv & & + ucx_b * ucx_b * t2cs4inv - usq_b ) ! update incomping link state( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_b + ( fTmpAll_2f(varPos(iDir)) - fEq_2f ) end if end do end do end associate end subroutine spc_mole_fraction_noneq_expol ! **************************************************************************** ! ! **************************************************************************** ! !> Open outflow boundary condition based on nonequilibrium extrapolation !! method. Default qVal = 0.0 !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_outflow', !! } !!} !!``` !! !! 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 spc_outflow( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmpAll_2f(nScalars), fTmpAll_3f(nScalars) real(kind=rk) :: fEq_b, fEq_2f real(kind=rk) :: velAvg_2f(3), vel_b(3), eqVel_2f(3) real(kind=rk) :: usq_2f, ucx_2f, usq_b, ucx_b real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_2f, mass_dens_3f real(kind=rk), dimension(varSys%nStateVars) :: num_dens_2f real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars) real(kind=rk) :: omega_fac, paramBInv, mass_dens_b integer :: iElem, iDir, QQ, iFld, vPos integer :: elemPos, neighPos_2f, neighPos_3f type(mus_varSys_data_type), pointer :: fPtr integer :: nFields ! ------------------------------------------------------------------------- nFields = varSys%nStateVars call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr ) ! resitivity coefficient of all species do iFld = 1, nFields !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:) resi_coeff(iFld, :) = fPtr%solverData%scheme & & %field(iFld)%fieldProp & & %species%resi_coeff(:) end do ! omega factor of LSE omega_fac = mixture%omega_diff * 0.5_rk paramBInv = 1.0_rk / mixture%paramB associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val, & & species => fPtr%solverData%scheme%field(:)%fieldProp%species, & & neighBufferPre_nNext => me%neigh(iLevel)%neighBufferPre_nNext, & & stencil => layout%fStencil ) QQ = layout%fStencil%QQ ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) ! Position of second neighbor element in state array neighPos_2f = me%neigh(iLevel)%posInState(1, iElem) ! Position of third neighbor element in state array neighPos_3f = me%neigh(iLevel)%posInState(2, iElem) ! Compute only mass density from third neighbor mass_dens_3f = 0.0_rk do iFld = 1, nFields do iDir = 1, QQ ! Position of current field variable current dir in state array vPos = varSys%method%val(iFld)%state_varPos(iDir) ! Pre-Collision values of second fluid element fTmpAll_2f(vPos) = state( & & neigh((idir-1)* nsize+ neighpos_2f)+( ifld-1)* qq+ nscalars*0 ) ! Pre-Collision values of third fluid element fTmpAll_3f(vPos) = state( & & neigh((idir-1)* nsize+ neighpos_3f)+( ifld-1)* qq+ nscalars*0 ) ! mass density of third neighbor mass_dens_3f(iFld) = mass_dens_3f(iFld) + fTmpAll_3f(vPos) end do end do ! calculate spc density and velAvg from PDF of second fluid call calcDensAndVelsFromPDF( nFields = nFields, & & iField = iField, & & mass_dens = mass_dens_2f, & & num_dens = num_dens_2f, & & velAvg = velAvg_2f, & & eqVel = eqVel_2f, & & varSys = varSys, & & pdf = fTmpAll_2f, & & stencil = layout%fStencil, & & species = species, & & resi_coeff = resi_coeff, & & omega_fac = omega_fac, & & paramBInv = paramBInv ) ! fluid node velocity square term of equilibrium function usq_2f = dot_product( velAvg_2f(1:3), velAvg_2f(1:3) ) * t2cs2inv ! second-order extrapolation of density of current species at boundary ! node from second and third fluid element mass_dens_b = ( 4.0_rk * mass_dens_2f(iField) - mass_dens_3f(iField) ) & & / 3.0_rk ! first-order extrapolation of velocity at boundary node from second ! fluid element vel_b = velAvg_2f ! second-order extrapolation of velocity at boundary node ! vel_b = ( 4.0_rk * velAvg_2f - velAvg_3f ) / 3.0_rk ! boundary node velocity square term of equilibrium function usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv ! Position of this boundary element in state array elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) do iDir = 1, layout%fStencil%QQN if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then ! Compute equilibrium on second fluid node ucx_2f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_2f ) fEq_2f = layout%weight(iDir) * mass_dens_2f(iField) & & * ( species(iField)%molWeigRatio + ucx_2f * cs2inv & & + ucx_2f * ucx_2f * t2cs4inv - usq_2f ) ! Compute equilibrium on boundary node ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b ) fEq_b = layout%weight(iDir) * mass_dens_b & & * ( species(iField)%molWeigRatio + ucx_b * cs2inv & & + ucx_b * ucx_b * t2cs4inv - usq_b ) ! update incomping direction state( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_b + ( fTmpAll_2f(varPos(iDir)) - fEq_2f ) end if end do end do end associate end subroutine spc_outflow ! **************************************************************************** ! ! **************************************************************************** ! !> Open outflow boundary condition for solvent based on nonequilibrium !! extrapolation. total moledens at boundary is enforced. !! method. Default qVal = 0.0 !! Usage !! ----- !!```lua !!boundary_condition = { !! { label = 'outlet', !! kind = 'spc_solvent_outflow', !! } !!} !!``` !! !! 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 spc_solvent_outflow( 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 ! -------------------------------------------------------------------- ! real(kind=rk) :: fTmpAll_2f(nScalars), fTmpAll_3f(nScalars) real(kind=rk) :: fEq_b, fEq_2f real(kind=rk) :: velAvg_2f(3), vel_b(3), eqVel_2f(3) real(kind=rk) :: usq_2f, ucx_2f, usq_b, ucx_b real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_2f, mass_dens_3f real(kind=rk), dimension(varSys%nStateVars) :: num_dens_2f, num_dens_3f real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars) real(kind=rk) :: omega_fac, paramBInv, mass_dens_b, num_dens_b, & & tot_NumDens_b integer :: iElem, iDir, QQ, iFld, vPos integer :: elemPos, neighPos_2f, neighPos_3f type(mus_varSys_data_type), pointer :: fPtr integer :: nFields ! ------------------------------------------------------------------------- nFields = varSys%nStateVars call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr ) ! resitivity coefficient of all species do iFld = 1, nFields !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:) resi_coeff(iFld, :) = fPtr%solverData%scheme & & %field(iFld)%fieldProp & & %species%resi_coeff(:) end do ! omega factor of LSE omega_fac = mixture%omega_diff * 0.5_rk paramBInv = 1.0_rk / mixture%paramB associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val, & & species => fPtr%solverData%scheme%field(:)%fieldProp%species, & & neighBufferPre_nNext => me%neigh(iLevel)%neighBufferPre_nNext, & & stencil => layout%fStencil ) QQ = layout%fStencil%QQ ! Calculate the density of current element do iElem = 1, globBC%nElems(iLevel) ! Position of second neighbor element in state array neighPos_2f = me%neigh(iLevel)%posInState(1, iElem) ! Position of third neighbor element in state array neighPos_3f = me%neigh(iLevel)%posInState(2, iElem) ! Compute only mass density from third neighbor mass_dens_3f = 0.0_rk do iFld = 1, nFields do iDir = 1, QQ ! Position of current field variable current dir in state array vPos = varSys%method%val(iFld)%state_varPos(iDir) ! Pre-Collision values of second fluid element fTmpAll_2f(vPos) = state( & & neigh((idir-1)* nsize+ neighpos_2f)+( ifld-1)* qq+ nscalars*0 ) ! Pre-Collision values of third fluid element fTmpAll_3f(vPos) = state( & & neigh((idir-1)* nsize+ neighpos_3f)+( ifld-1)* qq+ nscalars*0 )