! Copyright (c) 2013, 2016, 2020 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2013-2014 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2013-2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2013-2014 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2013-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2015-2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de> ! Copyright (c) 2016 Philipp Otte <otte@mathcces.rwth-aachen.de> ! Copyright (c) 2016-2017 Raphael Haupt <raphael.haupt@uni-siegen.de> ! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 2021-2022 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ****************************************************************************** ! !> author: Kannan Masilamani !! author: Jiaxing Qi !! This module provides the MUSUBI specific functions for calculating !! macroscopic quantities from the state variables. !! !! The depending common interface between MUSUBI and ATELES is defined in the !! tem_derived_module. The functionality for accessing a variable from the state !! and evaluating a lua function are also provided in the tem_derived module. !! !! Do not use get_Element or get_Point routines to update the state ! !! ! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de> ! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2014-2015, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2015 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2015 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2013-2014 Nikhil Anand <nikhil.anand@uni-siegen.de> ! Copyright (c) 2014, 2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2015, 2018, 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 2016 Verena Krupp <verena.krupp@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 COPYRIGHT HOLDERS AND CONTRIBUTORS "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 THE COPYRIGHT HOLDER 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. !-------------------------------------------- ! A O S - Array of structures layout new !------------------------------------------- ! Access to get_point value output ! Access to get_element value output module mus_derQuan_module use iso_c_binding, only: c_loc, c_ptr, c_f_pointer ! include treelm modules use tem_param_module, only: div1_2, div1_3, div1_54, div1_9, div3_4, & & sqrt3, cs2inv, cs2, t2cs2inv, t2cs4inv, & & cs4inv, q000, rho0 use env_module, only: rk, long_k, labelLen use tem_float_module, only: operator(.feq.), operator(.fge.), & & operator(.fle.) use tem_variable_module, only: tem_variable_type use tem_stencil_module, only: tem_stencilHeader_type use tem_topology_module, only: tem_levelOf use tem_time_module, only: tem_time_type use treelmesh_module, only: treelmesh_type use tem_logging_module, only: logUnit use tem_varSys_module, only: tem_varSys_type, tem_varSys_op_type, & & tem_varSys_append_derVar, & & tem_varSys_proc_point, & & tem_varSys_proc_element, & & tem_varSys_proc_setParams, & & tem_varSys_proc_getParams, & & tem_varSys_proc_setupIndices, & & tem_varSys_proc_getValOfIndex, & & tem_varSys_getPoint_dummy, & & tem_varSys_getElement_dummy, & & tem_varSys_setupIndices_dummy, & & tem_varSys_getValOfIndex_dummy, & & tem_varSys_setParams_dummy, & & tem_varSys_getParams_dummy use tem_aux_module, only: tem_abort use tem_operation_var_module, only: tem_evalMag_forElement, & & tem_evalMag_forPoint, & & tem_evalMag_fromIndex, & & tem_opVar_setupIndices, & & tem_get_new_varSys_data_ptr, & & tem_opVar_setParams, & & tem_opVar_getParams use tem_property_module, only: prp_hasBnd, prp_hasQval use tem_tools_module, only: tem_PositionInSorted use tem_debug_module, only: dbgUnit use tem_grow_array_module, only: grw_labelarray_type, append use tem_dyn_array_module, only: PositionOfVal ! include musubi modules use mus_source_type_module, only: mus_source_op_type use mus_pdf_module, only: pdf_data_type use mus_scheme_header_module, only: mus_scheme_header_type use mus_scheme_layout_module, only: mus_scheme_layout_type use mus_scheme_type_module, only: mus_scheme_type use mus_derivedQuantities_module2, only: secondMom, getShearRate use mus_varSys_module, only: mus_varSys_data_type, & & mus_varSys_solverData_type, & & mus_get_new_solver_ptr, & & mus_deriveVar_forPoint, & & mus_generic_varFromPDF_fromIndex, & & mus_generic_fromPDF_forElement, & & mus_derive_fromPDF use mus_stateVar_module, only: mus_accessVar_setupIndices, & & mus_stateVar_Fetch_fromIndex, & & mus_stateVar_Fetch_now_fromIndex, & & mus_access_stateFetch_ForElement, & & mus_access_stateFetch_now_ForElement use mus_operation_var_module, only: mus_opVar_setupIndices, & & mus_opVar_gradU_forElement, & & mus_opVar_vorticity_forElement, & & mus_opVar_QCriterion_forElement use mus_stateVar_module, only: mus_accessVar_setupIndices use mus_derVarPos_module, only: mus_derVarPos_type use mus_physics_module, only: mus_convertFac_type implicit none private public :: mus_append_derVar_fluid public :: derivePressure, derivePressure_fromIndex public :: deriveDensity public :: deriveDensity_fromIndex public :: deriveShearStress public :: deriveShearMag public :: deriveWSS3D public :: deriveWSS2D public :: deriveTemp public :: deriveShearRate public :: deriveBndForce ! equilbrium from macro uses different interface defined in ! mus_variable_module public :: deriveEquil_FromMacro public :: deriveRho_fromState public :: deriveVel_fromState public :: deriveVel_FromPreColState public :: deriveEq_fromState public :: deriveEquil_fromAux public :: deriveAux_fromState public :: deriveAux_fromState_d2q9 public :: deriveAux_fromState_d3q19 public :: deriveAux_fromState_d3q27 ! source variable public :: derive_absorbLayer public :: derive_force public :: derive_force1stOrd public :: derive_HRRCorrection ! Apply source add source term to state public :: applySrc_absorbLayer public :: applySrc_absorbLayer_MRT public :: applySrc_absorbLayerDyn public :: applySrc_absorbLayerDyn_MRT public :: applySrc_force public :: applySrc_force_MRT public :: applySrc_force_MRT_d2q9 public :: applySrc_force_MRT_d3q19 public :: applySrc_force_MRT_d3q27 public :: applySrc_turbChanForce public :: applySrc_turbChanForce_MRT public :: applySrc_turbChanForce_MRT_d2q9 public :: applySrc_turbChanForce_MRT_d3q19 public :: applySrc_turbChanForce_MRT_d3q27 public :: applySrc_force1stOrd public :: applySrc_HRRCorrection contains ! **************************************************************************** ! !> subroutine to add derive variables for weakly compressible LBM !! (schemekind = 'fluid') to the varsys. subroutine mus_append_derVar_fluid( varSys, solverData, schemeHeader, & & stencil, fldLabel, derVarName ) ! --------------------------------------------------------------------------- !> global variable system type(tem_varSys_type), intent(inout) :: varSys !> Contains pointer to solver data types type(mus_varSys_solverData_type), target, intent(in) :: solverData !> identifier of the scheme type(mus_scheme_header_type), intent(in) :: schemeHeader !> compute stencil defintion type(tem_stencilHeader_type), intent(in) :: stencil !> array of field label prefix. Size=nFields character(len=*), intent(in) :: fldLabel !> array of derive physical variables type(grw_labelarray_type), intent(inout) :: derVarName ! --------------------------------------------------------------------------- ! number of derive variables integer :: nDerVars, iVar, nComponents, addedPos, iIn logical :: wasAdded character(len=labelLen), allocatable :: input_varname(:) character(len=labelLen) :: varName procedure(tem_varSys_proc_point), pointer :: get_point => NULL() procedure(tem_varSys_proc_element), pointer :: get_element => NULL() procedure(tem_varSys_proc_setParams), pointer :: set_params => null() procedure(tem_varSys_proc_getParams), pointer :: get_params => null() procedure(tem_varSys_proc_setupIndices), pointer :: & & setup_indices => null() procedure(tem_varSys_proc_getValOfIndex), pointer :: & & get_valOfIndex => null() type(c_ptr) :: method_data character(len=labelLen), allocatable :: derVarName_loc(:) ! --------------------------------------------------------------------------- nullify(get_point, get_element, set_params, get_params, setup_indices, & & get_valOfIndex) nDerVars = 19 allocate(derVarName_loc(nDerVars)) derVarName_loc = [ 'fetch_pdf ', 'fetch_pdf_now ', & & 'pressure ', 'equilibrium ', & & 'non_equilibrium', 'kinetic_energy ', & & 'shear_stress ', 'strain_rate ', & & 'shear_rate ', 'wss ', & & 'momentum ', 'vel_mag ', & & 'bnd_force ', 'moment ', & & 'shear_mag ', 'temperature ', & & 'grad_velocity ', 'vorticity ', & & 'q_criterion ' ] do iVar = 1, nDerVars call append(derVarName, derVarName_loc(iVar)) ! set default pointers, overwrite if neccessary get_element => tem_varSys_getElement_dummy get_point => mus_deriveVar_forPoint setup_indices => mus_opVar_setupIndices get_valOfIndex => tem_varSys_getValOfIndex_dummy method_data = mus_get_new_solver_ptr(solverData) set_params => tem_varSys_setParams_dummy get_params => tem_varSys_getParams_dummy select case(trim(adjustl(derVarName_loc(iVar)))) case ('fetch_pdf') get_element => mus_access_stateFetch_ForElement get_valOfIndex => mus_stateVar_Fetch_fromIndex setup_indices => mus_accessVar_setupIndices nComponents = stencil%QQ allocate(input_varname(1)) input_varname(1) = 'pdf' case ('fetch_pdf_now') get_element => mus_access_stateFetch_now_ForElement get_valOfIndex => mus_stateVar_Fetch_now_fromIndex setup_indices => mus_accessVar_setupIndices nComponents = stencil%QQ allocate(input_varname(1)) input_varname(1) = 'pdf' case ('pressure') get_element => derivePressure get_valOfIndex => derivePressure_fromIndex nComponents = 1 allocate(input_varname(1)) input_varname(1) = 'density' case ('bnd_force') get_element => deriveBndForce nComponents = 3 allocate(input_varname(1)) input_varname(1) = 'pdf' case ('equilibrium') get_element => deriveEquil get_valOfIndex => deriveEquil_fromIndex nComponents = stencil%QQ allocate(input_varname(1)) input_varname(1) = 'pdf' case ('non_equilibrium') get_element => deriveNonEquil get_valOfIndex => deriveNonEquil_fromIndex nComponents = stencil%QQ allocate(input_varname(1)) input_varname(1) = 'fetch_pdf_now' case ('kinetic_energy') get_element => deriveKE get_valOfIndex => deriveKe_fromIndex nComponents = 1 allocate(input_varname(2)) input_varname(1) = 'density' input_varname(2) = 'velocity' case ('temperature') get_element => deriveTemp nComponents = 1 allocate(input_varname(0)) case ('shear_stress') nComponents = 6 get_element => deriveShearStress allocate(input_varname(1)) input_varname(1) = 'non_equilibrium' case ('strain_rate') nComponents = 6 get_element => deriveStrainRate get_valOfIndex => deriveStrainRate_fromIndex allocate(input_varname(1)) input_varname(1) = 'fetch_pdf_now' case ('shear_rate') get_element => deriveShearRate nComponents = 1 allocate(input_varname(1)) input_varname(1) = 'strain_rate' case ('wss') nComponents = 1 allocate(input_varname(1)) input_varname(1) = 'shear_stress' if (stencil%nDims == 2) then get_element => deriveWSS2D else if (stencil%nDims == 3) then get_element => deriveWSS3D else write(logUnit(1),*) 'WARNING: WSS does not support 1D' end if case ('momentum') get_element => deriveMomentum get_valOfIndex => deriveMomentum_fromIndex nComponents = 3 allocate(input_varname(2)) input_varname(1) = 'density' input_varname(2) = 'velocity' case ('grad_velocity') get_element => mus_opVar_gradU_forElement nComponents = 9 allocate(input_varname(1)) input_varname(1) = 'velocity' case ('vorticity') get_element => mus_opVar_vorticity_forElement nComponents = 3 allocate(input_varname(1)) input_varname(1) = 'velocity' case ('q_criterion') get_element => mus_opVar_QCriterion_forElement nComponents = 1 allocate(input_varname(1)) input_varname(1) = 'velocity' case ('vel_mag') get_element => tem_evalMag_forElement get_point => tem_evalMag_forPoint get_valOfIndex => tem_evalMag_fromIndex setup_indices => tem_opVar_setupIndices set_params => tem_opVar_setParams get_params => tem_opVar_getParams method_data = tem_get_new_varSys_data_ptr(method_data) nComponents = 1 allocate(input_varname(1)) input_varname(1) = 'velocity' case ('shear_mag') get_element => deriveShearMag nComponents = 1 allocate(input_varname(1)) input_varname(1) = 'shear_stress' case ('moment') !ONLY FOR 2D get_element => deriveMoment nComponents = 9 allocate(input_varname(1)) input_varname(1) = 'pdf' case default write(logUnit(1),*) 'WARNING: Unknown variable: '//& & trim(derVarName_loc(iVar)) cycle !go to next variable end select ! update variable names with field label varname = trim(fldLabel)//trim(adjustl(derVarName_loc(iVar))) do iIn = 1, size(input_varname) input_varname(iIn) = trim(fldLabel)//trim(input_varname(iIn)) end do ! append variable to varSys call tem_varSys_append_derVar( me = varSys, & & varName = trim(varname), & & nComponents = nComponents, & & input_varname = input_varname, & & method_data = method_data, & & get_point = get_point, & & get_element = get_element, & & set_params = set_params, & & get_params = get_params, & & setup_indices = setup_indices, & & get_valOfIndex = get_valOfIndex, & & pos = addedPos, & & wasAdded = wasAdded ) if (wasAdded) then write(logUnit(10),*) ' Appended variable:'//trim(varname) else if (addedpos < 1) then write(logUnit(1),*) 'Error: variable '//trim(varname)// & & ' is not added to variable system' end if deallocate(input_varname) end do end subroutine mus_append_derVar_fluid ! **************************************************************************** ! ! ****************************************************************************** ! ! Subroutines with common interface for the function pointers ! ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of density !! This routine sets the function Pointer for density calcualtion and calls !! the generice get Element from PDF routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveDensity(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! -------------------------------------------------------------------- ! fnCalcPtr => mus_derivedensity call mus_generic_fromPDF_forElement( & & fun = fun, & & varSys = varSys, & & elempos = elempos, & & tree = tree, & & time = time, & & nVals = nElems, & & fnCalcPtr = fnCalcPtr, & & nDofs = nDofs, & & res = res ) end subroutine deriveDensity ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate the temperature of a given set of elements (sum up all links). !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveTemp(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! res( 1:nElems ) = div1_3 end subroutine deriveTemp ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate the pressure of a given set of elements (sum up all links). !! !! Pressure calculation according to the isentropic equation of state for !! the LBM \( p = \rho c_s^2 \) !! with the calculation of density as in deriveDensity !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine derivePressure(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! integer :: dens_pos ! -------------------------------------------------------------------- ! ! position of density in glob system dens_pos = fun%input_varPos(1) ! get density variable from auxField call varSys%method%val(dens_pos)%get_element( varSys = varSys, & & elemPos = elemPos, & & time = time, & & tree = tree, & & nElems = nElems, & & nDofs = nDofs, & & res = res ) ! convert density to pressure res(1:nElems) = res(1:nElems) * cs2 end subroutine derivePressure ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate the force on the boundary of a given set of elements !! !! Force calculation in LBM: !! \[ F = \sum_\alpha{\alpha \ne 0} \vec{e}_{\alphabar} (2f_\alphabar !! + 6\omega_\alpha \vec{e}_\alpha \rho \vec{u_b}) \] !! \n where, \n !! \alpha = direction towards boundary !! \alphabar = direction opposite to boundary !! \vec{u_b} = velocity of the obstacle cell. !! For stationary object, this should be zero. !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveBndForce(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr integer :: iElem, posInBndID ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) res = 0.0_rk ! compute bnd_force only if shape is defined as boundary kind do iElem = 1, nElems ! compute force only if current element has boundary property if ( btest( tree%elemPropertyBits( elemPos(iElem) ), & & prp_hasBnd ) ) then ! position of current element in boundary_ID list posInBndID = fPtr%solverData%geometry%posInBndID( elempos(iElem) ) if (posInBndID > 0) then res( (iElem-1)*fun%nComponents+1: iElem*fun%nComponents ) & & = fPtr%solverData%geometry%bndForce(posInBndID, :) end if ! posInBndID > 0 end if ! element has_boundary end do ! iElem end subroutine deriveBndForce ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate momentum from density and velocity stored in auxField !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveMomentum(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! !> Function pointer to perform specific operation. integer :: statePos, iElem, iLevel, elemOff type(mus_varSys_data_type), pointer :: fPtr integer :: dens_pos, vel_pos(3) real(kind=rk) :: rho ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) res = 0.0_rk associate( scheme => fPtr%solverData%scheme, & & levelPointer => fPtr%solverData%geometry%levelPointer, & & auxField => fPtr%solverData%scheme%auxField ) do iElem = 1, nElems ! if state array is defined level wise then use levelPointer(pos) ! to access state array statePos = levelPointer( elemPos(iElem) ) iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) ) ! element offset for auxField elemoff = (statePos-1)*varSys%nAuxScalars ! position of density in auxField array dens_pos = varSys%method%val( fun%input_varPos(1) ) & & %auxField_varPos(1) ! position of velocity in auxField array vel_pos = varSys%method%val( fun%input_varPos(2) ) & & %auxField_varPos(1:3) ! density rho = auxField(iLevel)%val( elemOff + dens_pos ) ! compute and store momentum res((iElem-1)*3 + 1) = rho * auxField(iLevel)%val(elemOff + vel_pos(1)) res((iElem-1)*3 + 2) = rho * auxField(iLevel)%val(elemOff + vel_pos(2)) res((iElem-1)*3 + 3) = rho * auxField(iLevel)%val(elemOff + vel_pos(3)) end do ! iElem end associate end subroutine deriveMomentum ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of equlibrium !! This routine sets the function Pointer for equlibrium calcualtion and calls !! the generice get Element from PDF routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveEquil(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! -------------------------------------------------------------------- ! fnCalcPtr => mus_deriveEquil call mus_generic_fromPDF_forElement( & & fun = fun, & & varSys = varSys, & & elempos = elempos, & & tree = tree, & & time = time, & & nVals = nElems, & & fnCalcPtr = fnCalcPtr, & & nDofs = nDofs, & & res = res ) end subroutine deriveEquil ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of NonEquil !! This routine sets the function Pointer for NonEquil calcualtion and calls !! the generice get Element from PDF routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveNonEquil(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! -------------------------------------------------------------------- ! fnCalcPtr => mus_deriveNonEquil call mus_generic_fromPDF_forElement( & & fun = fun, & & varSys = varSys, & & elempos = elempos, & & tree = tree, & & time = time, & & nVals = nElems, & & fnCalcPtr = fnCalcPtr, & & nDofs = nDofs, & & res = res ) end subroutine deriveNonEquil ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of StrainRate !! This routine sets the function Pointer for StrainRate calcualtion and calls !! the generice get Element from PDF routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveStrainRate(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! -------------------------------------------------------------------- ! fnCalcPtr => mus_deriveStrainRate call mus_generic_fromPDF_forElement( & & fun = fun, & & varSys = varSys, & & elempos = elempos, & & tree = tree, & & time = time, & & nVals = nElems, & & fnCalcPtr = fnCalcPtr, & & nDofs = nDofs, & & res = res ) end subroutine deriveStrainRate ! ****************************************************************************** ! ! ****************************************************************************** ! !> author: Jiaxing Qi !! Calculate the shear rate !! !! The Shear Rate is defined as !! \[ !! \dot{\gamma} = \sqrt{ 2D_{II} } !! \] !! where \( D_{II} \) is the second invariant of the strain rate tensor and !! defined as !! \[ !! D_{II} = \sum^{l}_{\alpha,\beta=l} S_{\alpha\beta} S_{\alpha\beta} !! \] !! where \( S_{\alpha\beta} \) is the strain rate tensor. !! recursive subroutine deriveShearRate(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! integer :: posStrain, iElem real(kind=rk) :: strain(6*nElems) ! -------------------------------------------------------------------- ! posStrain = fun%input_varPos(1) call varSys%method%val(posStrain)%get_element( varSys = varSys, & & elemPos = elemPos, & & time = time, & & tree = tree, & & nElems = nElems, & & nDofs = nDofs, & & res = strain ) do iElem = 1, nElems res( iElem ) = getShearRate( strain( (iElem-1)*6+1:iElem*6 ) ) end do ! iElem end subroutine deriveShearRate ! ****************************************************************************** ! ! ****************************************************************************** ! !> author: Jiaxing Qi !! Calculate the deviatoric shear stress for Newtonian fluid !! (exclude pressure) (no mixtures).\n !! Shear Stress depends on variable: nonEquilibirium !! !! The formula is !! \[ !! \tau_{\alpha \beta}= !! -(1-\frac{\omega}{2}) \sum_{i} f^{neq}_{i} c_{i\alpha} c_{i\beta} !! \] !! where \( \tau_{\alpha \beta}\) is the stress !! in the \(\beta\)-direction on a face normal to the \(\alpha\)-axis,\n !! \( f^{neq}_i = f_i - f^{eq}_i\) is the non-equilibirium density.\n !! For more information, please refer to:\n !! Krueger T, Varnik F, Raabe D. Shear stress in lattice Boltzmann !! simulations. Physical Review E. 2009;79(4):1-14.\n !! Thus this variable is dependent on another variable: nonEquilibirium. !! !! For multi-level mesh, Omega on finer level needs to be adjusted in order to !! get the correct shearstress calculation.\n !! First, we defines c as the dx ratio between finer and coarse level.\n !! \( c={ \Delta dx }_{ c }/{ \Delta dx }_{ f } \) !! Then the viscosity on the different levels must satisfy:\n !! \( \frac { { \nu }_{ f } }{ { \nu }_{ c } } =c \) !! This constrain leads to a relationship of omega on different levels:\n !! \( {\omega}_f = \frac {1}{ {\lambda}(\frac{1}{{\omega}_c}-0.5)+0.5 } \) !! For more information, please refer to:\n !! Manuel H, Harald K, Joerg B, Sabine R. Aeroacoustic validation of the !! lattice boltzmann method on non-uniform grids. ECCOMAS 2012 !! recursive subroutine deriveShearStress(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme integer :: posNonEq, nCompsNonEq, iLevel, iElem real(kind=rk) :: omega real(kind=rk), allocatable :: nonEq(:) real(kind=rk), allocatable :: tau(:) ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme posNonEq = fun%input_varPos(1) nCompsNonEq = varSys%method%val( posNonEq )%nComponents allocate( nonEq( nCompsNonEq * nElems ) ) allocate( tau( fun%nComponents ) ) ! calculate nonEq for all elements call varSys%method%val(posNonEq)%get_element( varSys = varSys, & & elemPos = elemPos, & & time = time, & & tree = tree, & & nElems = nElems, & & nDofs = nDofs, & & res = nonEq ) do iElem = 1, nElems iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) ) ! get the correct omega value omega = scheme%field(1)%fieldProp%fluid%viscKine%omLvl(iLevel)%val(iElem) ! compute second moment of nonEq tau(:) = secondMom( cxcx = scheme%layout%fStencil%cxcx(:,:), & & f = nonEq( (iElem-1)*nCompsNonEq+1 & & :iElem*nCompsNonEq ), & & QQ = scheme%layout%fStencil%QQ ) res( (iElem-1)*fun%nComponents+1: iElem*fun%nComponents ) = & & tau(:) * ( 0.5_rk * omega - 1._rk ) end do ! iElem deallocate( nonEq ) deallocate( tau ) end subroutine deriveShearStress ! ****************************************************************************** ! ! ****************************************************************************** ! !> author: Jiaxing Qi !! Calculate the wall shear stress (WSS) of a given element with the given !! input !! !! The wall shear stress is shear stress exerted on the wall. Since it is well !! known that shear stress gets its maximum value on the wall, here we can !! directly calculate the principle stress of the stress tensor instead of !! multiplying the stress tensor by normal vector of the plane. !! The principle stresses are just the eigenvalues of the stress tensor. To !! get those eigenvalues, we need solve the characteristic equation:\n !! \( {\tau}^3 + a_{2}{\tau}^2 + a_1{\tau} + a_0 = 0 \) !! where \( a_2 = \tau_{ii} \) !! \( a_1 = (\tau_{ii}\tau_{jj} - \tau_{ij}\tau_{ji}) / 2 \) !! \( a_0 = det(\tau_{ij}) \). Here Einstein notation is used. !! recursive subroutine deriveWSS3D(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! ! local variables integer :: iElem, shear_pos real(kind=rk) :: shear(6*nElems) ! shear stress of all ! -------------------------------------------------------------------- ! ! wall shear stress is depend on shear stress. get its position in globSys shear_pos = fun%input_varPos(1) ! first compute shear stress, save it to array tau ! the shear tensor has the following order: x, y, z, xy, xz, yz call varSys%method%val(shear_pos)%get_element( varSys = varSys, & & elemPos = elemPos, & & time = time, & & tree = tree, & & nElems = nElems, & & nDofs = nDofs, & & res = shear ) do iElem = 1, nElems res( iElem ) = getWSS( shear( (iElem-1)*6+1 : iElem*6 ) ) end do ! iElem end subroutine deriveWSS3D ! ****************************************************************************** ! ! ****************************************************************************** ! recursive subroutine deriveWSS2D(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! ! local variables integer :: iElem, shear_pos real(kind=rk) :: a0, a1 real(kind=rk) :: tau(6*nElems) ! shear stress of all ! -------------------------------------------------------------------- ! ! wall shear stress is depend on shear stress. get its position in globSys shear_pos = fun%input_varPos(1) ! first compute shear stress, save it to array tau ! the shear tensor has the following order: x, y, z, xy, xz, yz call varSys%method%val(shear_pos)%get_element( varSys = varSys, & & elemPos = elemPos, & & time = time, & & tree = tree, & & nElems = nElems, & & nDofs = nDofs, & & res = tau ) do iElem = 1, nElems a1 = - tau((iElem-1)*6+1) - tau((iElem-1)*6+2) a0 = tau((iElem-1)*6+1) * tau((iElem-1)*6+2) & & - tau((iElem-1)*6+3) * tau((iElem-1)*6+3) res(iElem) = sqrt( a1 * a1 - 4._rk * a0 ) * div1_2 end do ! iElem end subroutine deriveWSS2D ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate the shear stress magnitude of a given element number with the !! given !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveShearMag(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! integer :: shear_pos, iElem real(kind=rk) :: shear(6*nElems) ! shear stress of all real(kind=rk) :: tau(6) ! -------------------------------------------------------------------- ! shear_pos = fun%input_varPos(1) ! first compute shear stress, save it to array tau ! the shear tensor has the following order: x, y, z, xy, xz, yz call varSys%method%val(shear_pos)%get_element( varSys = varSys, & & elemPos = elemPos, & & time = time, & & tree = tree, & & nElems = nElems, & & nDofs = nDofs, & & res = shear ) do iElem = 1, nElems tau = shear( (iElem-1)*6+1 : iElem*6 ) ! Von mises-criterion (see wikipedia article on von Mises yield criterion) ! wss = (tau_xx-tau_yy)**2 + (tau_yy-tau_zz)**2 + (tau_xx-tau_zz)**2 & ! & + 6._rk*( tau_xy*tau_xy + tau_yz*tau_yz + tau_xz*tau_xz) ! wss = sqrt(0.5_rk*wss ) res( iElem ) = sqrt( 0.5_rk*( (tau(1)-tau(2))**2 & & + (tau(2)-tau(3))**2 & & + (tau(1)-tau(3))**2 & & + 6._rk * ( tau(4)*tau(4) & & + tau(5)*tau(5) & & + tau(6)*tau(6) ))) end do ! iElem end subroutine deriveShearMag ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate kinetic energy from density and velocity in auxField !! This routine sets the function Pointer for kinetic energy calcualtion and !! calls the generice get Element from PDF routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveKE(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! !> Function pointer to perform specific operation. integer :: statePos, iElem, iLevel, elemOff type(mus_varSys_data_type), pointer :: fPtr integer :: dens_pos, vel_pos(3) real(kind=rk) :: rho, vel(3) ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) res = 0.0_rk associate( scheme => fPtr%solverData%scheme, & & levelPointer => fPtr%solverData%geometry%levelPointer, & & auxField => fPtr%solverData%scheme%auxField ) do iElem = 1, nElems ! if state array is defined level wise then use levelPointer(pos) ! to access state array statePos = levelPointer( elemPos(iElem) ) iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) ) ! element offset for auxField elemoff = (statePos-1)*varSys%nAuxScalars ! position of density in auxField array dens_pos = varSys%method%val( fun%input_varPos(1) ) & & %auxField_varPos(1) ! position of velocity in auxField array vel_pos = varSys%method%val( fun%input_varPos(2) ) & & %auxField_varPos(1:3) ! density rho = auxField(iLevel)%val( elemOff + dens_pos ) ! velocity vel(1) = auxField(iLevel)%val(elemOff + vel_pos(1)) vel(2) = auxField(iLevel)%val(elemOff + vel_pos(2)) vel(3) = auxField(iLevel)%val(elemOff + vel_pos(3)) ! compute and store kinetic energy res(iElem) = sum( vel(:)*vel(:) ) * 0.5_rk * rho end do ! iElem end associate end subroutine deriveKE ! ****************************************************************************** ! ! ****************************************************************************** ! !> Derive absorb layer variable defined as a source term. recursive subroutine derive_absorbLayer(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! call tem_abort('Not implemented') end subroutine derive_absorbLayer ! ****************************************************************************** ! ! ****************************************************************************** ! !> Derive external force variable defined as a source term. !! It evaluates spacetime function defined in lua file for force variable !! and convert it to state value which is to be added to the state recursive subroutine derive_force(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: forceField(nElems*3) real(kind=rk) :: G(3), velocity(3), ucx, uMinusCX(3) integer :: iElem, iDir, QQ, nScalars, posInTotal, elemOff integer :: vel_pos(3), iLevel real(kind=rk) :: omegaKine, omegaBulk, discForce real(kind=rk) :: forceTerm(27), momForce(27), s_mrt(27) real(kind=rk) :: mInvXOmega(27,27) ! -------------------------------------------------------------------- ! !call tem_abort('Not implemented') ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%input_varPos(1) )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%input_varPos(2))%get_element( & & varSys = varSys, & & elemPos = elemPos, & & time = time, & & tree = tree, & & nElems = nElems, & & nDofs = nDofs, & & res = forceField ) ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of velocity variable in auxField vel_pos = varSys%method%val(scheme%derVarPos(1)%velocity)%auxField_varPos(1:3) do iElem = 1, nElems ! get iLevel for element iLevel = tem_levelOf( tree%treeID( elemPos(iElem ) ) ) posInTotal = fPtr%solverData%geometry%levelPointer( elemPos(iElem) ) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain velocity from auxField velocity(1) = scheme%auxField(iLevel)%val(elemOff + vel_pos(1)) velocity(2) = scheme%auxField(iLevel)%val(elemOff + vel_pos(2)) velocity(3) = scheme%auxField(iLevel)%val(elemOff + vel_pos(3)) ! convert physical to lattice. ! force field on current element ! For incompressible model: this forceField should be divided by rho0. ! Since rho0 =1, this term is also valid for incompressible model G = forceField((iElem-1)*3+1 : iElem*3) & & / fPtr%solverData%physics%fac(iLevel)%body_force ! get the correct omega value omegaKine = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(1:QQ) = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=omegaKine, omegaBulk=omegaBulk, QQ=QQ) ! M^-1 * (I-0.5 S) s_mrt(1:QQ) = 1.0_rk - 0.5_rk * s_mrt(1:QQ) do iDir = 1, QQ mInvXOmega(1:QQ,iDir) = scheme%layout%moment%toPDF%A(1:QQ,iDir) & & * s_mrt(iDir) end do ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} do iDir = 1, QQ ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), & & velocity ) uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - velocity forceTerm(iDir) = scheme%layout%weight(iDir) & & * dot_product( uMinusCx * cs2inv & & + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, G ) end do ! Force moments: M * F !do iDir = 1, QQ ! momForce(iDir) = sum(scheme%layout%moment%toMoments%A(iDir,:) * forceTerm) !end do momForce(1:QQ) = matmul( scheme%layout%moment%toMoments%A(1:QQ, 1:QQ), & & forceTerm(1:QQ) ) !discForce = matmul( omegaTerm, forceTerm ) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = dot_product(mInvXOmega(iDir,1:QQ), momForce(1:QQ)) res((iElem-1)*fun%nComponents + iDir) = discForce end do end do !iElem end subroutine derive_force ! ****************************************************************************** ! ! ************************************************************************** ! !> Derive external force variable defined as a source term. !! It evaluates spacetime function defined in lua file for force variable !! and convert it to state value which is to be added to the state !! @todo KM: Not use we need seperate force for incompressible fluid model recursive subroutine derive_force1stOrd(fun, varsys, elempos, time, tree, & & nElems, nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! call tem_abort('Not implemented') end subroutine derive_force1stOrd ! ************************************************************************** ! ! ****************************************************************************** ! !> Derive external force variable defined as a source term. !! It evaluates spacetime function defined in lua file for force variable !! and convert it to state value which is to be added to the state recursive subroutine derive_HRRCorrection(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! integer :: statePos, iElem, iLevel, iComp, QQ, iVar, HRRVar type(mus_varSys_data_type), pointer :: fPtr ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) associate (source => fPtr%solverData%scheme%field(1)%source) HRRVar = -1 ! res is always AOS layout !write(*,*) 'Source nVals = ', source%varDict%nVals do iVar = 1, source%varDict%nVals !write(*,*) 'Source varDictKey = ', trim(source%varDict%val(iVar)%key) if (trim(source%varDict%val(iVar)%key) == 'hrr_correction') then HRRVar = iVar end if end do if (HRRVar == -1) then call tem_abort ('Source variable HRR_Correction not found') end if res = 0.0_rk QQ = fun%nComponents do iElem = 1, nElems ! if state array is defined level wise then use levelPointer(pos) ! to access state array statePos = fPtr%solverData%geometry%levelPointer( elemPos(iElem) ) iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) ) do iComp = 1, QQ res( (iElem-1)*QQ+iComp ) = source%method(HRRVar)%elemLvl(iLevel) & & %HRR_Corr%phi((statePos-1) * QQ + iComp) end do !iComp end do !iElem end associate end subroutine derive_HRRCorrection ! ****************************************************************************** ! ! **************************************************************************** ! !> Update state with source variable "absorb_layer". !! absorb_layer is used to absorb the flow and gradually reduce the flow !! quantities like pressure and velocity to a fixed value. !! It is based on: !! Xu, H., & Sagaut, P. (2013). Analysis of the absorbing layers for the !! weakly-compressible lattice Boltzmann methods. Journal of Computational !! Physics, 245(x), 14–42. !! Jacob, J.; Sagaut, P. (2019): Solid wall and open boundary conditions in !! hybrid recursive regularized lattice Boltzmann method for compressible !! flows. In Physics of Fluids 31 (12), p. 126103. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_absorbLayer( fun, inState, outState, neigh, auxField, & & nPdfSize, iLevel, varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: spongeField(fun%elemLvl(iLevel)%nElems*5) real(kind=rk) :: dens, vel(3), ucx, uMinusCX(3), sponge_velTerm integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: sigma, dens_ref, vel_ref(3), sponge_vel(3), sponge_dens real(kind=rk) :: inv_rho_phy, inv_vel_phy real(kind=rk) :: omega, omega_fac ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- !call tem_abort('Error: Absorb layer is not yet implemented') ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = spongeField ) inv_rho_phy = 1.0_rk / fPtr%solverData%physics%fac(iLevel)%press * cs2inv inv_vel_phy = 1.0_rk / fPtr%solverData%physics%fac(iLevel)%vel ! target pressure and velocity in lattice unit dens_ref = fun%absLayer%config%target_pressure * inv_rho_phy vel_ref(1:3) = fun%absLayer%config%target_velocity(1:3) * inv_vel_phy ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField vel(1) = auxField(elemOff + vel_pos(1)) vel(2) = auxField(elemOff + vel_pos(2)) vel(3) = auxField(elemOff + vel_pos(3)) ! get the correct omega value omega = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) omega_fac = 1.0_rk - omega * 0.5_rk ! SpongeField contains: spongeStrength sigma = spongeField(iElem) sponge_vel(:) = -sigma*(vel - vel_ref) sponge_dens = -sigma*(dens - dens_ref) ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} do iDir = 1, QQ ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), & & vel ) uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - vel sponge_velTerm = dot_product( uMinusCx * cs2inv & & + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, sponge_vel ) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) & & + omega_fac * scheme%layout%weight( iDir ) & & * (sponge_dens + sponge_velTerm) end do end do !iElem end subroutine applySrc_absorbLayer ! **************************************************************************** ! ! **************************************************************************** ! !> Update state with source variable "absorb_layer". !! absorb_layer is used to absorb the flow and gradually reduce the flow !! quantities like pressure and velocity to a fixed value. !! It is based on: !! Xu, H., & Sagaut, P. (2013). Analysis of the absorbing layers for the !! weakly-compressible lattice Boltzmann methods. Journal of Computational !! Physics, 245(x), 14–42. !! Jacob, J.; Sagaut, P. (2019): Solid wall and open boundary conditions in !! hybrid recursive regularized lattice Boltzmann method for compressible !! flows. In Physics of Fluids 31 (12), p. 126103. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_absorbLayer_MRT( fun, inState, outState, neigh, auxField,& & nPdfSize, iLevel, varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: spongeField(fun%elemLvl(iLevel)%nElems*5) real(kind=rk) :: dens, vel(3), ucx, uMinusCX(3) integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: sigma, dens_ref, vel_ref(3), sponge_vel(3), sponge_dens real(kind=rk) :: inv_rho_phy, inv_vel_phy real(kind=rk) :: omegaKine, omegaBulk, discForce, sponge_velTerm real(kind=rk) :: momForce(27), s_mrt(27) real(kind=rk) :: mInvXOmega(27,27), sponge_Term(27) ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- !call tem_abort('Error: Absorb layer is not yet implemented') ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = spongeField ) inv_rho_phy = 1.0_rk / fPtr%solverData%physics%fac(iLevel)%press * cs2inv inv_vel_phy = 1.0_rk / fPtr%solverData%physics%fac(iLevel)%vel ! target pressure and velocity in lattice unit dens_ref = fun%absLayer%config%target_pressure * inv_rho_phy vel_ref(1:3) = fun%absLayer%config%target_velocity(1:3) * inv_vel_phy ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField vel(1) = auxField(elemOff + vel_pos(1)) vel(2) = auxField(elemOff + vel_pos(2)) vel(3) = auxField(elemOff + vel_pos(3)) ! get the correct omega value omegaKine = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(1:QQ) = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=omegaKine, omegaBulk=omegaBulk, QQ=QQ) ! M^-1 * (I-0.5 S) s_mrt(1:QQ) = 1.0_rk - 0.5_rk * s_mrt(1:QQ) do iDir = 1, QQ mInvXOmega(1:QQ,iDir) = scheme%layout%moment%toPDF%A(1:QQ,iDir) & & * s_mrt(iDir) end do ! SpongeField contains: spongeStrength sigma = spongeField(iElem) ! Sponge factor for density and velocity field sponge_dens = -sigma*(dens - dens_ref) sponge_vel(:) = -sigma*(vel - vel_ref) ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} do iDir = 1, QQ ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), vel ) uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - vel sponge_velTerm = dot_product( uMinusCx * cs2inv & & + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, sponge_vel ) sponge_Term(iDir) = scheme%layout%weight(iDir) & & * ( sponge_dens + sponge_velTerm ) end do ! Force moments: M * F momForce(1:QQ) = matmul( scheme%layout%moment%toMoments%A(1:QQ, 1:QQ), & & sponge_Term(1:QQ) ) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = dot_product(mInvXOmega(iDir,1:QQ), momForce(1:QQ)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_absorbLayer_MRT ! **************************************************************************** ! ! **************************************************************************** ! !> Update state with source variable "absorb_layer". !! absorb_layer is used to absorb the flow and gradually reduce the flow !! quantities like pressure and velocity to a fixed value. !! It is based on: !! Xu, H., & Sagaut, P. (2013). Analysis of the absorbing layers for the !! weakly-compressible lattice Boltzmann methods. Journal of Computational !! Physics, 245(x), 14–42. !! Jacob, J.; Sagaut, P. (2019): Solid wall and open boundary conditions in !! hybrid recursive regularized lattice Boltzmann method for compressible !! flows. In Physics of Fluids 31 (12), p. 126103. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_absorbLayerDyn( fun, inState, outState, neigh, auxField, & & nPdfSize, iLevel, varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr real(kind=rk) :: spongeField(fun%elemLvl(iLevel)%nElems*5) real(kind=rk) :: dens, vel(3), ucx, uMinusCX(3), sponge_velTerm integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: sigma, sponge_vel(3), sponge_dens real(kind=rk) :: omega, omega_fac ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- !call tem_abort('Error: Absorb layer is not yet implemented') ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = spongeField ) ! constant parameter nScalars = varSys%nScalars ! Position of velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) associate( layout => fPtr%solverData%scheme%layout, & & fluid => fPtr%solverData%scheme%field(1)%fieldProp%fluid, & & dynAvg => fun%elemLvl(iLevel)%dynAvg ) QQ = layout%fStencil%QQ !$omp parallel do schedule(static), private( posInTotal, vel, ucx, uMinusCx, omega ) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField vel(1) = auxField(elemOff + vel_pos(1)) vel(2) = auxField(elemOff + vel_pos(2)) vel(3) = auxField(elemOff + vel_pos(3)) ! get the correct omega value omega = fluid%viscKine%omLvl(iLevel)%val(posInTotal) omega_fac = 1.0_rk - omega * 0.5_rk ! SpongeField contains: spongeStrength sigma = spongeField(iElem) sponge_dens = -sigma*(dens - dynAvg%dens(iElem)) sponge_vel(1) = -sigma*(vel(1) - dynAvg%velX(iElem)) sponge_vel(2) = -sigma*(vel(2) - dynAvg%velY(iElem)) sponge_vel(3) = -sigma*(vel(3) - dynAvg%velZ(iElem)) ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} do iDir = 1, QQ ucx = dot_product( layout%fStencil%cxDirRK(:, iDir), vel ) uMinusCx = layout%fStencil%cxDirRK(:, iDir) - vel sponge_velTerm = dot_product( uMinusCx * cs2inv & & + ucx * layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, sponge_vel ) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) & & + omega_fac * layout%weight( iDir ) & & * (sponge_dens + sponge_velTerm) end do end do !iElem !$omp end parallel do end associate end subroutine applySrc_absorbLayerDyn ! **************************************************************************** ! ! **************************************************************************** ! !> Update state with source variable "absorb_layer". !! absorb_layer is used to absorb the flow and gradually reduce the flow !! quantities like pressure and velocity to a fixed value. !! It is based on: !! Xu, H., & Sagaut, P. (2013). Analysis of the absorbing layers for the !! weakly-compressible lattice Boltzmann methods. Journal of Computational !! Physics, 245(x), 14–42. !! Jacob, J.; Sagaut, P. (2019): Solid wall and open boundary conditions in !! hybrid recursive regularized lattice Boltzmann method for compressible !! flows. In Physics of Fluids 31 (12), p. 126103. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_absorbLayerDyn_MRT( fun, inState, outState, neigh, & & auxField, nPdfSize, iLevel, varSys, & & time, phyConvFac, derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr real(kind=rk) :: spongeField(fun%elemLvl(iLevel)%nElems*5) real(kind=rk) :: dens, vel(3), ucx, uMinusCX(3) integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: sigma, sponge_vel(3), sponge_dens real(kind=rk) :: omegaKine, omegaBulk, discForce, sponge_velTerm real(kind=rk) :: momForce(27), s_mrt(27) real(kind=rk) :: mInvXOmega(27,27), sponge_Term(27) ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- !call tem_abort('Error: Absorb layer is not yet implemented') ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = spongeField ) ! constant parameter nScalars = varSys%nScalars ! Position of velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) associate( layout => fPtr%solverData%scheme%layout, & & fluid => fPtr%solverData%scheme%field(1)%fieldProp%fluid, & & dynAvg => fun%elemLvl(iLevel)%dynAvg ) QQ = layout%fStencil%QQ omegaBulk = fluid%omegaBulkLvl(iLevel) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField vel(1) = auxField(elemOff + vel_pos(1)) vel(2) = auxField(elemOff + vel_pos(2)) vel(3) = auxField(elemOff + vel_pos(3)) ! get the correct omega value omegaKine = fluid%viscKine%omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(1:QQ) = fluid%mrtPtr(omegaKine=omegaKine, omegaBulk=omegaBulk, QQ=QQ) ! M^-1 * (I-0.5 S) s_mrt(1:QQ) = 1.0_rk - 0.5_rk * s_mrt(1:QQ) do iDir = 1, QQ mInvXOmega(1:QQ,iDir) = layout%moment%toPDF%A(1:QQ,iDir) & & * s_mrt(iDir) end do ! SpongeField contains: spongeStrength sigma = spongeField(iElem) ! Sponge factor for density and velocity field sponge_dens = -sigma*(dens - dynAvg%dens(iElem)) sponge_vel(1) = -sigma*(vel(1) - dynAvg%velX(iElem)) sponge_vel(2) = -sigma*(vel(2) - dynAvg%velY(iElem)) sponge_vel(3) = -sigma*(vel(3) - dynAvg%velZ(iElem)) ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} do iDir = 1, QQ ucx = dot_product( layout%fStencil%cxDirRK(:, iDir), vel ) uMinusCx = layout%fStencil%cxDirRK(:, iDir) - vel sponge_velTerm = dot_product( uMinusCx * cs2inv & & + ucx * layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, sponge_vel ) sponge_Term(iDir) = layout%weight(iDir) & & * ( sponge_dens + sponge_velTerm ) end do ! Force moments: M * F momForce(1:QQ) = matmul( layout%moment%toMoments%A(1:QQ, 1:QQ), & & sponge_Term(1:QQ) ) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = dot_product(mInvXOmega(iDir,1:QQ), momForce(1:QQ)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end associate end subroutine applySrc_absorbLayerDyn_MRT ! **************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force". !! Force term used here is from: !! "Discrete lattice effects on the forcing term in the lattice Boltzmann !! method", Zhaoli Guo, Chugung Zheng and Baochang Shi. !! In the paper, use force term is referred as Method 2 as: !! \[ F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + !! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} \] !! Force must be defined as body force per unit volume !! KM: If this force formula is used then velocity needs to be !! computed as u = \sum c_i f_i + \vec{F}/2 !! !! Simuilar to derive routine but it updates the state whereas derive !! is used for tracking. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_force( fun, inState, outState, neigh, auxField, & & nPdfSize, iLevel, varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3) real(kind=rk) :: G(3), velocity(3), ucx, uMinusCX(3), forceTerm integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3) real(kind=rk) :: omega, omega_fac ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = forceField ) ! convert physical to lattice forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of velocity variable in auxField vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) velocity(3) = auxField(elemOff + vel_pos(3)) ! force field on current element ! For incompressible model: this forceField should be divided by rho0. ! Since rho0 =1, this term is also valid for incompressible model G = forceField((iElem-1)*3+1 : iElem*3) ! get the correct omega value omega = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) omega_fac = 1.0_rk - omega * 0.5_rk ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} do iDir = 1, QQ ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), & & velocity ) uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - velocity forceTerm = dot_product( uMinusCx * cs2inv & & + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, G ) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) & & + omega_fac * scheme%layout%weight( iDir ) * forceTerm end do end do !iElem end subroutine applySrc_force ! ****************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force" for MRT collision model. !! Force term used here is from: !! Chai, Z., & Zhao, T. (2012). Effect of the forcing term in the !! multiple-relaxation-time lattice Boltzmann equation on the shear stress !! or the strain rate tensor. Physical Review E, 86(1), 1–11. !! Force term for MRT is !! \[ \bar{F} = M^-1 (I-0.5 S) M F' \] and !! \[ F'_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + !! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} \] !! !! \vec{F} is the force that must be defined as body force per unit volume !! !! Simuilar to derive routine but it updates the state whereas derive !! is used for tracking. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_force_MRT( fun, inState, outState, neigh, auxField, & & nPdfSize, iLevel, varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3) real(kind=rk) :: G(3), velocity(3), ucx, uMinusCX(3) integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3) real(kind=rk) :: omegaKine, omegaBulk, discForce real(kind=rk) :: forceTerm(27), momForce(27), s_mrt(27) real(kind=rk) :: mInvXOmega(27,27) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = forceField ) ! convert physical to lattice forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of velocity variable in auxField vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) velocity(3) = auxField(elemOff + vel_pos(3)) ! force field on current element ! For incompressible model: this forceField should be divided by rho0. ! Since rho0 =1, this term is also valid for incompressible model G = forceField((iElem-1)*3+1 : iElem*3) ! get the correct omega value omegaKine = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(1:QQ) = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=omegaKine, omegaBulk=omegaBulk, QQ=QQ) ! M^-1 * (I-0.5 S) s_mrt(1:QQ) = 1.0_rk - 0.5_rk * s_mrt(1:QQ) do iDir = 1, QQ mInvXOmega(1:QQ,iDir) = scheme%layout%moment%toPDF%A(1:QQ,iDir) & & * s_mrt(iDir) end do ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} do iDir = 1, QQ ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), & & velocity ) uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - velocity forceTerm(iDir) = scheme%layout%weight(iDir) & & * dot_product( uMinusCx * cs2inv & & + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, G ) end do ! Force moments: M * F !do iDir = 1, QQ ! momForce(iDir) = sum(scheme%layout%moment%toMoments%A(iDir,:) * forceTerm) !end do momForce(1:QQ) = matmul( scheme%layout%moment%toMoments%A(1:QQ, 1:QQ), & & forceTerm(1:QQ) ) !discForce = matmul( omegaTerm, forceTerm ) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = dot_product(mInvXOmega(iDir,1:QQ), momForce(1:QQ)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_force_MRT ! ****************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force" for MRT collision model. !! Force term used here is from: !! Chai, Z., & Zhao, T. (2012). Effect of the forcing term in the !! multiple-relaxation-time lattice Boltzmann equation on the shear stress !! or the strain rate tensor. Physical Review E, 86(1), 1–11. !! Force term for MRT is !! \[ \bar{F} = M^-1 (I-0.5 S) M F' \] and !! \[ F'_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + !! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} \] !! !! \vec{F} is the force that must be defined as body force per unit volume !! !! Simuilar to derive routine but it updates the state whereas derive !! is used for tracking. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_force_MRT_d3q27( fun, inState, outState, neigh, & & auxField, nPdfSize, iLevel, varSys, time, & & phyConvFac,derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3) real(kind=rk) :: F(3), velocity(3) integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3) real(kind=rk) :: omegaBulk, discForce real(kind=rk) :: momForce(27), s_mrt(27) real(kind=rk) :: mInvXOmega(27,27) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = forceField ) ! convert physical to lattice forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of velocity variable in auxField vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) ! MRT omegas s_mrt(1:QQ) = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=1._rk, omegaBulk=omegaBulk, QQ=QQ) ! M^-1 * (I-0.5 S) s_mrt(2:4) = 1.0_rk - 0.5_rk * s_mrt(2:4) s_mrt(10) = 1.0_rk - 0.5_rk * s_mrt(10) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) velocity(3) = auxField(elemOff + vel_pos(3)) ! force field on current element ! For incompressible model: this forceField should be divided by rho0. ! Since rho0 =1, this term is also valid for incompressible model F = forceField((iElem-1)*3+1 : iElem*3) ! MRT omegas ! overwrite omegaKine term in the element loop ! get the correct omega value s_mrt(5:9) = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! M^-1 * (I-0.5 S) s_mrt(5:9) = 1.0_rk - 0.5_rk * s_mrt(5:9) do iDir = 2, 10 mInvXOmega(1:QQ,iDir) = scheme%layout%moment%toPDF%A(1:QQ,iDir) & & * s_mrt(iDir) end do ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} ! Force moments: M * F momForce(1:QQ) = 0._rk momForce(2:4) = F(1:3) momForce(5) = F(1) * velocity(2) + F(2) * velocity(1) momForce(6) = F(2) * velocity(3) + F(3) * velocity(2) momForce(7) = F(1) * velocity(3) + F(3) * velocity(1) momForce(8) = -2._rk * ( F(2) * velocity(2) - 2._rk * F(1) * velocity(1) & & + F(3) * velocity(3) ) momForce(9) = 2._rk * ( F(2) * velocity(2) - F(3) * velocity(3) ) momForce(10) = 2._rk * ( F(1) * velocity(1) + F(2) * velocity(2) & & + F(3) * velocity(3) ) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = dot_product(mInvXOmega(iDir,2:10), momForce(2:10)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_force_MRT_d3q27 ! ****************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force" for d3q19 MRT collision model. !! Force term used here is from: !! Chai, Z., & Zhao, T. (2012). Effect of the forcing term in the !! multiple-relaxation-time lattice Boltzmann equation on the shear stress !! or the strain rate tensor. Physical Review E, 86(1), 1–11. !! Force term for MRT is !! \[ \bar{F} = M^-1 (I-0.5 S) M F' \] and !! \[ F'_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + !! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} \] !! !! Force must be defined as body force per unit volume !! !! Simuilar to derive routine but it updates the state whereas derive !! is used for tracking. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_force_MRT_d3q19( fun, inState, outState, neigh, & & auxField, nPdfSize, iLevel, varSys, & & time, phyConvFac, derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3) real(kind=rk) :: F(3), velocity(3) integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3) real(kind=rk) :: omegaKine, omegaBulk, discForce real(kind=rk) :: momForce(19), s_mrt(19) real(kind=rk) :: mInvXOmega(19,19) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = forceField ) ! convert physical to lattice forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force ! constant parameter QQ = 19 nScalars = varSys%nScalars ! Position of velocity variable in auxField vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) ! MRT omegas ! overwrite omegaKine term in the element loop ! KM: For incompressible model: omegaBulk is unused in mrtPtr s_mrt = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ) ! F = M^-1 (I-0.5 S) M F ! (I-0.5 S) - omega for force term s_mrt = 1.0_rk - 0.5_rk * s_mrt do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) velocity(3) = auxField(elemOff + vel_pos(3)) ! force field on current element ! For incompressible model: this forceField should be divided by rho0. ! Since rho0 =1, this term is also valid for incompressible model F = forceField((iElem-1)*3+1 : iElem*3) ! get the correct omega value omegaKine = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(10) = 1.0_rk - 0.5_rk * omegaKine s_mrt(12) = s_mrt(10) s_mrt(14) = s_mrt(10) s_mrt(15) = s_mrt(10) s_mrt(16) = s_mrt(10) ! M^-1 (1-0.5 S) do iDir = 1, QQ mInvXOmega(:,iDir) = scheme%layout%moment%toPDF%A(:,iDir) * s_mrt(iDir) end do ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} ! Force moments: M * F momForce(1:QQ) = 0._rk momForce(2) = 2._rk * ( F(1) * velocity(1) + F(2) * velocity(2) & & + F(3) * velocity(3) ) momForce(4) = F(1) momForce(6) = F(2) momForce(8) = F(3) momForce(10) = -2._rk * ( F(2) * velocity(2) - 2._rk * F(1) * velocity(1) & & + F(3) * velocity(3) ) momForce(12) = 2._rk * ( F(2) * velocity(2) - F(3) * velocity(3) ) momForce(14) = F(1) * velocity(2) + F(2) * velocity(1) momForce(15) = F(2) * velocity(3) + F(3) * velocity(2) momForce(16) = F(1) * velocity(3) + F(3) * velocity(1) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = sum(mInvXOmega(iDir,1:QQ) * momForce(1:QQ)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_force_MRT_d3q19 ! ****************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force" for d3q19 MRT collision model. !! Force term used here is from: !! Chai, Z., & Zhao, T. (2012). Effect of the forcing term in the !! multiple-relaxation-time lattice Boltzmann equation on the shear stress !! or the strain rate tensor. Physical Review E, 86(1), 1–11. !! Force term for MRT is !! \[ \bar{F} = M^-1 (I-0.5 S) M F' \] and !! \[ F'_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + !! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} \] !! !! Force must be defined as body force per unit volume !! !! Simuilar to derive routine but it updates the state whereas derive !! is used for tracking. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_force_MRT_d2q9( fun, inState, outState, neigh, & & auxField, nPdfSize, iLevel, varSys, & & time, phyConvFac, derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3) real(kind=rk) :: F(3), velocity(3) integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3) real(kind=rk) :: omegaKine, omegaBulk, discForce real(kind=rk) :: momForce(9), s_mrt(9) real(kind=rk) :: mInvXOmega(9,9) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = forceField ) ! convert physical to lattice forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force ! constant parameter QQ = 9 nScalars = varSys%nScalars ! Position of velocity variable in auxField vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) ! MRT omegas ! overwrite omegaKine term in the element loop ! KM: For incompressible model: omegaBulk is unused in mrtPtr s_mrt = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ) ! F = M^-1 (I-0.5 S) M F ! (I-0.5 S) - omega for force term s_mrt = 1.0_rk - 0.5_rk * s_mrt do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) ! force field on current element ! For incompressible model: this forceField should be divided by rho0. ! Since rho0 =1, this term is also valid for incompressible model F = forceField((iElem-1)*3+1 : iElem*3) ! get the correct omega value omegaKine = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(8:9) = 1.0_rk - 0.5_rk * omegaKine ! M^-1 (1-0.5 S) do iDir = 1, QQ mInvXOmega(:,iDir) = scheme%layout%moment%toPDF%A(:,iDir) * s_mrt(iDir) end do ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} ! Force moments: M * F momForce(1) = 0._rk momForce(2) = 6._rk * ( F(1) * velocity(1) + F(2) * velocity(2) ) momForce(3) = -momForce(2) momForce(4) = F(1) momForce(5) = -F(1) momForce(6) = F(2) momForce(7) = -F(2) momForce(8) = 2._rk * ( F(1) * velocity(1) - F(2) * velocity(2) ) momForce(9) = F(1) * velocity(2) + F(2) * velocity(1) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = sum(mInvXOmega(iDir,1:QQ) * momForce(1:QQ)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_force_MRT_d2q9 ! ****************************************************************************** ! ! ************************************************************************** ! !> Update state with source variable "force_1stOrd" !! Force term used here is from: !! "A D3Q27 multiple-relaxation-time lattice Boltzmann method for !! turbulent flows", K. Suga, Y. Kuwata, K. Takashima, R. Chikasue !! !! \[ F_i = w_i/c_s^2 ( \vec{e}_i \cdot \vec{F} ) \] !! Force must be defined as body force per unit volume !! This force term can be applied for both compressible and incompressible !! LBM models !! !! Similar to derive routine but it updates the state whereas derive !! is used for tracking !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_force1stOrd( fun, inState, outState, neigh, auxField, & & nPdfSize, iLevel, varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3) real(kind=rk) :: FF_elem(3) real(kind=rk) :: forceTerm integer :: nElems, iElem, iDir, QQ, nScalars integer :: posInTotal ! ---------------------------------------------------------------------- ! ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Get force which is refered in config file either its ! spacetime variable or operation variable call varSys%method%val(fun%data_varPos)%get_valOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fun%elemLvl(iLevel)%idx(1:nElems), & & nVals = nElems, & & res = forceField ) !write(dbgUnit(1),*) 'ApplySrc_force1stOrdIncomp' ! do iElem = 1, nElems ! posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! FF_elem = forceField((iElem-1)*3+1 : iElem*3) !write(dbgUnit(1),*) 'treeID ', scheme%levelDesc(iLevel)%total(posInTotal) & ! & , 'force ', FF_elem ! end do ! convert physical to lattice forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars do iElem = 1, nElems ! to access level wise state array posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! force field on current element FF_elem = forceField((iElem-1)*3+1 : iElem*3) ! force term: ! F_i = w_i/cs2 (\vec{e}_i \cdot \vec{F} do iDir = 1, QQ forceTerm = dot_product( scheme%layout%fStencil & & %cxDirRK(:,iDir), FF_elem ) outState( (posintotal-1)*nscalars+idir+(1-1)*qq ) & & = outState( & & (posintotal-1)*nscalars+idir+(1-1)*qq ) & & + scheme%layout%weight( iDir ) * cs2inv * forceTerm end do end do !iElem end subroutine applySrc_force1stOrd ! ************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "turb_channel_force" for BGK. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_turbChanForce( fun, inState, outState, neigh, auxField, & & nPdfSize, iLevel, varSys, time, & & phyConvFac, derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: velocity(3), ucx, uMinusCX(3), forceTerm, dens integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: omega, omega_fac real(kind=rk) :: forceDynL(3) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Convert Dynamic force in m/s^2 to lattice forceDynL = fun%turbChanForce%forceDyn & & / fPtr%solverData%physics%fac(iLevel)%accel ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of density and velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain density from auxField dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) velocity(3) = auxField(elemOff + vel_pos(3)) ! get the correct omega value omega = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) omega_fac = 1.0_rk - omega * 0.5_rk ! force term: ! F_i = w_i*(1-omega/2)*rho*( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{g} do iDir = 1, QQ ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), & & velocity ) uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - velocity forceTerm = dot_product( uMinusCx * cs2inv & & + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, forceDynL ) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) & & + omega_fac * scheme%layout%weight( iDir ) * dens * forceTerm end do end do !iElem end subroutine applySrc_turbChanForce ! ****************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force" for generic MRT collision model !! for turb_channel_force. It uses velocityX average in bulk to adapt the !! driving force for turbulent channel. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_turbChanForce_MRT( fun, inState, outState, neigh, & & auxField, nPdfSize, iLevel, & & varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: velocity(3), ucx, uMinusCX(3), dens integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: omegaKine, omegaBulk, discForce real(kind=rk) :: forceTerm(27), momForce(27), s_mrt(27) real(kind=rk) :: mInvXOmega(27,27) real(kind=rk) :: forceDynL(3) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! convert physical to lattice, force is defined in acceleration m/s^2 ! Convert Dynamic force in m/s^2 to lattice forceDynL = fun%turbChanForce%forceDyn & & / fPtr%solverData%physics%fac(iLevel)%accel ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of density and velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain density from auxField dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) velocity(3) = auxField(elemOff + vel_pos(3)) ! get the correct omega value omegaKine = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(1:QQ) = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=omegaKine, omegaBulk=omegaBulk, QQ=QQ) ! M^-1 * (I-0.5 S) s_mrt(1:QQ) = 1.0_rk - 0.5_rk * s_mrt(1:QQ) do iDir = 1, QQ mInvXOmega(1:QQ,iDir) = scheme%layout%moment%toPDF%A(1:QQ,iDir) & & * s_mrt(iDir) end do ! force term: ! F_i = w_i*(1-omega/2)*rho*( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} do iDir = 1, QQ ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), & & velocity ) uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - velocity forceTerm(iDir) = scheme%layout%weight(iDir) * dens & & * dot_product( uMinusCx * cs2inv & & + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) & & * cs4inv, forceDynL ) end do ! Force moments: M * F !do iDir = 1, QQ ! momForce(iDir) = sum(scheme%layout%moment%toMoments%A(iDir,:) * forceTerm) !end do momForce(1:QQ) = matmul(scheme%layout%moment%toMoments%A(1:QQ,1:QQ), & & forceTerm(1:QQ)) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = sum(mInvXOmega(iDir,1:QQ) * momForce(1:QQ)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_turbChanForce_MRT ! ****************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force" for generic MRT collision model !! for turb_channel_force. It uses velocityX average in bulk to adapt the !! driving force for turbulent channel. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_turbChanForce_MRT_d3q27( fun, inState, outState, neigh, & & auxField, nPdfSize, iLevel, & & varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: F(3), velocity(3), dens integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: omegaBulk, discForce real(kind=rk) :: momForce(27), s_mrt(27) real(kind=rk) :: mInvXOmega(27,27) real(kind=rk) :: forceDynL(3) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! convert physical to lattice, force is defined in acceleration m/s^2 ! Convert Dynamic force in m/s^2 to lattice forceDynL = fun%turbChanForce%forceDyn & & / fPtr%solverData%physics%fac(iLevel)%accel ! constant parameter QQ = scheme%layout%fStencil%QQ nScalars = varSys%nScalars ! Position of density and velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) ! MRT omegas s_mrt(1:QQ) = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=1._rk, omegaBulk=omegaBulk, QQ=QQ) ! M^-1 * (I-0.5 S) s_mrt(2:4) = 1.0_rk - 0.5_rk * s_mrt(2:4) s_mrt(10) = 1.0_rk - 0.5_rk * s_mrt(10) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain density from auxField dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) velocity(3) = auxField(elemOff + vel_pos(3)) ! MRT omegas ! overwrite omegaKine term in the element loop ! get the correct omega value s_mrt(5:9) = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! M^-1 * (I-0.5 S) s_mrt(5:9) = 1.0_rk - 0.5_rk * s_mrt(5:9) do iDir = 2, 10 mInvXOmega(1:QQ,iDir) = scheme%layout%moment%toPDF%A(1:QQ,iDir) & & * s_mrt(iDir) end do ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} ! Force moments: M * F F = dens * forceDynL momForce(1:QQ) = 0._rk momForce(2:4) = F(1:3) momForce(5) = F(1) * velocity(2) + F(2) * velocity(1) momForce(6) = F(2) * velocity(3) + F(3) * velocity(2) momForce(7) = F(1) * velocity(3) + F(3) * velocity(1) momForce(8) = -2._rk * ( F(2) * velocity(2) - 2._rk * F(1) * velocity(1) & & + F(3) * velocity(3) ) momForce(9) = 2._rk * ( F(2) * velocity(2) - F(3) * velocity(3) ) momForce(10) = 2._rk * ( F(1) * velocity(1) + F(2) * velocity(2) & & + F(3) * velocity(3) ) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = sum(mInvXOmega(iDir,2:10) * momForce(2:10)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_turbChanForce_MRT_d3q27 ! ****************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force" for d3q19 MRT collision model !! for turb_channel_force. It uses velocityX average in bulk to adapt the !! driving force for turbulent channel. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_turbChanForce_MRT_d3q19( fun, inState, outState, neigh, & & auxField, nPdfSize, iLevel, & & varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: F(3), velocity(3), dens integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: omegaKine, omegaBulk, discForce real(kind=rk) :: momForce(19), s_mrt(19) real(kind=rk) :: mInvXOmega(19,19) real(kind=rk) :: forceDynL(3) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! convert physical to lattice, force is defined in acceleration m/s^2 ! Convert Dynamic force in m/s^2 to lattice forceDynL = fun%turbChanForce%forceDyn & & / fPtr%solverData%physics%fac(iLevel)%accel ! constant parameter QQ = 19 nScalars = varSys%nScalars ! Position of density and velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) ! MRT omegas ! overwrite omegaKine term in the element loop ! KM: For incompressible model: omegaBulk is unused in mrtPtr s_mrt = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ) ! F = M^-1 (I-0.5 S) M F ! (I-0.5 S) - omega for force term s_mrt = 1.0_rk - 0.5_rk * s_mrt do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain density from auxField dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) velocity(3) = auxField(elemOff + vel_pos(3)) ! get the correct omega value omegaKine = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(10) = 1.0_rk - 0.5_rk * omegaKine s_mrt(12) = s_mrt(10) s_mrt(14) = s_mrt(10) s_mrt(15) = s_mrt(10) s_mrt(16) = s_mrt(10) ! M^-1 (1-0.5 S) do iDir = 1, QQ mInvXOmega(:,iDir) = scheme%layout%moment%toPDF%A(:,iDir) * s_mrt(iDir) end do ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} ! Force moments: M * F F = dens * forceDynL momForce(1:QQ) = 0._rk momForce(2) = 2._rk * ( F(1) * velocity(1) + F(2) * velocity(2) & & + F(3) * velocity(3) ) momForce(4) = F(1) momForce(6) = F(2) momForce(8) = F(3) momForce(10) = -2._rk * ( F(2) * velocity(2) - 2._rk * F(1) * velocity(1) & & + F(3) * velocity(3) ) momForce(12) = 2._rk * ( F(2) * velocity(2) - F(3) * velocity(3) ) momForce(14) = F(1) * velocity(2) + F(2) * velocity(1) momForce(15) = F(2) * velocity(3) + F(3) * velocity(2) momForce(16) = F(1) * velocity(3) + F(3) * velocity(1) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = sum(mInvXOmega(iDir,1:QQ) * momForce(1:QQ)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_turbChanForce_MRT_d3q19 ! ****************************************************************************** ! ! ****************************************************************************** ! !> Update state with source variable "force" for d3q19 MRT collision model !! for turb_channel_force. It uses velocityX average in bulk to adapt the !! driving force for turbulent channel. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_turbChanForce_MRT_d2q9( fun, inState, outState, neigh, & & auxField, nPdfSize, iLevel, & & varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk) :: F(3), velocity(3), dens integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff integer :: vel_pos(3), dens_pos real(kind=rk) :: omegaKine, omegaBulk, discForce real(kind=rk) :: momForce(9), s_mrt(9) real(kind=rk) :: mInvXOmega(9,9) real(kind=rk) :: forceDynL(3) ! --------------------------------------------------------------------------- ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) scheme => fPtr%solverData%scheme ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! convert physical to lattice, force is defined in acceleration m/s^2 ! Convert Dynamic force in m/s^2 to lattice forceDynL = fun%turbChanForce%forceDyn & & / fPtr%solverData%physics%fac(iLevel)%accel ! constant parameter QQ = 9 nScalars = varSys%nScalars ! Position of density and velocity variable in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel) ! MRT omegas ! overwrite omegaKine term in the element loop ! KM: For incompressible model: omegaBulk is unused in mrtPtr s_mrt = scheme%field(1)%fieldProp%fluid & & %mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ) ! F = M^-1 (I-0.5 S) M F ! (I-0.5 S) - omega for force term s_mrt = 1.0_rk - 0.5_rk * s_mrt do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! element offset elemoff = (posInTotal-1)*varSys%nAuxScalars ! obtain density from auxField dens = auxField(elemOff + dens_pos) ! obtain velocity from auxField velocity(1) = auxField(elemOff + vel_pos(1)) velocity(2) = auxField(elemOff + vel_pos(2)) ! get the correct omega value omegaKine = scheme%field(1)%fieldProp%fluid%viscKine & & %omLvl(iLevel)%val(posInTotal) ! MRT omegas ! overwrite omegaKine term in the element loop s_mrt(8:9) = 1.0_rk - 0.5_rk * omegaKine ! M^-1 (1-0.5 S) do iDir = 1, QQ mInvXOmega(:,iDir) = scheme%layout%moment%toPDF%A(:,iDir) * s_mrt(iDir) end do ! force term: ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 + ! (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} ! Force moments: M * F F = dens * forceDynL momForce(1) = 0._rk momForce(2) = 6._rk * ( F(1) * velocity(1) + F(2) * velocity(2) ) momForce(3) = -momForce(2) momForce(4) = F(1) momForce(5) = -F(1) momForce(6) = F(2) momForce(7) = -F(2) momForce(8) = 2._rk * ( F(1) * velocity(1) - F(2) * velocity(2) ) momForce(9) = F(1) * velocity(2) + F(2) * velocity(1) do iDir = 1, QQ ! discrete force ! \bar{F} = M^-1 (I-S/2) M F discForce = sum(mInvXOmega(iDir,1:QQ) * momForce(1:QQ)) ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + discForce end do end do !iElem end subroutine applySrc_turbChanForce_MRT_d2q9 ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate wss from shear stress (tau) !! tau: x, y, z, xy, yz, xz pure function getWSS( tau ) result( wss ) ! --------------------------------------------------------------------------- real(kind=rk), intent(in) :: tau(6) real(kind=rk) :: a0, a1, a2 real(kind=rk) :: q_stress, q_stress_cube, q_term, r_stress, d_stress real(kind=rk) :: theta_stress, cos_term, shear_1, shear_2, shear_3 real(kind=rk) :: wss ! --------------------------------------------------------------------------- ! tau_x + tau_y + tau_z a2 = -1._rk * ( tau(1) + tau(2) + tau(3) ) ! tau_x * tau_y + tau_y * tau_z + tau_x * tau_z - tau_xy^2 - tau_xz^2 - tau_yz^2 a1 = tau(1) * tau(2) + tau(2) * tau(3) + tau(1) * tau(3) & & - tau(4) * tau(4) - tau(5) * tau(5) - tau(6) * tau(6) a0 = -1._rk * ( tau(1) * tau(2) * tau(3) & & + 2._rk * tau(4) * tau(5) * tau(6) & & - tau(1) * tau(6) * tau(6) & & - tau(2) * tau(5) * tau(5) & & - tau(3) * tau(4) * tau(4) ) q_stress = div1_9 * ( 3._rk * a1 - a2 * a2 ) q_stress_cube = q_stress * q_stress * q_stress if ( q_stress < 0._rk ) then q_term = 2._rk * sqrt( -q_stress ) else q_term = 0._rk end if r_stress = div1_54 * ( 9._rk*a2*a1 - 27._rk*a0 - 2._rk*a2*a2*a2 ) d_stress = q_stress_cube + r_stress * r_stress if ( d_stress < 0._rk ) then ! the solutions are real and unequal theta_stress = div1_3 * acos( r_stress / sqrt( -q_stress_cube )) cos_term = div3_4 * q_term * cos( theta_stress ) ! three values of maximum shear stress shear_1 = abs( sqrt3 * div1_2 * q_term * sin( theta_stress )) shear_2 = abs( cos_term - shear_1*div1_2 ) shear_3 = abs( cos_term + shear_1*div1_2 ) ! choose the max shear as Wall Shear Stress(wss) wss = max(shear_1, shear_2, shear_3) else wss = 0._rk end if end function getWSS ! ****************************************************************************** ! ! ****************************************************************************** ! !> This routine computes equilbrium from density and velocity !! This must comply with interface in mus_variable_module derive_FromMacro !! !! This subroutine's interface must match the abstract interface definition !! [[derive_FromMacro]] in derived/[[mus_derVarPos_module]].f90 in order to be !! callable via [[mus_derVarPos_type:equilFromMacro]] function pointer. subroutine deriveEquil_FromMacro( density, velocity, iField, nElems, varSys, & & layout, res ) ! -------------------------------------------------------------------- ! !> Array of density. !! Single species: dens_1, dens_2 .. dens_n !! multi-species: dens_1_sp1, dens_1_sp2, dens_2_sp1, dens_2_sp2 ... !! dens_n_sp1, dens_n_sp2 real(kind=rk), intent(in) :: density(:) !> Array of velocity. !! Size: dimension 1: n*nFields. dimension 2: 3 (nComp) !! 1st dimension arrangement for multi-species is same as density real(kind=rk), intent(in) :: velocity(:, :) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> scheme layout contains stencil definition and lattice weights type(mus_scheme_layout_type), intent(in) :: layout !> Output of this routine !! Dimension: n*nComponents of res real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! real(kind=rk) :: fEq(layout%fStencil%QQ), vel(3), usq, ucx integer :: QQ, iElem, iDir ! --------------------------------------------------------------------------- QQ = layout%fStencil%QQ do iElem = 1, nElems vel = velocity(:, iElem) usq = dot_product(vel, vel)*t2cs2inv do iDir = 1, QQ ucx = dot_product( layout%fStencil%cxDirRK(:, iDir), vel ) ! calculate equilibrium density fEq( iDir ) = layout%weight( iDir ) * density(iElem) & & * ( 1._rk + cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) enddo res( (iElem-1)*QQ+1: iElem*QQ ) = fEq end do end subroutine deriveEquil_FromMacro ! ****************************************************************************** ! ! ************************************************************************** ! !> This routine computes equilbrium from auxField !! !! This subroutine's interface must match the abstract interface definition !! [[derive_equilFromAux]] in derived/[[mus_derVarPos_module]].f90 in order to !! be callable via [[mus_derVarPos_type:equilFromAux]] function pointer. subroutine deriveEquil_fromAux( derVarPos, auxField, iField, nElems, varSys, & & layout, fEq ) ! -------------------------------------------------------------------- ! !> Position of derive variable in variable system class(mus_derVarPos_type), intent(in) :: derVarPos !> Array of auxField. !! Single species: dens_1, vel_1, dens_2, vel_2, .. dens_n, vel_n !! multi-species: dens_1_sp1, vel_1_spc1, dens_1_sp2, vel_1_spc2, !! dens_2_sp1, vel_2_spc2, dens_2_sp2, vel_2_spc2 ... !! dens_n_sp1, vel_n_sp1, dens_n_sp2, vel_n_spc2 !! Access: (iElem-1)*nAuxScalars + auxField_varPos real(kind=rk), intent(in) :: auxField(:) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> scheme layout contains stencil definition and lattice weights type(mus_scheme_layout_type), intent(in) :: layout !> Output of this routine !! Dimension: n*QQ of res real(kind=rk), intent(out) :: fEq(:) ! -------------------------------------------------------------------- ! real(kind=rk) :: rho, vel(3), usq, ucx integer :: QQ, iElem, iDir, elemOff integer :: dens_pos, vel_pos(3) ! ---------------------------------------------------------------------- ! dens_pos = varSys%method%val(derVarPos%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3) QQ = layout%fStencil%QQ !NEC$ ivdep do iElem = 1, nElems ! element offset elemoff = (iElem-1)*varSys%nAuxScalars ! density rho = auxField(elemOff + dens_pos) ! velocity vel(1) = auxField(elemOff + vel_pos(1)) vel(2) = auxField(elemOff + vel_pos(2)) vel(3) = auxField(elemOff + vel_pos(3)) usq = ( vel(1)*vel(1) + vel(2)*vel(2) + vel(3)*vel(3) )*t2cs2inv !NEC$ shortloop do iDir = 1, QQ ucx = layout%fStencil%cxDirRK(1, iDir) * vel(1) & & + layout%fStencil%cxDirRK(2, iDir) * vel(2) & & + layout%fStencil%cxDirRK(3, iDir) * vel(3) ! calculate equilibrium density fEq( (iElem-1)*QQ + iDir ) = layout%weight( iDir ) * rho * ( 1.0_rk & & + cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) enddo end do end subroutine deriveEquil_fromAux ! ************************************************************************** ! ! ****************************************************************************** ! !> This routine computes density from state array !! !! This subroutine's interface must match the abstract interface definition !! [[derive_FromState]] in derived/[[mus_derVarPos_module]].f90 in order to be !! callable via [[mus_derVarPos_type:velFromState]], !! mus_derVarPos_type:equilFromState, mus_derVarPos_type:momFromState, !! mus_derVarPos_type:velocitiesFromState, and !! mus_derVarPos_type:momentaFromState function pointers. subroutine deriveRho_FromState( state, iField, nElems, varSys, layout, res ) ! -------------------------------------------------------------------- ! !> Array of state !! n * layout%fStencil%QQ * nFields real(kind=rk), intent(in) :: state(:) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> scheme layout contains stencil definition and lattice weights type(mus_scheme_layout_type), intent(in) :: layout !> Output of this routine !! Dimension: n * nComponents of res real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! integer :: iElem, iDir ! --------------------------------------------------------------------------- res = 0.0_rk do iElem = 1, nElems do iDir = 1, layout%fStencil%QQ res( iElem ) = res( iElem ) + state( iDir+(iElem-1)*varSys%nScalars ) end do end do end subroutine deriveRho_FromState ! ****************************************************************************** ! ! ****************************************************************************** ! !> This routine computes velocity from state array !! !! This subroutine's interface must match the abstract interface definition !! [[derive_FromState]] in derived/[[mus_derVarPos_module]].f90 in order to be !! callable via [[mus_derVarPos_type:velFromState]], !! [[mus_derVarPos_type:equilFromState]], [[mus_derVarPos_type:momFromState], !! [[mus_derVarPos_type:velocitiesFromState]], and !! [[mus_derVarPos_type:momentaFromState] function pointers. subroutine deriveVel_FromState( state, iField, nElems, varSys, layout, res ) ! -------------------------------------------------------------------- ! !> Array of state !! n * layout%fStencil%QQ * nFields real(kind=rk), intent(in) :: state(:) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> scheme layout contains stencil definition and lattice weights type(mus_scheme_layout_type), intent(in) :: layout !> Output of this routine !! Dimension: n * nComponents of res real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! integer :: iElem, iDir real(kind=rk) :: rho, pdf( layout%fStencil%QQ ) ! --------------------------------------------------------------------------- do iElem = 1, nElems do iDir = 1, layout%fStencil%QQ pdf(iDir) = state( iDir+(iElem-1)*varSys%nScalars ) end do rho = sum( pdf ) res( (iElem-1)*3+1 ) = sum( pdf * layout%fStencil%cxDirRK(1,:) ) / rho res( (iElem-1)*3+2 ) = sum( pdf * layout%fStencil%cxDirRK(2,:) ) / rho res( (iElem-1)*3+3 ) = sum( pdf * layout%fStencil%cxDirRK(3,:) ) / rho end do end subroutine deriveVel_FromState ! ****************************************************************************** ! ! ****************************************************************************** ! !> This routine computes velocity from pre collision state array using Fetch !! !! This subroutine's interface must match the abstract interface definition !! [[derive_FromPreColState]] in derived/[[mus_derVarPos_module]].f90 in order !! to be callable via [[mus_derVarPos_type:velFromPreColState]] function !! pointer. subroutine deriveVel_FromPreColState( state, neigh, iField, nSize, nElems, & varSys, layout, res ) ! -------------------------------------------------------------------- ! !> Array of state !! n * layout%fStencil%QQ * nFields real(kind=rk), intent(in) :: state(:) !> connectivity array integer, intent(in) :: neigh(:) !> Current field integer, intent(in) :: iField !> number of elements in state array integer, intent(in) :: nSize !> number of elements integer, intent(in) :: nElems !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> scheme layout contains stencil definition and lattice weights type(mus_scheme_layout_type), intent(in) :: layout !> Output of this routine !! Dimension: n * nComponents of res real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! integer :: iElem, iDir real(kind=rk) :: rho, pdf( layout%fStencil%QQ ), vel(3) integer :: QQ, nScalars, nDims ! --------------------------------------------------------------------------- QQ = layout%fStencil%QQ nScalars = varSys%nScalars nDims = layout%fStencil%nDims do iElem = 1, nElems do iDir = 1, QQ pdf(iDir) = state( & & neigh((idir-1)* nsize+ ielem)+( ifield-1)* qq+ nscalars*0) end do ! density rho = sum( pdf ) ! momentum vel( 1 ) = sum( pdf * layout%fStencil%cxDirRK(1,:) ) vel( 2 ) = sum( pdf * layout%fStencil%cxDirRK(2,:) ) vel( 3 ) = sum( pdf * layout%fStencil%cxDirRK(3,:) ) ! return velocity field according on stencil dimensions res( (iElem-1)*nDims+1:iElem*nDims) = vel(1:nDims) / rho end do end subroutine deriveVel_FromPreColState ! ****************************************************************************** ! ! **************************************************************************** ! !> This routine computes auxField from state array !! !! This subroutine's interface must match the abstract interface definition !! [[derive_auxFromState]] in derived/[[mus_derVarPos_module]].f90 in order to !! be callable via [[mus_derVarPos_type:auxFieldFromState]] function pointer. subroutine deriveAux_fromState( derVarPos, state, neigh, iField, nElems, & & nSize, iLevel, stencil, varSys, auxField ) ! -------------------------------------------------------------------- ! !> Position of derive variable in variable system class(mus_derVarPos_type), intent(in) :: derVarPos !> Array of state !! n * layout%stencil(1)%QQ * nFields real(kind=rk), intent(in) :: state(:) !> connectivity vector integer, intent(in) :: neigh(:) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> number of elements in state array integer, intent(in) :: nSize !> current level integer, intent(in) :: iLevel !> stencil header contains discrete velocity vectors type(tem_stencilHeader_type), intent(in) :: stencil !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> Output of this routine !! Size: nElems*nAuxScalars real(kind=rk), intent(inout) :: auxField(:) ! -------------------------------------------------------------------- ! integer :: dens_pos, vel_pos(3) integer :: iElem, iDir, elemOff real(kind=rk) :: pdf( stencil%QQ ), rho, inv_rho ! ------------------------------------------------------------------------ ! dens_pos = varSys%method%val(derVarPos%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3) !NEC$ ivdep do iElem = 1, nElems !NEC$ shortloop do iDir = 1, stencil%QQ pdf(iDir) = state( & & ( ielem-1)* varsys%nscalars+idir+( 1-1)* stencil%qq) end do ! element offset elemoff = (iElem-1)*varSys%nAuxScalars ! density rho = sum( pdf ) auxField(elemOff+dens_pos) = rho ! velocity inv_rho = 1.0_rk/rho auxField(elemOff+vel_pos(1)) = sum( pdf * stencil%cxDirRK(1,:) ) & & * inv_rho auxField(elemOff+vel_pos(2)) = sum( pdf * stencil%cxDirRK(2,:) ) & & * inv_rho auxField(elemOff+vel_pos(3)) = sum( pdf * stencil%cxDirRK(3,:) ) & & * inv_rho end do end subroutine deriveAux_fromState ! **************************************************************************** ! ! **************************************************************************** ! !> This routine computes auxField from state array !! !! This subroutine's interface must match the abstract interface definition !! [[derive_auxFromState]] in derived/[[mus_derVarPos_module]].f90 in order to !! be callable via [[mus_derVarPos_type:auxFieldFromState]] function pointer. subroutine deriveAux_fromState_d2q9( derVarPos, state, neigh, iField, & & nElems, nSize, iLevel, stencil, & & varSys, auxField ) ! -------------------------------------------------------------------- ! !> Position of derive variable in variable system class(mus_derVarPos_type), intent(in) :: derVarPos !> Array of state !! n * layout%stencil(1)%QQ * nFields real(kind=rk), intent(in) :: state(:) !> connectivity vector integer, intent(in) :: neigh(:) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> number of elements in state array integer, intent(in) :: nSize !> current level integer, intent(in) :: iLevel !> stencil header contains discrete velocity vectors type(tem_stencilHeader_type), intent(in) :: stencil !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> Output of this routine !! Size: nElems*nAuxScalars real(kind=rk), intent(inout) :: auxField(:) ! -------------------------------------------------------------------- ! integer :: dens_pos, vel_pos(3) integer :: iElem, QQ, nScalars, elemOff real(kind=rk) :: pdf(9), rho, inv_rho, u_x, u_y ! ------------------------------------------------------------------------ ! dens_pos = varSys%method%val(derVarPos%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3) QQ = 9 nScalars = varSys%nScalars !NEC$ ivdep do iElem = 1, nElems pdf(1) = state(( ielem-1)* nscalars+1+( 1-1)* qq) pdf(2) = state(( ielem-1)* nscalars+2+( 1-1)* qq) pdf(3) = state(( ielem-1)* nscalars+3+( 1-1)* qq) pdf(4) = state(( ielem-1)* nscalars+4+( 1-1)* qq) pdf(5) = state(( ielem-1)* nscalars+5+( 1-1)* qq) pdf(6) = state(( ielem-1)* nscalars+6+( 1-1)* qq) pdf(7) = state(( ielem-1)* nscalars+7+( 1-1)* qq) pdf(8) = state(( ielem-1)* nscalars+8+( 1-1)* qq) pdf(9) = state(( ielem-1)* nscalars+9+( 1-1)* qq) ! element offset elemoff = (iElem-1)*varSys%nAuxScalars ! density rho = sum( pdf ) auxField(elemOff+dens_pos) = rho ! velocity inv_rho = 1.0_rk/rho u_x = pdf(3) + pdf(7) + pdf(8) - pdf(1) - pdf(5) - pdf(6) u_y = pdf(4) + pdf(6) + pdf(8) - pdf(2) - pdf(5) - pdf(7) u_x = u_x * inv_rho u_y = u_y * inv_rho auxField(elemOff+vel_pos(1)) = u_x auxField(elemOff+vel_pos(2)) = u_y auxField(elemOff+vel_pos(3)) = 0.0_rk end do end subroutine deriveAux_fromState_d2q9 ! **************************************************************************** ! ! **************************************************************************** ! !> This routine computes auxField from state array !! !! This subroutine's interface must match the abstract interface definition !! [[derive_auxFromState]] in derived/[[mus_derVarPos_module]].f90 in order to !! be callable via [[mus_derVarPos_type:auxFieldFromState]] function pointer. subroutine deriveAux_fromState_d3q19( derVarPos, state, neigh, iField, & & nElems, nSize, iLevel, stencil, & & varSys, auxField ) ! -------------------------------------------------------------------- ! !> Position of derive variable in variable system class(mus_derVarPos_type), intent(in) :: derVarPos !> Array of state !! n * layout%stencil(1)%QQ * nFields real(kind=rk), intent(in) :: state(:) !> connectivity vector integer, intent(in) :: neigh(:) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> number of elements in state array integer, intent(in) :: nSize !> current level integer, intent(in) :: iLevel !> stencil header contains discrete velocity vectors type(tem_stencilHeader_type), intent(in) :: stencil !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> Output of this routine !! Size: nElems*nAuxScalars real(kind=rk), intent(inout) :: auxField(:) ! -------------------------------------------------------------------- ! integer :: dens_pos, vel_pos(3) integer :: iElem, QQ, nScalars, elemOff real(kind=rk) :: pdf(19), rho, inv_rho, u_x, u_y, u_z ! ------------------------------------------------------------------------ ! dens_pos = varSys%method%val(derVarPos%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3) QQ = 19 nScalars = varSys%nScalars !NEC$ ivdep do iElem = 1, nElems pdf( 1) = state(( ielem-1)* nscalars+ 1+( 1-1)* qq) pdf( 2) = state(( ielem-1)* nscalars+ 2+( 1-1)* qq) pdf( 3) = state(( ielem-1)* nscalars+ 3+( 1-1)* qq) pdf( 4) = state(( ielem-1)* nscalars+ 4+( 1-1)* qq) pdf( 5) = state(( ielem-1)* nscalars+ 5+( 1-1)* qq) pdf( 6) = state(( ielem-1)* nscalars+ 6+( 1-1)* qq) pdf( 7) = state(( ielem-1)* nscalars+ 7+( 1-1)* qq) pdf( 8) = state(( ielem-1)* nscalars+ 8+( 1-1)* qq) pdf( 9) = state(( ielem-1)* nscalars+ 9+( 1-1)* qq) pdf(10) = state(( ielem-1)* nscalars+10+( 1-1)* qq) pdf(11) = state(( ielem-1)* nscalars+11+( 1-1)* qq) pdf(12) = state(( ielem-1)* nscalars+12+( 1-1)* qq) pdf(13) = state(( ielem-1)* nscalars+13+( 1-1)* qq) pdf(14) = state(( ielem-1)* nscalars+14+( 1-1)* qq) pdf(15) = state(( ielem-1)* nscalars+15+( 1-1)* qq) pdf(16) = state(( ielem-1)* nscalars+16+( 1-1)* qq) pdf(17) = state(( ielem-1)* nscalars+17+( 1-1)* qq) pdf(18) = state(( ielem-1)* nscalars+18+( 1-1)* qq) pdf(19) = state(( ielem-1)* nscalars+19+( 1-1)* qq) ! element offset elemoff = (iElem-1)*varSys%nAuxScalars ! density rho = sum( pdf ) auxField(elemOff+dens_pos) = rho ! velocity inv_rho = 1.0_rk/rho u_x = (pdf(4) - pdf(1)) & & + (pdf(12) - pdf(11)) & & + (pdf(14) - pdf(13)) & & + (pdf(17) - pdf(15)) & & + (pdf(18) - pdf(16)) u_y = (pdf(5) - pdf(2)) & & + (pdf(9) - pdf(7)) & & + (pdf(10) - pdf(8)) & & + (pdf(16) - pdf(15)) & & + (pdf(18) - pdf(17)) u_z = (pdf(6) - pdf(3)) & & + (pdf(8) - pdf(7)) & & + (pdf(10) - pdf(9)) & & + (pdf(13) - pdf(11)) & & + (pdf(14) - pdf(12)) u_x = u_x * inv_rho u_y = u_y * inv_rho u_z = u_z * inv_rho auxField(elemOff+vel_pos(1)) = u_x auxField(elemOff+vel_pos(2)) = u_y auxField(elemOff+vel_pos(3)) = u_z end do end subroutine deriveAux_fromState_d3q19 ! **************************************************************************** ! ! **************************************************************************** ! !> This routine computes auxField from state array !! !! This subroutine's interface must match the abstract interface definition !! [[derive_auxFromState]] in derived/[[mus_derVarPos_module]].f90 in order to !! be callable via [[mus_derVarPos_type:auxFieldFromState]] function pointer. subroutine deriveAux_fromState_d3q27( derVarPos, state, neigh, iField, & & nElems, nSize, iLevel, stencil, & & varSys, auxField ) ! -------------------------------------------------------------------- ! !> Position of derive variable in variable system class(mus_derVarPos_type), intent(in) :: derVarPos !> Array of state !! n * layout%stencil(1)%QQ * nFields real(kind=rk), intent(in) :: state(:) !> connectivity vector integer, intent(in) :: neigh(:) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> number of elements in state array integer, intent(in) :: nSize !> current level integer, intent(in) :: iLevel !> stencil header contains discrete velocity vectors type(tem_stencilHeader_type), intent(in) :: stencil !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> Output of this routine !! Size: nElems*nAuxScalars real(kind=rk), intent(inout) :: auxField(:) ! -------------------------------------------------------------------- ! integer :: dens_pos, vel_pos(3) integer :: iElem, QQ, nScalars, elemOff real(kind=rk) :: pdf(27), rho, inv_rho, u_x, u_y, u_z ! ------------------------------------------------------------------------ ! dens_pos = varSys%method%val(derVarPos%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3) QQ = 27 nScalars = varSys%nScalars !NEC$ ivdep do iElem = 1, nElems pdf(1) = state(( ielem-1)* nscalars+1+( 1-1)* qq) pdf(2) = state(( ielem-1)* nscalars+2+( 1-1)* qq) pdf(3) = state(( ielem-1)* nscalars+3+( 1-1)* qq) pdf(4) = state(( ielem-1)* nscalars+4+( 1-1)* qq) pdf(5) = state(( ielem-1)* nscalars+5+( 1-1)* qq) pdf(6) = state(( ielem-1)* nscalars+6+( 1-1)* qq) pdf(7) = state(( ielem-1)* nscalars+7+( 1-1)* qq) pdf(8) = state(( ielem-1)* nscalars+8+( 1-1)* qq) pdf(9) = state(( ielem-1)* nscalars+9+( 1-1)* qq) pdf(10) = state(( ielem-1)* nscalars+10+( 1-1)* qq) pdf(11) = state(( ielem-1)* nscalars+11+( 1-1)* qq) pdf(12) = state(( ielem-1)* nscalars+12+( 1-1)* qq) pdf(13) = state(( ielem-1)* nscalars+13+( 1-1)* qq) pdf(14) = state(( ielem-1)* nscalars+14+( 1-1)* qq) pdf(15) = state(( ielem-1)* nscalars+15+( 1-1)* qq) pdf(16) = state(( ielem-1)* nscalars+16+( 1-1)* qq) pdf(17) = state(( ielem-1)* nscalars+17+( 1-1)* qq) pdf(18) = state(( ielem-1)* nscalars+18+( 1-1)* qq) pdf(19) = state(( ielem-1)* nscalars+19+( 1-1)* qq) pdf(20) = state(( ielem-1)* nscalars+20+( 1-1)* qq) pdf(21) = state(( ielem-1)* nscalars+21+( 1-1)* qq) pdf(22) = state(( ielem-1)* nscalars+22+( 1-1)* qq) pdf(23) = state(( ielem-1)* nscalars+23+( 1-1)* qq) pdf(24) = state(( ielem-1)* nscalars+24+( 1-1)* qq) pdf(25) = state(( ielem-1)* nscalars+25+( 1-1)* qq) pdf(26) = state(( ielem-1)* nscalars+26+( 1-1)* qq) pdf(27) = state(( ielem-1)* nscalars+27+( 1-1)* qq) ! element offset elemoff = (iElem-1)*varSys%nAuxScalars ! density rho = sum( pdf ) auxField(elemOff+dens_pos) = rho ! velocity inv_rho = 1.0_rk/rho u_x = (pdf(4) - pdf(1)) & & + (pdf(12) - pdf(11)) & & + (pdf(14) - pdf(13)) & & + (pdf(17) - pdf(15)) & & + (pdf(18) - pdf(16)) & & + (pdf(23) - pdf(19)) & & + (pdf(24) - pdf(20)) & & + (pdf(25) - pdf(21)) & & + (pdf(26) - pdf(22)) u_y = (pdf(5) - pdf(2)) & & + (pdf(9) - pdf(7)) & & + (pdf(10) - pdf(8)) & & + (pdf(16) - pdf(15)) & & + (pdf(18) - pdf(17)) & & + (pdf(21) - pdf(19)) & & + (pdf(22) - pdf(20)) & & + (pdf(25) - pdf(23)) & & + (pdf(26) - pdf(24)) u_z = (pdf(6) - pdf(3)) & & + (pdf(8) - pdf(7)) & & + (pdf(10) - pdf(9)) & & + (pdf(13) - pdf(11)) & & + (pdf(14) - pdf(12)) & & + (pdf(20) - pdf(19)) & & + (pdf(22) - pdf(21)) & & + (pdf(24) - pdf(23)) & & + (pdf(26) - pdf(25)) u_x = u_x * inv_rho u_y = u_y * inv_rho u_z = u_z * inv_rho auxField(elemOff+vel_pos(1)) = u_x auxField(elemOff+vel_pos(2)) = u_y auxField(elemOff+vel_pos(3)) = u_z end do end subroutine deriveAux_fromState_d3q27 ! **************************************************************************** ! ! ****************************************************************************** ! !> This routine computes velocity from state array !! !! This subroutine's interface must match the abstract interface definition !! [[derive_FromState]] in derived/[[mus_derVarPos_module]].f90 in order to be !! callable via [[mus_derVarPos_type:velFromState]], !! [[mus_derVarPos_type:equilFromState]], !! [[mus_derVarPos_type:momFromState]], !! [[mus_derVarPos_type:velocitiesFromState]], and !! [[mus_derVarPos_type:momentaFromState]] function pointers. subroutine deriveEq_FromState( state, iField, nElems, varSys, layout, res ) ! -------------------------------------------------------------------- ! !> Array of state !! n * layout%fStencil%QQ * nFields real(kind=rk), intent(in) :: state(:) !> Current field integer, intent(in) :: iField !> number of elements integer, intent(in) :: nElems !> variable system which is required to access fieldProp !! information via variable method data c_ptr type(tem_varSys_type), intent(in) :: varSys !> scheme layout contains stencil definition and lattice weights type(mus_scheme_layout_type), intent(in) :: layout !> Output of this routine !! Dimension: n * nComponents of res real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! integer :: QQ, iElem, iDir real(kind=rk) :: pdf( layout%fStencil%QQ ) real(kind=rk) :: rho, vel(3), usq, ucx ! --------------------------------------------------------------------------- QQ = layout%fStencil%QQ do iElem = 1, nElems do iDir = 1, QQ pdf(iDir) = state( iDir+(iElem-1)*varSys%nScalars ) end do rho = sum( pdf ) vel(1) = sum( pdf * layout%fStencil%cxDirRK(1,:) ) / rho vel(2) = sum( pdf * layout%fStencil%cxDirRK(2,:) ) / rho vel(3) = sum( pdf * layout%fStencil%cxDirRK(3,:) ) / rho usq = dot_product(vel, vel) * t2cs2inv do iDir = 1, QQ ucx = dot_product( layout%fStencil%cxDirRK(:, iDir), vel ) ! calculate equilibrium density res( (iElem-1)*QQ+iDir ) = layout%weight( iDir ) * rho & & * ( 1._rk + cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) enddo end do end subroutine deriveEq_FromState ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of moment for 2D !! This routine sets the function Pointer for moment for 2D calcualtion and !! calls the generice get Element from PDF routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! !> For 2D only! recursive subroutine deriveMoment(fun, varsys, elempos, time, tree, nElems, & & nDofs, res ) ! -------------------------------------------------------------------- ! !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Position of the TreeID of the element to get the variable for in the !! global treeID list. integer, intent(in) :: elempos(:) !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nElems !> Number of degrees of freedom within an element. integer, intent(in) :: nDofs !> Resulting values for the requested variable. !! !! Linearized array dimension: !! (n requested entries) x (nComponents of this variable) !! x (nDegrees of freedom) !! Access: (iElem-1)*fun%nComponents*nDofs + !! (iDof-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! -------------------------------------------------------------------- ! !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! -------------------------------------------------------------------- ! fnCalcPtr => mus_deriveMoment call mus_generic_fromPDF_forElement( & & fun = fun, & & varSys = varSys, & & elempos = elempos, & & tree = tree, & & time = time, & & nVals = nElems, & & fnCalcPtr = fnCalcPtr, & & nDofs = nDofs, & & res = res ) end subroutine deriveMoment ! ****************************************************************************** ! ! ************************************************************************* ! ! Subroutines with common interface for values from index ! ! ************************************************************************* ! ! ****************************************************************************** ! !> Initiates the calculation of density. !! This routine sets the function Pointer for density calcualtion and calls !! the generice get Value of Index routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]]. !! recursive subroutine deriveDensity_fromIndex( fun, varSys, time, iLevel, & & idx, idxLen, nVals, res ) !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Level on which values are requested integer, intent(in) :: iLevel !> Index of points in the growing array and variable val array to !! return. !! Size: most times nVals, if contiguous arrays are used it depends !! on the number of first indices integer, intent(in) :: idx(:) !> With idx as start index in contiguous memory, !! idxLength defines length of each contiguous memory !! Size: dependes on number of first index for contiguous array, !! but the sum of all idxLen is equal to nVals integer, optional, intent(in) :: idxLen(:) !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nVals !> Resulting values for the requested variable. !! !! Dimension: n requested entries x nComponents of this variable !! Access: (iElem-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! --------------------------------------------------------------------------- !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! --------------------------------------------------------------------------- fnCalcPtr => mus_derivedensity call mus_generic_varFromPDF_fromIndex( & & fun = fun, & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = idx, & & nVals = nVals, & & fnCalcPtr = fnCalcPtr, & & res = res ) end subroutine deriveDensity_fromIndex ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of pressure. !! This routine sets the function Pointer for pressure calcualtion and calls !! the generice get Value of Index routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]]. !! recursive subroutine derivePressure_fromIndex( fun, varSys, time, iLevel, & & idx, idxLen, nVals, res ) !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Level on which values are requested integer, intent(in) :: iLevel !> Index of points in the growing array and variable val array to !! return. !! Size: most times nVals, if contiguous arrays are used it depends !! on the number of first indices integer, intent(in) :: idx(:) !> With idx as start index in contiguous memory, !! idxLength defines length of each contiguous memory !! Size: dependes on number of first index for contiguous array, !! but the sum of all idxLen is equal to nVals integer, optional, intent(in) :: idxLen(:) !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nVals !> Resulting values for the requested variable. !! !! Dimension: n requested entries x nComponents of this variable !! Access: (iElem-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! --------------------------------------------------------------------------- integer :: dens_pos type(mus_varSys_data_type), pointer :: fPtr ! --------------------------------------------------------------------------- !convert pointer from C to Fotran call C_F_POINTER( fun%method_Data, fPtr ) ! position of density in glob system dens_pos = fun%input_varPos(1) ! get density variable from auxField call varSys%method%val( dens_pos )%get_ValOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fPtr%opData%input_pntIndex(1) & & %indexLvl(iLevel)%val( idx(:) ), & & nVals = nVals, & & res = res ) ! convert density to pressure res(1:nVals) = res(1:nVals) * cs2 end subroutine derivePressure_fromIndex ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of equilibrium. !! This routine sets the function Pointer for equilibrium calcualtion and calls !! the generice get Value of Index routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]]. !! recursive subroutine deriveEquil_fromIndex( fun, varSys, time, iLevel, idx, & & idxLen, nVals, res ) !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Level on which values are requested integer, intent(in) :: iLevel !> Index of points in the growing array and variable val array to !! return. !! Size: most times nVals, if contiguous arrays are used it depends !! on the number of first indices integer, intent(in) :: idx(:) !> With idx as start index in contiguous memory, !! idxLength defines length of each contiguous memory !! Size: dependes on number of first index for contiguous array, !! but the sum of all idxLen is equal to nVals integer, optional, intent(in) :: idxLen(:) !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nVals !> Resulting values for the requested variable. !! !! Dimension: n requested entries x nComponents of this variable !! Access: (iElem-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! --------------------------------------------------------------------------- !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! --------------------------------------------------------------------------- fnCalcPtr => mus_deriveEquil call mus_generic_varFromPDF_fromIndex( & & fun = fun, & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = idx, & & nVals = nVals, & & fnCalcPtr = fnCalcPtr, & & res = res ) end subroutine deriveEquil_fromIndex ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of non_equilibrium. !! This routine sets the function Pointer for non_equilibrium calcualtion and !! calls the generice get Value of Index routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]]. !! recursive subroutine deriveNonEquil_fromIndex( fun, varSys, time, iLevel, idx, & & idxLen, nVals, res ) !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Level on which values are requested integer, intent(in) :: iLevel !> Index of points in the growing array and variable val array to !! return. !! Size: most times nVals, if contiguous arrays are used it depends !! on the number of first indices integer, intent(in) :: idx(:) !> With idx as start index in contiguous memory, !! idxLength defines length of each contiguous memory !! Size: dependes on number of first index for contiguous array, !! but the sum of all idxLen is equal to nVals integer, optional, intent(in) :: idxLen(:) !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nVals !> Resulting values for the requested variable. !! !! Dimension: n requested entries x nComponents of this variable !! Access: (iElem-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! --------------------------------------------------------------------------- !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! --------------------------------------------------------------------------- fnCalcPtr => mus_deriveNonEquil call mus_generic_varFromPDF_fromIndex( & & fun = fun, & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = idx, & & nVals = nVals, & & fnCalcPtr = fnCalcPtr, & & res = res ) end subroutine deriveNonEquil_fromIndex ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of kinetic_energy. !! This routine sets the function Pointer for kinetic_energy calcualtion and !! calls the generice get Value of Index routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]]. !! recursive subroutine deriveKe_fromIndex( fun, varSys, time, iLevel, idx, & & idxLen, nVals, res ) ! --------------------------------------------------------------------------- !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Level on which values are requested integer, intent(in) :: iLevel !> Index of points in the growing array and variable val array to !! return. !! Size: most times nVals, if contiguous arrays are used it depends !! on the number of first indices integer, intent(in) :: idx(:) !> With idx as start index in contiguous memory, !! idxLength defines length of each contiguous memory !! Size: dependes on number of first index for contiguous array, !! but the sum of all idxLen is equal to nVals integer, optional, intent(in) :: idxLen(:) !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nVals !> Resulting values for the requested variable. !! !! Dimension: n requested entries x nComponents of this variable !! Access: (iElem-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! --------------------------------------------------------------------------- type(mus_varSys_data_type), pointer :: fPtr integer :: dens_pos, vel_pos, iVal real(kind=rk) :: mass_dens(nVals), vel(3) ! --------------------------------------------------------------------------- call C_F_POINTER( fun%method_Data, fPtr ) res = 0.0_rk ! get mass density values for IDX dens_pos = fun%input_varPos(1) call varSys%method%val( dens_pos )%get_ValOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fPtr%opData%input_pntIndex(1) & & %indexLvl(iLevel)%val( idx(:) ), & & nVals = nVals, & & res = mass_dens ) ! get velocity values for IDX vel_pos = fun%input_varPos(2) call varSys%method%val( vel_pos )%get_ValOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fPtr%opData%input_pntIndex(2) & & %indexLvl(iLevel)%val( idx(:) ), & & nVals = nVals, & & res = res(:) ) ! convert velocity to momentum do iVal = 1, nVals vel(1) = res((iVal-1)*3 + 1) vel(2) = res((iVal-1)*3 + 2) vel(3) = res((iVal-1)*3 + 3) res(iVal) = sum( vel(:)*vel(:) ) * 0.5_rk * mass_dens(iVal) end do end subroutine deriveKe_fromIndex ! ****************************************************************************** ! ! ****************************************************************************** ! !> Initiates the calculation of kinetic_energy. !! This routine sets the function Pointer for kinetic_energy calcualtion and !! calls the generice get Value of Index routine !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]]. !! recursive subroutine deriveStrainRate_fromIndex( fun, varSys, time, iLevel, & & idx, idxLen, nVals, res ) !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Level on which values are requested integer, intent(in) :: iLevel !> Index of points in the growing array and variable val array to !! return. !! Size: most times nVals, if contiguous arrays are used it depends !! on the number of first indices integer, intent(in) :: idx(:) !> With idx as start index in contiguous memory, !! idxLength defines length of each contiguous memory !! Size: dependes on number of first index for contiguous array, !! but the sum of all idxLen is equal to nVals integer, optional, intent(in) :: idxLen(:) !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nVals !> Resulting values for the requested variable. !! !! Dimension: n requested entries x nComponents of this variable !! Access: (iElem-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! --------------------------------------------------------------------------- !> Function pointer to perform specific operation. procedure(mus_derive_fromPDF), pointer :: fnCalcPtr ! --------------------------------------------------------------------------- fnCalcPtr => mus_deriveStrainRate call mus_generic_varFromPDF_fromIndex( & & fun = fun, & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = idx, & & nVals = nVals, & & fnCalcPtr = fnCalcPtr, & & res = res ) end subroutine deriveStrainRate_fromIndex ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate Momentum from density and velocity in auxField. !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]]. !! recursive subroutine deriveMomentum_fromIndex( fun, varSys, time, iLevel, & & idx, idxLen, nVals, res ) ! --------------------------------------------------------------------------- !> Description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Level on which values are requested integer, intent(in) :: iLevel !> Index !! return. !! Size: most times nVals, if contiguous arrays are used it depends !! on the number of first indices integer, intent(in) :: idx(:) !> With idx as start index in contiguous memory, !! idxLength defines length of each contiguous memory !! Size: dependes on number of first index for contiguous array, !! but the sum of all idxLen is equal to nVals integer, optional, intent(in) :: idxLen(:) !> Number of values to obtain for this variable (vectorized access). integer, intent(in) :: nVals !> Resulting values for the requested variable. !! !! Dimension: n requested entries x nComponents of this variable !! Access: (iElem-1)*fun%nComponents + iComp real(kind=rk), intent(out) :: res(:) ! --------------------------------------------------------------------------- type(mus_varSys_data_type), pointer :: fPtr integer :: dens_pos, vel_pos, iVal real(kind=rk) :: mass_dens(nVals) ! --------------------------------------------------------------------------- call C_F_POINTER( fun%method_Data, fPtr ) res = 0.0_rk ! get mass density values for IDX dens_pos = fun%input_varPos(1) call varSys%method%val( dens_pos )%get_ValOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fPtr%opData%input_pntIndex(1) & & %indexLvl(iLevel)%val( idx(:) ), & & nVals = nVals, & & res = mass_dens ) ! get velocity values for IDX vel_pos = fun%input_varPos(2) call varSys%method%val( vel_pos )%get_ValOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fPtr%opData%input_pntIndex(2) & & %indexLvl(iLevel)%val( idx(:) ), & & nVals = nVals, & & res = res(:) ) ! convert velocity to momentum do iVal = 1, nVals res( (iVal-1)*3 + 1 : iVal*3 ) = res( (iVal-1)*3 + 1 : iVal*3 ) & & * mass_dens(iVal) end do end subroutine deriveMomentum_fromIndex ! ****************************************************************************** ! ! ************************************************************************* ! ! Subroutines with common interface called by fnCalcPtr ! ! ************************************************************************* ! ! ****************************************************************************** ! !> Calculate the density of a given set of elements (sum up all links). !! This routine is used to compute density for all scheme kinds !! For multispecies, it can compute both species density and mixture density !! !! The interface has to comply to the abstract interface !! [[mus_varSys_module:mus_derive_fromPDF]]. !! recursive subroutine mus_derivedensity(fun, varsys, stencil, iLevel, & & posInState, pdf, res, nVals ) !> description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varsys_op_type), intent(in) :: fun !> the variable system to obtain the variable from. type(tem_varsys_type), intent(in) :: varsys !> fluid stencil defintion type(tem_stencilHeader_type), intent(in) :: stencil !> current Level integer, intent(in) :: iLevel !> Position of element in levelwise state array integer, intent(in) :: posInState(:) !> pdf array real(kind=rk), intent(in) :: pdf(:) !> results real(kind=rk), intent(out) :: res(:) !> nVals to get integer, intent(in) :: nVals ! --------------------------------------------------------------------------- real(kind=rk), allocatable :: tmpPDF(:) real(kind=rk) :: dens integer :: iDep, iVal, iComp integer :: pdfPos, nCompPDF ! --------------------------------------------------------------------------- pdfPos = fun%input_varPos(1) nCompPDF = varSys%method%val(pdfPos)%nComponents allocate( tmpPDF( nCompPDF ) ) res = 0.0_rk do iVal = 1, nVals dens = 0.0_rk do iDep = 1, fun%nInputs tmpPDF = pdf( & & (iDep-1)*nVals + (iVal-1)*nCompPDF + 1 & & : & & (iDep-1)*nVals + iVal*nCompPDF ) do iComp = 1, nCompPDF dens = dens + tmpPDF( iComp ) end do !iComp end do !iDep res( iVal ) = dens end do !iVal deallocate( tmpPDF ) end subroutine mus_deriveDensity ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate the equlibrium of given elements with the given input state !! array. !! !! The interface has to comply to the abstract interface !! [[mus_varSys_module:mus_derive_fromPDF]]. !! recursive subroutine mus_deriveEquil(fun, varsys, stencil, iLevel, & & posInState, pdf, res, nVals ) !> description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varsys_op_type), intent(in) :: fun !> the variable system to obtain the variable from. type(tem_varsys_type), intent(in) :: varsys !> fluid stencil defintion type(tem_stencilHeader_type), intent(in) :: stencil !> current Level integer, intent(in) :: iLevel !> Position of element in levelwise state array integer, intent(in) :: posInState(:) !> pdf array real(kind=rk), intent(in) :: pdf(:) !> results real(kind=rk), intent(out) :: res(:) !> nVals to get integer, intent(in) :: nVals ! --------------------------------------------------------------------------- type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk), allocatable :: tmpPDF(:) real(kind=rk), allocatable :: fEq(:) real(kind=rk) :: dens, vel(3), densInv integer :: pdfPos, nCompsPDF, iVal, iDir ! --------------------------------------------------------------------------- call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme pdfPos = fun%input_varPos(1) nCompsPDF = varSys%method%val( pdfPos )%nComponents allocate( tmpPDF( nCompsPDF ) ) allocate( fEq( fun%nComponents ) ) res = 0.0_rk do iVal = 1 , nVals tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF ) ! computes density and velocity dens = sum(tmpPDF) densInv = 1.0_rk / dens vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:)) & & * densInv vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) & & * densInv vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) & & * densInv do iDir = 1, scheme%layout%fStencil%QQ ! calculate equilibrium density feq(idir) = scheme%layout%weight( idir ) & & * dens & & * ( 1._rk & & + ( cs2inv & & * sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) end do res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = fEq end do !iVal deallocate( tmpPDF ) deallocate( fEq ) end subroutine mus_deriveEquil ! ****************************************************************************** ! ! ****************************************************************************** ! !> For 2D only! !! !! The interface has to comply to the abstract interface !! [[mus_varSys_module:mus_derive_fromPDF]]. !! recursive subroutine mus_deriveMoment(fun, varsys, stencil, iLevel, & & posInState, pdf, res, nVals ) !> description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varsys_op_type), intent(in) :: fun !> the variable system to obtain the variable from. type(tem_varsys_type), intent(in) :: varsys !> fluid stencil defintion type(tem_stencilHeader_type), intent(in) :: stencil !> current Level integer, intent(in) :: iLevel !> Position of element in levelwise state array integer, intent(in) :: posInState(:) !> pdf array real(kind=rk), intent(in) :: pdf(:) !> results real(kind=rk), intent(out) :: res(:) !> nVals to get integer, intent(in) :: nVals ! --------------------------------------------------------------------------- real(kind=rk), allocatable :: tmpPDF(:) integer :: pdfPos, nCompsPDF integer :: iVal ! --------------------------------------------------------------------------- pdfPos = fun%input_varPos(1) nCompsPDF = varSys%method%val( pdfPos )%nComponents allocate( tmpPDF( nCompsPDF ) ) res = 0.0_rk do iVal = 1, nVals tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF ) res( (iVal-1)*9+1 ) = sum( tmpPDF ) - 1.0_rk res( (iVal-1)*9+2 ) = - sum(tmpPDF(1:4)) + 2.0_rk & & * sum(tmpPDF(5:8)) - 4.0_rk & & * tmpPDF(9) + 2.0_rk res( (iVal-1)*9+3 ) = - 2.0_rk * sum(tmpPDF(1:4)) & & + sum(tmpPDF(5:8)) + 4.0_rk & & * tmpPDF(9) - 1.0_rk res( (iVal-1)*9+4 ) = tmpPDF(3) + tmpPDF(7) + tmpPDF(8) & & - tmpPDF(1) - tmpPDF(5) - tmpPDF(6) res( (iVal-1)*9+6 ) = tmpPDF(4) + tmpPDF(6) + tmpPDF(8) & & - tmpPDF(2) - tmpPDF(5) - tmpPDF(7) res( (iVal-1)*9+5 ) = 2.0_rk * ( tmpPDF(1) - tmpPDf(3) ) - tmpPDF(5) & & - tmpPDF(6) + tmpPDF(7) + tmpPDF(8) res( (iVal-1)*9+7 ) = 2.0_rk * ( tmpPDF(2) - tmpPDF(4) ) - tmpPDf(5) & & + tmpPDF(6) - tmpPDF(7) + tmpPDF(8) res( (iVal-1)*9+8 ) = tmpPDF(1) - tmpPDF(2) + tmpPDF(3) - tmpPDF(4) res( (iVal-1)*9+9 ) = tmpPDF(5) - tmpPDF(6) - tmpPDF(7) + tmpPDF(8) end do ! iVal deallocate( tmpPDF ) end subroutine mus_deriveMoment ! ****************************************************************************** ! ! ****************************************************************************** ! !> Calculate the Non-Equlibrium !! recursive subroutine mus_deriveNonEquil(fun, varsys, stencil, iLevel, & & posInState, pdf, res, nVals ) !> description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varsys_op_type), intent(in) :: fun !> the variable system to obtain the variable from. type(tem_varsys_type), intent(in) :: varsys !> fluid stencil defintion type(tem_stencilHeader_type), intent(in) :: stencil !> current Level integer, intent(in) :: iLevel !> Position of element in levelwise state array integer, intent(in) :: posInState(:) !> pdf array real(kind=rk), intent(in) :: pdf(:) !> results real(kind=rk), intent(out) :: res(:) !> nVals to get integer, intent(in) :: nVals ! --------------------------------------------------------------------------- type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk), allocatable :: tmpPDF(:) real(kind=rk), allocatable :: fEq(:) real(kind=rk) :: dens, vel(3) integer :: iDir, iVal integer :: pdfPos, nCompsPDF ! --------------------------------------------------------------------------- call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme pdfPos = fun%input_varPos(1) nCompsPDF = varSys%method%val( pdfPos )%nComponents allocate( fEq( fun%nComponents ) ) allocate( tmpPDF( nCompsPDF ) ) res = 0.0_rk do iVal = 1 , nVals tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF ) ! computes density and velocity dens = sum(tmpPDF) vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:)) / dens vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) / dens vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) / dens ! computes equilibrium do iDir = 1, scheme%layout%fStencil%QQ ! calculate equilibrium density feq(idir) = scheme%layout%weight( idir ) & & * dens & & * ( 1._rk & & + ( cs2inv & & * sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) end do res( (iVal-1)*fun%nComponents+1 : ival*fun%nComponents ) = tmpPDF - fEq end do !iVal deallocate( tmpPDF ) deallocate( fEq ) end subroutine mus_deriveNonEquil ! ****************************************************************************** ! ! ****************************************************************************** ! !> author: Jiaxing Qi !! Calculate the strain rate ( or rate of strain, or rate of deformation) !! !! The equation is: !! \[ !! \tau_{\alpha \beta}= !! -\frac{3\omega}{2\rho} \sum_{i} f^{neq}_{i} c_{i\alpha} c_{i\beta} !! \] !! where \( \tau_{\alpha \beta}\) is the stress !! in the \(\beta\)-direction on a face normal to the \(\alpha\)-axis,\n !! \( f^{neq}_i = f_i - f^{eq}_i\) is the non-equilibrium pdf.\n !! For more information, please refer to: equation 45 in\n !! Krueger T, Varnik F, Raabe D. Shear stress in lattice Boltzmann !! simulations. Physical Review E. 2009;79(4):1-14.\n !! !! For multi-level mesh, Omega on finer level needs to be adjusted in order to !! get the correct shearstress calculation.\n !! First, we defines c as the dx ratio between finer and coarse level.\n !! \( c={ \Delta dx }_{ c }/{ \Delta dx }_{ f } \) !! Then the viscosity on the different levels must satisfy:\n !! \( \frac { { \nu }_{ f } }{ { \nu }_{ c } } =c \) !! This constrain leads to a relationship of omega on different levels:\n !! \( {\omega}_f = \frac {1}{ {\lambda}(\frac{1}{{\omega}_c}-0.5)+0.5 } \) !! For more information, please refer to:\n !! Manuel H, Harald K, Joerg B, Sabine R. Aeroacoustic validation of the !! lattice boltzmann method on non-uniform grids. ECCOMAS 2012 !! recursive subroutine mus_deriveStrainRate(fun, varsys, stencil, iLevel, & & posInState, pdf, res, nVals ) !> description of the method to obtain the variables, here some preset !! values might be stored, like the space time function to use or the !! required variables. class(tem_varsys_op_type), intent(in) :: fun !> the variable system to obtain the variable from. type(tem_varsys_type), intent(in) :: varsys !> fluid stencil defintion type(tem_stencilHeader_type), intent(in) :: stencil !> current level integer, intent(in) :: iLevel !> Position of element in levelwise state array integer, intent(in) :: posInState(:) !> pdf array real(kind=rk), intent(in) :: pdf(:) !> results real(kind=rk), intent(out) :: res(:) !> nVals to get integer, intent(in) :: nVals ! --------------------------------------------------------------------------- type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk), allocatable :: tmpPDF(:) real(kind=rk), allocatable :: nonEq(:) real(kind=rk), allocatable :: tau(:) real(kind=rk), allocatable :: fEq(:) real(kind=rk) :: dens, vel(3), omega integer :: pdfPos, nCompsPDF, iVal, iDir ! --------------------------------------------------------------------------- call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme pdfPos = fun%input_varPos(1) nCompsPDF = varSys%method%val( pdfPos )%nComponents allocate( tmpPDF( nCompsPDF ) ) allocate( nonEq( nCompsPDF ) ) allocate( fEq( nCompsPDF ) ) allocate( tau( fun%nComponents ) ) do iVal = 1, nVals tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF ) ! computes density and velocity dens = sum(tmpPDF) vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:)) / dens vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) / dens vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) / dens ! computes equilibrium do iDir = 1, scheme%layout%fStencil%QQ ! calculate equilibrium density feq(idir) = scheme%layout%weight( idir ) & & * dens & & * ( 1._rk & & + ( cs2inv & & * sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(scheme%layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) end do ! Non-Eq nonEq = tmpPDF - fEq ! get the correct omega value omega = scheme%field(1)%fieldProp%fluid%viscKine%omLvl(iLevel) & & %val(posInState(iVal)) ! compute shear stress tau(:) = secondMom( scheme%layout%fStencil%cxcx(:,:), & & nonEq(:), & & scheme%layout%fStencil%QQ ) ! S_{ij} = - (omega/(2*rho*cs2))* sum_k(c_{ki} c{kj} f^neq_k) res( (iVal-1)*fun%nComponents + 1 : iVal*fun%nComponents ) & & = tau(:) * (-1.5_rk) * omega / dens end do !iVal deallocate( nonEq ) deallocate( tau ) deallocate( fEq ) deallocate( tmpPDF ) end subroutine mus_deriveStrainRate ! ****************************************************************************** ! ! **************************************************************************** ! !> Update state with source variable "absorb_layer". !! absorb_layer is used to absorb the flow and gradually reduce the flow !! quantities like pressure and velocity to a fixed value. !! It is based on: !! Xu, H., & Sagaut, P. (2013). Analysis of the absorbing layers for the !! weakly-compressible lattice Boltzmann methods. Journal of Computational !! Physics, 245(x), 14–42. !! Jacob, J.; Sagaut, P. (2019): Solid wall and open boundary conditions in !! hybrid recursive regularized lattice Boltzmann method for compressible !! flows. In Physics of Fluids 31 (12), p. 126103. !! !! This subroutine's interface must match the abstract interface definition !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to !! be callable via [[mus_source_op_type:applySrc]] function pointer. subroutine applySrc_HRRCorrection( fun, inState, outState, neigh, auxField, & & nPdfSize, iLevel, varSys, time, phyConvFac, & & derVarPos ) ! -------------------------------------------------------------------- ! !> Description of method to apply source terms class(mus_source_op_type), intent(in) :: fun !> input pdf vector real(kind=rk), intent(in) :: inState(:) !> output pdf vector real(kind=rk), intent(inout) :: outState(:) !> connectivity Array corresponding to state vector integer,intent(in) :: neigh(:) !> auxField array real(kind=rk), intent(in) :: auxField(:) !> number of elements in state Array integer, intent(in) :: nPdfSize !> current level integer, intent(in) :: iLevel !> variable system type(tem_varSys_type), intent(in) :: varSys !> Point in time at which to evaluate the variable. type(tem_time_type), intent(in) :: time !> Physics conversion factor for current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------- ! type(mus_varSys_data_type), pointer :: fPtr integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- !call tem_abort('Error: Absorb layer is not yet implemented') ! convert c pointer to solver type fortran pointer call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, & & fPtr ) ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! constant parameter nScalars = varSys%nScalars associate( layout => fPtr%solverData%scheme%layout, & & HRR_Corr => fun%elemLvl(iLevel)%HRR_Corr ) QQ = layout%fStencil%QQ !$omp parallel do schedule(static), private( posInTotal, iElem, iDir, statepos ) do iElem = 1, nElems posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem) ! correction term ! f + dt/2 * HRR_Corr do iDir = 1, QQ ! position in state array statePos = ( posintotal-1)* nscalars+idir+( 1-1)* qq ! update outstate outState(statePos) = outState(statePos) + 0.5_rk & & * HRR_Corr%phi((iElem-1)*QQ+iDir) end do end do !iElem !$omp end parallel do end associate end subroutine applySrc_HRRCorrection ! **************************************************************************** ! end module mus_derQuan_module