! Copyright (c) 2013-2015 Nikhil Anand <nikhil.anand@uni-siegen.de> ! Copyright (c) 2013,2015-2016,2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2014 Jens Zudrop <j.zudrop@grs-sim.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de> ! Copyright (c) 2018 Neda Ebrahimi Pour <neda.epour@uni-siegen.de> ! Copyright (c) 2019 Harald Klimach <harald.klimach@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. ! ****************************************************************************** ! !> This module gathers the various spatial spongeLayer definitions !! module tem_spongeLayer_module ! include treelm modules use env_module, only: rk, long_k, labelLen use tem_param_module, only: PI use tem_aux_module, only: tem_abort use tem_logging_module, only: logUnit use treelmesh_module, only: treelmesh_type use tem_geometry_module, only: tem_BaryOfId ! include aotus modules use aotus_module, only: flu_State, aot_get_val, & & aoterr_NonExistent, aoterr_Fatal use aot_table_module, only: aot_table_open, aot_table_close, & & aot_table_length, aot_get_val implicit none private public :: load_spongeLayer, spongeLayer_type public :: tem_spongelayer_for, tem_spongelayer_for_treeids public :: tem_radialSpongeLayer_type public :: tem_load_radialSpongeLayer public :: tem_radialSpongeLayer2d_for, tem_radialSpongeLayer_for public :: tem_radialViscSpongeLayer2d_for, tem_radialViscSpongeLayer_for !> This type contains data to define spongeLayer type spongeLayer_type !> Sponge Plane origin real(kind=rk) :: plane_origin(3) !> Sponge Plane normal real(kind=rk) :: plane_normal(3) !> Damp factor for the sponge Layer real(kind=rk) :: dampFactor !> damping exponent for the sponge layer real(kind=rk) :: dampExponent !> target state real(kind=rk), allocatable :: targetState(:) end type spongeLayer_type !> This type contains data to define radial spongeLayer type tem_radialSpongeLayer_type !> Sponge radial origin real(kind=rk) :: origin(3) !> Sponge inner radius i.e. sponge start real(kind=rk) :: innerRadius !> Sponge outer radius i.e. sponge end real(kind=rk) :: outerRadius !> Damp factor for the sponge Layer real(kind=rk) :: dampFactor !> damping exponent for the sponge layer real(kind=rk) :: dampExponent !> target state. !! For viscous sponge, viscosity is stored and multiplied with sponge !! strength real(kind=rk), allocatable :: targetState(:) end type tem_radialSpongeLayer_type contains ! -------------------------------------------------------------------------- ! subroutine load_spongeLayer(conf, thandle, me, ndim) ! -------------------------------------------------------------------------- !> lua state type type(flu_State) :: conf !> aotus parent handle integer, intent(in) :: thandle !> Global spongeLayer data type type(spongeLayer_type), intent(out) :: me !> number of spatial dimensions integer, intent(in) :: ndim ! -------------------------------------------------------------------------- integer :: iError integer :: vError(3), errfatal(3) integer :: iState, nState character(len=labelLen), allocatable :: stateName(:) ! -------------------------------------------------------------------------- errfatal = aotErr_Fatal ! Plane_origin call aot_get_val( L = conf, & & thandle = thandle, & & key = 'plane_origin', & & val = me%plane_origin, & & ErrCode = vError ) if (any(btest( vError, errFatal )) ) then write(logUnit(1),*) 'ERROR reading the plane_origin of sponge layer. ' & & // 'It should have 3 entries for each coordinate.' call tem_abort() end if write(logUnit(1),*) ' * plane_origin =', me%plane_origin ! Plane_normal call aot_get_val( L = conf, & & thandle = thandle, & & key = 'plane_normal', & & val = me%plane_normal, & & ErrCode = vError ) if (any(btest( vError, errFatal )) ) then write(logUnit(1),*) 'ERROR reading the plane_normal of sponge layer. ' & & // 'It should have 3 entries for each coordinate.' call tem_abort() end if write(logUnit(1),*) ' * plane_normal =', me%plane_normal !damp_factor call aot_get_val( L = conf, & & thandle = thandle, & & key = 'damp_factor', & & val = me%dampFactor, & & ErrCode = iError ) write(logUnit(1),*) ' * Damping_factor =', me%dampFactor !damp_exponent call aot_get_val( L = conf, & & thandle = thandle, & & key = 'damp_exponent', & & val = me%dampExponent, & & ErrCode = iError ) if (iError .ne. 0) then me%dampExponent = 1.0 end if write(logUnit(1),*) ' * Damping_exponent =', me%dampExponent ! Load target states depending on nDim call load_defaultTargetState( conf = conf, & & parent = thandle, & & nDim = nDim, & & targetState = me%targetState ) end subroutine load_spongeLayer ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This subroutine load data for radial sponge layer !! Example: !! !!```lua !! spatial = { !! predefined = 'radial_viscous_spongelayer', !! origin = {0.0,0.0,0.0}, !! inner_radius = 1.0, -- Sponge start !! outer_radius = 1.3, -- Sponge end !! damp_factor = 0.5, !! damp_exponent = 1.0, !! target_state = { !! Default: density, velocityX, velocityY, velocityZ and pressure !! viscosity = 1e-3 !! } !!``` subroutine tem_load_radialSpongeLayer(me, conf, thandle, nDim, stateName) ! -------------------------------------------------------------------------- !> Radial spongeLayer data type type(tem_radialSpongeLayer_type), intent(out) :: me !> lua state type type(flu_State) :: conf !> aotus parent handle integer, intent(in) :: thandle !> number of Dimension for nonViscous sponges integer, intent(in), optional :: nDim !> Load stateName from target_state table character(len=*), intent(in), optional :: stateName ! -------------------------------------------------------------------------- integer :: iError, ts_handle integer :: vError(3), errfatal(3) integer :: nDim_loc ! -------------------------------------------------------------------------- if (present(nDim)) then nDim_loc = nDim else nDim_loc = 1 end if errfatal = aotErr_Fatal ! Sponge origin call aot_get_val( L = conf, & & thandle = thandle, & & key = 'origin', & & val = me%origin, & & ErrCode = vError ) if (any(btest( vError, errFatal )) ) then write(logUnit(1),*) 'ERROR reading the sponge origin, ' & & // 'origin is not well defined. ' & & // 'It should have 3 entries for each coordinate.' call tem_abort() end if write(logUnit(1),*) ' * origin =', me%origin ! Sponge inner radius call aot_get_val( L = conf, & & thandle = thandle, & & key = 'inner_radius', & & val = me%innerRadius, & & ErrCode = iError ) if (btest(iError, aotErr_Fatal)) then write(logUnit(1),*) 'ERROR reading the sponge inner_radius. ' call tem_abort() end if write(logUnit(1),*) ' * inner_radius =', me%innerRadius ! Sponge outer radius call aot_get_val( L = conf, & & thandle = thandle, & & key = 'outer_radius', & & val = me%outerRadius, & & ErrCode = iError ) if (btest(iError, aotErr_Fatal)) then write(logUnit(1),*) 'ERROR reading the sponge outer_radius. ' call tem_abort() end if write(logUnit(1),*) ' * outer_radius =', me%outerRadius ! damp_factor call aot_get_val( L = conf, & & thandle = thandle, & & key = 'damp_factor', & & val = me%dampFactor, & & ErrCode = iError ) if (btest(iError, aotErr_Fatal)) then write(logUnit(1),*) 'ERROR reading the sponge outer_radius. ' call tem_abort() end if write(logUnit(1),*) ' * Damping_factor =', me%dampFactor !damp_exponent call aot_get_val( L = conf, & & thandle = thandle, & & key = 'damp_exponent', & & val = me%dampExponent, & & ErrCode = iError ) if (iError .ne. 0) then me%dampExponent = 1.0 end if write(logUnit(1),*) ' * Damping_exponent =', me%dampExponent ! Load stateName provided by caller function if ( present(stateName) ) then allocate(me%targetState(1)) call aot_table_open( L = conf, & & parent = thandle, & & thandle = ts_handle, & & key = 'target_state' ) call aot_get_val( L = conf, & & thandle = ts_handle, & & key = trim(stateName), & & val = me%targetState(1), & & ErrCode = iError ) if (btest(iError, aotErr_Fatal)) then write(logUnit(1),*) 'ERROR reading the target state: '//trim(stateName) call tem_abort() end if call aot_table_close(conf, thandle) write(logUnit(1),*) ' * Target state:' write(logUnit(1),*) ' '//trim(stateName)//' =', me%targetState(1) else call load_defaultTargetState( conf = conf, & & parent = thandle, & & nDim = nDim_loc, & & targetState = me%targetState ) end if end subroutine tem_load_radialSpongeLayer ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine loads state names from target_state table subroutine load_defaultTargetState(conf, parent, nDim, targetState) ! -------------------------------------------------------------------------- !> lua state type type(flu_State) :: conf !> aotus parent handle integer, intent(in) :: parent !> Number of dimension integer, intent(in) :: nDim !> Target state value real(kind=rk), allocatable, intent(out) :: targetState(:) ! -------------------------------------------------------------------------- integer :: thandle, iState, nState integer :: iError !> List of stateNames character(len=labelLen), allocatable :: stateName(:) ! -------------------------------------------------------------------------- nState = 0 select case (nDim) case (3) nState = 5 allocate(stateName(nState)) stateName = [ 'density ', 'velocityX', & & 'velocityY', 'velocityZ', & & 'pressure ' ] case (2) nState = 4 allocate(stateName(nState)) stateName = [ 'density ', 'velocityX', & & 'velocityY', 'pressure ' ] case (1) nState = 3 allocate(stateName(nState)) stateName = [ 'density ', 'velocityX', & & 'pressure ' ] end select ! target_state allocate(targetState(nState)) call aot_table_open( L = conf, & & parent = parent, & & thandle = thandle, & & key = 'target_state' ) write(logUnit(1),*) ' * Target state:' do iState = 1, nState call aot_get_val( L = conf, & & thandle = thandle, & & key = trim(stateName(iState)), & & val = targetState(iState), & & ErrCode = iError ) if (btest(iError, aoterr_Fatal)) then write(*,*) 'FATAL Error occured, when loading target state: ' & & // trim(stateName(iState))// '! Aborting' call tem_abort() end if write(logUnit(1),*) ' '//trim(stateName(iState))//' =', & & targetState(iState) end do call aot_table_close(conf, thandle) end subroutine load_defaultTargetState ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the spongelayer and fills up !! the res with the target state function tem_spongelayer_for(me, nComp, coord, n) & & result(res) !> Spacetime function to evaluate type(spongeLayer_type) :: me !> Number of arrays to return integer, intent(in) :: n !> Number of entrys in each array integer, intent(in) :: ncomp !> barycentric Ids of an elements. !! 1st index goes over number of elements and !! 2nd index goes over x,y,z coordinates real(kind=rk), intent( in ) :: coord(n,3) !> return value real(kind=rk) :: res(n,ncomp) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: sigma, origin(3), normal(3), vec(3), proj_len ! -------------------------------------------------------------------------- origin(:) = me%plane_origin normal(:) = me%plane_normal do i = 1,n vec (:) = coord(i,:) - origin(:) proj_len = (vec(1)*normal(1)+ vec(2)*normal(2)+vec(3)*normal(3))/ & (normal(1)**2 + normal(2)**2 + normal(3)**2) sigma = me%dampFactor*((proj_len)**me%dampExponent) if (proj_len .gt. 0) then res(i,1) = sigma else res(i,1) = 0.0 end if end do if (ncomp > 1) then do i = 1,n res(i,2:) = me%targetState(:) end do end if end function tem_spongelayer_for ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the spongelayer and fills up !! the res with the target state function tem_spongelayer_for_treeids(me, ncomp, treeids, tree, n) & & result(res) !> Spacetime function to evaluate type(spongeLayer_type) :: me !> Number of arrays to return integer, intent(in) :: n !> Number of entrys in each array integer, intent(in) :: ncomp !> global treelm mesh type( treelmesh_type ), intent(in) ::tree !> treeIds of elements in given level integer(kind=long_k), intent(in) :: treeIds(n) !> return value real(kind=rk) :: res(n,ncomp) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: sigma, origin(3), normal(3), vec(3), coord(3), proj_len ! -------------------------------------------------------------------------- origin(:) = me%plane_origin normal(:) = me%plane_normal do i = 1,n !barycentric coordinate coord = tem_BaryOfId( tree, treeIds(i) ) vec (:) = coord(:) - origin(:) proj_len = (vec(1)*normal(1)+ vec(2)*normal(2)+vec(3)*normal(3))/ & (normal(1)**2 + normal(2)**2 + normal(3)**2) sigma = me%dampFactor*((proj_len)**me%dampExponent) if (proj_len .gt. 0) then res(i,1) = sigma else res(i,1) = 0.0 end if end do if (ncomp > 1) then do i = 1,n res(i,2:) = me%targetState(:) end do end if end function tem_spongelayer_for_treeids ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the 2d radial viscosity spongelayer !! and fills up rest with target_state. !! This function is currectly used to define viscosity sponge in musubi. function tem_radialSpongelayer2d_for(me, nComp, coord, n) & & result(res) !> Spacetime function to evaluate type(tem_radialSpongeLayer_type) :: me !> Number of arrays to return integer, intent(in) :: n !> Number of entrys in each array integer, intent(in) :: nComp !> barycentric Ids of an elements. !! 1st index goes over number of elements and !! 2nd index goes over x,y,z coordinates real(kind=rk), intent( in ) :: coord(n,3) !> return value real(kind=rk) :: res(n,nComp) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: sigma, rad, origin(2), length ! -------------------------------------------------------------------------- origin(:) = me%origin(1:2) length = me%outerRadius - me%innerRadius do i = 1, n rad = sqrt( (coord(i,1) - origin(1))**2 & & + (coord(i,2) - origin(2))**2 ) if ( rad > me%innerRadius .and. rad < me%outerRadius ) then sigma = ( (rad-me%innerRadius)/length ) * me%dampFactor else sigma = 0.0_rk end if res(i,1) = sigma end do if(ncomp > 1) then do i = 1,n res(i,2:) = me%targetState(:) enddo end if end function tem_radialSpongelayer2d_for ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the radial viscosity spongelayer !! and fills up rest with target_state. !! This function is currectly used to define viscosity sponge in musubi. function tem_radialSpongelayer_for(me, nComp, coord, n) & & result(res) !> Spacetime function to evaluate type(tem_radialSpongeLayer_type) :: me !> Number of arrays to return integer, intent(in) :: n !> Number of entrys in each array integer, intent(in) :: nComp !> barycentric Ids of an elements. !! 1st index goes over number of elements and !! 2nd index goes over x,y,z coordinates real(kind=rk), intent( in ) :: coord(n,3) !> return value real(kind=rk) :: res(n, nComp) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: sigma, rad, origin(3), length ! -------------------------------------------------------------------------- origin(:) = me%origin length = me%outerRadius - me%innerRadius do i = 1, n rad = sqrt( (coord(i,1) - origin(1))**2 & & + (coord(i,2) - origin(2))**2 & & + (coord(i,3) - origin(3))**2 ) if ( rad > me%innerRadius .and. rad < me%outerRadius ) then sigma = ( (rad-me%innerRadius)/length ) * me%dampFactor else sigma = 0.0_rk end if res(i,1) = sigma end do if (ncomp > 1) then do i = 1,n res(i,2:) = me%targetState(:) end do end if end function tem_radialSpongelayer_for ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the 2d radial viscosity spongelayer !! and multiply with targetState. !! This function is currectly used to define viscosity sponge in musubi. function tem_radialViscSpongelayer2d_for(me, coord, n) & & result(res) !> Spacetime function to evaluate type(tem_radialSpongeLayer_type) :: me !> Number of arrays to return integer, intent(in) :: n !> barycentric Ids of an elements. !! 1st index goes over number of elements and !! 2nd index goes over x,y,z coordinates real(kind=rk), intent( in ) :: coord(n,3) !> return value real(kind=rk) :: res(n) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: sigma, rad, origin(2), length ! -------------------------------------------------------------------------- origin(:) = me%origin(1:2) length = me%outerRadius - me%innerRadius do i = 1, n rad = sqrt( (coord(i,1) - origin(1))**2 & & + (coord(i,2) - origin(2))**2 ) if ( rad > me%innerRadius .and. rad < me%outerRadius ) then sigma = 1.0_rk + ( (rad-me%innerRadius)/length ) & & * (me%dampFactor-1.0_rk) else if (rad > me%outerRadius) then sigma = me%dampFactor else sigma = 1.0_rk end if res(i) = sigma * me%targetState(1) enddo end function tem_radialViscSpongelayer2d_for ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the radial viscosity spongelayer !! and multiply with targetState. !! This function is currectly used to define viscosity sponge in musubi. function tem_radialViscSpongelayer_for(me, coord, n) & & result(res) !> Spacetime function to evaluate type(tem_radialSpongeLayer_type) :: me !> Number of arrays to return integer, intent(in) :: n !> barycentric Ids of an elements. !! 1st index goes over number of elements and !! 2nd index goes over x,y,z coordinates real(kind=rk), intent( in ) :: coord(n,3) !> return value real(kind=rk) :: res(n) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: sigma, rad, origin(3), length ! -------------------------------------------------------------------------- origin(:) = me%origin length = me%outerRadius - me%innerRadius do i = 1, n rad = sqrt( (coord(i,1) - origin(1))**2 & & + (coord(i,2) - origin(2))**2 & & + (coord(i,3) - origin(3))**2 ) if ( rad > me%innerRadius .and. rad < me%outerRadius ) then sigma = 1.0_rk + ( (rad-me%innerRadius)/length ) & & * (me%dampFactor-1.0_rk) else if (rad > me%outerRadius) then sigma = me%dampFactor else sigma = 1.0_rk end if res(i) = sigma * me%targetState(1) enddo end function tem_radialViscSpongelayer_for ! -------------------------------------------------------------------------- ! end module tem_spongeLayer_module