PPM_functions.F90

!> Provides functions used with the Piecewise-Parabolic-Method in the vertical ALE algorithm.
module ppm_functions

! This file is part of MOM6. See LICENSE.md for the license.

! First version was created by Laurent White, June 2008.
! Substantially re-factored January 2016.

!! @todo Re-factor PPM_boundary_extrapolation to give round-off safe and
!!       optimization independent results.

use regrid_edge_values, only : bound_edge_values, check_discontinuous_edge_values

implicit none ; private

public ppm_reconstruction, ppm_boundary_extrapolation

!> A tiny width that is so small that adding it to cell widths does not
!! change the value due to a computational representation. It is used
!! to avoid division by zero.
!! @note This is a dimensional parameter and should really include a unit
!!       conversion.
real, parameter :: hneglect_dflt = 1.e-30

contains

!> Builds quadratic polynomials coefficients from cell mean and edge values.
subroutine ppm_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answers_2018)
  integer,              intent(in)    :: n !< Number of cells
  real, dimension(N),   intent(in)    :: h !< Cell widths [H]
  real, dimension(N),   intent(in)    :: u !< Cell averages [A]
  real, dimension(N,2), intent(inout) :: edge_values !< Edge values [A]
  real, dimension(N,3), intent(inout) :: ppoly_coef !< Polynomial coefficients, mainly [A]
  real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
  logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.

  ! Local variables
  integer   :: k              ! Loop index
  real      :: edge_l, edge_r ! Edge values (left and right)

  ! PPM limiter
  call ppm_limiter_standard( n, h, u, edge_values, h_neglect, answers_2018 )

  ! Loop over all cells
  do k = 1,n

    edge_l = edge_values(k,1)
    edge_r = edge_values(k,2)

    ! Store polynomial coefficients
    ppoly_coef(k,1) = edge_l
    ppoly_coef(k,2) = 4.0 * ( u(k) - edge_l ) + 2.0 * ( u(k) - edge_r )
    ppoly_coef(k,3) = 3.0 * ( ( edge_r - u(k) ) + ( edge_l - u(k) ) )

  enddo

end subroutine ppm_reconstruction

!> Adjusts edge values using the standard PPM limiter (Colella & Woodward, JCP 1984)
!! after first checking that the edge values are bounded by neighbors cell averages
!! and that the edge values are monotonic between cell averages.
subroutine ppm_limiter_standard( N, h, u, edge_values, h_neglect, answers_2018 )
  integer,              intent(in)    :: N !< Number of cells
  real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
  real, dimension(:),   intent(in)    :: u !< cell average properties (size N) [A]
  real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A]
  real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
  logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.

  ! Local variables
  integer   :: k              ! Loop index
  real      :: u_l, u_c, u_r  ! Cell averages (left, center and right)
  real      :: edge_l, edge_r ! Edge values (left and right)
  real      :: expr1, expr2

  ! Bound edge values
  call bound_edge_values( n, h, u, edge_values, h_neglect, answers_2018 )

  ! Make discontinuous edge values monotonic
  call check_discontinuous_edge_values( n, u, edge_values )

  ! Loop on interior cells to apply the standard
  ! PPM limiter (Colella & Woodward, JCP 84)
  do k = 2,n-1

    ! Get cell averages
    u_l = u(k-1)
    u_c = u(k)
    u_r = u(k+1)

    edge_l = edge_values(k,1)
    edge_r = edge_values(k,2)

    if ( (u_r - u_c)*(u_c - u_l) <= 0.0) then
      ! Flatten extremum
      edge_l = u_c
      edge_r = u_c
    else
      expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r))
      expr2 = (edge_r - edge_l) * (edge_r - edge_l)
      if ( expr1 > expr2 ) then
        ! Place extremum at right edge of cell by adjusting left edge value
        edge_l = u_c + 2.0 * ( u_c - edge_r )
        edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) ) ! In case of round off
      elseif ( expr1 < -expr2 ) then
        ! Place extremum at left edge of cell by adjusting right edge value
        edge_r = u_c + 2.0 * ( u_c - edge_l )
        edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) ) ! In case of round off
      endif
    endif
    ! This checks that the difference in edge values is representable
    ! and avoids overshoot problems due to round off.
    !### The 1.e-60 needs to have units of [A], so this dimensionally inconsisent.
    if ( abs( edge_r - edge_l )<max(1.e-60,epsilon(u_c)*abs(u_c)) ) then
      edge_l = u_c
      edge_r = u_c
    endif

    edge_values(k,1) = edge_l
    edge_values(k,2) = edge_r

  enddo ! end loop on interior cells

  ! PCM within boundary cells
  edge_values(1,:) = u(1)
  edge_values(n,:) = u(n)

end subroutine ppm_limiter_standard


