regrid_interp.F90

Go to the documentation of this file.

module regrid_interp

use mom_error_handler, only : mom_error, fatal
use mom_string_functions, only : uppercase

use regrid_edge_values, only : edge_values_explicit_h2, edge_values_explicit_h4
use regrid_edge_values, only : edge_values_implicit_h4, edge_values_implicit_h6
use regrid_edge_slopes, only : edge_slopes_implicit_h3, edge_slopes_implicit_h5

use plm_functions, only : plm_reconstruction, plm_boundary_extrapolation
use ppm_functions, only : ppm_reconstruction, ppm_boundary_extrapolation
use pqm_functions, only : pqm_reconstruction, pqm_boundary_extrapolation_v1

use p1m_functions, only : p1m_interpolation, p1m_boundary_extrapolation
use p3m_functions, only : p3m_interpolation, p3m_boundary_extrapolation

implicit none ; private

type, public :: interp_cs_type
  private

  !> The following parameter is only relevant when used with the target
  !! interface densities regridding scheme. It indicates which interpolation
  !! to use to determine the grid.
  integer :: interpolation_scheme = -1

  !> Indicate whether high-order boundary extrapolation should be used within
  !! boundary cells
  logical :: boundary_extrapolation
end type interp_cs_type

public regridding_set_ppolys, interpolate_grid
public build_and_interpolate_grid
public set_interp_scheme, set_interp_extrap

! List of interpolation schemes
integer, parameter :: interpolation_p1m_h2     = 0 !< O(h^2)
integer, parameter :: interpolation_p1m_h4     = 1 !< O(h^2)
integer, parameter :: interpolation_p1m_ih4    = 2 !< O(h^2)
integer, parameter :: interpolation_plm        = 3 !< O(h^2)
integer, parameter :: interpolation_ppm_h4     = 4 !< O(h^3)
integer, parameter :: interpolation_ppm_ih4    = 5 !< O(h^3)
integer, parameter :: interpolation_p3m_ih4ih3 = 6 !< O(h^4)
integer, parameter :: interpolation_p3m_ih6ih5 = 7 !< O(h^4)
integer, parameter :: interpolation_pqm_ih4ih3 = 8 !< O(h^4)
integer, parameter :: interpolation_pqm_ih6ih5 = 9 !< O(h^5)

!> List of interpolant degrees
integer, parameter :: degree_1 = 1, degree_2 = 2, degree_3 = 3, degree_4 = 4
integer, public, parameter :: degree_max = 5

!> When the N-R algorithm produces an estimate that lies outside [0,1], the
!! estimate is set to be equal to the boundary location, 0 or 1, plus or minus
!! an offset, respectively, when the derivative is zero at the boundary.
real, public, parameter    :: nr_offset = 1e-6
!> Maximum number of Newton-Raphson iterations. Newton-Raphson iterations are
!! used to build the new grid by finding the coordinates associated with
!! target densities and interpolations of degree larger than 1.
integer, public, parameter :: nr_iterations = 8
!> Tolerance for Newton-Raphson iterations (stop when increment falls below this)
real, public, parameter    :: nr_tolerance = 1e-12

contains

