This function calculates the sigma for the spongeLayer from coord for the polynomial order n5. Sponge profile: where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
Profile is taken from: 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.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(tem_spongeLayer_plane_type) | :: | me |
Spatial sponge layer to evaluate |
|||
real(kind=rk), | intent(in) | :: | coord(n,3) |
barycentric Ids of an elements. 1st index goes over number of elements and 2nd index goes over x,y,z coordinates |
||
integer, | intent(in) | :: | n |
Number of arrays to return |
return value
function spongeLayer_plane_polyn5_for_coord(me, coord, n) result(res)
! --------------------------------------------------------------------------
!> Spatial sponge layer to evaluate
type(tem_spongeLayer_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%normal
const_fac = 3125_rk/(256_rk*me%thickness**5)
do i = 1,n
vec1(:) = coord(i,:) - origin(:)
vec2(:) = me%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%dampFactor
else if (proj_len2 < 0) then ! If coord is beyond thickness
res(i) = me%dampFactor
else
res(i) = 0.0_rk
end if
end do
end function spongeLayer_plane_polyn5_for_coord