! Copyright (c) 2017, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! **************************************************************************** ! !> author: Kannan Masilamani, Verena Krupp !! This module contains routine to retrieve state variables for getElement, !! getPoint, setupIndices and getValOfIndex !! ! 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) 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_stateVar_module use, intrinsic :: iso_c_binding, only: c_ptr, c_f_pointer, c_loc use env_module, only: rk, long_k use tem_logging_module, only: logUnit use tem_debug_module, only: dbgUnit use tem_aux_module, only: tem_abort use tem_varSys_module, only: tem_varSys_type, & & tem_varSys_op_type use tem_grow_array_module, only: grw_realArray_type, & & grw_intArray_type, & & append, init, truncate use tem_pointData_module, only: tem_pointData_list_type, & & init, append, truncate use treelmesh_module, only: treelmesh_type use tem_geometry_module, only: tem_CoordOfReal, & & tem_PosofId, tem_BaryOfId use tem_time_module, only: tem_time_type use tem_topology_module, only: tem_IdOfCoord use tem_topology_module, only: tem_levelOf use tem_element_module, only: eT_minRelevant, eT_maxRelevant, & & eT_fluid ! include musubi modules use mus_scheme_type_module, only: mus_scheme_type use mus_varSys_module, only: mus_varSys_data_type use mus_connectivity_module, only: mus_getSrcElemPosForIntp implicit none private public :: mus_access_state_forElement public :: mus_access_stateFetch_forElement public :: mus_access_stateFetch_now_forElement public :: mus_stateVar_forPoint public :: mus_accessVar_setupIndices public :: mus_stateVar_fromIndex public :: mus_stateVar_Fetch_fromIndex public :: mus_stateVar_Fetch_now_fromIndex contains ! ************************************************************************* ! !> Return the solver state variable for a given set of elements !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine mus_access_state_forElement(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, iComp, iLevel type(mus_varSys_data_type), pointer :: fPtr integer :: nSize, iFld, QQ, nScalars ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) ! myPos can be only used to access state variable since state variables are ! added to varSys first. iFld = fun%myPos QQ = fPtr%solverData%scheme%layout%fStencil%QQ nScalars = varSys%nSCalars ! res is always AOS layout res = 0.0_rk associate( state => fPtr%solverData%scheme%state, & & pdf => fPtr%solverData%scheme%pdf, & & levelPointer => fPtr%solverData%geometry%levelPointer ) 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) ) ) nSize = pdf( iLevel )%nSize do iComp = 1, fun%nComponents res( (iElem-1)*fun%nComponents+iComp ) = & & state( iLevel )%val( & ! position of this state variable in the state array & ( statepos-1)* nscalars+icomp+( ifld-1)* qq, & & pdf( iLevel )%nNext ) end do !iComp end do !iElem end associate end subroutine mus_access_state_forElement ! ************************************************************************* ! ! ************************************************************************* ! !> Return the solver state variable for a given set of elements by using !! FETCH macro for nNow !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine mus_access_stateFetch_now_forElement( & & 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, iComp, iLevel type(mus_varSys_data_type), pointer :: fPtr integer :: nSize, iFld, QQ, nScalars ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) QQ = fPtr%solverData%scheme%layout%fStencil%QQ nScalars = varSys%nSCalars iFld = fun%input_varPos(1) ! res is always AOS layout res = 0.0_rk associate( state => fPtr%solverData%scheme%state, & & pdf => fPtr%solverData%scheme%pdf, & & levelPointer => fPtr%solverData%geometry%levelPointer ) 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) ) ) nSize = pdf( iLevel )%nSize do iComp = 1, fun%nComponents res( (iElem-1)*fun%nComponents+iComp ) = & & state( iLevel )%val( & ! position of this state variable in the state array & pdf(ilevel)%neigh((icomp-1)* nsize+ statepos)+( ifld-1)* qq+ nscalars*0, & & pdf( iLevel )%nNow ) end do !iComp end do !iElem end associate end subroutine mus_access_stateFetch_now_forElement ! ************************************************************************* ! ! ************************************************************************* ! !> Return the solver state variable for a given set of elements by using !! FETCH macro for nNext !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_element]]. !! recursive subroutine mus_access_stateFetch_forElement(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, iComp, iLevel type(mus_varSys_data_type), pointer :: fPtr integer :: nSize, iFld, QQ, nScalars ! -------------------------------------------------------------------- ! call C_F_POINTER( fun%method_Data, fPtr ) QQ = fPtr%solverData%scheme%layout%fStencil%QQ nScalars = varSys%nSCalars iFld = fun%input_varPos(1) ! res is always AOS layout res = 0.0_rk associate( state => fPtr%solverData%scheme%state, & & pdf => fPtr%solverData%scheme%pdf, & & levelPointer => fPtr%solverData%geometry%levelPointer ) 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) ) ) nSize = pdf( iLevel )%nSize do iComp = 1, fun%nComponents res( (iElem-1)*fun%nComponents+iComp ) = & & state( iLevel )%val( & ! position of this state variable in the state array & pdf(ilevel)%neigh((icomp-1)* nsize+ statepos)+( ifld-1)* qq+ nscalars*0, & & pdf( iLevel )%nNext ) end do !iComp end do !iElem end associate end subroutine mus_access_stateFetch_forElement ! ************************************************************************* ! ! ************************************************************************* ! !> State variable for a given set of points using linear interpolation. !! !! The interface has to comply to the abstract interface !! [[tem_varSys_module:tem_varSys_proc_point]]. !! recursive subroutine mus_stateVar_forPoint(fun, varsys, point, time, tree, & & nPnts, 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 !> Three-dimensional coordinates at which the variable should be !! evaluated. Only useful for variables provided as space-time functions. real(kind=rk), intent(in) :: point(:,:) !> 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) :: nPnts !> 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 type(mus_scheme_type), pointer :: scheme integer :: statePos, iPnt, iLevel, elemPos integer :: coord(4) integer, allocatable :: srcElemPos(:) real(kind=rk), allocatable :: weights(:), srcRes(:), pntVal(:) integer :: iSrc, iComp, nSrcElems, nSize, iFld, QQ, nScalars ! -------------------------------------------------------------------- ! !write(dbgUnit(1),*) 'Derive for point :'//trim(varSys%varname%val(fun%myPos)) call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme ! myPos can be only used to access state variable since state variables are ! added to varSys first. iFld = fun%myPos nScalars = varSys%nScalars QQ = scheme%layout%fStencil%QQ allocate(srcElemPos(QQ)) allocate(weights(QQ)) allocate(srcRes(QQ*fun%nComponents)) allocate(pntVal(fun%nComponents)) res = 0.0_rk do iPnt = 1, nPnts weights = 0.0_rk srcRes = 0.0_rk coord = tem_CoordOfReal(tree, point(iPnt,:), tree%global%maxLevel ) ! returns position of existing element in tree which contains point ! 0 if no corresponding node is found, ! or the negative of the found ID, if it is a virtual node. elemPos = abs(tem_PosofId( tem_IdOfCoord(coord), tree%treeID )) ! skip this point if no corresponding element is found if (elemPos == 0) cycle ! level of existing element iLevel = tem_levelOf( tree%treeID( elemPos ) ) ! position of element in levelDesc total list statePos = fPtr%solverData%geometry%levelPointer( elemPos ) ! get position of source element position in state array for interpolation ! using neighbor connectivity array. ! Also calculate weights for interpolation using distance between the ! point and source element barycenter call mus_getSrcElemPosForIntp( & & srcElemPos = srcElemPos, & & weights = weights, & & nSrcElems = nSrcElems, & & point = point(iPnt,:), & & elemPos = statePos, & & neigh = fPtr%solverData%scheme%pdf(iLevel)%neigh, & & baryOfTotal = scheme%levelDesc(iLevel)%baryOfTotal, & & nElems = scheme%pdf(iLevel)%nSize, & & nSolve = scheme%pdf(iLevel)%nElems_computed, & & stencil = scheme%layout%fStencil, & & excludeHalo = .true., & & nScalars = varSys%nScalars ) ! get source element values do iSrc = 1, nSrcElems ! position of element in levelDesc total list statePos = srcElemPos(iSrc) nSize = scheme%pdf( iLevel )%nSize do iComp = 1, fun%nComponents srcRes( (iSrc-1)*fun%nComponents + iComp ) & & = scheme%state( iLevel )%val( & & ( statepos-1)* nscalars+icomp+( ifld-1)* qq,& & scheme%pdf( iLevel )%nNext ) end do !iComp end do !iSrc ! Linear interpolation res = sum(weight_i*phi_i) pntVal = 0.0_rk do iSrc = 1, nSrcElems pntVal(:) = pntVal(:) + weights(iSrc) & & * srcRes( (iSrc-1)*fun%nComponents+1 : iSrc*fun%nComponents ) end do res( (iPnt-1)*fun%nComponents+1 : iPnt*fun%nComponents ) = pntVal end do !iPnt end subroutine mus_stateVar_forPoint ! ************************************************************************* ! ! ************************************************************************* ! !> This routine takes points coordinates, stores them in the method_data and !! return indices where points are located in the growing array of points or !! values ( sometimes we do not need to store the points ) !! It is need to setup points for every variable. Points will be provided by !! boundaries or sources depends on what uses the variable. This points do not !! change with time . This indices will be stored in corresponding boundary !! or source to evaluate a value on that point later using !! tem_varSys_proc_getValOfIndex. subroutine mus_accessVar_setupIndices( fun, varSys, point, offset_bit, & & iLevel, tree, nPnts, idx ) ! -------------------------------------------------------------------------- ! !> Description of the method to obtain the variables, for this routine !! we need the location where to store the points. class(tem_varSys_op_type), intent(in) :: fun !> The variable system to obtain the variable from. type(tem_varSys_type), intent(in) :: varSys !> arrays of points for which the indices are returned real(kind=rk), intent(in) :: point(:,:) !> Offset bit encoded as character for every point. !! If not present default is to center i.e offset_bit = achar(1+4+16) character, optional, intent(in) :: offset_bit(:) !> the point data need to be loaded levelwise, we need the current level integer, intent(in) :: iLevel !> global treelm mesh info type(treelmesh_type), intent(in) :: tree !> number of points integer, intent(in) :: nPnts integer, intent(out) :: idx(:) ! -------------------------------------------------------------------------- ! integer :: iPoint integer :: coord(4) integer :: elemPos, loc_level, statePos type(mus_varSys_data_type), pointer :: fPtr type(mus_scheme_type), pointer :: scheme real(kind=rk), allocatable :: weights(:) integer, allocatable :: srcElemPos(:) integer :: nSrcElems integer(kind=long_k) :: treeID ! -------------------------------------------------------------------------- ! ! write(dbgUnit(4),*) 'setup indices for the points of variable ', & ! & trim(varSys%varname%val(fun%myPos)) call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme allocate(srcElemPos(scheme%layout%fStencil%QQ)) allocate(weights(scheme%layout%fStencil%QQ)) idx = 0 ! store elemPos of each point in global treeID list do iPoint = 1, nPnts weights = 0.0_rk ! Get the treeID of the element and then get the element position ! in the levelwise array coord = tem_CoordOfReal( mesh = tree, & & point = point(iPoint,:), & & level = tree%global%maxLevel ) treeID = tem_IdOfCoord(coord) ! returns position of existing element in tree which contains point ! 0 if no corresponding node is found, ! or the negative of the found ID, if it is a virtual node. elemPos = abs(tem_PosofId( treeID, tree%treeID )) ! write(dbgUnit(1),*) 'iPnt ', iPoint, ' elemPos ', elemPos, & ! & 'treeID ', treeID, 'point ', point(iPoint,:) ! skip this point if no corresponding element is found if (elemPos == 0) cycle ! store point data only if element position is found in tree ! level of existing element loc_level = tem_levelOf( tree%treeID( elemPos ) ) ! append grw array for this point's element position in tree call append(me = fPtr%pointData%pntLvl(iLevel)%elemPos, & & val = elemPos ) ! Actual level on which this point was found call append(me = fPtr%pointData%pntLvl(iLevel)%pntLevel, & & val = loc_level ) ! Append points to growing array of points call append(me = fPtr%pointData%pntLvl(iLevel)%grwPnt, & & val = point(iPoint,:) ) ! position of element in levelDesc total list statePos = fPtr%solverData%geometry%levelPointer( elemPos ) ! get position of source element position in state array for interpolation ! using neighbor connectivity array. ! Also calculate weights for interpolation using distance between the ! point and source element barycenter ! Excluded halo elements as they require communication of all links call mus_getSrcElemPosForIntp( & & srcElemPos = srcElemPos, & & weights = weights, & & nSrcElems = nSrcElems, & & point = point(iPoint,:), & & elemPos = statePos, & & neigh = fPtr%solverData%scheme%pdf(loc_level)%neigh, & & baryOfTotal = scheme%levelDesc(loc_level)%baryOfTotal, & & nElems = scheme%pdf(loc_level)%nSize, & & nSolve = scheme%pdf(loc_level)%nElems_computed, & & stencil = scheme%layout%fStencil, & & excludeHalo = .true., & & nScalars = varSys%nScalars ) ! First position of source elements call append(me = fPtr%pointData%pntLvl(iLevel)%srcElem%first, & & val = fPtr%pointData%pntLvl(iLevel)%srcElem%elemPos & & %nVals + 1 ) ! Append all src elemPos call append(me = fPtr%pointData%pntLvl(iLevel)%srcElem%elemPos, & & val = srcElemPos(1:nSrcElems) ) ! last position of source elements call append(me = fPtr%pointData%pntLvl(iLevel)%srcElem%last, & & val = fPtr%pointData%pntLvl(iLevel)%srcElem%elemPos & & %nVals ) ! weights for srcElements call append(me = fPtr%pointData%pntLvl(iLevel)%srcElem%weight, & & val = weights(1:nSrcElems) ) ! store the index, idx(iPoint)= fPtr%pointData%pntLvl(iLevel)%elemPos%nVals end do ! write(dbgUnit(1),*) 'varName ', trim(varSys%varName%val(fun%myPos)) ! write(dbgUnit(1),*) 'idx ', idx fPtr%pointData%pntLvl(iLevel)%nPnts = fPtr%pointData%pntLvl(iLevel)%grwPnt & & %coordX%nVals call truncate(fPtr%pointData%pntLvl(iLevel)%grwPnt) call truncate(fPtr%pointData%pntLvl(iLevel)%elemPos) call truncate(fPtr%pointData%pntLvl(iLevel)%srcElem%first) call truncate(fPtr%pointData%pntLvl(iLevel)%srcElem%last) call truncate(fPtr%pointData%pntLvl(iLevel)%srcElem%elemPos) call truncate(fPtr%pointData%pntLvl(iLevel)%srcElem%weight) end subroutine mus_accessVar_setupIndices ! ************************************************************************* ! ! ************************************************************************* ! !> Routine to get the actual value for a given array of indices. !! The indices belong to the grwarray of points storing levelwise in !! Pointdata%pntLvl(iLevel). !! Hence this routines takes the indeices as input, can refer to the pointData !! and evaluate the variable and returns the values subroutine mus_stateVar_fromIndex( fun, varSys, time, iLevel, idx, & & idxLen, nVals, res ) ! -------------------------------------------------------------------------- ! !> Description of the method to obtain the 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: n 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 type(mus_scheme_type), pointer :: scheme integer :: iComp, iSrc, iVal, elemPos, first, last, nSize, loc_level integer :: iSrcElems real(kind=rk), allocatable :: srcRes(:), pntVal(:) real(kind=rk) :: weight integer :: iFld, QQ, nScalars ! -------------------------------------------------------------------------- ! ! write(dbgUnit(4),*) 'get the values of indices for variable ', & ! & trim(varSys%varname%val(fun%myPos)) call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme ! myPos can be only used to access state variable since state variables are ! added to varSys first. iFld = fun%myPos nScalars = varSys%nScalars QQ = scheme%layout%fStencil%QQ allocate(srcRes(QQ*fun%nComponents)) allocate(pntVal(fun%nComponents)) res = 0.0_rk ! distinguish if we have an array of index or we have contingous memory ! access where index are always first entries! if (present(idxLen)) then call tem_abort('Error: idxLen is not supported in get_valOfIndex for ' & & //'state variable') else ! Get the state value at the specific point do iVal = 1, nVals if (idx(iVal)>0) then ! elemPos in tree ! elemPos = fPtr%pointData%pntLvl(iLevel)%elemPos%val(idx(iVal)) first = fPtr%pointData%pntLvl(iLevel)%srcElem%first%val(idx(iVal)) last = fPtr%pointData%pntLvl(iLevel)%srcElem%last%val(idx(iVal)) loc_level = fPtr%pointData%pntLvl(iLevel)%pntLevel%val(idx(iVal)) nSize = scheme%pdf( loc_level )%nSize ! get pdf's of source elements srcRes = 0.0_rk iSrcElems = 0 do iSrc = first, last iSrcElems = iSrcElems + 1 ! position of element in levelDesc total list elemPos = fPtr%pointData%pntLvl(iLevel)%srcElem%elemPos%val(iSrc) do iComp = 1, fun%nComponents srcRes( (iSrcElems-1)*fun%nComponents + iComp ) & & = scheme%state( loc_level )%val( & & ( elempos-1)* nscalars+icomp+( ifld-1)* qq, & & scheme%pdf( loc_level )%nNext ) end do !iComp end do !iSrc ! Linear interpolation res = sum(weight_i*phi_i) pntVal = 0.0_rk iSrcElems = 0 do iSrc = first, last weight = fPtr%pointData%pntLvl(iLevel)%srcElem%weight%val(iSrc) iSrcElems = iSrcElems + 1 pntVal(:) = pntVal(:) + weight & & * srcRes( (iSrcElems-1)*fun%nComponents+1 & & : iSrcElems*fun%nComponents ) end do ! get the state value for each component of this res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = pntVal end if !idx>0 end do !iVal end if ! idx_len end subroutine mus_stateVar_fromIndex ! ************************************************************************* ! ! ************************************************************************* ! !> Routine to get the actual value for a given array of indices. !! The indices belong to the grwarray of points storing levelwise in !! Pointdata%pntLvl(iLevel). !! Hence this routines takes the indeices as input, can refer to the pointData !! and evaluate the variable and returns the values subroutine mus_stateVar_Fetch_fromIndex( fun, varSys, time, iLevel, idx, & & idxLen, nVals, res ) ! -------------------------------------------------------------------------- ! !> Description of the method to obtain the 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: n 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 type(mus_scheme_type), pointer :: scheme integer :: iComp, iSrc, iVal, elemPos, first, last, nSize, loc_level integer :: nSrcElems, iFld, QQ, nScalars real(kind=rk), allocatable :: srcRes(:), pntVal(:) real(kind=rk) :: weight ! -------------------------------------------------------------------------- ! ! write(dbgUnit(4),*) 'get the values of indices for variable ', & ! & trim(varSys%varname%val(fun%myPos)) call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme QQ = scheme%layout%fStencil%QQ nScalars = varSys%nSCalars allocate(srcRes(QQ*fun%nComponents)) allocate(pntVal(fun%nComponents)) iFld = fun%input_varPos(1) res = 0.0_rk ! distinguish if we have an array of index or we have contingous memory ! access where index are always first entries! if (present(idxLen)) then call tem_abort('Error: idxLen is not supported in get_valOfIndex for ' & & //'state variable') else ! Get the state value at the specific point do iVal = 1, nVals if (idx(iVal)>0) then ! elemPos in tree ! elemPos = fPtr%pointData%pntLvl(iLevel)%elemPos%val(idx(iVal)) first = fPtr%pointData%pntLvl(iLevel)%srcElem%first%val(idx(iVal)) last = fPtr%pointData%pntLvl(iLevel)%srcElem%last%val(idx(iVal)) loc_level = fPtr%pointData%pntLvl(iLevel)%pntLevel%val(idx(iVal)) nSize = scheme%pdf( loc_level )%nSize ! get pdf's of source elements srcRes = 0.0_rk nSrcElems = 0 do iSrc = first, last nSrcElems = nSrcElems + 1 ! position of element in levelDesc total list elemPos = fPtr%pointData%pntLvl(iLevel)%srcElem%elemPos%val(iSrc) do iComp = 1, fun%nComponents srcRes( (nSrcElems-1)*fun%nComponents + iComp ) & & = scheme%state( loc_level )%val( & & scheme%pdf(loc_level)%neigh((icomp-1)*nsize+elempos)+(ifld-1)*qq+nscalars*0, & & scheme%pdf( loc_level )%nNext ) end do !iComp end do !iSrc ! Linear interpolation res = sum(weight_i*phi_i) pntVal = 0.0_rk nSrcElems = 0 do iSrc = first, last weight = fPtr%pointData%pntLvl(iLevel)%srcElem%weight%val(iSrc) nSrcElems = nSrcElems + 1 pntVal(:) = pntVal(:) + weight & & * srcRes( (nSrcElems-1)*fun%nComponents+1 & & : nSrcElems*fun%nComponents ) end do ! get the state value for each component of this ! get the state value for each component of this res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = pntVal end if !idx>0 end do !iVal end if ! idx_len end subroutine mus_stateVar_Fetch_fromIndex ! ************************************************************************* ! ! ************************************************************************* ! !> Routine to get the actual value for a given array of indices. !! The indices belong to the grwarray of points storing levelwise in !! Pointdata%pntLvl(iLevel). !! Hence this routines takes the indeices as input, can refer to the pointData !! and evaluate the variable and returns the values subroutine mus_stateVar_Fetch_now_fromIndex( fun, varSys, time, iLevel, & & idx, idxLen, nVals, res ) ! -------------------------------------------------------------------------- ! !> Description of the method to obtain the 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: n 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 type(mus_scheme_type), pointer :: scheme integer :: iComp, iSrc, iVal, elemPos, first, last, nSize, loc_level integer :: nSrcElems, iFld, QQ, nScalars real(kind=rk), allocatable :: srcRes(:), pntVal(:) real(kind=rk) :: weight ! -------------------------------------------------------------------------- ! ! write(dbgUnit(4),*) 'get the values of indices for variable ', & ! & trim(varSys%varname%val(fun%myPos)) call C_F_POINTER( fun%method_Data, fPtr ) scheme => fPtr%solverData%scheme QQ = scheme%layout%fStencil%QQ nScalars = varSys%nSCalars allocate(srcRes(QQ*fun%nComponents)) allocate(pntVal(fun%nComponents)) iFld = fun%input_varPos(1) res = 0.0_rk ! distinguish if we have an array of index or we have contingous memory ! access where index are always first entries! if (present(idxLen)) then call tem_abort('Error: idxLen is not supported in get_valOfIndex for ' & & //'state variable') else ! Get the state value at the specific point do iVal = 1, nVals if (idx(iVal)>0) then ! elemPos in tree ! elemPos = fPtr%pointData%pntLvl(iLevel)%elemPos%val(idx(iVal)) first = fPtr%pointData%pntLvl(iLevel)%srcElem%first%val(idx(iVal)) last = fPtr%pointData%pntLvl(iLevel)%srcElem%last%val(idx(iVal)) loc_level = fPtr%pointData%pntLvl(iLevel)%pntLevel%val(idx(iVal)) nSize = scheme%pdf( loc_level )%nSize ! get pdf's of source elements srcRes = 0.0_rk nSrcElems = 0 do iSrc = first, last nSrcElems = nSrcElems + 1 ! position of element in levelDesc total list elemPos = fPtr%pointData%pntLvl(iLevel)%srcElem%elemPos%val(iSrc) do iComp = 1, fun%nComponents srcRes( (nSrcElems-1)*fun%nComponents + iComp ) & & = scheme%state( loc_level )%val( & & scheme%pdf(loc_level)%neigh((icomp-1)*nsize+elempos)+(ifld-1)*qq+nscalars*0, & & scheme%pdf( loc_level )%nNow ) end do !iComp end do !iSrc ! Linear interpolation res = sum(weight_i*phi_i) pntVal = 0.0_rk nSrcElems = 0 do iSrc = first, last weight = fPtr%pointData%pntLvl(iLevel)%srcElem%weight%val(iSrc) nSrcElems = nSrcElems + 1 pntVal(:) = pntVal(:) + weight & & * srcRes( (nSrcElems-1)*fun%nComponents+1 & & : nSrcElems*fun%nComponents ) end do ! get the state value for each component of this ! get the state value for each component of this res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = pntVal end if !idx>0 end do !iVal end if ! idx_len end subroutine mus_stateVar_Fetch_now_fromIndex ! ************************************************************************* ! end module mus_stateVar_module ! **************************************************************************** !