!> Given the set of target values and cell densities, this routine
!! builds an interpolated profile for the densities within each grid cell.
!! It may happen that, given a high-order interpolator, the number of
!! available layers is insufficient (e.g., there are two available layers for
!! a third-order PPM ih4 scheme). In these cases, we resort to the simplest
!! continuous linear scheme (P1M h2).
subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, &
     ppoly0_coefficients, degree)
  type(interp_cs_type),intent(in)    :: CS !< Interpolation control structure
  real, dimension(:),  intent(in)    :: densities !< Actual cell densities
  integer,             intent(in)    :: n0 !< Number of cells on source grid
  real, dimension(:),  intent(in)    :: h0 !< cell widths on source grid
  real, dimension(:,:),intent(inout) :: ppoly0_E            !< Edge value of polynomial
  real, dimension(:,:),intent(inout) :: ppoly0_S            !< Edge slope of polynomial
  real, dimension(:,:),intent(inout) :: ppoly0_coefficients !< Coefficients of polynomial
  integer,             intent(inout) :: degree  !< The degree of the polynomials

  logical :: extrapolate

  ! Reset piecewise polynomials
  ppoly0_e(:,:) = 0.0
  ppoly0_s(:,:) = 0.0
  ppoly0_coefficients(:,:) = 0.0

  extrapolate = cs%boundary_extrapolation

  ! Compute the interpolated profile of the density field and build grid
  select case (cs%interpolation_scheme)

    case ( interpolation_p1m_h2 )
      degree = degree_1
      call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
      call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
      if (extrapolate) then
        call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
      end if

    case ( interpolation_p1m_h4 )
      degree = degree_1
      if ( n0 >= 4 ) then
        call edge_values_explicit_h4( n0, h0, densities, ppoly0_e )
      else
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
      end if
      call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
      if (extrapolate) then
        call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
      end if

    case ( interpolation_p1m_ih4 )
      degree = degree_1
      if ( n0 >= 4 ) then
        call edge_values_implicit_h4( n0, h0, densities, ppoly0_e )
      else
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
      end if
      call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
      if (extrapolate) then
        call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
      end if

    case ( interpolation_plm )
      degree = degree_1
      call plm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
      if (extrapolate) then
        call plm_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
      end if

    case ( interpolation_ppm_h4 )
      if ( n0 >= 4 ) then
        degree = degree_2
        call edge_values_explicit_h4( n0, h0, densities, ppoly0_e )
        call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        if (extrapolate) then
          call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        end if
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        end if
      end if

    case ( interpolation_ppm_ih4 )
      if ( n0 >= 4 ) then
        degree = degree_2
        call edge_values_implicit_h4( n0, h0, densities, ppoly0_e )
        call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        if (extrapolate) then
          call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        end if
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        end if
      end if

    case ( interpolation_p3m_ih4ih3 )
      if ( n0 >= 4 ) then
        degree = degree_3
        call edge_values_implicit_h4( n0, h0, densities, ppoly0_e )
        call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s )
        call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, ppoly0_coefficients )
        if (extrapolate) then
          call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, ppoly0_coefficients )
        end if
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        end if
      end if

    case ( interpolation_p3m_ih6ih5 )
      if ( n0 >= 6 ) then
        degree = degree_3
        call edge_values_implicit_h6( n0, h0, densities, ppoly0_e )
        call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s )
        call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, ppoly0_coefficients )
        if (extrapolate) then
          call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, ppoly0_coefficients )
        end if
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        end if
      end if

    case ( interpolation_pqm_ih4ih3 )
      if ( n0 >= 4 ) then
        degree = degree_4
        call edge_values_implicit_h4( n0, h0, densities, ppoly0_e )
        call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s )
        call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, ppoly0_coefficients )
        if (extrapolate) then
          call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, ppoly0_coefficients )
        end if
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        end if
      end if

    case ( interpolation_pqm_ih6ih5 )
      if ( n0 >= 6 ) then
        degree = degree_4
        call edge_values_implicit_h6( n0, h0, densities, ppoly0_e )
        call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s )
        call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, ppoly0_coefficients )
        if (extrapolate) then
          call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, ppoly0_coefficients )
        end if
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefficients )
        end if
      end if
  end select
end subroutine regridding_set_ppolys

!> Given target values (e.g., density), build new grid based on polynomial
!!
!! Given the grid 'grid0' and the piecewise polynomial interpolant
!! 'ppoly0' (possibly discontinuous), the coordinates of the new grid 'grid1'
!! are determined by finding the corresponding target interface densities.
subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefficients, target_values, degree, n1, h1, x1 )
  ! Arguments
  integer,            intent(in)    :: n0                  !< Number of points on source grid
  real, dimension(:), intent(in)    :: h0                  !< Thicknesses of source grid cells
  real, dimension(:), intent(in)    :: x0                  !< Source interface positions
  real, dimension(:,:), intent(in)  :: ppoly0_E            !< Edge values of interpolating polynomials
  real, dimension(:,:), intent(in)  :: ppoly0_coefficients !< Coefficients of interpolating polynomials
  real, dimension(:), intent(in)    :: target_values       !< Target values of interfaces
  integer,            intent(in)    :: degree              !< Degree of interpolating polynomials
  integer,            intent(in)    :: n1                  !< Number of points on target grid
  real, dimension(:), intent(inout) :: h1                  !< Thicknesses of target grid cells
  real, dimension(:), intent(inout) :: x1                  !< Target interface positions

  ! Local variables
  integer        :: k ! loop index
  real           :: t ! current interface target density

  ! Make sure boundary coordinates of new grid coincide with boundary
  ! coordinates of previous grid
  x1(1) = x0(1)
  x1(n1+1) = x0(n0+1)

  ! Find coordinates for interior target values
  do k = 2,n1
    t = target_values(k)
    x1(k) = get_polynomial_coordinate( n0, h0, x0, ppoly0_e, ppoly0_coefficients, t, degree )
    h1(k-1) = x1(k) - x1(k-1)
  end do
  h1(n1) = x1(n1+1) - x1(n1)

