! Copyright (c) 2016-2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2019 Seyfettin Bilgi <seyfettin.bilgi@student.uni-siegen.de> ! Copyright (c) 2021 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.de> ! Copyright (c) 2022 Kannan Masilamani <kannan.masilamani@dlr.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ***************************************************************************** ! !> author: Kannan Masilamani !! Module containing subroutines for building MUSUBI specific source !! variables !! module mus_source_var_module use, intrinsic :: iso_c_binding, only: c_ptr, c_f_pointer ! include treelm modules use mpi use env_module, only: rk, rk_mpi use tem_param_module, only: rho0, rho0Inv, cs2, cs2inv, cs4inv use tem_aux_module, only: tem_abort use tem_varSys_module, only: tem_varSys_type use tem_time_module, only: tem_time_type use treelmesh_module, only: treelmesh_type use tem_topology_module, only: tem_levelOf use tem_float_module, only: operator(.feq.) use tem_logging_module, only: logUnit ! include musubi modules use mus_physics_module, only: mus_convertFac_type use mus_derVarPos_module, only: mus_derVarPos_type use mus_source_type_module, only: mus_source_op_type use mus_varSys_module, only: mus_varSys_data_type implicit none private ! update auxField dependent source variables public :: mus_updateSrcVar_dynSponFld public :: mus_updateSrcVar_turbChanForce contains ! ************************************************************************** ! !> Compute density and velocity in sponge layer for dynamic sponge subroutine mus_updateSrcVar_dynSponFld(fun, auxField, iLevel, varSys, & & phyConvFac, derVarPos) ! ------------------------------------------------------------------------ ! !> Description of method to update source class(mus_source_op_type), intent(inout) :: fun !> input auxField array on current level real(kind=rk), intent(in) :: auxField(:) !> current level integer, intent(in) :: iLevel !> variable system definition type(tem_varSys_type), intent(in) :: varSys !> Physics conversion factor on current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------------- integer :: dens_pos, vel_pos(3) integer :: iElem, nElems, elemOff real(kind=rk) :: dens, vel(3) real(kind=rk) :: smoothFac real(kind=rk) :: inv_rho_phy, inv_vel_phy real(kind=rk) :: dens_ref, vel_ref(3) ! -------------------------------------------------------------------------- ! position of density and velocity field in auxField dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1) vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) inv_rho_phy = 1.0_rk / phyConvFac%press * cs2inv inv_vel_phy = 1.0_rk / phyConvFac%vel ! Number of elements to apply source terms nElems = fun%elemLvl(iLevel)%nElems ! Calculate time average value ! Set initial density and velocity if no initialized associate( dynAvg => fun%elemLvl(iLevel)%dynAvg, & & posInTotal => fun%elemLvl(iLevel)%posInTotal ) ! Do Exponential moving average over nRecord window ! P_n = P_{n-1} + smoothFac * (P - P_{n-1}) smoothFac = fun%absLayer%smoothFac ! DO time average for density if (fun%absLayer%config%isPressDyn) then if (dynAvg%isInitDens) then dynAvg%isInitDens = .false. ! Initialize dynAvg%dens with initialize density do iElem = 1, nElems ! element offset elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars dynAvg%dens(iElem) = auxField(elemOff+dens_pos) end do else do iElem = 1, nElems ! element offset elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars ! Local density and velocity dens = auxField(elemOff+dens_pos) ! New average dynAvg%dens(iElem) = dynAvg%dens(iElem) & & + smoothFac * (dens - dynAvg%dens(iElem) ) end do end if else ! target pressure in lattice unit dens_ref = fun%absLayer%config%target_pressure * inv_rho_phy dynAvg%dens(:) = dens_ref end if ! DO time average for velocity if (fun%absLayer%config%isVelDyn) then if (dynAvg%isInitVel) then dynAvg%isInitVel = .false. ! Initialize dynAvg%vel with initialize velocity do iElem = 1, nElems ! element offset elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars dynAvg%velX(iElem) = auxField(elemOff+vel_pos(1)) dynAvg%velY(iElem) = auxField(elemOff+vel_pos(2)) dynAvg%velZ(iElem) = auxField(elemOff+vel_pos(3)) end do else do iElem = 1, nElems ! element offset elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars ! Local velocity vel(1) = auxField(elemOff+vel_pos(1)) vel(2) = auxField(elemOff+vel_pos(2)) vel(3) = auxField(elemOff+vel_pos(3)) ! New average dynAvg%velX(iElem) = dynAvg%velX(iElem) & & + smoothFac * (vel(1) - dynAvg%velX(iElem) ) dynAvg%velY(iElem) = dynAvg%velY(iElem) & & + smoothFac * (vel(2) - dynAvg%velY(iElem) ) dynAvg%velZ(iElem) = dynAvg%velZ(iElem) & & + smoothFac * (vel(3) - dynAvg%velZ(iElem) ) end do end if else ! target velocity in lattice unit vel_ref(1:3) = fun%absLayer%config%target_velocity(1:3) * inv_vel_phy dynAvg%velX(:) = vel_ref(1) dynAvg%velY(:) = vel_ref(2) dynAvg%velZ(:) = vel_ref(3) end if end associate end subroutine mus_updateSrcVar_dynSponFld ! ************************************************************************** ! ! ************************************************************************** ! !> Compute dynamic force term using auxField for turbulent channel !! force. subroutine mus_updateSrcVar_turbChanForce(fun, auxField, iLevel, varSys, & & phyConvFac, derVarPos) ! ------------------------------------------------------------------------ ! !> Description of method to update source class(mus_source_op_type), intent(inout) :: fun !> input auxField array on current level real(kind=rk), intent(in) :: auxField(:) !> current level integer, intent(in) :: iLevel !> variable system definition type(tem_varSys_type), intent(in) :: varSys !> Physics conversion factor on current level type(mus_convertFac_type), intent(in) :: phyConvFac !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos(:) ! -------------------------------------------------------------------------- integer :: vel_pos(3) integer :: iElem, nElemsGlobal_uTau, posInTotal, elemOff type(mus_varSys_data_type), pointer :: fPtr integer :: iErr, iBC, iLvlLoc real(kind=rk) :: vel_bulk, vel_tau real(kind=rk) :: avgVelBulk, avgVelTau ! index 1 for vel_tau and 2 for vel_bulk to get global avg with one mpi call real(kind=rk) :: avgVel(2), avgVelGlobal(2) real(kind=rk) :: forceDyn, gradU(3,3,1) logical :: velTauFromBnd ! -------------------------------------------------------------------------- ! position of velocity field in auxField vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3) ! First variable is always pdf so just use to access method_data call C_F_POINTER( varSys%method%val(1)%method_Data, fPtr ) ! Calculate average bulk velocityX over defined shape. associate( levelPointer => fPtr%solverData%geometry%levelPointer, & & auxFieldTot => fPtr%solverData%scheme%auxField, & & gradDataTot => fPtr%solverData%scheme%gradData, & & viscKineTot => fPtr%solverData%scheme%field(1)%fieldProp & & %fluid%viscKine%dataOnLvl, & & map2Global_umean => fun%turbChanForce%subTree_umean%map2global, & & map2Global_utau => fun%turbChanForce%subTree_utau%map2global, & & physics => fPtr%solverData%physics, & & tree => fPtr%solverData%geometry%tree, & & grad => fPtr%solverData%scheme%Grad ) ! Do average only at every multilevel cycle i.e. when iLevel=minLevel if (iLevel == tree%global%minLevel) then ! Initialize variables vel_bulk = 0.0_rk vel_tau = 0.0_rk velTauFromBnd = .false. nElemsGlobal_uTau = 0 !> If wall model BC is applied than compute the friction velocity ! that is obtained directly from the wall model and perform the spatial ! averaging. do iBC = 1, fPtr%solverData%geometry%boundary%nBCtypes select case(trim(fPtr%solverData%scheme%field(1)%BC(iBC)%BC_kind)) case('turbulent_wall', 'turbulent_wall_noneq_expol', & & 'turbulent_wall_eq') if (allocated(fPtr%solverData%scheme%field(1)%BC(iBC) & & %turbWallFunc%dataOnLvl)) then velTauFromBnd = .true. do iLvlLoc = tree%global%minLevel, tree%global%maxLevel ! vel_tau in turbWallFunc is in lattice unit so convert it vel_tau = vel_tau + sum( fPtr%solverData%scheme%field(1) & & %BC(iBC)%turbWallFunc & & %dataOnLvl(iLvlLoc)%velTau(:) ) & & * physics%fac(iLvlLoc)%vel end do nElemsGlobal_uTau = nElemsGlobal_uTau & & + fPtr%solverData%scheme%globBC(iBC)%nElems_Total end if end select end do !> If wall model BC is not used than compute the friction velocity from !! the single sided finite difference's and perform the spatial averaging. !! Friction velocity is computed only on elements intersected by !! shape_utau defined in musubi.lua. if (.not. velTauFromBnd) then nElemsGlobal_uTau = fun%turbChanForce%nElemsGlobal_utau select case(fun%turbChanForce%flow_direction) case (1) ! x-direction do iElem = 1, fun%turbChanForce%subTree_utau%nElems ! map2Global refers to position in treeid list ! levelPointer refers to position in level wise total list posInTotal = levelPointer( map2Global_utau(iElem) ) iLvlLoc = tem_levelOf( tree%treeID( map2Global_utau(iElem) ) ) ! elemOffset refers to the previous number of elements processed gradU(:,:,:) = grad%U_ptr( & & auxField = auxFieldTot(iLvlLoc)%val(:), & & gradData = gradDataTot(iLvlLoc), & & velPos = vel_Pos, & & nAuxScalars = varSys%nAuxScalars, & & nDims = 3, & & nSolve = 1, & & elemOffset = posInTotal - 1 ) ! compute uTau using the dudy component of velocity gradient tensor vel_tau = vel_tau + sqrt( viscKineTot(iLvlLoc)%val(posInTotal) & & * abs(gradU(1,2,1)) ) & & * physics%fac(iLvlLoc)%vel !write(dbgunit(1), *) 'posInTotal: ', posInTotal, 'gradU: ', gradU(1,2,1) end do case (2) ! y-direction do iElem = 1, fun%turbChanForce%subTree_utau%nElems ! map2Global refers to position in treeid list ! levelPointer refers to position in level wise total list posInTotal = levelPointer( map2Global_utau(iElem) ) iLvlLoc = tem_levelOf( tree%treeID( map2Global_utau(iElem) ) ) gradU(:,:,:) = grad%U_ptr( & & auxField = auxFieldTot(iLvlLoc)%val(:), & & gradData = gradDataTot(iLvlLoc), & & velPos = vel_Pos, & & nAuxScalars = varSys%nAuxScalars, & & nDims = 3, & & nSolve = 1, & & elemOffset = posInTotal - 1 ) ! compute uTau using the dvdx component of velocity gradient tensor vel_tau = vel_tau + sqrt( viscKineTot(iLvlLoc)%val(posInTotal) & & * abs(gradU(2,1,1)) ) & & * physics%fac(iLvlLoc)%vel end do case (3) ! z-direction do iElem = 1, fun%turbChanForce%subTree_utau%nElems ! map2Global refers to position in treeid list ! levelPointer refers to position in level wise total list posInTotal = levelPointer( map2Global_utau(iElem) ) iLvlLoc = tem_levelOf( tree%treeID( map2Global_utau(iElem) ) ) gradU(:,:,:) = grad%U_ptr( & & auxField = auxFieldTot(iLvlLoc)%val(:), & & gradData = gradDataTot(iLvlLoc), & & velPos = vel_Pos, & & nAuxScalars = varSys%nAuxScalars, & & nDims = 3, & & nSolve = 1, & & elemOffset = posInTotal - 1 ) ! compute uTau using the dwdy component of velocity gradient tensor vel_tau = vel_tau + sqrt( viscKineTot(iLvlLoc)%val(posInTotal) & & * abs(gradU(3,2,1)) ) & & * physics%fac(iLvlLoc)%vel end do end select end if !> Bulk mean velocity part of forcing is independent whether a wall ! modelBC is used or not select case(fun%turbChanForce%flow_direction) case (1) ! x-direction do iElem = 1, fun%turbChanForce%subTree_umean%nElems ! map2Global refers to position in treeid list ! levelPointer refers to position in level wise total list posInTotal = levelPointer( map2Global_umean(iElem) ) ! elemoffset for auxField elemoff = (posInTotal-1)*varSys%nAuxScalars ! velocity X in defined shape to compute average iLvlLoc = tem_levelOf( tree%treeID( map2Global_umean(iElem) ) ) vel_bulk = vel_bulk + auxFieldTot(iLvlLoc)%val(elemOff+vel_pos(1)) & & * physics%fac(iLvlLoc)%vel end do case (2) ! y-direction do iElem = 1, fun%turbChanForce%subTree_umean%nElems ! map2Global refers to position in treeid list ! levelPointer refers to position in level wise total list posInTotal = levelPointer( map2Global_umean(iElem) ) ! elemoffset for auxField elemoff = (posInTotal-1)*varSys%nAuxScalars ! velocity X in defined shape to compute average iLvlLoc = tem_levelOf( tree%treeID( map2Global_umean(iElem) ) ) vel_bulk = vel_bulk + auxFieldTot(iLvlLoc)%val(elemOff+vel_pos(2)) & & * physics%fac(iLvlLoc)%vel end do case (3) ! z-direction do iElem = 1, fun%turbChanForce%subTree_umean%nElems ! map2Global refers to position in treeid list ! levelPointer refers to position in level wise total list posInTotal = levelPointer( map2Global_umean(iElem) ) ! elemoffset for auxField elemoff = (posInTotal-1)*varSys%nAuxScalars ! velocity X in defined shape to compute average iLvlLoc = tem_levelOf( tree%treeID( map2Global_umean(iElem) ) ) vel_bulk = vel_bulk + auxFieldTot(iLvlLoc)%val(elemOff+vel_pos(3)) & & * physics%fac(iLvlLoc)%vel end do end select avgVel(1) = vel_tau avgVel(2) = vel_bulk ! Calculate total friction and bulk velocity call mpi_allreduce( avgVel, avgVelGlobal, & & 2, rk_mpi, mpi_sum, tree%global%comm, iErr ) ! compute average friction velocity avgVelTau = avgVelGlobal(1) / nElemsGlobal_uTau ! compute average bulk velocity in physical unit avgVelBulk = avgVelGlobal(2) / fun%turbChanForce%nElemsGlobal_umean ! Dynamic force term for turbulent channel ! F_dyn = (refVelBulk-avgVelBulk) * refVelBulk / refHeight forceDyn = avgVelTau**2 / fun%turbChanForce%refHeight & & + ( fun%turbChanForce%refVelBulk - avgVelBulk ) & & * fun%turbChanForce%refVelBulk / fun%turbChanForce%refHeight fun%turbChanForce%forceDyn = 0.0_rk !write(dbgunit(1), *) 'avgVelBulk: ', avgVelBulk , 'avgVelTau:', avgVelTau !write(dbgunit(1), *) 'vel_tau: ', vel_tau !flush(dbgunit(1)) !write(dbgunit(1), *) 'forceDyn_phy: ', forceDyn select case(fun%turbChanForce%flow_direction) case(1) ! x-direction fun%turbChanForce%forceDyn(1) = forceDyn case(2) ! y-direction fun%turbChanForce%forceDyn(2) = forceDyn case(3) ! z-direction fun%turbChanForce%forceDyn(3) = forceDyn end select end if end associate end subroutine mus_updateSrcVar_turbChanForce ! ************************************************************************** ! end module mus_source_var_module ! **************************************************************************** !