Merge branch 'master' into reconnection

sources/interpolations.F90

@@ -381,6 +381,12 @@ module interpolations
       reconstruct_states => reconstruct_ocmp5
       order              = 5
       nghosts            = max(nghosts, 4)
+    case ("ocmp7", "OCMP7")
+      name_rec           =  "7th order Optimized Compact Monotonicity Preserving"
+      interfaces         => interfaces_dir
+      reconstruct_states => reconstruct_ocmp7
+      order              = 7
+      nghosts            = max(nghosts, 4)
     case ("gp", "GP")
       write(stmp, '(f16.1)') sgp
       write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp      &
@@ -4489,30 +4495,30 @@ module interpolations
     real(kind=8), dimension(size(fc)) :: r
     real(kind=8), dimension(size(fc)) :: u
 
+    real(kind=8), parameter :: a1 = 5.0163016d-01
+    real(kind=8), parameter :: a2 = 2.5394716d-01
     real(kind=8), dimension(3), parameter ::                                   &
-            di  = [ 5.0163016d-01, 1.0d+00, 2.5394716d-01 ]
+            di5 = [ a1, 1.0d+00, a2 ]
     real(kind=8), dimension(5), parameter ::                                   &
-            ci5 = [- 3.0d+00 * di(1) - 3.0d+00 * di(3) + 2.0d+00,              &
+            ci5 = [- 3.0d+00 * a1 - 3.0d+00 * a2 + 2.0d+00,                    &
-                     2.7d+01 * di(1) + 1.7d+01 * di(3) - 1.3d+01,              &
+                     2.7d+01 * a1 + 1.7d+01 * a2 - 1.3d+01,                    &
-                     4.7d+01 * di(1) - 4.3d+01 * di(3) + 4.7d+01,              &
+                     4.7d+01 * a1 - 4.3d+01 * a2 + 4.7d+01,                    &
-                   - 1.3d+01 * di(1) + 7.7d+01 * di(3) + 2.7d+01,              &
+                   - 1.3d+01 * a1 + 7.7d+01 * a2 + 2.7d+01,                    &
-                     2.0d+00 * di(1) + 1.2d+01 * di(3) - 3.0d+00 ] / 6.0d+01
+                     2.0d+00 * a1 + 1.2d+01 * a2 - 3.0d+00 ] / 6.0d+01
 !
 !-------------------------------------------------------------------------------
 !
     n = size(fc)
 
-! prepare the diagonals of the tridiagonal matrix
-!
     do i = 1, ng
       a(i) = 0.0d+00
       b(i) = 1.0d+00
       c(i) = 0.0d+00
     end do
     do i = ng + 1, n - ng - 1
-      a(i) = di(1)
+      a(i) = di5(1)
-      b(i) = di(2)
+      b(i) = di5(2)
-      c(i) = di(3)
+      c(i) = di5(3)
     end do
     do i = n - ng, n
       a(i) = 0.0d+00
@@ -4522,14 +4528,10 @@ module interpolations
 
 !! === left-side interpolation ===
 !!
-! prepare the right-hand side of the linear system
-!
     do i = ng, n - ng + 1
       r(i) = sum(ci5(:) * fc(i-2:i+2))
     end do
 
-! use explicit methods for ghost zones
-!
     r(  1) = sum(ce2(:) * fc(  1:  2))
     r(  2) = sum(ce3(:) * fc(  1:  3))
     do i = 3, ng
@@ -4541,32 +4543,20 @@ module interpolations
     r(n-1) = sum(ce3(:) * fc(n-2:  n))
     r(n  ) =              fc(n      )
 
-! solve the tridiagonal system of equations
-!
     call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
 
-! apply the monotonicity preserving limiter
-!
     call mp_limiting(fc(:), u(:))
 
-! return the interpolated values of the left state
-!
     fl(1:n) = u(1:n)
 
 !! === right-side interpolation ===
 !!
-! invert the cell-centered integrals
-!
     fi(1:n) = fc(n:1:-1)
 
-! prepare the right-hand side of the linear system
-!
     do i = ng, n - ng + 1
       r(i) = sum(ci5(:) * fi(i-2:i+2))
     end do ! i = ng, n - ng + 1
 
-! use explicit methods for ghost zones
-!
     r(  1) = sum(ce2(:) * fi(  1:  2))
     r(  2) = sum(ce3(:) * fi(  1:  3))
     do i = 3, ng
@@ -4578,20 +4568,12 @@ module interpolations
     r(n-1) = sum(ce3(:) * fi(n-2:  n))
     r(n  ) =              fi(n      )
 
-! solve the tridiagonal system of equations
-!
     call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
 
-! apply the monotonicity preserving limiter
-!
     call mp_limiting(fi(:), u(:))
 
-! return the interpolated values of the right state
-!
     fr(1:n-1) = u(n-1:1:-1)
 
-! update the extremum points
-!
     i     = n - 1
     fl(1) = 0.5d+00 * (fc(1) + fc(2))
     fr(i) = 0.5d+00 * (fc(i) + fc(n))
@@ -4604,6 +4586,147 @@ module interpolations
 !
 !===============================================================================
 !