end subroutine interpolate_grid

subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, n1, h1, x1)
  type(interp_cs_type), intent(in) :: CS
  real, dimension(:), intent(in) :: densities, target_values
  integer, intent(in) :: n0, n1
  real, dimension(:), intent(in) :: h0, x0
  real, dimension(:), intent(inout) :: h1, x1

  real, dimension(n0,2) :: ppoly0_E, ppoly0_S
  real, dimension(n0,DEGREE_MAX+1) :: ppoly0_C
  integer :: degree

  call regridding_set_ppolys(cs, densities, n0, h0, ppoly0_e, ppoly0_s, ppoly0_c, &
       degree)
  call interpolate_grid(n0, h0, x0, ppoly0_e, ppoly0_c, target_values, degree, &
       n1, h1, x1)
end subroutine build_and_interpolate_grid

!> Given a target value, find corresponding coordinate for given polynomial
!!
!! Here, 'ppoly' is assumed to be a piecewise discontinuous polynomial of degree
!! 'degree' throughout the domain defined by 'grid'. A target value is given
!! and we need to determine the corresponding grid coordinate to define the
!! new grid.
!!
!! If the target value is out of range, the grid coordinate is simply set to
!! be equal to one of the boundary coordinates, which results in vanished layers
!! near the boundaries.
!!
!! IT IS ASSUMED THAT THE PIECEWISE POLYNOMIAL IS MONOTONICALLY INCREASING.
!! IF THIS IS NOT THE CASE, THE NEW GRID MAY BE ILL-DEFINED.
!!
!! It is assumed that the number of cells defining 'grid' and 'ppoly' are the
!! same.
function get_polynomial_coordinate ( N, h, x_g, ppoly_E, ppoly_coefficients, &
                                     target_value, degree ) result ( x_tgt )
  ! Arguments
  integer,              intent(in) :: N                  !< Number of grid cells
  real, dimension(:),   intent(in) :: h                  !< Grid cell thicknesses
  real, dimension(:),   intent(in) :: x_g                !< Grid interface locations
  real, dimension(:,:), intent(in) :: ppoly_E            !< Edge values of interpolating polynomials
  real, dimension(:,:), intent(in) :: ppoly_coefficients !< Coefficients of interpolating polynomials
  real,                 intent(in) :: target_value       !< Target value to find position for
  integer,              intent(in) :: degree             !< Degree of the interpolating polynomials

  real :: x_tgt !< The position of x_g at which target_value is found.

  ! Local variables
  integer                     :: i, k        ! loop indices
  integer                     :: k_found     ! index of target cell
  integer                     :: iter
  real                        :: xi0         ! normalized target coordinate
  real, dimension(DEGREE_MAX) :: a           ! polynomial coefficients
  real                        :: numerator
  real                        :: denominator
  real                        :: delta       ! Newton-Raphson increment
  real                        :: x           ! global target coordinate
  real                        :: eps         ! offset used to get away from
                                             ! boundaries
  real                        :: grad        ! gradient during N-R iterations

  eps = nr_offset
  k_found = -1

  ! If the target value is outside the range of all values, we
  ! force the target coordinate to be equal to the lowest or
  ! largest value, depending on which bound is overtaken
  if ( target_value <= ppoly_e(1,1) ) then
    x_tgt = x_g(1)
    return  ! return because there is no need to look further
  end if

  ! Since discontinuous edge values are allowed, we check whether the target
  ! value lies between two discontinuous edge values at interior interfaces
  do k = 2,n
    if ( ( target_value >= ppoly_e(k-1,2) ) .AND. &
      ( target_value <= ppoly_e(k,1) ) ) then
      x_tgt = x_g(k)
      return   ! return because there is no need to look further
      exit
    end if
  end do

  ! If the target value is outside the range of all values, we
  ! force the target coordinate to be equal to the lowest or
  ! largest value, depending on which bound is overtaken
  if ( target_value >= ppoly_e(n,2) ) then
    x_tgt = x_g(n+1)
    return  ! return because there is no need to look further
  end if

  ! At this point, we know that the target value is bounded and does not
  ! lie between discontinuous, monotonic edge values. Therefore,
  ! there is a unique solution. We loop on all cells and find which one
  ! contains the target value. The variable k_found holds the index value
  ! of the cell where the taregt value lies.
  do k = 1,n
    if ( ( target_value > ppoly_e(k,1) ) .AND. &
         ( target_value < ppoly_e(k,2) ) ) then
      k_found = k
      exit
    end if
  end do

  ! At this point, 'k_found' should be strictly positive. If not, this is
  ! a major failure because it means we could not find any target cell
  ! despite the fact that the target value lies between the extremes. It
  ! means there is a major problem with the interpolant. This needs to be
  ! reported.
  if ( k_found == -1 ) then
      write(*,*) target_value, ppoly_e(1,1), ppoly_e(n,2)
      write(*,*) 'Could not find target coordinate in ' //&
                 '"get_polynomial_coordinate". This is caused by an '//&
                 'inconsistent interpolant (perhaps not monotonically '//&
                 'increasing)'
      call mom_error( fatal, 'Aborting execution' )
  end if

  ! Reset all polynomial coefficients to 0 and copy those pertaining to
  ! the found cell
  a(:) = 0.0
  do i = 1,degree+1
    a(i) = ppoly_coefficients(k_found,i)
  end do

  ! Guess value to start Newton-Raphson iterations (middle of cell)
  xi0 = 0.5
  iter = 1
  delta = 1e10

  ! Newton-Raphson iterations
  do
    ! break if converged or too many iterations taken
    if ( ( iter > nr_iterations ) .OR. &
         ( abs(delta) < nr_tolerance ) ) then
      exit
    end if

    numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + &
                a(5)*xi0*xi0*xi0*xi0 - target_value

    denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0

    delta = -numerator / denominator

    xi0 = xi0 + delta

    ! Check whether new estimate is out of bounds. If the new estimate is
    ! indeed out of bounds, we manually set it to be equal to the overtaken
    ! bound with a small offset towards the interior when the gradient of
    ! the function at the boundary is zero (in which case, the Newton-Raphson
    ! algorithm does not converge).
    if ( xi0 < 0.0 ) then
      xi0 = 0.0
      grad = a(2)
      if ( grad == 0.0 ) xi0 = xi0 + eps
    end if

    if ( xi0 > 1.0 ) then
      xi0 = 1.0
      grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5)
      if ( grad == 0.0 ) xi0 = xi0 - eps
    end if

    iter = iter + 1
  end do ! end Newton-Raphson iterations

  x_tgt = x_g(k_found) + xi0 * h(k_found)