!------------------------------------------------------------------------------
!> Reconstruction by parabolas within boundary cells
subroutine ppm_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_neglect)
!------------------------------------------------------------------------------
! Reconstruction by parabolas within boundary cells.
!
! The following explanations apply to the left boundary cell. The same
! reasoning holds for the right boundary cell.
!
! A parabola needs to be built in the cell and requires three degrees of
! freedom, which are the right edge value and slope and the cell average.
! The right edge values and slopes are taken to be that of the neighboring
! cell (i.e., the left edge value and slope of the neighboring cell).
! The resulting parabola is not necessarily monotonic and the traditional
! PPM limiter is used to modify one of the edge values in order to yield
! a monotonic parabola.
!
! N:     number of cells in grid
! h:     thicknesses of grid cells
! u:     cell averages to use in constructing piecewise polynomials
! edge_values : edge values of piecewise polynomials
! ppoly_coef : coefficients of piecewise polynomials
!
! It is assumed that the size of the array 'u' is equal to the number of cells
! defining 'grid' and 'ppoly'. No consistency check is performed here.
!------------------------------------------------------------------------------

  ! Arguments
  integer,              intent(in)    :: n !< Number of cells
  real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
  real, dimension(:),   intent(in)    :: u !< cell averages (size N) [A]
  real, dimension(:,:), intent(inout) :: edge_values    !< edge values of piecewise polynomials [A]
  real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly [A]
  real,       optional, intent(in)    :: h_neglect  !< A negligibly small width for
                                           !! the purpose of cell reconstructions [H]

  ! Local variables
  integer :: i0, i1
  real    :: u0, u1
  real    :: h0, h1
  real    :: a, b, c
  real    :: u0_l, u0_r
  real    :: u1_l, u1_r
  real    :: slope
  real    :: exp1, exp2
  real    :: hneglect

  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect

  ! ----- Left boundary -----
  i0 = 1
  i1 = 2
  h0 = h(i0)
  h1 = h(i1)
  u0 = u(i0)
  u1 = u(i1)

  ! Compute the left edge slope in neighboring cell and express it in
  ! the global coordinate system
  b = ppoly_coef(i1,2)
  u1_r = b *((h0+hneglect)/(h1+hneglect))     ! derivative evaluated at xi = 0.0,
                        ! expressed w.r.t. xi (local coord. system)

  ! Limit the right slope by the PLM limited slope
  slope = 2.0 * ( u1 - u0 )
  if ( abs(u1_r) > abs(slope) ) then
    u1_r = slope
  endif

  ! The right edge value in the boundary cell is taken to be the left
  ! edge value in the neighboring cell
  u0_r = edge_values(i1,1)

  ! Given the right edge value and slope, we determine the left
  ! edge value and slope by computing the parabola as determined by
  ! the right edge value and slope and the boundary cell average
  u0_l = 3.0 * u0 + 0.5 * u1_r - 2.0 * u0_r

  ! Apply the traditional PPM limiter
  exp1 = (u0_r - u0_l) * (u0 - 0.5*(u0_l+u0_r))
  exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0

  if ( exp1 > exp2 ) then
    u0_l = 3.0 * u0 - 2.0 * u0_r
  endif

  if ( exp1 < -exp2 ) then
    u0_r = 3.0 * u0 - 2.0 * u0_l
  endif

  edge_values(i0,1) = u0_l
  edge_values(i0,2) = u0_r

  a = u0_l
  b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r
  c = 3.0 * ( u0_r + u0_l - 2.0 * u0 )

  ppoly_coef(i0,1) = a
  ppoly_coef(i0,2) = b
  ppoly_coef(i0,3) = c

  ! ----- Right boundary -----
  i0 = n-1
  i1 = n
  h0 = h(i0)
  h1 = h(i1)
  u0 = u(i0)
  u1 = u(i1)

  ! Compute the right edge slope in neighboring cell and express it in
  ! the global coordinate system
  b = ppoly_coef(i0,2)
  c = ppoly_coef(i0,3)
  u1_l = (b + 2*c)                  ! derivative evaluated at xi = 1.0
  u1_l = u1_l * ((h1+hneglect)/(h0+hneglect))

  ! Limit the left slope by the PLM limited slope
  slope = 2.0 * ( u1 - u0 )
  if ( abs(u1_l) > abs(slope) ) then
    u1_l = slope
  endif

  ! The left edge value in the boundary cell is taken to be the right
  ! edge value in the neighboring cell
  u0_l = edge_values(i0,2)

  ! Given the left edge value and slope, we determine the right
  ! edge value and slope by computing the parabola as determined by
  ! the left edge value and slope and the boundary cell average
  u0_r = 3.0 * u1 - 0.5 * u1_l - 2.0 * u0_l

  ! Apply the traditional PPM limiter
  exp1 = (u0_r - u0_l) * (u1 - 0.5*(u0_l+u0_r))
  exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0

  if ( exp1 > exp2 ) then
    u0_l = 3.0 * u1 - 2.0 * u0_r
  endif

  if ( exp1 < -exp2 ) then
    u0_r = 3.0 * u1 - 2.0 * u0_l
  endif

  edge_values(i1,1) = u0_l
  edge_values(i1,2) = u0_r

  a = u0_l
  b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r
  c = 3.0 * ( u0_r + u0_l - 2.0 * u1 )

  ppoly_coef(i1,1) = a
  ppoly_coef(i1,2) = b
  ppoly_coef(i1,3) = c

end subroutine ppm_boundary_extrapolation

end module ppm_functions