+! subroutine RECONSTRUCT_OCMP7:
+! -----------------------------
+!
+!   Subroutine reconstructs the interface states using the 7th order Optimized
+!   Compact Reconstruction Monotonicity Preserving (CRMP) method.
+!
+!   Arguments are described in subroutine reconstruct().
+!
+!   References:
+!
+!     [1] Myeong-Hwan Ahn, Duck-Joo Lee,
+!         "Modified Monotonicity Preserving Constraints for High-Resolution
+!          Optimized Compact Scheme",
+!         Journal of Scientific Computing,
+!         2020, vol. 83, p. 34
+!         https://doi.org/10.1007/s10915-020-01221-0
+!
+!===============================================================================
+!
+  subroutine reconstruct_ocmp7(h, fc, fl, fr)
+
+    use algebra, only : pentadiag
+
+    implicit none
+
+    real(kind=8)              , intent(in)  :: h
+    real(kind=8), dimension(:), intent(in)  :: fc
+    real(kind=8), dimension(:), intent(out) :: fl, fr
+
+    integer :: n, i
+
+    real(kind=8), dimension(size(fc)) :: fi
+    real(kind=8), dimension(size(fc)) :: e, c, d, a, b
+    real(kind=8), dimension(size(fc)) :: r
+    real(kind=8), dimension(size(fc)) :: u
+
+    real(kind=8), parameter :: a1 = 6.6850691831375709863684029643567d-01
+    real(kind=8), parameter :: a2 = 3.3644225201902153852210572056440d-01
+    real(kind=8), dimension(5), parameter ::                                   &
+            di7 = [ (3.0d+00 * a1 + 2.0d+00 * a2 - 2.0d+00) / 8.0d+00,         &
+                    a1, 1.0d+00, a2,                                           &
+                    (          a1 + 6.0d+00 * a2 - 2.0d+00) / 4.0d+01 ]
+    real(kind=8), dimension(5), parameter ::                                   &
+            ci7 = [ 1.80d+01 * a1 + 1.80d+01 * a2 - 1.60d+01,                  &
+                    5.43d+02 * a1 + 2.68d+02 * a2 - 2.91d+02,                  &
+                    3.43d+02 * a1 - 3.32d+02 * a2 + 5.09d+02,                  &
+                   -1.07d+02 * a1 + 5.68d+02 * a2 + 3.09d+02,                  &
+                    4.30d+01 * a1 + 3.18d+02 * a2 - 9.10d+01 ] / 6.0d+02
+!
+!-------------------------------------------------------------------------------
+!
+    n = size(fc)
+
+! prepare the diagonals of the tridiagonal matrix
+!
+    do i = 1, ng
+      e(i) = 0.0d+00
+      c(i) = 0.0d+00
+      d(i) = 1.0d+00
+      a(i) = 0.0d+00
+      b(i) = 0.0d+00
+    end do
+    do i = ng + 1, n - ng - 1
+      e(i) = di7(1)
+      c(i) = di7(2)
+      d(i) = di7(3)
+      a(i) = di7(4)
+      b(i) = di7(5)
+    end do
+    do i = n - ng, n
+      e(i) = 0.0d+00
+      c(i) = 0.0d+00
+      d(i) = 1.0d+00
+      a(i) = 0.0d+00
+      b(i) = 0.0d+00
+    end do
+
+!! === left-side interpolation ===
+!!
+    do i = ng, n - ng + 1
+      r(i) = sum(ci7(:) * fc(i-2:i+2))
+    end do
+
+    r(  1) = sum(ce2(:) * fc(  1:  2))
+    r(  2) = sum(ce3(:) * fc(  1:  3))
+    r(  3) = sum(ce5(:) * fc(  1:  5))
+    do i = 4, ng
+      r(i) = sum(ce7(:) * fc(i-3:i+3))
+    end do
+    do i = n - ng, n - 3
+      r(i) = sum(ce7(:) * fc(i-3:i+3))
+    end do
+    r(n-2) = sum(ce5(:) * fc(n-4:  n))
+    r(n-1) = sum(ce3(:) * fc(n-2:  n))
+    r(n  ) =              fc(n      )
+
+    call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
+
+    call mp_limiting(fc(:), u(:))
+
+    fl(1:n) = u(1:n)
+
+!! === right-side interpolation ===
+!!
+    fi(1:n) = fc(n:1:-1)
+
+    do i = ng, n - ng + 1
+      r(i) = sum(ci7(:) * fi(i-2:i+2))
+    end do ! i = ng, n - ng + 1
+
+    r(  1) = sum(ce2(:) * fi(  1:  2))
+    r(  2) = sum(ce3(:) * fi(  1:  3))
+    r(  3) = sum(ce5(:) * fi(  1:  5))
+    do i = 4, ng
+      r(i) = sum(ce7(:) * fi(i-3:i+3))
+    end do
+    do i = n - ng, n - 3
+      r(i) = sum(ce7(:) * fi(i-3:i+3))
+    end do
+    r(n-2) = sum(ce5(:) * fi(n-4:  n))
+    r(n-1) = sum(ce3(:) * fi(n-2:  n))
+    r(n  ) =              fi(n      )
+
+    call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
+
+    call mp_limiting(fi(:), u(:))
+
+    fr(1:n-1) = u(n-1:1:-1)
+
+    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_ocmp7
+!
+!===============================================================================
+!
 ! subroutine PREPARE_GP:
 ! ---------------------
 !