end function get_polynomial_coordinate

!> Numeric value of interpolation_scheme corresponding to scheme name
integer function interpolation_scheme(interp_scheme)
  character(len=*), intent(in) :: interp_scheme !< Name of interpolation scheme

  select case ( uppercase(trim(interp_scheme)) )
    case ("P1M_H2");     interpolation_scheme = interpolation_p1m_h2
    case ("P1M_H4");     interpolation_scheme = interpolation_p1m_h4
    case ("P1M_IH2");    interpolation_scheme = interpolation_p1m_ih4
    case ("PLM");        interpolation_scheme = interpolation_plm
    case ("PPM_H4");     interpolation_scheme = interpolation_ppm_h4
    case ("PPM_IH4");    interpolation_scheme = interpolation_ppm_ih4
    case ("P3M_IH4IH3"); interpolation_scheme = interpolation_p3m_ih4ih3
    case ("P3M_IH6IH5"); interpolation_scheme = interpolation_p3m_ih6ih5
    case ("PQM_IH4IH3"); interpolation_scheme = interpolation_pqm_ih4ih3
    case ("PQM_IH6IH5"); interpolation_scheme = interpolation_pqm_ih6ih5
    case default ; call mom_error(fatal, "regrid_interp: "//&
     "Unrecognized choice for INTERPOLATION_SCHEME ("//trim(interp_scheme)//").")
  end select
end function interpolation_scheme

subroutine set_interp_scheme(CS, interp_scheme)
  type(interp_cs_type),  intent(inout) :: CS
  character(len=*), intent(in)    :: interp_scheme

  cs%interpolation_scheme = interpolation_scheme(interp_scheme)
end subroutine set_interp_scheme

subroutine set_interp_extrap(CS, extrapolation)
  type(interp_cs_type), intent(inout) :: CS
  logical,         intent(in)    :: extrapolation

  cs%boundary_extrapolation = extrapolation
end subroutine set_interp_extrap

end module regrid_interp