! 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. ! **************************************************************************** ! !> This module contains routines to bnd_force type and routines to initialize !! bndForce array and compute bndForce on all boundary elements !! !! author: Kannan Masilamani ! 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. module mus_bndForce_module ! include treelm modules use env_module, only: rk use tem_debug_module, only: dbgUnit use tem_logging_module, only: logUnit use tem_aux_module, only: tem_abort use tem_bc_prop_module, only: tem_bc_prop_type use tem_construction_module, only: tem_levelDesc_type use tem_varSys_module, only: tem_varSys_type ! include musubi modules use mus_scheme_layout_module, only: mus_scheme_layout_type use mus_scheme_header_module, only: mus_scheme_header_type use mus_scheme_type_module, only: array2D_type use mus_pdf_module, only: pdf_data_type use mus_field_module, only: mus_field_type use mus_bc_header_module, only: boundary_type, glob_boundary_type implicit none private public :: mus_init_bndForce public :: mus_calcBndForce contains ! ************************************************************************** ! !> This routine initialize bndForce array assign function pointer to !! calculate bndForce subroutine mus_init_bndForce(bndForce, bc_prop, schemeHeader, bc) !--------------------------------------------------------------------------- !> Bnd force to allocate real(kind=rk), allocatable, intent(out) :: bndForce(:,:) !> Boundary property type(tem_bc_prop_type), intent(in) :: bc_prop !> scheme header info type(mus_scheme_header_type), intent(in) :: schemeHeader !> global array boundary type type(boundary_type), intent(inout) :: bc(:) !--------------------------------------------------------------------------- integer :: iBnd !--------------------------------------------------------------------------- write(logUnit(1),'(A)') 'Initialize BndForce' select case( trim(schemeHeader%kind) ) case('fluid', 'fluid_incompressible', 'isotherm_acEq') ! deallocated for dynamic load balancing if (allocated(bndForce)) deallocate(bndForce) allocate(bndForce(bc_prop%property%nElems, 3)) bndForce = 0.0_rk do iBnd = 1, bc_prop%nBCtypes select case (trim(bc(iBnd)%BC_kind)) case('wall') bc(iBnd)%calcBndForce => mus_calcBndForce_wall case('wall_libb') bc(iBnd)%calcBndForce => mus_calcBndForce_wall_libb case default bc(iBnd)%calcBndForce => mus_calcBndForce_dummy end select end do case default write(logUnit(1),'(A)') 'WARNING: BndForce calculation is not supported' & & //' for '//trim(schemeHeader%kind) do iBnd = 1, bc_prop%nBCtypes bc(iBnd)%calcBndForce => mus_calcBndForce_dummy end do end select end subroutine mus_init_bndForce ! ************************************************************************** ! ! -------------------------------------------------------------------------- ! !> This routine computes force on boundary elements which are used to !! compute lift and drag coefficient on boundary subroutine mus_calcBndForce( bndForce, posInBndID, nBCs, field, globBC, & & minLevel, maxLevel, state, pdf, levelDesc, & & layout, varSys ) ! -------------------------------------------------------------------- ! !> Boundary force on wall boundary real(kind=rk), intent(inout) :: bndForce(:,:) !> Mapping from global tree%treeid to global boundary%boundaryID integer, intent(in) :: posInBndID(:) !> number of BC integer, intent(in) :: nBCs !> fluid parameters and properties type(mus_field_type), intent(in) :: field(:) !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC(:) !> minlevel and maxLevel integer, intent(in) :: minLevel, maxLevel !> contains global state vector type(pdf_data_type), intent(in) :: pdf(minLevel: maxLevel) !> state arrays fo current iLevel both now and next type(array2D_type), intent(in) :: state(minLevel: maxLevel) !> Level Descriptor type(tem_leveldesc_type), intent(in) :: levelDesc(minLevel: maxLevel) !> scheme layout type type(mus_scheme_layout_type), intent(in) ::layout !> scheme variable system type(tem_varSys_type), intent(in) :: varSys ! -------------------------------------------------------------------- ! integer :: iLevel, iBnd, nFields ! -------------------------------------------------------------------- ! nFields = size( field ) ! Boundary force calculation is valid only for single field schemes ! like fluid and fluid_incompressible if ( nBCs > 0 .and. nFields == 1) then do iLevel = minLevel, maxLevel ! Treat all boundary conditions do iBnd = 1, nBCs ! Calculate force on wall and wall_libb boundaries call field( nFields )%bc( iBnd )%calcBndForce( & & bndForce = bndForce, & & posInBndID = posInBndID, & & globBC = globBC( iBnd ), & & currState = state( iLevel )%val( :, pdf( iLevel )%nNext ), & & levelDesc = levelDesc( iLevel ), & & nSize = pdf( iLevel )%nSize, & & iLevel = iLevel, & & neigh = pdf( iLevel )%neigh, & & layout = layout, & & nScalars = varSys%nScalars ) end do ! iBnd end do ! iLevel end if ! nBCs>0 end subroutine mus_calcBndForce ! -------------------------------------------------------------------------- ! ! ************************************************************************** ! !> Dummy routine for calcBndForce subroutine mus_calcBndForce_dummy( me, bndForce, posInBndID, globBC, & & currState, levelDesc, nSize, iLevel, & & neigh, layout, nScalars ) ! --------------------------------------------------------------------- ! !> field boundary type class( boundary_type ), intent(in) :: me !> bndForce to fill real(kind=rk), intent(inout) :: bndForce(:,:) ! position of boundary element in boundary%bcID integer, intent(in) :: posInBndID(:) !> scheme global boundary type type( glob_boundary_type ), intent(in) :: globBC !> current state array to access post-collision values real(kind=rk), intent(in) :: currState(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> level which invokes boundary integer,intent(in) :: iLevel !> global parameters ! type(mus_param_type),intent(in) :: params integer,intent(in) :: neigh(:) !< connectivity array !> scheme layout type( mus_scheme_layout_type ),intent(in) :: layout !> number of Scalars in the scheme var system integer, intent(in) :: nScalars ! --------------------------------------------------------------------- ! !call tem_abort('Dummy routine for calcBndForce') end subroutine mus_calcBndForce_dummy ! ************************************************************************** ! ! ************************************************************************** ! !> This routine computes bndForce on wall boundary elements subroutine mus_calcBndForce_wall( me, bndForce, posInBndID, globBC, & & currState, levelDesc, nSize, iLevel, & & neigh, layout, nScalars ) ! --------------------------------------------------------------------- ! !> field boundary type class( boundary_type ), intent(in) :: me !> bndForce to fill real(kind=rk), intent(inout) :: bndForce(:,:) ! position of boundary element in boundary%bcID integer, intent(in) :: posInBndID(:) !> scheme global boundary type type( glob_boundary_type ), intent(in) :: globBC !> current state array to access post-collision values real(kind=rk), intent(in) :: currState(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> level which invokes boundary integer,intent(in) :: iLevel !> global parameters ! type(mus_param_type),intent(in) :: params integer,intent(in) :: neigh(:) !< connectivity array !> scheme layout type( mus_scheme_layout_type ),intent(in) :: layout !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !--------------------------------------------------------------------------- integer :: elemPos, iElem, iDir, QQN, QQ integer :: invDir real(kind=rk) :: force(3), fOut !--------------------------------------------------------------------------- QQ = layout%fStencil%QQ QQN = layout%fStencil%QQN do iElem = 1, globBC%nElems_Fluid(iLevel) force = 0.0_rk elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) do iDir = 1, QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then ! direction towards boundary invDir = layout%fStencil%cxDirInv( iDir ) ! Post-collision PDF fOut = currState( & & ( elempos-1)* nscalars+invdir+( 1-1)* qq ) ! For wall, qVal=0.5 so fIn = fOut force = force + layout%fStencil%cxDirRK(:,invDir) * 2.0_rk * fOut end if end do !iComp ! store computed force on bndForce val bndForce(posInBndID( levelDesc%pntTID(elemPos) ), :) = force end do end subroutine mus_calcBndForce_wall ! ************************************************************************** ! ! ************************************************************************** ! !> This routine computes bndForce on wall_libb boundary elements subroutine mus_calcBndForce_wall_libb( me, bndForce, posInBndID, globBC, & & currState, levelDesc, nSize, iLevel, & & neigh, layout, nScalars ) ! --------------------------------------------------------------------- ! !> field boundary type class( boundary_type ), intent(in) :: me !> bndForce to fill real(kind=rk), intent(inout) :: bndForce(:,:) ! position of boundary element in boundary%boundaryID integer, intent(in) :: posInBndID(:) !> scheme global boundary type type( glob_boundary_type ), intent(in) :: globBC !> current state array to access post-collision values real(kind=rk), intent(in) :: currState(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> level which invokes boundary integer,intent(in) :: iLevel !> global parameters ! type(mus_param_type),intent(in) :: params integer,intent(in) :: neigh(:) !< connectivity array !> scheme layout type( mus_scheme_layout_type ),intent(in) :: layout !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !--------------------------------------------------------------------------- integer :: elemPos, posInBuffer, iElem, iDir, QQN, QQ integer :: invDir, iLink real(kind=rk) :: force(3), fIn, fOut, fOutNeigh, fIntp real(kind=rk) :: qVal, cIn, cOut, cNgh !--------------------------------------------------------------------------- QQ = layout%fStencil%QQ QQN = layout%fStencil%QQN iLink = 0 do iElem = 1, globBC%nElems_Fluid(iLevel) force = 0.0_rk elemPos = globBC%elemLvl( iLevel )%elem%val( iElem ) posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ) do iDir = 1, QQN if( globBC%elemLvl( iLevel )%bitmask%val( iDir, iElem )) then iLink = iLink + 1 ! direction towards boundary invDir = layout%fStencil%cxDirInv( iDir ) ! qValues qVal = globBC%elemLvl( iLevel )%qVal%val( invDir, iElem ) cIn = me%bouzidi(iLevel)% cIn( iLink ) cOut = me%bouzidi(iLevel)%cOut( iLink ) cNgh = me%bouzidi(iLevel)%cNgh( iLink ) ! Incoming state after collision fIn = currState( & & ( elempos-1)* nscalars+idir+( 1-1)* qq ) ! Outgoing state after collision fOut = currState( & & ( elempos-1)* nscalars+invdir+( 1-1)* qq ) ! Outgoing direction of neighbor after collision can be obtained ! using FETCH in invDir fOutNeigh = currState( & & neigh((invdir-1)* nsize+ elempos)+( 1-1)* qq+ nscalars*0 ) ! Incoming interpolated state according to bouzidi rule fIntp = cIn*fIn + cOut*fOut + cNgh*fOutNeigh ! ! For outgoing direction, post-collision is extracted from ! ! currState using SAVE macro in invDir since this value is ! ! not modified by any BC. So no need to use bcBuffer. ! fOut = currState( & ! & ( elempos-1)* nscalars+invdir+( 1-1)* qq ) ! force = force + layout%fStencil%cxDirRK(:, invDir) * (fOut + fIntp) end if end do !iComp ! store computed force on bndForce val bndForce(posInBndID( levelDesc%pntTID(elemPos) ), :) = force end do end subroutine mus_calcBndForce_wall_libb ! ************************************************************************** ! end module mus_bndForce_module