! Copyright (c) 2021 Kannan Masilamani <kannan.masilamani@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 absorbLayer definitions !! !! The definition of sponge profile in the absorb layer is based on: !! 1) Xu, Hui; Sagaut, Pierre (2013): Analysis of the absorbing layers for the !! weakly-compressible lattice Boltzmann methods. In Journal of Computational !! Physics 245, pp. 14–42. DOI: 10.1016/j.jcp.2013.02.051. !! 2) Jacob, J.; Sagaut, P. (2019): Solid wall and open boundary conditions in !! hybrid recursive regularized lattice Boltzmann method for compressible flows. !! In Physics of Fluids 31 (12), p. 126103. DOI: 10.1063/1.5129138. !! module tem_absorbLayer_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 :: tem_load_absorbLayer_plane, tem_load_absorbLayer_radial public :: tem_load_absorbLayer_box public :: tem_absorbLayer_plane_type, tem_absorbLayer_plane_for public :: tem_absorbLayer_radial_type, tem_absorbLayer_radial_for public :: tem_absorbLayer_box_type, tem_absorbLayer_box_for !> This type contains data to define base absorbLayer type absorbLayer_base_type !> Thickness of absorb layer real(kind=rk) :: thickness !> Strength factor for the absorb Layer real(kind=rk) :: strength !> target state: pressure, velocityX, velocityY and velocityZ real(kind=rk) :: targetState(4) end type absorbLayer_base_type !> This type contains data to define planar absorbLayer type tem_absorbLayer_plane_type !> Plane origin, start of absorb layer real(kind=rk) :: origin(3) !> Plane normal should point outwards the domain and it should be unitNormal real(kind=rk) :: unitNormal(3) !> Basic common data for all absorb layers type(absorbLayer_base_type) :: base end type tem_absorbLayer_plane_type !> This type contains data to define box absorbLayer. type tem_absorbLayer_box_type !> Box origin, bottom left corner of absorb layer real(kind=rk) :: origin(3) !> Length of box in each dimension real(kind=rk) :: extent(3) !> Basic common data for all absorb layers type(absorbLayer_base_type) :: base end type tem_absorbLayer_box_type !> This type contains data to define radial absorbLayer. type tem_absorbLayer_radial_type !> origin for radial absorblayer real(kind=rk) :: origin(3) !> inner radius, start of radial absorblayer real(kind=rk) :: radius !> Basic common data for all absorb layers type(absorbLayer_base_type) :: base end type tem_absorbLayer_radial_type !> Interface for planar absorbLayer interface tem_absorbLayer_plane_for module procedure absorbLayer_plane_for_coord module procedure absorbLayer_plane_for_treeIDs module procedure absorbLayer_plane_vector_for_coord module procedure absorbLayer_plane_vector_for_treeIDs end interface tem_absorbLayer_plane_for !> Interface for radial absorbLayer interface tem_absorbLayer_radial_for module procedure absorbLayer_radial_for_coord module procedure absorbLayer_radial_for_treeIDs module procedure absorbLayer_radial_vector_for_coord module procedure absorbLayer_radial_vector_for_treeIDs end interface tem_absorbLayer_radial_for !> Interface for box absorbLayer interface tem_absorbLayer_box_for module procedure absorbLayer_box_for_coord module procedure absorbLayer_box_for_treeIDs module procedure absorbLayer_box_vector_for_coord module procedure absorbLayer_box_vector_for_treeIDs end interface tem_absorbLayer_box_for contains ! -------------------------------------------------------------------------- ! !> This routine load base info for absorbLayer subroutine load_absorbLayer_base(me, nComp, conf, thandle) ! -------------------------------------------------------------------------- !> Base info for absorb layer type(absorbLayer_base_type), intent(out) :: me !> Number of components, load target_state only if nComp > 1 integer, intent(in) :: nComp !> lua state type type(flu_State) :: conf !> aotus parent handle integer, intent(in) :: thandle ! -------------------------------------------------------------------------- integer :: iError, state_handle integer :: vError(3), errfatal(3) ! -------------------------------------------------------------------------- errfatal = aotErr_Fatal ! sponge thickness call aot_get_val( L = conf, & & thandle = thandle, & & key = 'thickness', & & val = me%thickness, & & ErrCode = iError ) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*) 'ERROR reading the thickness of absorb layer. ' call tem_abort() end if write(logUnit(1),*) ' * thickness =', me%thickness ! sponge strength call aot_get_val( L = conf, & & thandle = thandle, & & key = 'strength', & & val = me%strength, & & ErrCode = iError ) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*) 'ERROR reading the strength of absorb layer. ' call tem_abort() end if write(logUnit(1),*) ' * Strength =', me%strength ! Load target_state only if nComp > 1, otherwise treat tem_absorblayer_for ! as scalar and return only sponge sigma if (nComp > 1) then write(logUnit(1),*) ' * Target state:' call aot_table_open( L = conf, & & parent = thandle, & & thandle = state_handle, & & key = 'target_state' ) ! target state: pressure call aot_get_val( L = conf, & & thandle = state_handle, & & key = 'pressure', & & val = me%targetState(1), & & ErrCode = iError ) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*) 'FATAL Error occured, when loading target state: '& & //'pressure' call tem_abort() end if write(logUnit(1),*) ' pressure =', me%targetState(1) ! target state: velocity call aot_get_val( L = conf, & & thandle = state_handle, & & key = 'velocity', & & val = me%targetState(2:4), & & ErrCode = vError ) if (any(btest( vError, errFatal )) ) then write(logUnit(1),*) 'FATAL Error occured, when loading target state: ' & & //'velocity' call tem_abort() end if write(logUnit(1),*) ' velocity =', me%targetState(2:4) call aot_table_close(conf, state_handle) end if end subroutine load_absorbLayer_base ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine load planar absorbLayer. !! Example: !!```lua !! spatial = { !! predefined = 'absorblayer_plane', !! origin = {0.0, 0.0, 0.0}, !! unit_normal = {1.0, 0.0, 0.0}, -- absorb layer on east boundary !! thickness = 2.0, !! strength = 1.0, !! target_state = { !! pressure = 1.0, !! velocity = {0.0, 0.0, 0.0} !! } !! } !!``` subroutine tem_load_absorbLayer_plane(me, conf, thandle, nComp) ! -------------------------------------------------------------------------- !> Plane absorb layer type(tem_absorbLayer_plane_type), intent(out) :: me !> lua state type type(flu_State) :: conf !> aotus parent handle integer, intent(in) :: thandle !> Number of components, load target_state only if nComp > 1 integer, intent(in) :: nComp ! -------------------------------------------------------------------------- integer :: vError(3), errfatal(3) ! -------------------------------------------------------------------------- errfatal = aotErr_Fatal ! Plane_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 origin of absorb layer. ' & & // 'It should have 3 entries for each coordinate.' call tem_abort() end if write(logUnit(1),*) ' * origin =', me%origin ! unit_normal call aot_get_val( L = conf, & & thandle = thandle, & & key = 'unit_normal', & & val = me%unitNormal, & & ErrCode = vError ) if (any(btest( vError, errFatal )) ) then write(logUnit(1),*) 'ERROR reading the unit_normal of absorb layer. ' & & // 'It should have 3 entries for each coordinate.' call tem_abort() end if write(logUnit(1),*) ' * unit_normal =', me%unitNormal ! Load base info for absorbLayer: thickness, strength and target_state call load_absorbLayer_base( me = me%base, & & nComp = nComp, & & conf = conf, & & thandle = thandle ) end subroutine tem_load_absorbLayer_plane ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine load box absorbLayer. !! Box must be axis-aligned so it is defined by an origin and extent in !! 3 dimensions. !! Example: !!```lua !! spatial = { !! predefined = 'absorblayer_box', !! origin = {0.0, 0.0, 0.0}, !! extent = {1.0, 3.0, 0.2}, !! }, !! thickness = 2.0, !! strength = 1.0, !! target_state = { !! pressure = 1.0, !! velocity = {0.0, 0.0, 0.0} !! } !! } !!``` subroutine tem_load_absorbLayer_box(me, conf, thandle, nComp) ! -------------------------------------------------------------------------- !> Plane absorb layer type(tem_absorbLayer_box_type), intent(out) :: me !> lua state type type(flu_State) :: conf !> aotus parent handle integer, intent(in) :: thandle !> Number of components, load target_state only if nComp > 1 integer, intent(in) :: nComp ! -------------------------------------------------------------------------- integer :: vError(3), errfatal(3) ! -------------------------------------------------------------------------- errfatal = aotErr_Fatal ! Box 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 origin of box absorb layer. ' & & // 'It should have 3 entries for each coordinate.' call tem_abort() end if write(logUnit(1),*) ' * origin =', me%origin ! Box extent call aot_get_val( L = conf, & & thandle = thandle, & & key = 'extent', & & val = me%extent, & & ErrCode = vError ) if (any(btest( vError, errFatal )) ) then write(logUnit(1),*) 'ERROR reading the extent of box absorb layer. ' & & // 'It should have 3 entries for each dimension.' call tem_abort() end if write(logUnit(1),*) ' * extent =', me%extent ! Load base info for absorbLayer: thickness, strength and target_state call load_absorbLayer_base( me = me%base, & & nComp = nComp, & & conf = conf, & & thandle = thandle ) end subroutine tem_load_absorbLayer_box ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine load radial absorbLayer. !! Example: !!```lua !! spatial = { !! predefined = 'absorblayer_radial', --or 'absorblayer_radial_2d' !! origin = {0.0, 0.0, 0.0}, -- Center of sponge !! radius = 1.0, -- Start of sponge !! thickness = 2.0, !! strength = 1.0, !! target_state = { !! pressure = 1.0, !! velocity = {0.0, 0.0, 0.0} !! } !! } !!``` subroutine tem_load_absorbLayer_radial(me, conf, thandle, nComp) ! -------------------------------------------------------------------------- !> Plane absorb layer type(tem_absorbLayer_radial_type), intent(out) :: me !> lua state type type(flu_State) :: conf !> aotus parent handle integer, intent(in) :: thandle !> Number of components, load target_state only if nComp > 1 integer, intent(in) :: nComp ! -------------------------------------------------------------------------- integer :: iError integer :: vError(3), errfatal(3) ! -------------------------------------------------------------------------- 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 origin for radial absorblayer, ' & & // '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 = 'radius', & & val = me%radius, & & ErrCode = iError ) if (btest(iError, aotErr_Fatal)) then write(logUnit(1),*) 'ERROR reading the radius for radial absorblayer' call tem_abort() end if write(logUnit(1),*) ' * radius =', me%radius ! Load base info for absorbLayer: thickness, strength and target_state call load_absorbLayer_base( me = me%base, & & nComp = nComp, & & conf = conf, & & thandle = thandle ) end subroutine tem_load_absorbLayer_radial ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from coord. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_plane_for_coord(me, coord, n) result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_plane_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, origin(3), normal(3), vec1(3), vec2(3) real(kind=rk) :: proj_len1, proj_len2, const_fac ! -------------------------------------------------------------------------- origin(:) = me%origin normal(:) = me%unitNormal const_fac = 3125_rk/(256_rk*me%base%thickness**5) do i = 1,n vec1(:) = coord(i,:) - origin(:) vec2(:) = me%base%thickness*normal(:) + origin(:) - coord(i,:) proj_len1 = vec1(1)*normal(1)+ vec1(2)*normal(2)+vec1(3)*normal(3) proj_len2 = vec2(1)*normal(1)+ vec2(2)*normal(2)+vec2(3)*normal(3) if (proj_len1 > 0 .and. proj_len2 > 0) then sigma = const_fac * proj_len2 * (proj_len1**4) res(i) = sigma*me%base%strength else res(i) = 0.0_rk end if end do end function absorbLayer_plane_for_coord ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from coord and !! fills up the res with the target state. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_plane_vector_for_coord(me, nComp, coord, n) & & result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_plane_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 ! -------------------------------------------------------------------------- res(:,1) = absorbLayer_plane_for_coord(me, coord, n) if (ncomp > 1) then do i = 1,n res(i,2:) = me%base%targetState(:) end do end if end function absorbLayer_plane_vector_for_coord ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from treeids. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_plane_for_treeids(me, treeids, tree, n) result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_plane_type) :: me !> Number of arrays to return integer, intent(in) :: n !> 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) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: sigma, origin(3), normal(3), vec1(3), vec2(3), coord(3) real(kind=rk) :: proj_len1, proj_len2, const_fac ! -------------------------------------------------------------------------- origin(:) = me%origin normal(:) = me%unitNormal const_fac = 3125_rk/(256_rk*me%base%thickness**5) do i = 1,n !barycentric coordinate coord = tem_BaryOfId( tree, treeIds(i) ) vec1(:) = coord(:) - origin(:) vec2(:) = me%base%thickness*normal(:) + origin(:) - coord(:) proj_len1 = vec1(1)*normal(1)+ vec1(2)*normal(2)+vec1(3)*normal(3) proj_len2 = vec2(1)*normal(1)+ vec2(2)*normal(2)+vec2(3)*normal(3) if (proj_len1 > 0 .and. proj_len2 > 0) then sigma = const_fac * proj_len2 * (proj_len1**4) res(i) = sigma*me%base%strength else res(i) = 0.0_rk end if end do end function absorbLayer_plane_for_treeids ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from treeids and !! fills up the res with the target state. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_plane_vector_for_treeids(me, nComp, treeids, tree, n) & & result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_plane_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 ! -------------------------------------------------------------------------- res(:,1) = absorbLayer_plane_for_treeids(me, treeids, tree, n) if (ncomp > 1) then do i = 1,n res(i,2:) = me%base%targetState(:) end do end if end function absorbLayer_plane_vector_for_treeids ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from coord and !! fills up the res with the target state. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_box_for_coord(me, coord, n) result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_box_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, origin(3), extent(3), coordLoc(3), normal, thickness real(kind=rk) :: proj_len1, proj_len2 real(kind=rk) :: vec_min(3), vec_max(3), vec_minSqr(3), vec_maxSqr(3) real(kind=rk) :: rad, const_fac, box_max(3) ! -------------------------------------------------------------------------- origin(:) = me%origin extent(:) = me%extent box_max(:) = origin(:) + extent(:) thickness = me%base%thickness const_fac = 3125_rk/(256_rk*me%base%thickness**5) do i = 1,n coordLoc = coord(i,:) vec_min(:) = coordLoc(:) - origin(:) vec_max(:) = coordLoc(:) - box_max(:) vec_minSqr(:) = vec_min(:)**2 vec_maxSqr(:) = vec_max(:)**2 proj_len1 = 0_rk proj_len2 = 0_rk normal = 1_rk rad = 0_rk ! Bottom-South-West: -x,-y,-z if (all(coordLoc(:) < origin(:))) then rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_minSqr(3)) ) else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) & & .and. coordLoc(3) < origin(3)) then ! Bottom-South-East: x, -y, -z rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_minSqr(3)) ) else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) & & .and. coordLoc(3) < origin(3)) then ! Bottom-North-West: -x, y, -z rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_minSqr(3)) ) else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2) & & .and. coordLoc(3) < origin(3)) then ! Bottom-North-East: x, y, -z rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_minSqr(3)) ) else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2) & & .and. coordLoc(3) > box_max(3)) then ! Top-South-West: -x, -y, z rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_maxSqr(3)) ) else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) & & .and. coordLoc(3) > box_max(3)) then ! Top-South-East: x, -y, z rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_maxSqr(3)) ) else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) & & .and. coordLoc(3) > box_max(3)) then ! Top-North-West: -x, y, z rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_maxSqr(3)) ) else if (all(coordLoc(:) > box_max(:))) then ! Bottom-North-East: x, y, -z rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_maxSqr(3)) ) else if (coordLoc(2) < origin(2) .and. coordLoc(3) < origin(3)) then ! Botton-South: -y,-z rad = sqrt( max(vec_minSqr(2), vec_minSqr(3)) ) else if (coordLoc(2) < origin(2) .and. coordLoc(3) > box_max(3)) then ! Top-South: -y, z rad = sqrt( max(vec_minSqr(2), vec_maxSqr(3)) ) else if (coordLoc(2) > box_max(2) .and. coordLoc(3) < origin(3)) then ! Botom-North: -y, z rad = sqrt( max(vec_maxSqr(2), vec_minSqr(3)) ) else if (coordLoc(2) > box_max(2) .and. coordLoc(3) > box_max(3)) then ! Top-North: y, z rad = sqrt( max(vec_maxSqr(2), vec_maxSqr(3)) ) else if (coordLoc(1) < origin(1) .and. coordLoc(3) < origin(3)) then ! Botton-West: -x,-z rad = sqrt( max(vec_minSqr(1), vec_minSqr(3)) ) else if (coordLoc(1) > box_max(1) .and. coordLoc(3) < origin(3)) then ! Bottom-East: x,-z rad = sqrt( max(vec_maxSqr(1), vec_minSqr(3)) ) else if (coordLoc(1) < origin(1) .and. coordLoc(3) > box_max(3)) then ! Top-West: -x, z rad = sqrt( max(vec_minSqr(1), vec_maxSqr(3)) ) else if (coordLoc(1) > box_max(1) .and. coordLoc(3) > box_max(3)) then ! Top-East: x, z rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(3)) ) else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2)) then ! South-West: -x,-y rad = sqrt( max(vec_minSqr(1), vec_minSqr(2)) ) else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2)) then ! North-West: -x, y rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2)) ) else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2)) then ! South-East: x, -y rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2)) ) else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2)) then ! North-East: x, y rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2)) ) else if (coordLoc(1) < origin(1)) then ! West: -x normal = -1_rk rad = vec_min(1) else if (coordLoc(2) < origin(2)) then ! South: -y normal = -1_rk rad = vec_min(2) else if (coordLoc(3) < origin(3)) then ! Bottom: -z normal = -1_rk rad = vec_min(3) else if (coordLoc(1) > box_max(1)) then ! East: x normal = 1_rk rad = vec_max(1) else if (coordLoc(2) > box_max(2)) then ! North: y normal = 1_rk rad = vec_max(2) else if (coordLoc(3) > box_max(3)) then ! Top: z normal = 1_rk rad = vec_max(3) end if proj_len1 = rad*normal proj_len2 = (thickness*normal - rad)*normal if (proj_len1 > 0 .and. proj_len2 > 0) then sigma = const_fac * proj_len2 * (proj_len1**4) res(i) = sigma*me%base%strength else res(i) = 0.0_rk end if end do end function absorbLayer_box_for_coord ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from coord and !! fills up the res with the target state. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_box_vector_for_coord(me, nComp, coord, n) & & result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_box_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 ! -------------------------------------------------------------------------- res(:,1) = absorbLayer_box_for_coord(me, coord, n) if (ncomp > 1) then do i = 1,n res(i,2:) = me%base%targetState(:) end do end if end function absorbLayer_box_vector_for_coord ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from treeids and !! fills up the res with the target state. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_box_for_treeids(me, treeids, tree, n) & & result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_box_type) :: me !> Number of arrays to return integer, intent(in) :: n !> 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) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: res_i(1), coord(3,1) ! -------------------------------------------------------------------------- do i = 1,n !barycentric coordinate coord(:,1) = tem_BaryOfId( tree, treeIds(i) ) res_i = absorbLayer_box_for_coord(me, coord, 1) res(i) = res_i(1) end do end function absorbLayer_box_for_treeids ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from treeids and !! fills up the res with the target state. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_box_vector_for_treeids(me, nComp, treeids, tree, n) & & result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_box_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 ! -------------------------------------------------------------------------- res(:,1) = absorbLayer_box_for_treeids(me, treeids, tree, n) if (ncomp > 1) then do i = 1,n res(i,2:) = me%base%targetState(:) end do end if end function absorbLayer_box_vector_for_treeids ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from coord. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_radial_for_coord(me, coord, n, nDim) & & result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_radial_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) !> Dimension integer, intent(in) :: nDim !> return value real(kind=rk) :: res(n) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: rad, outer_radius, sigma, vec(3) real(kind=rk) :: const_fac ! -------------------------------------------------------------------------- outer_radius = me%radius + me%base%thickness vec = 0.0_rk const_fac = 3125_rk/(256_rk*me%base%thickness**5) do i = 1,n vec(1:nDim) = coord(i,1:nDim) - me%origin(1:nDim) rad = sqrt( vec(1)**2 + vec(2)**2 + vec(3)**2 ) if (rad > me%radius .and. rad < outer_radius) then sigma = const_fac * (outer_radius - rad) * (rad - me%radius)**4 res(i) = sigma * me%base%strength else res(i) = 0.0_rk end if end do end function absorbLayer_radial_for_coord ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from coord and !! fills up the res with the target state. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_radial_vector_for_coord(me, nComp, coord, n, nDim) & & result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_radial_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) !> Dimension integer, intent(in) :: nDim !> return value real(kind=rk) :: res(n,ncomp) ! -------------------------------------------------------------------------- integer :: i ! -------------------------------------------------------------------------- res(:,1) = absorbLayer_radial_for_coord(me, coord, n, nDim) if (ncomp > 1) then do i = 1,n res(i,2:) = me%base%targetState(:) end do end if end function absorbLayer_radial_vector_for_coord ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from treeids. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_radial_for_treeids(me, treeids, tree, n, nDim) & & result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_radial_type) :: me !> Number of arrays to return integer, intent(in) :: n !> global treelm mesh type( treelmesh_type ), intent(in) ::tree !> treeIds of elements in given level integer(kind=long_k), intent(in) :: treeIds(n) !> Dimension integer, intent(in) :: nDim !> return value real(kind=rk) :: res(n) ! -------------------------------------------------------------------------- integer :: i real(kind=rk) :: rad, outer_radius, sigma, vec(3), coord(3) real(kind=rk) :: const_fac ! -------------------------------------------------------------------------- outer_radius = me%radius + me%base%thickness vec = 0.0_rk const_fac = 3125_rk/(256_rk*me%base%thickness**5) do i = 1,n !barycentric coordinate coord = tem_BaryOfId( tree, treeIds(i) ) vec(1:nDim) = coord(1:nDim) - me%origin(1:nDim) rad = sqrt( vec(1)**2 + vec(2)**2 + vec(3)**2 ) if (rad > me%radius .and. rad < outer_radius) then sigma = const_fac * (outer_radius - rad) * (rad - me%radius)**4 res(i) = sigma * me%base%strength else res(i) = 0.0_rk end if end do end function absorbLayer_radial_for_treeids ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This function calculates the sigma for the absorblayer from treeids and !! fills up the res with the target state. !! Sponge profile: !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$ !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge. function absorbLayer_radial_vector_for_treeids(me, nComp, treeids, tree, n, & & nDim) result(res) ! -------------------------------------------------------------------------- !> Spatial absorb layer to evaluate type(tem_absorbLayer_radial_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) !> Dimension integer, intent(in) :: nDim !> return value real(kind=rk) :: res(n,ncomp) ! -------------------------------------------------------------------------- integer :: i ! -------------------------------------------------------------------------- res(:,1) = absorbLayer_radial_for_treeids(me, treeids, tree, n, nDim) if (ncomp > 1) then do i = 1,n res(i,2:) = me%base%targetState(:) end do end if end function absorbLayer_radial_vector_for_treeids ! -------------------------------------------------------------------------- ! end module tem_absorbLayer_module