spongeLayer_box_roundCornerPolyn5_for_coord Function

private function spongeLayer_box_roundCornerPolyn5_for_coord(me, coord, n) result(res)

This function calculates the sigma for the box shape sponge layer from coord for polynomial 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.

Arguments

Type IntentOptional Attributes Name
type(tem_spongeLayer_box_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 real(kind=rk), (n)

return value


Called by

proc~~spongelayer_box_roundcornerpolyn5_for_coord~~CalledByGraph proc~spongelayer_box_roundcornerpolyn5_for_coord spongeLayer_box_roundCornerPolyn5_for_coord proc~spongelayer_box_roundcornerpolyn5_for_treeids spongelayer_box_roundCornerPolyn5_for_treeIDs proc~spongelayer_box_roundcornerpolyn5_for_treeids->proc~spongelayer_box_roundcornerpolyn5_for_coord proc~spongelayer_box_scalar_for_coord spongelayer_box_scalar_for_coord proc~spongelayer_box_scalar_for_coord->proc~spongelayer_box_roundcornerpolyn5_for_coord interface~tem_spongelayer_box_for tem_spongeLayer_box_for interface~tem_spongelayer_box_for->proc~spongelayer_box_scalar_for_coord proc~spongelayer_box_scalar_for_treeids spongelayer_box_scalar_for_treeIDs interface~tem_spongelayer_box_for->proc~spongelayer_box_scalar_for_treeids proc~spongelayer_box_vector_for_coord spongelayer_box_vector_for_coord interface~tem_spongelayer_box_for->proc~spongelayer_box_vector_for_coord proc~spongelayer_box_vector_for_treeids spongelayer_box_vector_for_treeIDs interface~tem_spongelayer_box_for->proc~spongelayer_box_vector_for_treeids proc~spongelayer_box_scalar_for_treeids->proc~spongelayer_box_roundcornerpolyn5_for_treeids proc~spongelayer_box_vector_for_coord->proc~spongelayer_box_scalar_for_coord proc~spongelayer_box_vector_for_treeids->proc~spongelayer_box_scalar_for_treeids proc~tem_spatial_for_coord tem_spatial_for_coord proc~tem_spatial_for_coord->interface~tem_spongelayer_box_for proc~tem_spatial_for_treeids tem_spatial_for_treeIDs proc~tem_spatial_for_treeids->interface~tem_spongelayer_box_for proc~tem_spatial_vector_for_coord tem_spatial_vector_for_coord proc~tem_spatial_vector_for_coord->interface~tem_spongelayer_box_for proc~tem_spatial_vector_for_treeids tem_spatial_vector_for_treeIDs proc~tem_spatial_vector_for_treeids->interface~tem_spongelayer_box_for interface~tem_spatial_for tem_spatial_for interface~tem_spatial_for->proc~tem_spatial_for_coord interface~tem_spatial_for->proc~tem_spatial_for_treeids interface~tem_spatial_for->proc~tem_spatial_vector_for_coord interface~tem_spatial_for->proc~tem_spatial_vector_for_treeids proc~tem_spatial_scalar_for_index tem_spatial_scalar_for_index proc~tem_spatial_scalar_for_index->proc~tem_spatial_for_coord proc~tem_spatial_vector_for_index tem_spatial_vector_for_index proc~tem_spatial_vector_for_index->proc~tem_spatial_vector_for_coord

Source Code

  function spongeLayer_box_roundCornerPolyn5_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_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
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_corMin(3), vec_corMax(3)
    real(kind=rk) :: vec_corMinSqr(3), vec_corMaxSqr(3) 
    real(kind=rk) :: rad, const_fac, box_max(3)
    real(kind=rk) :: corInRad, corOutRad, corMin(3), corMax(3)
    logical :: isCorner
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    extent(:) = me%extent
    box_max(:) = origin(:) + extent(:)
    corInRad = me%corner_radius
    corOutRad = corInRad + me%thickness
    corMin(:) = origin(:) + corInRad
    corMax(:) = box_max(:) - corInRad

    const_fac = 3125_rk/(256_rk*me%thickness**5)

    do i = 1,n
      coordLoc = coord(i,:)
      vec_corMin(:) = coordLoc(:) - corMin(:)
      vec_corMax(:) = coordLoc(:) - corMax(:)
      vec_corMinSqr(:) = vec_corMin(:)**2
      vec_corMaxSqr(:) = vec_corMax(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk
      isCorner = .true.

      ! Bottom-South-West: -x,-y,-z
      if (all(coordLoc(:) < corMin(:))) then
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-South-East: x, -y, -z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-North-West: -x, y, -z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-North-East: x, y, -z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-South-West: -x, -y, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-South-East: x, -y, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-North-West: -x, y, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (all(coordLoc(:) > corMax(:))) then ! Top-North-East: x, y, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(2) < corMin(2) .and. coordLoc(3) < corMin(3)) then
        ! Botton-South: -y,-z
        rad = sqrt( vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(2) < corMin(2) .and. coordLoc(3) > corMax(3)) then
        ! Top-South: -y, z
        rad = sqrt( vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(2) > corMax(2) .and. coordLoc(3) < corMin(3)) then
        ! Bottom-North: y, -z
        rad = sqrt( vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(2) > corMax(2) .and. coordLoc(3) > corMax(3)) then
        ! Top-North: y, z
        rad = sqrt( vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(3) < corMin(3)) then
        ! Botton-West: -x,-z
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(3) < corMin(3)) then
        ! Bottom-East: x,-z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(3) > corMax(3)) then
        ! Top-West: -x, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(3) > corMax(3)) then
        ! Top-East: x, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2)) then
        ! South-West: -x,-y
        rad = sqrt(vec_corMinSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2)) then
        ! North-West: -x, y
        rad = sqrt(vec_corMinSqr(1) + vec_corMaxSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2)) then
        ! South-East: x, -y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2)) then
        ! North-East: x, y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMaxSqr(2))

      else
        isCorner = .false.
        if (coordLoc(1) < origin(1)) then ! West: -x
          normal = -1_rk
          rad = coordLoc(1) - origin(1)

        else if (coordLoc(2) < origin(2)) then ! South: -y
          normal = -1_rk
          rad = coordLoc(2) - origin(2)

        else if (coordLoc(3) < origin(3)) then ! Bottom: -z
          normal = -1_rk
          rad = coordLoc(3) - origin(3)

        else if (coordLoc(1) > box_max(1)) then ! East: x
          normal = 1_rk
          rad = coordLoc(1) - box_max(1)

        else if (coordLoc(2) > box_max(2)) then ! North: y
          normal = 1_rk
          rad = coordLoc(2) - box_max(2)

        else if (coordLoc(3) > box_max(3)) then ! Top: z
          normal = 1_rk
          rad = coordLoc(3) - box_max(3)
        end if
      end if

      if (isCorner) then
        proj_len1 = rad - corInRad
        proj_len2 = corOutRad - rad
      else
        proj_len1 = rad*normal
        proj_len2 = (me%thickness*normal - rad)*normal
      end if

      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_box_roundCornerPolyn5_for_coord