! Copyright (c) 2013, 2016, 2019 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2013-2020 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) 2016 Verena Krupp <verena.krupp@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2016-2018 Raphael Haupt <raphael.haupt@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. ! ************************************************************************** ! !> author: Kannan Masilamani !! author: Jiaxing Qi !! This module provides the MUSUBI specific functions for calculating !! macroscopic quantities from the state variables for incompressible LBM !! models.\n !! Notice that only those quantities that related to density should have a !! formula which differs from normal LBM model. !! !! 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-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) 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_derQuanIncomp_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, & & div1_36, div3_4h, & & sqrt3, cs2inv, cs2, t2cs4inv, t2cs2inv, & & cs4inv, rho0, rho0Inv, q000 use env_module, only: rk, long_k, labelLen use tem_float_module, only: operator(.feq.), operator(.fge.), & & operator(.fle.) 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_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_subTree_type_module, only: tem_subTree_type, tem_treeIDfrom_subTree use tem_aux_module, only: tem_abort use tem_logging_module, only: logUnit 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_debug_module, only: dbgUnit use tem_grow_array_module, only: grw_labelarray_type, append ! include musubi modules use mus_source_type_module, only: mus_source_op_type use mus_pdf_module, only: pdf_data_type 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_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 use mus_derQuan_module, only: deriveDensity, & & deriveDensity_fromIndex, & & derivePressure, & & derivePressure_fromIndex, & & deriveShearStress, & & deriveShearMag, & & deriveWSS2D, & & deriveWSS3D, & & deriveTemp, & & deriveShearRate, & & deriveBndForce use mus_operation_var_module, only: mus_opVar_setupIndices, & & mus_opVar_gradU_forElement, & & mus_opVar_vorticity_forElement, & & mus_opVar_QCriterion_forElement use mus_derVarPos_module, only: mus_derVarPos_type use mus_physics_module, only: mus_convertFac_type implicit none private public :: mus_append_derVar_fluidIncomp ! equilbrium from macro uses different interface defined in ! mus_variable_module public :: deriveEquilIncomp_FromMacro public :: deriveEquilIncomp_FromMacro_d3q19 public :: deriveVelIncomp_FromState public :: deriveVelIncomp_FromState_d3q19 public :: deriveVelIncomp_FromPreColState public :: deriveEqIncomp_FromState public :: deriveAuxIncomp_fromState public :: deriveEquilIncomp_fromAux ! source variables public :: derive_absorbLayerIncomp public :: applySrc_absorbLayerIncomp contains ! ************************************************************************ ! !> subroutine to add derive variables for incompressible LBM !! (schemekind = 'fluid_incompressible') to the varsys. subroutine mus_append_derVar_fluidIncomp( 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 ', 'fetch_pdf ', & & '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 set_params => tem_varSys_setParams_dummy get_params => tem_varSys_getParams_dummy method_data = mus_get_new_solver_ptr(solverData) 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 => deriveEquilIncomp get_valOfIndex => deriveEquilIncomp_fromIndex nComponents = stencil%QQ allocate(input_varname(1)) input_varname(1) = 'pdf' case ('non_equilibrium') get_element => deriveNonEquilIncomp get_valOfIndex => deriveNonEquilIncomp_fromIndex nComponents = stencil%QQ allocate(input_varname(1)) input_varname(1) = 'fetch_pdf_now' case ('kinetic_energy') get_element => deriveKeIncomp get_ValOfIndex => deriveKEIncomp_fromIndex nComponents = 1 allocate(input_varname(1)) input_varname(1) = '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 => deriveStrainRateIncomp get_ValOfIndex => deriveStrainRateIncomp_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 => deriveMomentumIncomp get_valOfIndex => deriveMomentumIncomp_fromIndex nComponents = 3 allocate(input_varname(1)) input_varname(1) = '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 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_fluidIncomp ! ************************************************************************ ! ! **************************************************************************** ! ! Subroutines with common interface for the function pointers ! ! **************************************************************************** ! ! **************************************************************************** ! !> Calculate momentum from velocity stored in auxField !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveMomentumIncomp(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 :: vel_pos(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 velocity in auxField array vel_pos = varSys%method%val( fun%input_varPos(1) ) & & %auxField_varPos(1:3) ! compute and store momentum res((iElem-1)*3 + 1) = rho0 * auxField(iLevel)%val(elemOff + vel_pos(1)) res((iElem-1)*3 + 2) = rho0 * auxField(iLevel)%val(elemOff + vel_pos(2)) res((iElem-1)*3 + 3) = rho0 * auxField(iLevel)%val(elemOff + vel_pos(3)) end do ! iElem end associate end subroutine deriveMomentumIncomp ! **************************************************************************** ! ! **************************************************************************** ! !> 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 deriveEquilIncomp(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_deriveEquilIncomp 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 deriveEquilIncomp ! **************************************************************************** ! ! **************************************************************************** ! !> 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 deriveNonEquilIncomp(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_deriveNonEquilIncomp 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 deriveNonEquilIncomp ! **************************************************************************** ! ! **************************************************************************** ! !> 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 deriveStrainRateIncomp(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_deriveStrainRateIncomp 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 deriveStrainRateIncomp ! **************************************************************************** ! ! **************************************************************************** ! ! Subroutines with common interface for the function pointers ! ! getValOfIndex ! ! **************************************************************************** ! ! **************************************************************************** ! !> 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 deriveMomentumIncomp_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 :: vel_pos, iVal ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) res = 0.0_rk ! get velocity values for IDX vel_pos = fun%input_varPos(1) call varSys%method%val( vel_pos )%get_ValOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fPtr%opData%input_pntIndex(1) & & %indexLvl(iLevel)%val( idx(:) ), & & nVals = nVals, & & res = res(:) ) ! convert velocity to momentum do iVal = 1, nVals res( (iVal-1)*3 + 1 : iVal*3 ) = rho0 * res( (iVal-1)*3 + 1 : iVal*3 ) end do end subroutine deriveMomentumIncomp_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 deriveEquilIncomp_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_deriveEquilIncomp call mus_generic_varFromPDF_fromIndex( & & fun = fun, & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = idx, & & nVals = nVals, & & fnCalcPtr = fnCalcPtr, & & res = res ) end subroutine deriveEquilIncomp_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 deriveNonEquilIncomp_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_deriveNonEquilIncomp call mus_generic_varFromPDF_fromIndex( & & fun = fun, & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = idx, & & nVals = nVals, & & fnCalcPtr = fnCalcPtr, & & res = res ) end subroutine deriveNonEquilIncomp_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 deriveKeIncomp_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 :: vel_pos, iVal real(kind=rk) :: vel(3) ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) res = 0.0_rk ! get velocity values for IDX vel_pos = fun%input_varPos(1) call varSys%method%val( vel_pos )%get_ValOfIndex( & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = fPtr%opData%input_pntIndex(1) & & %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 * rho0 end do end subroutine deriveKeIncomp_fromIndex ! **************************************************************************** ! ! **************************************************************************** ! !> Initiates the calculation of StrainRate. !! This routine sets the function Pointer for StrainRate 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 deriveStrainRateIncomp_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_deriveStrainRateIncomp call mus_generic_varFromPDF_fromIndex( & & fun = fun, & & varSys = varSys, & & time = time, & & iLevel = iLevel, & & idx = idx, & & nVals = nVals, & & fnCalcPtr = fnCalcPtr, & & res = res ) end subroutine deriveStrainRateIncomp_fromIndex ! **************************************************************************** ! ! **************************************************************************** ! ! Calculation routines ! ! **************************************************************************** ! ! **************************************************************************** ! !> Calculate kinetic energy from velocity in auxField !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine deriveKeIncomp(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 :: vel_pos(3) real(kind=rk) :: 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 velocity in auxField array vel_pos = varSys%method%val( fun%input_varPos(1) ) & & %auxField_varPos(1:3) ! 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 * rho0 end do ! iElem end associate end subroutine deriveKeIncomp ! ************************************************************************** ! ! **************************************************************************** ! !> Derive absorb layer variable defined as a source term. recursive subroutine derive_absorbLayerIncomp(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_absorbLayerIncomp ! **************************************************************************** ! ! **************************************************************************** ! !> 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 for incompressible !! model. 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. !! !! 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_absorbLayerIncomp( 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(:) ! -------------------------------------------------------------------- ! call tem_abort('Error: Absorb layer is not yet implemented') end subroutine applySrc_absorbLayerIncomp ! **************************************************************************** ! ! **************************************************************************** ! !> This routine computes equilbrium from density and velocity !! !! 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 deriveEquilIncomp_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) & & + ( cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) ) enddo res( (iElem-1)*QQ+1: iElem*QQ ) = fEq end do end subroutine deriveEquilIncomp_FromMacro ! **************************************************************************** ! ! **************************************************************************** ! !> This routine computes equilbrium from density and velocity !! !! 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 deriveEquilIncomp_FromMacro_d3q19( 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(:) ! -------------------------------------------------------------------- ! integer :: iElem real(kind=rk) :: fEq(19), u_x, u_y, u_z real(kind=rk) :: usq, usqn, usqn_o2 real(kind=rk) :: coeff_1, coeff_2 real(kind=rk) :: ui1, ui3, ui10, ui11, ui12, ui13 real(kind=rk) :: fac_1, fac_2, fac_3, fac_4, fac_9, fac_10, fac_11, fac_12,& & fac_13 real(kind=rk) :: sum1_1, sum1_2, sum2_1, sum2_2, sum3_1, sum3_2, sum4_1, & & sum4_2, sum9_1, sum9_2, sum10_1, sum10_2, sum11_1, & & sum11_2, sum12_1, sum12_2, sum13_1, sum13_2 ! -------------------------------------------------------------------- ! do iElem = 1, nElems u_x = velocity(1,iElem) u_y = velocity(2,iElem) u_z = velocity(3,iElem) usq = u_x*u_x + u_y*u_y + u_z*u_z usqn = div1_36 * (density(ielem) - 1.5_rk * usq * rho0) feq(19) = 12.0_rk * usqn coeff_1 = 0.125_rk * rho0 ui1 = u_x + u_y fac_1 = coeff_1 * ui1 sum1_1 = fac_1 * div3_4h sum1_2 = fac_1 * ui1 + usqn feq(18) = sum1_1 +sum1_2 feq(15) = -sum1_1 +sum1_2 ui3 = -u_x + u_y fac_3 = coeff_1 * ui3 sum3_1 = fac_3 * div3_4h sum3_2 = fac_3 * ui3 + usqn feq(16) = sum3_1 +sum3_2 feq(17) = -sum3_1 +sum3_2 ui10 = u_x + u_z fac_10 = coeff_1 * ui10 sum10_1 = fac_10 * div3_4h sum10_2 = fac_10 * ui10 + usqn feq(14) = sum10_1+sum10_2 feq(11) = -sum10_1+sum10_2 ui12 = -u_x + u_z fac_12 = coeff_1 * ui12 sum12_1 = fac_12 * div3_4h sum12_2 = fac_12 * ui12 + usqn feq(13) = sum12_1+sum12_2 feq(12) = -sum12_1+sum12_2 ui11 = u_y + u_z fac_11 = coeff_1 * ui11 sum11_1 = fac_11 * div3_4h sum11_2 = fac_11 * ui11 + usqn feq(10) = sum11_1+sum11_2 feq( 7) = -sum11_1+sum11_2 ui13 = -u_y + u_z fac_13 = coeff_1 * ui13 sum13_1 = fac_13 * div3_4h sum13_2 = fac_13 * ui13 + usqn feq( 8) = sum13_1+sum13_2 feq( 9) = -sum13_1+sum13_2 coeff_2 = 0.25_rk * rho0 usqn_o2 = 2.0_rk * usqn fac_2 = coeff_2 * u_y sum2_1 = fac_2 * div3_4h sum2_2 = fac_2 * u_y + usqn_o2 feq( 5) = sum2_1 +sum2_2 feq( 2) = -sum2_1 +sum2_2 fac_4 = coeff_2 * u_x sum4_1 = fac_4 * div3_4h sum4_2 = fac_4 * u_x + usqn_o2 feq( 1) = -sum4_1 +sum4_2 feq( 4) = sum4_1 +sum4_2 fac_9 = coeff_2 * u_z sum9_1 = fac_9 * div3_4h sum9_2 = fac_9 * u_z + usqn_o2 feq( 6) = sum9_1 +sum9_2 feq( 3) = -sum9_1 +sum9_2 res( (iElem-1)*19 + 1) = fEq( 1) res( (iElem-1)*19 + 2) = fEq( 2) res( (iElem-1)*19 + 3) = fEq( 3) res( (iElem-1)*19 + 4) = fEq( 4) res( (iElem-1)*19 + 5) = fEq( 5) res( (iElem-1)*19 + 6) = fEq( 6) res( (iElem-1)*19 + 7) = fEq( 7) res( (iElem-1)*19 + 8) = fEq( 8) res( (iElem-1)*19 + 9) = fEq( 9) res( (iElem-1)*19 + 10) = fEq(10) res( (iElem-1)*19 + 11) = fEq(11) res( (iElem-1)*19 + 12) = fEq(12) res( (iElem-1)*19 + 13) = fEq(13) res( (iElem-1)*19 + 14) = fEq(14) res( (iElem-1)*19 + 15) = fEq(15) res( (iElem-1)*19 + 16) = fEq(16) res( (iElem-1)*19 + 17) = fEq(17) res( (iElem-1)*19 + 18) = fEq(18) res( (iElem-1)*19 + 19) = fEq(19) end do end subroutine deriveEquilIncomp_FromMacro_d3q19 ! **************************************************************************** ! ! **************************************************************************** ! !> 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 deriveEquilIncomp_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 & & + rho0*( cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) ) enddo end do end subroutine deriveEquilIncomp_fromAux ! **************************************************************************** ! ! **************************************************************************** ! !> 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 deriveVelIncomp_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) :: f( layout%fStencil%QQ ) ! -------------------------------------------------------------------- ! do iElem = 1, nElems do iDir = 1, layout%fStencil%QQ f(iDir) = state( iDir+(iElem-1)* varSys%nScalars ) end do res( (iElem-1)*3+1 ) = sum( f * layout%fStencil%cxDirRK(1,:) ) res( (iElem-1)*3+2 ) = sum( f * layout%fStencil%cxDirRK(2,:) ) res( (iElem-1)*3+3 ) = sum( f * layout%fStencil%cxDirRK(3,:) ) end do end subroutine deriveVelIncomp_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 deriveVelIncomp_FromState_d3q19( 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 real(kind=rk) :: f(18), u_x, u_y, u_z ! -------------------------------------------------------------------- ! !NEC$ ivdep do iElem = 1, nElems f( 1) = state( 1+(iElem-1)*varSys%nScalars ) f( 2) = state( 2+(iElem-1)*varSys%nScalars ) f( 3) = state( 3+(iElem-1)*varSys%nScalars ) f( 4) = state( 4+(iElem-1)*varSys%nScalars ) f( 5) = state( 5+(iElem-1)*varSys%nScalars ) f( 6) = state( 6+(iElem-1)*varSys%nScalars ) f( 7) = state( 7+(iElem-1)*varSys%nScalars ) f( 8) = state( 8+(iElem-1)*varSys%nScalars ) f( 9) = state( 9+(iElem-1)*varSys%nScalars ) f(10) = state( 10+(iElem-1)*varSys%nScalars ) f(11) = state( 11+(iElem-1)*varSys%nScalars ) f(12) = state( 12+(iElem-1)*varSys%nScalars ) f(13) = state( 13+(iElem-1)*varSys%nScalars ) f(14) = state( 14+(iElem-1)*varSys%nScalars ) f(15) = state( 15+(iElem-1)*varSys%nScalars ) f(16) = state( 16+(iElem-1)*varSys%nScalars ) f(17) = state( 17+(iElem-1)*varSys%nScalars ) f(18) = state( 18+(iElem-1)*varSys%nScalars ) u_x = (f(4) - f(1)) & & + (f(12) - f(11)) & & + (f(14) - f(13)) & & + (f(17) - f(15)) & & + (f(18) - f(16)) u_y = (f(5) - f(2)) & & + (f(9) - f(7)) & & + (f(10) - f(8)) & & + (f(16) - f(15)) & & + (f(18) - f(17)) u_z = (f(6) - f(3)) & & + (f(8) - f(7)) & & + (f(10) - f(9)) & & + (f(13) - f(11)) & & + (f(14) - f(12)) u_x = u_x * 1.0_rk u_y = u_y * 1.0_rk u_z = u_z * 1.0_rk res( (iElem-1)*3+1 ) = u_x res( (iElem-1)*3+2 ) = u_y res( (iElem-1)*3+3 ) = u_z end do end subroutine deriveVelIncomp_FromState_d3q19 ! **************************************************************************** ! ! **************************************************************************** ! !> This routine computes velocity from precollision state array using FETCH !! macro. !! !! 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 deriveVelIncomp_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) :: 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 ! 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) end do ! convert to velocity res = res / rho0 end subroutine deriveVelIncomp_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 deriveAuxIncomp_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 ) ! -------------------------------------------------------------------- ! 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 auxField(elemOff+dens_pos) = sum( pdf ) ! velocity auxField(elemOff+vel_pos(1)) = sum( pdf * stencil%cxDirRK(1,:) ) * rho0Inv auxField(elemOff+vel_pos(2)) = sum( pdf * stencil%cxDirRK(2,:) ) * rho0Inv auxField(elemOff+vel_pos(3)) = sum( pdf * stencil%cxDirRK(3,:) ) * rho0Inv end do end subroutine deriveAuxIncomp_fromState ! **************************************************************************** ! ! **************************************************************************** ! !> This routine computes equilibirium from state array !! Here it is assumed that rho0 = 1.0 !! !! 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 deriveEqIncomp_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) :: f( layout%fStencil%QQ ) real(kind=rk) :: rho, vel(3), usq, ucx ! -------------------------------------------------------------------- ! QQ = layout%fStencil%QQ do iElem = 1, nElems do iDir = 1, QQ f(iDir) = state( iDir+(iElem-1)*varSys%nScalars ) end do rho = sum( f ) vel(1) = sum( f * layout%fStencil%cxDirRK(1,:) ) vel(2) = sum( f * layout%fStencil%cxDirRK(2,:) ) vel(3) = sum( f * layout%fStencil%cxDirRK(3,:) ) 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 + rho0 * ( cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) ) enddo end do end subroutine deriveEqIncomp_FromState ! **************************************************************************** ! ! **************************************************************************** ! !> Calculate the equlibrium of a given element number with the given input !! state vector for incompressible model !! \( \rho0 = 1.0 \) !! \( f_eq = w_\alpha*(\rho + \rho0( (u . c_i) + (u . c_i)^2 - u^2 ) ) \) !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine mus_deriveEquilIncomp(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 :: 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 ) dens = sum(tmpPDF) vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:)) vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) do iDir = 1, scheme%layout%fStencil%QQ feq(idir) = scheme%layout%weight( idir ) & & * ( dens + rho0 & & * ( 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_deriveEquilIncomp ! **************************************************************************** ! ! **************************************************************************** ! !> Calculate the Non-Equlibrium !! recursive subroutine mus_deriveNonEquilIncomp(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,:)) vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) ! computes equilibrium do iDir = 1, scheme%layout%fStencil%QQ feq(idir) = scheme%layout%weight( idir ) & & * ( dens + rho0 & & * ( 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_deriveNonEquilIncomp ! **************************************************************************** ! ! **************************************************************************** ! !> author: Jiaxing Qi !! Calculate the strain rate ( or rate of strain, or rate of deformation) !! !! The formula is: !! \[ !! \tau_{\alpha \beta}= !! -\frac{3\omega}{2{\rho}_0} \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_deriveStrainRateIncomp( 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 integer :: nCompsPDF, iVal, pdfPos, iDir real(kind=rk) :: dens, vel(3), omega real(kind=rk), allocatable :: nonEq(:) real(kind=rk), allocatable :: tmpPDF(:) real(kind=rk), allocatable :: fEq(:) real(kind=rk), allocatable :: tau(:) integer :: nScalars ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme pdfPos = fun%input_varPos(1) nCompsPDF = varSys%method%val( pdfPos )%nComponents nScalars = varSys%nScalars allocate( tmpPDF( nCompsPDF ) ) allocate( fEq( nCompsPDF ) ) allocate( nonEq( 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,:)) vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) ! computes equilibrium do iDir = 1, scheme%layout%fStencil%QQ feq(idir) = scheme%layout%weight( idir ) & & * ( dens + rho0 & & * ( 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( cxcx = scheme%layout%fStencil%cxcx(:,:), & & f = nonEq(:), & & QQ = scheme%layout%fStencil%QQ ) res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = & & tau(:) * (-1.5_rk) * omega / rho0 end do ! iVal deallocate( tmpPDF ) deallocate( fEq ) deallocate( nonEq ) deallocate( tau ) end subroutine mus_deriveStrainRateIncomp ! ****************************************************************************** ! end module mus_derQuanIncomp_module ! **************************************************************************** !