From 97b10008f0765405614cdf529c3a59f732978f02 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Fri, 5 May 2017 07:40:35 -0300
Subject: [PATCH] INTERPOLATIONS: Implement 7th order MP reconstruction.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
 src/interpolations.F90 | 138 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 138 insertions(+)

diff --git a/src/interpolations.F90 b/src/interpolations.F90
index f732fd7..e65884d 100644
--- a/src/interpolations.F90
+++ b/src/interpolations.F90
@@ -264,6 +264,13 @@ module interpolations
       if (verbose .and. ng < 4)                                                &
                   call print_warning("interpolations:initialize_interpolation" &
                          , "Increase the number of ghost cells (at least 4).")
+    case ("mp7", "MP7")
+      name_rec           =  "7th order Monotonicity Preserving"
+      interfaces         => interfaces_dir
+      reconstruct_states => reconstruct_mp7
+      if (verbose .and. ng < 4)                                                &
+                  call print_warning("interpolations:initialize_interpolation" &
+                         , "Increase the number of ghost cells (at least 4).")
     case ("crmp5", "CRMP5")
       name_rec           =  "5th order Compact Monotonicity Preserving"
       interfaces         => interfaces_dir
@@ -3496,6 +3503,137 @@ module interpolations
 !
 !===============================================================================
 !
+! subroutine RECONSTRUCT_MP7:
+! --------------------------
+!
+!   Subroutine reconstructs the interface states using the seventh order
+!   Monotonicity Preserving (MP) method.
+!
+!   Arguments are described in subroutine reconstruct().
+!
+!   References:
+!
+!     [1] Suresh, A. & Huynh, H. T.,
+!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
+!          Time Stepping"
+!         Journal on Computational Physics,
+!         1997, vol. 136, pp. 83-99,
+!         http://dx.doi.org/10.1006/jcph.1997.5745
+!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
+!         "A 5th order monotonicity-preserving upwind compact difference
+!          scheme",
+!         Science China Physics, Mechanics and Astronomy,
+!         Volume 54, Issue 3, pp. 511-522,
+!         http://dx.doi.org/10.1007/s11433-010-4220-x
+!
+!===============================================================================
+!
+  subroutine reconstruct_mp7(n, h, fc, fl, fr)
+
+! local variables are not implicit by default
+!
+    implicit none
+
+! subroutine arguments
+!
+    integer                   , intent(in)  :: n
+    real(kind=8)              , intent(in)  :: h
+    real(kind=8), dimension(n), intent(in)  :: fc
+    real(kind=8), dimension(n), intent(out) :: fl, fr
+
+! local variables
+!
+    integer :: i
+
+! local arrays for derivatives
+!
+    real(kind=8), dimension(n) :: fi
+    real(kind=8), dimension(n) :: u
+
+! local parameters
+!
+    real(kind=8), dimension(7), parameter ::                                   &
+                                ce7 = (/-3.0d+00, 2.5d+01,-1.01d+02, 3.19d+02  &
+                                      , 2.14d+02,-3.8d+01, 4.0d+00 /) / 4.2d+02
+    real(kind=8), dimension(5), parameter ::                                   &
+                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
+                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
+    real(kind=8), dimension(3), parameter :: &
+                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
+    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
+!
+!-------------------------------------------------------------------------------
+!
+!! === left-side interpolation ===
+!!
+! reconstruct the interface state using the 5th order interpolation
+!
+    do i = 4, n - 3
+      u(i) = sum(ce7(:) * fc(i-3:i+3))
+    end do
+
+! interpolate the interface state of the ghost zones using the interpolations
+! of lower orders
+!
+    u(  1) = sum(ce2(:) * fc(  1:  2))
+    u(  2) = sum(ce3(:) * fc(  1:  3))
+    u(  3) = sum(ce5(:) * fc(  1:  5))
+    u(n-2) = sum(ce5(:) * fc(n-4:  n))
+    u(n-1) = sum(ce3(:) * fc(n-2:  n))
+    u(n  ) =              fc(n      )
+
+! apply the monotonicity preserving limiting
+!
+    call mp_limiting(n, fc(1:n), u(1:n))
+
+! copy the interpolation to the respective vector
+!
+    fl(1:n) = u(1:n)
+
+!! === right-side interpolation ===
+!!
+! invert the cell-centered value vector
+!
+    fi(1:n) = fc(n:1:-1)
+
+! reconstruct the interface state using the 5th order interpolation
+!
+    do i = 4, n - 3
+      u(i) = sum(ce7(:) * fi(i-3:i+3))
+    end do
+
+! interpolate the interface state of the ghost zones using the interpolations
+! of lower orders
+!
+    u(  1) = sum(ce2(:) * fi(  1:  2))
+    u(  2) = sum(ce3(:) * fi(  1:  3))
+    u(  3) = sum(ce5(:) * fi(  1:  5))
+    u(n-2) = sum(ce5(:) * fi(n-4:  n))
+    u(n-1) = sum(ce3(:) * fi(n-2:  n))
+    u(n  ) =              fi(n      )
+
+! apply the monotonicity preserving limiting
+!
+    call mp_limiting(n, fi(1:n), u(1:n))
+
+! copy the interpolation to the respective vector
+!
+    fr(1:n-1) = u(n-1:1:-1)
+
+! update the interpolation of the first and last points
+!
+    i     = n - 1
+    fl(1) = 0.5d+00 * (fc(1) + fc(2))
+    fr(i) = 0.5d+00 * (fc(i) + fc(n))
+    fl(n) = fc(n)
+    fr(n) = fc(n)
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine reconstruct_mp7
+!
+!===============================================================================
+!
 ! subroutine RECONSTRUCT_CRMP5:
 ! ----------------------------
 !