From a93bf0daea363f81ba609cc550c544c50634e63a Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 08:46:42 -0300
Subject: [PATCH 01/12] INTERPOLATIONS: Determine ci5(:) coefficients using
di5(:) in OCMP5.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 7 +++++--
1 file changed, 5 insertions(+), 2 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index cf8b6e6..600690b 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -5360,8 +5360,11 @@ module interpolations
real(kind=8), dimension(3), parameter :: &
di = [ 5.0163016d-01, 1.0d+00, 2.5394716d-01 ]
real(kind=8), dimension(5), parameter :: &
- ci5 = [-4.44553d-03,8.101861d-02, 9.9428149d-01, &
- 6.6721233d-01, 1.751043d-02 ]
+ ci5 = [- 3.0d+00 * di(1) - 3.0d+00 * di(3) + 2.0d+00, &
+ 2.7d+01 * di(1) + 1.7d+01 * di(3) - 1.3d+01, &
+ 4.7d+01 * di(1) - 4.3d+01 * di(3) + 4.7d+01, &
+ - 1.3d+01 * di(1) + 7.7d+01 * di(3) + 2.7d+01, &
+ 2.0d+00 * di(1) + 1.2d+01 * di(3) - 3.0d+00 ] / 6.0d+01
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 :: &
From a8d85390b4e0a4ec4da9968c0d7257d1d23b0503 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 12:28:17 -0300
Subject: [PATCH 02/12] INTERPOLATIONS: Update coefficients of the MP7LD
method.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 57 ++++++++------------------------------
1 file changed, 11 insertions(+), 46 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 600690b..02bb288 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -4338,66 +4338,48 @@ module interpolations
! 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
+! [2] Jian Fang, Yufeng Yao, Zhaorui Li, Lipeng Lu,
+! "Investigation of low-dissipation monotonicity-preserving scheme
+! for direct numerical simulation of compressible turbulent flows",
+! Computers & Fluids,
+! 2014, vol. 104, pp. 55-72,
+! http://dx.doi.org/10.1016/j.compfluid.2014.07.024
!
!===============================================================================
!
subroutine reconstruct_mp7ld(h, fc, fl, fr)
-! local variables are not implicit by default
-!
implicit none
-! subroutine arguments
-!
real(kind=8) , intent(in) :: h
real(kind=8), dimension(:), intent(in) :: fc
real(kind=8), dimension(:), intent(out) :: fl, fr
-! local variables
-!
integer :: n, i
-! local arrays for derivatives
-!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: u
-! local parameters
-!
real(kind=8), dimension(8), parameter :: &
- ce7 = [-8.00d+00, 6.80d+01,-2.82d+02, 9.22d+02, &
- 6.77d+02,-1.35d+02, 1.9d+01, -1.0d+00 ] / 1.26d+03
+ ce7 = [-3.9d+01, 3.53d+02,-1.579d+03, 5.645d+03, 5.015d+03, &
+ -1.201d+03, 2.27d+02, -2.1d+01 ] / 8.4d+03
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
+ 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
+ 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 ]
+ ce2 = [ 5.0d-01, 5.0d-01 ]
!
!-------------------------------------------------------------------------------
-!
-! get the input vector length
!
n = size(fc)
!! === left-side interpolation ===
!!
-! reconstruct the interface state using the 5th order interpolation
-!
do i = 4, n - 4
u(i) = sum(ce7(:) * fc(i-3:i+4))
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))
@@ -4406,29 +4388,18 @@ module interpolations
u(n-1) = sum(ce3(:) * fc(n-2:n ))
u(n ) = fc(n )
-! apply the monotonicity preserving limiting
-!
call mp_limiting(fc(:), u(:))
-! 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 - 4
u(i) = sum(ce7(:) * fi(i-3:i+4))
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))
@@ -4437,16 +4408,10 @@ module interpolations
u(n-1) = sum(ce3(:) * fi(n-2:n ))
u(n ) = fi(n )
-! apply the monotonicity preserving limiting
-!
call mp_limiting(fi(:), u(:))
-! 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))
From 64cc9794e20e63a5b29da57ab98f7ca3bd281962 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 14:17:55 -0300
Subject: [PATCH 03/12] INTERPOLATIONS: Move explicit interpolation
coefficients.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 133 +++++++------------------------------
1 file changed, 24 insertions(+), 109 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 02bb288..f20508c 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -125,6 +125,21 @@ module interpolations
logical , save :: positivity = .false.
logical , save :: clip = .false.
+! interpolation coefficients
+!
+ real(kind=8), dimension(9), parameter :: &
+ ce9 = [ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03, &
+ 1.375d+03,-3.050d+02, 5.500d+01,-5.000d+00 ] / 2.520d+03
+ real(kind=8), dimension(7), parameter :: &
+ ce7 = [-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, &
+ 2.140d+02,-3.800d+01, 4.000d+00 ] / 4.200d+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 ]
+
! by default everything is private
!
private
@@ -3836,15 +3851,6 @@ module interpolations
!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- 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 /)
!
!-------------------------------------------------------------------------------
!
@@ -3963,18 +3969,6 @@ module interpolations
!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: 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 /)
!
!-------------------------------------------------------------------------------
!
@@ -4097,20 +4091,6 @@ module interpolations
!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- real(kind=8), dimension(9), parameter :: &
- ce9 = (/ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03, &
- 1.375d+03,-3.050d+02, 5.500d+01,-5.000d+00 /) / 2.520d+03
- real(kind=8), dimension(7), parameter :: &
- ce7 = (/-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, 2.140d+02, &
- -3.800d+01, 4.000d+00 /) / 4.200d+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 /)
!
!-------------------------------------------------------------------------------
!
@@ -4241,12 +4221,8 @@ module interpolations
! local parameters
!
real(kind=8), dimension(6), parameter :: &
- ce5 = [ 1.20d+01,-8.10d+01, 3.09d+02, &
+ ce5ld = [ 1.20d+01,-8.10d+01, 3.09d+02, &
2.09d+02,-3.10d+01, 2.00d+00 ] / 4.2d+02
- real(kind=8), dimension(3), parameter :: &
- ce3 = [-1.00d+00, 5.00d+00, 2.00d+00 ] / 6.0d+00
- real(kind=8), dimension(2), parameter :: &
- ce2 = [ 5.0d-01, 5.0d-01 ]
!
!-------------------------------------------------------------------------------
!
@@ -4259,7 +4235,7 @@ module interpolations
! reconstruct the interface state using the 5th order interpolation
!
do i = 3, n - 3
- u(i) = sum(ce5(:) * fc(i-2:i+3))
+ u(i) = sum(ce5ld(:) * fc(i-2:i+3))
end do
! interpolate the interface state of the ghost zones using the interpolations
@@ -4288,7 +4264,7 @@ module interpolations
! reconstruct the interface state using the 5th order interpolation
!
do i = 3, n - 3
- u(i) = sum(ce5(:) * fi(i-2:i+3))
+ u(i) = sum(ce5ld(:) * fi(i-2:i+3))
end do
! interpolate the interface state of the ghost zones using the interpolations
@@ -4361,14 +4337,8 @@ module interpolations
real(kind=8), dimension(size(fc)) :: u
real(kind=8), dimension(8), parameter :: &
- ce7 = [-3.9d+01, 3.53d+02,-1.579d+03, 5.645d+03, 5.015d+03, &
+ ce7ld = [-3.9d+01, 3.53d+02,-1.579d+03, 5.645d+03, 5.015d+03, &
-1.201d+03, 2.27d+02, -2.1d+01 ] / 8.4d+03
- 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 ]
!
!-------------------------------------------------------------------------------
!
@@ -4377,7 +4347,7 @@ module interpolations
!! === left-side interpolation ===
!!
do i = 4, n - 4
- u(i) = sum(ce7(:) * fc(i-3:i+4))
+ u(i) = sum(ce7ld(:) * fc(i-3:i+4))
end do
u( 1) = sum(ce2(:) * fc( 1: 2))
@@ -4397,7 +4367,7 @@ module interpolations
fi(1:n) = fc(n:1:-1)
do i = 4, n - 4
- u(i) = sum(ce7(:) * fi(i-3:i+4))
+ u(i) = sum(ce7ld(:) * fi(i-3:i+4))
end do
u( 1) = sum(ce2(:) * fi( 1: 2))
@@ -4473,18 +4443,9 @@ module interpolations
! local parameters
!
real(kind=8), dimension(10), parameter :: &
- ce9 = [ 4.0000d+01,-4.1500d+02, 2.0450d+03,-6.7150d+03, &
+ ce9ld = [ 4.0000d+01,-4.1500d+02, 2.0450d+03,-6.7150d+03, &
2.0165d+04, 1.5629d+04,-3.6910d+03, 7.4900d+02, &
-9.1000d+01, 4.0000d+00 ] / 2.7720d+04
- real(kind=8), dimension(7), parameter :: &
- ce7 = [-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, 2.140d+02, &
- -3.800d+01, 4.000d+00 ] / 4.200d+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 ]
!
!-------------------------------------------------------------------------------
!
@@ -4497,7 +4458,7 @@ module interpolations
! reconstruct the interface state using the 9th order interpolation
!
do i = 5, n - 5
- u(i) = sum(ce9(:) * fc(i-4:i+5))
+ u(i) = sum(ce9ld(:) * fc(i-4:i+5))
end do
! interpolate the interface state of the ghost zones using the interpolations
@@ -4530,7 +4491,7 @@ module interpolations
! reconstruct the interface state using the 9th order interpolation
!
do i = 5, n - 5
- u(i) = sum(ce9(:) * fi(i-4:i+5))
+ u(i) = sum(ce9ld(:) * fi(i-4:i+5))
end do
! interpolate the interface state of the ghost zones using the interpolations
@@ -4626,12 +4587,6 @@ module interpolations
di = (/ 3.0d-01, 6.0d-01, 1.0d-01 /)
real(kind=8), dimension(3), parameter :: &
ci5 = (/ 1.0d+00, 1.9d+01, 1.0d+01 /) / 3.0d+01
- 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 /)
!
!-------------------------------------------------------------------------------
!
@@ -4807,12 +4762,6 @@ module interpolations
real(kind=8), dimension(4), parameter :: &
ci5 = (/ 3.0d+00, 6.7d+01, 4.9d+01 &
, 1.0d+00 /) / 1.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 /)
!
!-------------------------------------------------------------------------------
!
@@ -4981,15 +4930,6 @@ module interpolations
real(kind=8), dimension(5), parameter :: &
ci = (/-1.0d+00, 1.9d+01, 2.39d+02 &
, 1.59d+02, 4.0d+00 /) / 4.2d+02
- 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.70d+01 &
- , 2.70d+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 /)
!
!-------------------------------------------------------------------------------
!
@@ -5163,16 +5103,6 @@ module interpolations
real(kind=8), dimension(6), parameter :: &
ci = [-2.00d+00, 4.20d+01, 6.37d+02, &
5.57d+02, 2.70d+01,-1.00d+00 ] / 1.26d+03
- 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.70d+01, &
- 2.70d+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 ]
!
!-------------------------------------------------------------------------------
!
@@ -5330,12 +5260,6 @@ module interpolations
4.7d+01 * di(1) - 4.3d+01 * di(3) + 4.7d+01, &
- 1.3d+01 * di(1) + 7.7d+01 * di(3) + 2.7d+01, &
2.0d+00 * di(1) + 1.2d+01 * di(3) - 3.0d+00 ] / 6.0d+01
- 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 ]
!
!-------------------------------------------------------------------------------
!
@@ -5593,15 +5517,6 @@ module interpolations
!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- real(kind=8), dimension(5), parameter :: &
- ce5 = (/ 2.0d+00,-1.3d+01, 4.70d+01 &
- , 2.70d+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 /)
!
!-------------------------------------------------------------------------------
!
From 2e2ff486937c7f95b71f1e72c95ba4658529897e Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 16:20:50 -0300
Subject: [PATCH 04/12] INTERPOLATIONS: Rewrite explicit low-dissipation
methods.
Introduce a parameter cweight to control the weight toward the central
scheme and so the amount of schemes' dissipation.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 233 +++++++++++++------------------------
1 file changed, 80 insertions(+), 153 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index f20508c..f49f022 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -114,6 +114,10 @@ module interpolations
!
real(kind=8), save :: ppm_const = 1.25d+00
+! weight toward central scheme for low-dissipation methods
+!
+ real(kind=8), save :: cweight = 7.0d-01
+
! Gaussian process reconstruction coefficients vector
!
real(kind=8), dimension(:) , allocatable, save :: cgp
@@ -127,6 +131,10 @@ module interpolations
! interpolation coefficients
!
+ real(kind=8), dimension(10), save :: ce9ld
+ real(kind=8), dimension( 8), save :: ce7ld
+ real(kind=8), dimension( 6), save :: ce5ld
+
real(kind=8), dimension(9), parameter :: &
ce9 = [ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03, &
1.375d+03,-3.050d+02, 5.500d+01,-5.000d+00 ] / 2.520d+03
@@ -240,12 +248,17 @@ module interpolations
call get_parameter("kappa" , kappa )
call get_parameter("kbeta" , kbeta )
call get_parameter("ppm_const" , ppm_const )
+ call get_parameter("central_weight" , cweight )
call get_parameter("cfl" , cfl )
! calculate κ = (1 - ν) / ν
!
kappa = min(kappa, (1.0d+00 - cfl) / cfl)
+! correct central weight
+!
+ cweight = max(0.0d+00, min(1.0d+00, cweight))
+
! check ngp
!
if (mod(ngp,2) == 0 .or. ngp < 3) then
@@ -542,6 +555,21 @@ module interpolations
!
ng = nghosts
+! calculate coefficients of the low-dissipation methods
+!
+ ce5ld = [ - cweight +2.000d+00, 5.00d+00 * cweight -1.300d+01, &
+ -1.00d+01 * cweight +4.700d+01, 1.00d+01 * cweight +2.700d+01, &
+ -5.00d+00 * cweight -3.000d+00, cweight ] / 6.00d+01
+ ce7ld = [ 3.00d+00 * cweight -6.000d+00, -2.10d+01 * cweight +5.000d+01, &
+ 6.30d+01 * cweight -2.020d+02, -1.05d+02 * cweight +6.380d+02, &
+ 1.05d+02 * cweight +4.280d+02, -6.30d+01 * cweight -7.600d+01, &
+ 2.10d+01 * cweight +8.000d+00, -3.00d+00 * cweight ] / 8.40d+02
+ ce9ld = [ -2.00d+00 * cweight +4.000d+00, 1.80d+01 * cweight -4.100d+01, &
+ -7.20d+01 * cweight +1.990d+02, 1.68d+02 * cweight -6.410d+02, &
+ -2.52d+02 * cweight +1.879d+03, 2.52d+02 * cweight +1.375d+03, &
+ -1.68d+02 * cweight -3.050d+02, 7.20d+01 * cweight +5.500d+01, &
+ -1.80d+01 * cweight -5.000d+00, 2.00d+00 * cweight ] / 2.52d+03
+
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
@@ -640,6 +668,7 @@ module interpolations
!-------------------------------------------------------------------------------
!
if (verbose) then
+
call print_section(verbose, "Interpolations")
call print_parameter(verbose, "reconstruction" , name_rec )
call print_parameter(verbose, "TVD limiter" , name_tlim)
@@ -660,6 +689,7 @@ module interpolations
else
call print_parameter(verbose, "clip extrema" , "off" )
end if
+ call print_parameter(verbose, "central weight", cweight)
end if
@@ -4188,104 +4218,66 @@ module interpolations
! 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
+! [2] Jian Fang, Yufeng Yao, Zhaorui Li, Lipeng Lu,
+! "Investigation of low-dissipation monotonicity-preserving scheme
+! for direct numerical simulation of compressible turbulent flows",
+! Computers & Fluids,
+! 2014, vol. 104, pp. 55-72,
+! http://dx.doi.org/10.1016/j.compfluid.2014.07.024
!
!===============================================================================
!
subroutine reconstruct_mp5ld(h, fc, fl, fr)
-! local variables are not implicit by default
-!
implicit none
-! subroutine arguments
-!
real(kind=8) , intent(in) :: h
real(kind=8), dimension(:), intent(in) :: fc
real(kind=8), dimension(:), intent(out) :: fl, fr
-! local variables
-!
integer :: n, i
-! local arrays for derivatives
-!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- real(kind=8), dimension(6), parameter :: &
- ce5ld = [ 1.20d+01,-8.10d+01, 3.09d+02, &
- 2.09d+02,-3.10d+01, 2.00d+00 ] / 4.2d+02
!
!-------------------------------------------------------------------------------
-!
-! get the input vector length
!
n = size(fc)
!! === left-side interpolation ===
!!
-! reconstruct the interface state using the 5th order interpolation
-!
do i = 3, n - 3
u(i) = sum(ce5ld(:) * fc(i-2: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(n-2) = sum(ce3(:) * fc(n-3:n-1))
- u(n-1) = sum(ce3(:) * fc(n-2:n ))
- u(n ) = fc(n )
+ u( 1) = sum(ce2(:) * fc( 1:2))
+ u( 2) = sum(ce3(:) * fc( 1:3))
+ 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(fc(:), u(:))
-! 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 = 3, n - 3
u(i) = sum(ce5ld(:) * fi(i-2: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(n-2) = sum(ce3(:) * fi(n-3:n-1))
- u(n-1) = sum(ce3(:) * fi(n-2:n ))
- u(n ) = fi(n )
+ u( 1) = sum(ce2(:) * fi( 1:2))
+ u( 2) = sum(ce3(:) * fi( 1:3))
+ 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(fi(:), u(:))
-! 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))
@@ -4308,18 +4300,7 @@ module interpolations
!
! 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] Jian Fang, Yufeng Yao, Zhaorui Li, Lipeng Lu,
-! "Investigation of low-dissipation monotonicity-preserving scheme
-! for direct numerical simulation of compressible turbulent flows",
-! Computers & Fluids,
-! 2014, vol. 104, pp. 55-72,
-! http://dx.doi.org/10.1016/j.compfluid.2014.07.024
+! - See RECONSTRUCT_MP5LD
!
!===============================================================================
!
@@ -4335,10 +4316,6 @@ module interpolations
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: u
-
- real(kind=8), dimension(8), parameter :: &
- ce7ld = [-3.9d+01, 3.53d+02,-1.579d+03, 5.645d+03, 5.015d+03, &
- -1.201d+03, 2.27d+02, -2.1d+01 ] / 8.4d+03
!
!-------------------------------------------------------------------------------
!
@@ -4350,13 +4327,13 @@ module interpolations
u(i) = sum(ce7ld(:) * fc(i-3:i+4))
end do
- u( 1) = sum(ce2(:) * fc( 1: 2))
- u( 2) = sum(ce3(:) * fc( 1: 3))
- u( 3) = sum(ce5(:) * fc( 1: 5))
- u(n-3) = sum(ce5(:) * fc(n-5:n-1))
- u(n-2) = sum(ce5(:) * fc(n-4:n ))
- u(n-1) = sum(ce3(:) * fc(n-2:n ))
- u(n ) = fc(n )
+ u( 1) = sum(ce2(:) * fc( 1:2))
+ u( 2) = sum(ce3(:) * fc( 1:3))
+ u( 3) = sum(ce5(:) * fc( 1:5))
+ u(n-3) = sum(ce7(:) * fc(n-6:n))
+ u(n-2) = sum(ce5(:) * fc(n-4:n))
+ u(n-1) = sum(ce3(:) * fc(n-2:n))
+ u(n ) = fc(n )
call mp_limiting(fc(:), u(:))
@@ -4370,13 +4347,13 @@ module interpolations
u(i) = sum(ce7ld(:) * fi(i-3:i+4))
end do
- u( 1) = sum(ce2(:) * fi( 1: 2))
- u( 2) = sum(ce3(:) * fi( 1: 3))
- u( 3) = sum(ce5(:) * fi( 1: 5))
- u(n-3) = sum(ce5(:) * fi(n-5:n-1))
- u(n-2) = sum(ce5(:) * fi(n-4:n ))
- u(n-1) = sum(ce3(:) * fi(n-2:n ))
- u(n ) = fi(n )
+ u( 1) = sum(ce2(:) * fi( 1:2))
+ u( 2) = sum(ce3(:) * fi( 1:3))
+ u( 3) = sum(ce5(:) * fi( 1:5))
+ u(n-3) = sum(ce7(:) * fi(n-6:n))
+ u(n-2) = sum(ce5(:) * fi(n-4:n))
+ u(n-1) = sum(ce3(:) * fi(n-2:n))
+ u(n ) = fi(n )
call mp_limiting(fi(:), u(:))
@@ -4404,119 +4381,69 @@ module interpolations
!
! 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
+! - See RECONSTRUCT_MP5LD
!
!===============================================================================
!
subroutine reconstruct_mp9ld(h, fc, fl, fr)
-! local variables are not implicit by default
-!
implicit none
-! subroutine arguments
-!
real(kind=8) , intent(in) :: h
real(kind=8), dimension(:), intent(in) :: fc
real(kind=8), dimension(:), intent(out) :: fl, fr
-! local variables
-!
integer :: n, i
-! local arrays for derivatives
-!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- real(kind=8), dimension(10), parameter :: &
- ce9ld = [ 4.0000d+01,-4.1500d+02, 2.0450d+03,-6.7150d+03, &
- 2.0165d+04, 1.5629d+04,-3.6910d+03, 7.4900d+02, &
- -9.1000d+01, 4.0000d+00 ] / 2.7720d+04
!
!-------------------------------------------------------------------------------
-!
-! get the input vector length
!
n = size(fc)
!! === left-side interpolation ===
!!
-! reconstruct the interface state using the 9th order interpolation
-!
do i = 5, n - 5
u(i) = sum(ce9ld(:) * fc(i-4:i+5))
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( 4) = sum(ce7(:) * fc( 1: 7))
- u(n-4) = sum(ce7(:) * fc(n-7:n-1))
- u(n-3) = sum(ce7(:) * fc(n-6:n ))
- u(n-2) = sum(ce5(:) * fc(n-4:n ))
- u(n-1) = sum(ce3(:) * fc(n-2:n ))
- u(n ) = fc(n )
+ u( 1) = sum(ce2(:) * fc( 1:2))
+ u( 2) = sum(ce3(:) * fc( 1:3))
+ u( 3) = sum(ce5(:) * fc( 1:5))
+ u( 4) = sum(ce7(:) * fc( 1:7))
+ u(n-4) = sum(ce9(:) * fc(n-8:n))
+ u(n-3) = sum(ce7(:) * fc(n-6:n))
+ 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(fc(:), u(:))
-! 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 9th order interpolation
-!
do i = 5, n - 5
u(i) = sum(ce9ld(:) * fi(i-4:i+5))
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( 4) = sum(ce7(:) * fi( 1: 7))
- u(n-4) = sum(ce7(:) * fi(n-7:n-1))
- u(n-3) = sum(ce7(:) * fi(n-6:n ))
- u(n-2) = sum(ce5(:) * fi(n-4:n ))
- u(n-1) = sum(ce3(:) * fi(n-2:n ))
- u(n ) = fi(n )
+ u( 1) = sum(ce2(:) * fi( 1:2))
+ u( 2) = sum(ce3(:) * fi( 1:3))
+ u( 3) = sum(ce5(:) * fi( 1:5))
+ u( 4) = sum(ce7(:) * fi( 1:7))
+ u(n-4) = sum(ce9(:) * fi(n-8:n))
+ u(n-3) = sum(ce7(:) * fi(n-6:n))
+ 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(fi(:), u(:))
-! 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))
From a19e8b885f03693dfe83c5da6065fad9978d7703 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 17:25:19 -0300
Subject: [PATCH 05/12] INTERPOLATIONS: Rewrite compact MP methods.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 346 ++++++++++++++-----------------------
1 file changed, 128 insertions(+), 218 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index f49f022..a352532 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -4483,55 +4483,39 @@ module interpolations
!
subroutine reconstruct_crmp5(h, fc, fl, fr)
-! include external procedures
-!
- use algebra , only : tridiag
+ use algebra, only : tridiag
-! local variables are not implicit by default
-!
implicit none
-! subroutine arguments
-!
real(kind=8) , intent(in) :: h
real(kind=8), dimension(:), intent(in) :: fc
real(kind=8), dimension(:), intent(out) :: fl, fr
-! local variables
-!
integer :: n, i
-! local arrays for derivatives
-!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: a, b, c
real(kind=8), dimension(size(fc)) :: r
real(kind=8), dimension(size(fc)) :: u
-! local parameters
-!
real(kind=8), dimension(3), parameter :: &
- di = (/ 3.0d-01, 6.0d-01, 1.0d-01 /)
+ di5 = [ 3.0d+00, 6.0d+00, 1.0d+00 ] / 6.0d+00
real(kind=8), dimension(3), parameter :: &
- ci5 = (/ 1.0d+00, 1.9d+01, 1.0d+01 /) / 3.0d+01
+ ci5 = [ 1.0d+00, 1.9d+01, 1.0d+01 ] / 1.8d+01
!
!-------------------------------------------------------------------------------
-!
-! get the input vector length
!
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)
- b(i) = di(2)
- c(i) = di(3)
+ a(i) = di5(1)
+ b(i) = di5(2)
+ c(i) = di5(3)
end do
do i = n - ng, n
a(i) = 0.0d+00
@@ -4541,14 +4525,10 @@ module interpolations
!! === left-side interpolation ===
!!
-! prepare the tridiagonal system coefficients for the interior part
-!
do i = ng, n - ng + 1
r(i) = sum(ci5(:) * fc(i-1:i+1))
end do
-! interpolate ghost zones using the explicit method
-!
r( 1) = sum(ce2(:) * fc( 1: 2))
r( 2) = sum(ce3(:) * fc( 1: 3))
do i = 3, ng
@@ -4560,32 +4540,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 limiting
-!
call mp_limiting(fc(:), u(:))
-! 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)
-! prepare the tridiagonal system coefficients for the interior part
-!
do i = ng, n - ng + 1
r(i) = sum(ci5(:) * fi(i-1:i+1))
end do ! i = ng, n - ng + 1
-! interpolate ghost zones using the explicit method
-!
r( 1) = sum(ce2(:) * fi( 1: 2))
r( 2) = sum(ce3(:) * fi( 1: 3))
do i = 3, ng
@@ -4597,20 +4565,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 limiting
-!
call mp_limiting(fi(:), u(:))
-! 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))
@@ -4623,6 +4583,128 @@ module interpolations
!
!===============================================================================
!
+! subroutine RECONSTRUCT_CRMP7:
+! ----------------------------
+!
+! Subroutine reconstructs the interface states using the seventh order
+! Compact Reconstruction Monotonicity Preserving (CRMP) method.
+!
+! Arguments are described in subroutine reconstruct().
+!
+! References:
+!
+! - See RECONSTRUCT_CRMP5
+!
+!===============================================================================
+!
+ subroutine reconstruct_crmp7(h, fc, fl, fr)
+
+ use algebra, only : tridiag
+
+ 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)) :: a, b, c
+ real(kind=8), dimension(size(fc)) :: r
+ real(kind=8), dimension(size(fc)) :: u
+
+! local parameters
+!
+ real(kind=8), dimension(3), parameter :: &
+ di7 = [ 5.0d-01, 1.0d+00, 2.5d-01 ]
+ real(kind=8), dimension(5), parameter :: &
+ ci7 = [-1.0d+00, 1.9d+01, 2.39d+02, 1.59d+02, 4.0d+00 ] / 2.4d+02
+!
+!-------------------------------------------------------------------------------
+!
+ n = size(fc)
+
+ 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) = di7(1)
+ b(i) = di7(2)
+ c(i) = di7(3)
+ end do
+ do i = n - ng, n
+ a(i) = 0.0d+00
+ b(i) = 1.0d+00
+ c(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 tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+
+ 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 tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+
+ 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_crmp7
+!
+!===============================================================================
+!
! subroutine RECONSTRUCT_CRMP5LD:
! ------------------------------
!
@@ -4798,178 +4880,6 @@ module interpolations
!
!===============================================================================
!
-! subroutine RECONSTRUCT_CRMP7:
-! ----------------------------
-!
-! Subroutine reconstructs the interface states using the seventh order
-! Compact Reconstruction Monotonicity Preserving (CRMP) 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_crmp7(h, fc, fl, fr)
-
-! include external procedures
-!
- use algebra , only : tridiag
-
-! local variables are not implicit by default
-!
- implicit none
-
-! subroutine arguments
-!
- real(kind=8) , intent(in) :: h
- real(kind=8), dimension(:), intent(in) :: fc
- real(kind=8), dimension(:), intent(out) :: fl, fr
-
-! local variables
-!
- integer :: n, i
-
-! local arrays for derivatives
-!
- real(kind=8), dimension(size(fc)) :: fi
- real(kind=8), dimension(size(fc)) :: a, b, c
- real(kind=8), dimension(size(fc)) :: r
- real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- real(kind=8), dimension(3), parameter :: &
- di = (/ 2.0d+00, 4.0d+00, 1.0d+00/) / 7.0d+00
- real(kind=8), dimension(5), parameter :: &
- ci = (/-1.0d+00, 1.9d+01, 2.39d+02 &
- , 1.59d+02, 4.0d+00 /) / 4.2d+02
-!
-!-------------------------------------------------------------------------------
-!
-! get the input vector length
-!
- 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)
- b(i) = di(2)
- c(i) = di(3)
- end do
- do i = n - ng, n
- a(i) = 0.0d+00
- b(i) = 1.0d+00
- c(i) = 0.0d+00
- end do
-
-!! === left-side interpolation ===
-!!
-! prepare the tridiagonal system coefficients for the interior part
-!
- do i = ng, n - ng + 1
- r(i) = sum( ci(:) * fc(i-2:i+2))
- end do
-
-! interpolate ghost zones using the explicit method
-!
- 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 )
-
-! 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 limiting
-!
- call mp_limiting(fc(:), u(:))
-
-! 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)
-
-! prepare the tridiagonal system coefficients for the interior part
-!
- do i = ng, n - ng + 1
- r(i) = sum( ci(:) * fi(i-2:i+2))
- end do ! i = ng, n - ng + 1
-
-! interpolate ghost zones using the explicit method
-!
- 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 )
-
-! 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 limiting
-!
- call mp_limiting(fi(:), u(:))
-
-! 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_crmp7
-!
-!===============================================================================
-!
! subroutine RECONSTRUCT_CRMP7LD:
! ------------------------------
!
From 352af9a1a68476dd2df5fd49449b5f9344b38c1d Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 18:19:05 -0300
Subject: [PATCH 06/12] INTERPOLATIONS: Rewrite compact low-dissipation
methods.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 143 +++++++------------------------------
1 file changed, 27 insertions(+), 116 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index a352532..937f11e 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -134,6 +134,9 @@ module interpolations
real(kind=8), dimension(10), save :: ce9ld
real(kind=8), dimension( 8), save :: ce7ld
real(kind=8), dimension( 6), save :: ce5ld
+ real(kind=8), dimension( 3), save :: di5ld, di7ld
+ real(kind=8), dimension( 4), save :: ci5ld
+ real(kind=8), dimension( 6), save :: ci7ld
real(kind=8), dimension(9), parameter :: &
ce9 = [ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03, &
@@ -570,6 +573,17 @@ module interpolations
-1.68d+02 * cweight -3.050d+02, 7.20d+01 * cweight +5.500d+01, &
-1.80d+01 * cweight -5.000d+00, 2.00d+00 * cweight ] / 2.52d+03
+ di5ld = [ -(cweight -3.0d+00) / 6.0d+00, 1.0d+00, &
+ (cweight +1.0d+00) / 6.0d+00 ]
+ ci5ld = [ - cweight +2.0d+00, -9.0d+00 * cweight +3.8d+01, &
+ 9.0d+00 * cweight +2.0d+01, cweight ] / 3.6d+01
+
+ di7ld = [ -(cweight -4.0d+00) / 8.0d+00, 1.0d+00, &
+ (cweight +2.0d+00) / 8.0d+00 ]
+ ci7ld = [ cweight -2.0d+00, -1.50d+01 * cweight +3.8d+01, &
+ -8.0d+01 * cweight +4.78d+02, 8.0d+01 * cweight +3.18d+02, &
+ 1.5d+01 * cweight +8.00d+00, -cweight ] / 4.8d+02
+
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
@@ -4739,56 +4753,34 @@ module interpolations
!
subroutine reconstruct_crmp5ld(h, fc, fl, fr)
-! include external procedures
-!
- use algebra , only : tridiag
+ use algebra, only : tridiag
-! local variables are not implicit by default
-!
implicit none
-! subroutine arguments
-!
real(kind=8) , intent(in) :: h
real(kind=8), dimension(:), intent(in) :: fc
real(kind=8), dimension(:), intent(out) :: fl, fr
-! local variables
-!
integer :: n, i
-! local arrays for derivatives
-!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: a, b, c
real(kind=8), dimension(size(fc)) :: r
real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- real(kind=8), dimension(3), parameter :: &
- di = (/ 2.5d-01, 6.0d-01, 1.5d-01 /)
- real(kind=8), dimension(4), parameter :: &
- ci5 = (/ 3.0d+00, 6.7d+01, 4.9d+01 &
- , 1.0d+00 /) / 1.2d+02
!
!-------------------------------------------------------------------------------
-!
-! get the input vector length
!
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)
- b(i) = di(2)
- c(i) = di(3)
+ a(i) = di5ld(1)
+ b(i) = di5ld(2)
+ c(i) = di5ld(3)
end do
do i = n - ng, n
a(i) = 0.0d+00
@@ -4798,14 +4790,10 @@ module interpolations
!! === left-side interpolation ===
!!
-! prepare the tridiagonal system coefficients for the interior part
-!
do i = ng, n - ng + 1
- r(i) = sum(ci5(:) * fc(i-1:i+2))
+ r(i) = sum(ci5ld(:) * fc(i-1:i+2))
end do
-! interpolate ghost zones using the explicit method
-!
r( 1) = sum(ce2(:) * fc( 1: 2))
r( 2) = sum(ce3(:) * fc( 1: 3))
do i = 3, ng
@@ -4817,32 +4805,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 limiting
-!
call mp_limiting(fc(:), u(:))
-! 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)
-! prepare the tridiagonal system coefficients for the interior part
-!
do i = ng, n - ng + 1
- r(i) = sum(ci5(:) * fi(i-1:i+2))
+ r(i) = sum(ci5ld(:) * fi(i-1:i+2))
end do ! i = ng, n - ng + 1
-! interpolate ghost zones using the explicit method
-!
r( 1) = sum(ce2(:) * fi( 1: 2))
r( 2) = sum(ce3(:) * fi( 1: 3))
do i = 3, ng
@@ -4854,20 +4830,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 limiting
-!
call mp_limiting(fi(:), u(:))
-! 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))
@@ -4891,73 +4859,40 @@ module interpolations
!
! 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
+! - See RECONSTRUCT_CRMP5LD
!
!===============================================================================
!
subroutine reconstruct_crmp7ld(h, fc, fl, fr)
-! include external procedures
-!
- use algebra , only : tridiag
+ use algebra, only : tridiag
-! local variables are not implicit by default
-!
implicit none
-! subroutine arguments
-!
real(kind=8) , intent(in) :: h
real(kind=8), dimension(:), intent(in) :: fc
real(kind=8), dimension(:), intent(out) :: fl, fr
-! local variables
-!
integer :: n, i
-! local arrays for derivatives
-!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: a, b, c
real(kind=8), dimension(size(fc)) :: r
real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- real(kind=8), dimension(3), parameter :: &
- di = [ 5.0d+00, 1.2d+01, 4.0d+00 ] / 2.1d+01
- real(kind=8), dimension(6), parameter :: &
- ci = [-2.00d+00, 4.20d+01, 6.37d+02, &
- 5.57d+02, 2.70d+01,-1.00d+00 ] / 1.26d+03
!
!-------------------------------------------------------------------------------
-!
-! get the input vector length
!
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)
- b(i) = di(2)
- c(i) = di(3)
+ a(i) = di7ld(1)
+ b(i) = di7ld(2)
+ c(i) = di7ld(3)
end do
do i = n - ng, n
a(i) = 0.0d+00
@@ -4967,14 +4902,10 @@ module interpolations
!! === left-side interpolation ===
!!
-! prepare the tridiagonal system coefficients for the interior part
-!
do i = ng, n - ng + 1
- r(i) = sum( ci(:) * fc(i-2:i+3))
+ r(i) = sum(ci7ld(:) * fc(i-2:i+3))
end do
-! interpolate ghost zones using the explicit method
-!
r( 1) = sum(ce2(:) * fc( 1: 2))
r( 2) = sum(ce3(:) * fc( 1: 3))
r( 3) = sum(ce5(:) * fc( 1: 5))
@@ -4988,32 +4919,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 limiting
-!
call mp_limiting(fc(:), u(:))
-! 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)
-! prepare the tridiagonal system coefficients for the interior part
-!
do i = ng, n - ng + 1
- r(i) = sum( ci(:) * fi(i-2:i+3))
+ r(i) = sum(ci7ld(:) * fi(i-2:i+3))
end do ! i = ng, n - ng + 1
-! interpolate ghost zones using the explicit method
-!
r( 1) = sum(ce2(:) * fi( 1: 2))
r( 2) = sum(ce3(:) * fi( 1: 3))
r( 3) = sum(ce5(:) * fi( 1: 5))
@@ -5027,20 +4946,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 limiting
-!
call mp_limiting(fi(:), u(:))
-! 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))
From 2d31bb2a4f8acea32b98e675dc6d6c9eea591880 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 18:31:05 -0300
Subject: [PATCH 07/12] ALGEBRA: Implement pentadiagona linear equations
solver.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/algebra.F90 | 83 ++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 82 insertions(+), 1 deletion(-)
diff --git a/sources/algebra.F90 b/sources/algebra.F90
index 7251e91..97bde8b 100644
--- a/sources/algebra.F90
+++ b/sources/algebra.F90
@@ -51,7 +51,7 @@ module algebra
public :: quadratic, quadratic_normalized
public :: cubic, cubic_normalized
public :: quartic
- public :: tridiag, invert
+ public :: tridiag, pentadiag, invert
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
@@ -1148,6 +1148,87 @@ module algebra
!
!===============================================================================
!
+! subroutine PENTADIAG:
+! --------------------
+!
+! Subroutine solves the system of pentadiagonal linear equations using
+! PTRANS-I algorithm.
+!
+! Arguments:
+!
+! n - the size of the matrix;
+! e(:) - the vector of the outer lower off-diagonal values;
+! c(:) - the vector of the inner lower off-diagonal values;
+! d(:) - the vector of the diagonal values;
+! a(:) - the vector of the inner upper off-diagonal values;
+! b(:) - the vector of the outer upper off-diagonal values;
+! y(:) - the vector of the right side values;
+! x(:) - the solution vector;
+!
+! References:
+!
+! [1] Askar, S. S., and Karawia, A. A.,
+! "On Solving Pentadiagonal Linear Systems via Transformations",
+! Mathematical Problems in Engineering,
+! 2015, Volume 2015, Article ID 232456, 9 pages
+! http://dx.doi.org/10.1155/2015/232456
+!
+!===============================================================================
+!
+ subroutine pentadiag(n, e, c, d, a, b, y, x)
+
+ implicit none
+
+ integer , intent(in) :: n
+ real(kind=8), dimension(n), intent(in) :: e, c, d, a, b, y
+ real(kind=8), dimension(n), intent(out) :: x
+
+ integer :: i
+ real(kind=8), dimension(n) :: al, bt, gm, mu, z
+!
+!-------------------------------------------------------------------------------
+!
+ mu(1) = d(1)
+ al(1) = a(1) / mu(1)
+ bt(1) = b(1) / mu(1)
+ z (1) = y(1) / mu(1)
+
+ gm(2) = c(2)
+ mu(2) = d(2) - al(1) * gm(2)
+ al(2) = (a(2) - bt(1) * gm(2)) / mu(2)
+ bt(2) = b(2) / mu(2)
+ z (2) = (y(2) - z(1) * gm(2)) / mu(2)
+
+ do i = 3, n-2
+ gm(i) = c(i) - al(i-2) * e(i)
+ mu(i) = d(i) - bt(i-2) * e(i) - al(i-1) * gm(i)
+ al(i) = (a(i) - bt(i-1) * gm(i)) / mu(i)
+ bt(i) = b(i) / mu(i)
+ z (i) = (y(i) - z(i-2) * e(i) - z(i-1) * gm(i)) / mu(i)
+ end do
+
+ gm(n-1) = c(n-1) - al(n-3) * e(n-1)
+ mu(n-1) = d(n-1) - bt(n-3) * e(n-1) - al(n-2) * gm(n-1)
+ al(n-1) = (a(n-1) - bt(n-2) * gm(n-1)) / mu(n-1)
+
+ gm(n) = c(n) - al(n-2) * e(n)
+ mu(n) = d(n) - bt(n-3) * e(n) - al(n-1) * gm(n)
+
+ z(n-1) = (y(n-1) - z(n-2) * e(n-1) - z(n-2) * gm(n-1)) / mu(n-1)
+ z(n ) = (y(n) - z(n-1) * e(n) - z(n-1) * gm(n)) / mu(n)
+
+ x(n ) = z(n)
+ x(n-1) = z(n-1) - al(n-1) * x(n)
+ do i = n - 2, 1, -1
+ x(i) = z(i) - al(i) * x(i+1) - bt(i) * x(i+2)
+ end do
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine pentadiag
+!
+!===============================================================================
+!
! subroutine INVERT:
! -----------------
!
From 20bae9dec1525d720c8be3676b10fac318ee7364 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 19:08:21 -0300
Subject: [PATCH 08/12] INTERPOLATIONS: Implement 9th order Compact MP method.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 150 +++++++++++++++++++++++++++++++++++--
1 file changed, 144 insertions(+), 6 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 937f11e..3ec3424 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -380,18 +380,24 @@ module interpolations
reconstruct_states => reconstruct_crmp5
order = 5
nghosts = max(nghosts, 4)
- case ("crmp5l", "crmp5ld", "CRMP5L", "CRMP5LD")
- name_rec = "5th order Low-Dissipation Compact Monotonicity Preserving"
- interfaces => interfaces_dir
- reconstruct_states => reconstruct_crmp5ld
- order = 5
- nghosts = max(nghosts, 4)
case ("crmp7", "CRMP7")
name_rec = "7th order Compact Monotonicity Preserving"
interfaces => interfaces_dir
reconstruct_states => reconstruct_crmp7
order = 7
nghosts = max(nghosts, 4)
+ case ("crmp9", "CRMP9")
+ name_rec = "9th order Compact Monotonicity Preserving"
+ interfaces => interfaces_dir
+ reconstruct_states => reconstruct_crmp9
+ order = 9
+ nghosts = max(nghosts, 6)
+ case ("crmp5l", "crmp5ld", "CRMP5L", "CRMP5LD")
+ name_rec = "5th order Low-Dissipation Compact Monotonicity Preserving"
+ interfaces => interfaces_dir
+ reconstruct_states => reconstruct_crmp5ld
+ order = 5
+ nghosts = max(nghosts, 4)
case ("crmp7l", "crmp7ld", "CRMP7L", "CRMP7LD")
name_rec = "7th order Low-Dissipation Compact Monotonicity Preserving"
interfaces => interfaces_dir
@@ -4719,6 +4725,138 @@ module interpolations
!
!===============================================================================
!
+! subroutine RECONSTRUCT_CRMP9:
+! ----------------------------
+!
+! Subroutine reconstructs the interface states using the 9th order
+! Compact Reconstruction Monotonicity Preserving (CRMP) method.
+!
+! Arguments are described in subroutine reconstruct().
+!
+! References:
+!
+! - See RECONSTRUCT_CRMP5
+!
+!===============================================================================
+!
+ subroutine reconstruct_crmp9(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
+
+! local parameters
+!
+ real(kind=8), dimension(5), parameter :: &
+ di9 = [ 5.0d+00, 4.0d+01, 6.0d+01, 2.0d+01, 1.0d+00 ] / 6.0d+01
+ real(kind=8), dimension(5), parameter :: &
+ ci9 = [ 6.0d+00, 4.81d+02, 1.881d+03, 1.281d+03, 1.31d+02 ] / 1.8d+03
+!
+!-------------------------------------------------------------------------------
+!
+ n = size(fc)
+
+ 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) = di9(1)
+ c(i) = di9(2)
+ d(i) = di9(3)
+ a(i) = di9(4)
+ b(i) = di9(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(ci9(:) * 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))
+ r( 4) = sum(ce7(:) * fc( 1: 7))
+ do i = 5, ng
+ r(i) = sum(ce9(:) * fc(i-4:i+4))
+ end do
+ do i = n - ng, n - 4
+ r(i) = sum(ce9(:) * fc(i-4:i+4))
+ end do
+ r(n-3) = sum(ce7(:) * fc(n-6: n))
+ 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(ci9(:) * 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))
+ r( 4) = sum(ce7(:) * fi( 1: 7))
+ do i = 5, ng
+ r(i) = sum(ce9(:) * fi(i-4:i+4))
+ end do
+ do i = n - ng, n - 4
+ r(i) = sum(ce9(:) * fi(i-4:i+4))
+ end do
+ r(n-3) = sum(ce7(:) * fi(n-6: n))
+ 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_crmp9
+!
+!===============================================================================
+!
! subroutine RECONSTRUCT_CRMP5LD:
! ------------------------------
!
From 662f35ff01f1570eb440ca37880e41b68614417a Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 22:08:55 -0300
Subject: [PATCH 09/12] INTERPOLATIONS: Implement 9th order low-dissipation
Compact MP method.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 146 ++++++++++++++++++++++++++++++++++++-
1 file changed, 143 insertions(+), 3 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 3ec3424..68d14a0 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -135,8 +135,9 @@ module interpolations
real(kind=8), dimension( 8), save :: ce7ld
real(kind=8), dimension( 6), save :: ce5ld
real(kind=8), dimension( 3), save :: di5ld, di7ld
+ real(kind=8), dimension( 5), save :: di9ld
real(kind=8), dimension( 4), save :: ci5ld
- real(kind=8), dimension( 6), save :: ci7ld
+ real(kind=8), dimension( 6), save :: ci7ld, ci9ld
real(kind=8), dimension(9), parameter :: &
ce9 = [ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03, &
@@ -404,6 +405,12 @@ module interpolations
reconstruct_states => reconstruct_crmp7ld
order = 7
nghosts = max(nghosts, 4)
+ case ("crmp9l", "crmp9ld", "CRMP9L", "CRMP9LD")
+ name_rec = "9th order Low-Dissipation Compact Monotonicity Preserving"
+ interfaces => interfaces_dir
+ reconstruct_states => reconstruct_crmp9ld
+ order = 9
+ nghosts = max(nghosts, 6)
case ("ocmp5", "OCMP5")
name_rec = "5th order Optimized Compact Monotonicity Preserving"
interfaces => interfaces_dir
@@ -582,13 +589,20 @@ module interpolations
di5ld = [ -(cweight -3.0d+00) / 6.0d+00, 1.0d+00, &
(cweight +1.0d+00) / 6.0d+00 ]
ci5ld = [ - cweight +2.0d+00, -9.0d+00 * cweight +3.8d+01, &
- 9.0d+00 * cweight +2.0d+01, cweight ] / 3.6d+01
+ 9.0d+00 * cweight +2.0d+01, cweight ] / 3.6d+01
di7ld = [ -(cweight -4.0d+00) / 8.0d+00, 1.0d+00, &
(cweight +2.0d+00) / 8.0d+00 ]
ci7ld = [ cweight -2.0d+00, -1.50d+01 * cweight +3.8d+01, &
-8.0d+01 * cweight +4.78d+02, 8.0d+01 * cweight +3.18d+02, &
- 1.5d+01 * cweight +8.00d+00, -cweight ] / 4.8d+02
+ 1.5d+01 * cweight +8.00d+00, -cweight ] / 4.8d+02
+
+ di9ld = [ -2.0d+00 * cweight +5.0d+00, -1.0d+01 * cweight +4.0d+01, &
+ 6.0d+01, 1.0d+01 * cweight +2.0d+01, &
+ 2.0d+00 * cweight +1.0d+00 ] / 6.0d+01
+ ci9ld = [ -3.00d+00 * cweight +6.000d+00,-1.75d+02 * cweight +4.810d+02, &
+ -3.00d+02 * cweight +1.881d+03, 3.00d+02 * cweight +1.281d+03, &
+ 1.75d+02 * cweight +1.310d+02, 3.00d+00 * cweight ] / 1.8d+03
#ifdef PROFILE
! stop accounting time for module initialization/finalization
@@ -5102,6 +5116,132 @@ module interpolations
!
!===============================================================================
!
+! subroutine RECONSTRUCT_CRMP9LD:
+! ------------------------------
+!
+! Subroutine reconstructs the interface states using the 9th order
+! low dissipation Compact Reconstruction Monotonicity Preserving (CRMP)
+! method.
+!
+! Arguments are described in subroutine reconstruct().
+!
+! References:
+!
+! - See RECONSTRUCT_CRMP5LD
+!
+!===============================================================================
+!
+ subroutine reconstruct_crmp9ld(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
+!
+!-------------------------------------------------------------------------------
+!
+ n = size(fc)
+
+ 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) = di9ld(1)
+ c(i) = di9ld(2)
+ d(i) = di9ld(3)
+ a(i) = di9ld(4)
+ b(i) = di9ld(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(ci9ld(:) * fc(i-2:i+3))
+ end do
+
+ r( 1) = sum(ce2(:) * fc( 1: 2))
+ r( 2) = sum(ce3(:) * fc( 1: 3))
+ r( 3) = sum(ce5(:) * fc( 1: 5))
+ r( 4) = sum(ce7(:) * fc( 1: 7))
+ do i = 5, ng
+ r(i) = sum(ce9(:) * fc(i-4:i+4))
+ end do
+ do i = n - ng, n - 4
+ r(i) = sum(ce9(:) * fc(i-4:i+4))
+ end do
+ r(n-3) = sum(ce7(:) * fc(n-6: n))
+ 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(ci9ld(:) * fi(i-2:i+3))
+ 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))
+ r( 4) = sum(ce7(:) * fi( 1: 7))
+ do i = 5, ng
+ r(i) = sum(ce9(:) * fi(i-4:i+4))
+ end do
+ do i = n - ng, n - 4
+ r(i) = sum(ce9(:) * fi(i-4:i+4))
+ end do
+ r(n-3) = sum(ce7(:) * fi(n-6: n))
+ 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_crmp9ld
+!
+!===============================================================================
+!
! subroutine RECONSTRUCT_OCMP5:
! -----------------------------
!
From 325dc36c44334fa7e6b5fc5d2b89a78526a6d89b Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 22:19:10 -0300
Subject: [PATCH 10/12] INTERPOLATIONS: Merge regular explicit MP methods.
They are identical to low-dissipation versions with center_weight = 0.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 460 +++----------------------------------
1 file changed, 38 insertions(+), 422 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 68d14a0..2172982 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -131,9 +131,9 @@ module interpolations
! interpolation coefficients
!
- real(kind=8), dimension(10), save :: ce9ld
- real(kind=8), dimension( 8), save :: ce7ld
- real(kind=8), dimension( 6), save :: ce5ld
+ real(kind=8), dimension(10), save :: ce9a
+ real(kind=8), dimension( 8), save :: ce7a
+ real(kind=8), dimension( 6), save :: ce5a
real(kind=8), dimension( 3), save :: di5ld, di7ld
real(kind=8), dimension( 5), save :: di9ld
real(kind=8), dimension( 4), save :: ci5ld
@@ -357,24 +357,6 @@ module interpolations
reconstruct_states => reconstruct_mp9
order = 9
nghosts = max(nghosts, 6)
- case ("mp5ld", "MP5LD")
- name_rec = "5th order Low-Dissipation Monotonicity Preserving"
- interfaces => interfaces_dir
- reconstruct_states => reconstruct_mp5ld
- order = 5
- nghosts = max(nghosts, 4)
- case ("mp7ld", "MP7LD")
- name_rec = "7th order Low-Dissipation Monotonicity Preserving"
- interfaces => interfaces_dir
- reconstruct_states => reconstruct_mp7ld
- order = 7
- nghosts = max(nghosts, 4)
- case ("mp9ld", "MP9LD")
- name_rec = "9th order Low-Dissipation Monotonicity Preserving"
- interfaces => interfaces_dir
- reconstruct_states => reconstruct_mp9ld
- order = 9
- nghosts = max(nghosts, 6)
case ("crmp5", "CRMP5")
name_rec = "5th order Compact Monotonicity Preserving"
interfaces => interfaces_dir
@@ -573,18 +555,18 @@ module interpolations
! calculate coefficients of the low-dissipation methods
!
- ce5ld = [ - cweight +2.000d+00, 5.00d+00 * cweight -1.300d+01, &
- -1.00d+01 * cweight +4.700d+01, 1.00d+01 * cweight +2.700d+01, &
- -5.00d+00 * cweight -3.000d+00, cweight ] / 6.00d+01
- ce7ld = [ 3.00d+00 * cweight -6.000d+00, -2.10d+01 * cweight +5.000d+01, &
- 6.30d+01 * cweight -2.020d+02, -1.05d+02 * cweight +6.380d+02, &
- 1.05d+02 * cweight +4.280d+02, -6.30d+01 * cweight -7.600d+01, &
- 2.10d+01 * cweight +8.000d+00, -3.00d+00 * cweight ] / 8.40d+02
- ce9ld = [ -2.00d+00 * cweight +4.000d+00, 1.80d+01 * cweight -4.100d+01, &
- -7.20d+01 * cweight +1.990d+02, 1.68d+02 * cweight -6.410d+02, &
- -2.52d+02 * cweight +1.879d+03, 2.52d+02 * cweight +1.375d+03, &
- -1.68d+02 * cweight -3.050d+02, 7.20d+01 * cweight +5.500d+01, &
- -1.80d+01 * cweight -5.000d+00, 2.00d+00 * cweight ] / 2.52d+03
+ ce5a = [ - cweight +2.000d+00, 5.00d+00 * cweight -1.300d+01, &
+ -1.00d+01 * cweight +4.700d+01, 1.00d+01 * cweight +2.700d+01, &
+ -5.00d+00 * cweight -3.000d+00, cweight ] / 6.00d+01
+ ce7a = [ 3.00d+00 * cweight -6.000d+00, -2.10d+01 * cweight +5.000d+01, &
+ 6.30d+01 * cweight -2.020d+02, -1.05d+02 * cweight +6.380d+02, &
+ 1.05d+02 * cweight +4.280d+02, -6.30d+01 * cweight -7.600d+01, &
+ 2.10d+01 * cweight +8.000d+00, -3.00d+00 * cweight ] / 8.40d+02
+ ce9a = [ -2.00d+00 * cweight +4.000d+00, 1.80d+01 * cweight -4.100d+01, &
+ -7.20d+01 * cweight +1.990d+02, 1.68d+02 * cweight -6.410d+02, &
+ -2.52d+02 * cweight +1.879d+03, 2.52d+02 * cweight +1.375d+03, &
+ -1.68d+02 * cweight -3.050d+02, 7.20d+01 * cweight +5.500d+01, &
+ -1.80d+01 * cweight -5.000d+00, 2.00d+00 * cweight ] / 2.52d+03
di5ld = [ -(cweight -3.0d+00) / 6.0d+00, 1.0d+00, &
(cweight +1.0d+00) / 6.0d+00 ]
@@ -3873,374 +3855,8 @@ module interpolations
! subroutine RECONSTRUCT_MP5:
! --------------------------
!
-! Subroutine reconstructs the interface states using the fifth 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_mp5(h, fc, fl, fr)
-
-! local variables are not implicit by default
-!
- implicit none
-
-! subroutine arguments
-!
- real(kind=8) , intent(in) :: h
- real(kind=8), dimension(:), intent(in) :: fc
- real(kind=8), dimension(:), intent(out) :: fl, fr
-
-! local variables
-!
- integer :: n, i
-
-! local arrays for derivatives
-!
- real(kind=8), dimension(size(fc)) :: fi
- real(kind=8), dimension(size(fc)) :: u
-!
-!-------------------------------------------------------------------------------
-!
-! get the input vector length
-!
- n = size(fc)
-
-!! === left-side interpolation ===
-!!
-! reconstruct the interface state using the 5th order interpolation
-!
- do i = 3, n - 2
- u(i) = sum(ce5(:) * fc(i-2:i+2))
- 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(n-1) = sum(ce3(:) * fc(n-2: n))
- u(n ) = fc(n )
-
-! apply the monotonicity preserving limiting
-!
- call mp_limiting(fc(:), u(:))
-
-! 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 = 3, n - 2
- u(i) = sum(ce5(:) * fi(i-2:i+2))
- 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(n-1) = sum(ce3(:) * fi(n-2: n))
- u(n ) = fi(n )
-
-! apply the monotonicity preserving limiting
-!
- call mp_limiting(fi(:), u(:))
-
-! 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_mp5
-!
-!===============================================================================
-!
-! 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(h, fc, fl, fr)
-
-! local variables are not implicit by default
-!
- implicit none
-
-! subroutine arguments
-!
- real(kind=8) , intent(in) :: h
- real(kind=8), dimension(:), intent(in) :: fc
- real(kind=8), dimension(:), intent(out) :: fl, fr
-
-! local variables
-!
- integer :: n, i
-
-! local arrays for derivatives
-!
- real(kind=8), dimension(size(fc)) :: fi
- real(kind=8), dimension(size(fc)) :: u
-!
-!-------------------------------------------------------------------------------
-!
-! get the input vector length
-!
- n = size(fc)
-
-!! === 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(fc(:), u(:))
-
-! 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(fi(:), u(:))
-
-! 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_MP9:
-! --------------------------
-!
-! Subroutine reconstructs the interface states using the ninth 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_mp9(h, fc, fl, fr)
-
-! local variables are not implicit by default
-!
- implicit none
-
-! subroutine arguments
-!
- real(kind=8) , intent(in) :: h
- real(kind=8), dimension(:), intent(in) :: fc
- real(kind=8), dimension(:), intent(out) :: fl, fr
-
-! local variables
-!
- integer :: n, i
-
-! local arrays for derivatives
-!
- real(kind=8), dimension(size(fc)) :: fi
- real(kind=8), dimension(size(fc)) :: u
-!
-!-------------------------------------------------------------------------------
-!
-! get the input vector length
-!
- n = size(fc)
-
-!! === left-side interpolation ===
-!!
-! reconstruct the interface state using the 9th order interpolation
-!
- do i = 5, n - 4
- u(i) = sum(ce9(:) * fc(i-4:i+4))
- 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( 4) = sum(ce7(:) * fc( 1: 7))
- u(n-3) = sum(ce7(:) * fc(n-6: n))
- 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(fc(:), u(:))
-
-! 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 9th order interpolation
-!
- do i = 5, n - 4
- u(i) = sum(ce9(:) * fi(i-4:i+4))
- 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( 4) = sum(ce7(:) * fi( 1: 7))
- u(n-3) = sum(ce7(:) * fi(n-6: n))
- 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(fi(:), u(:))
-
-! 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_mp9
-!
-!===============================================================================
-!
-! subroutine RECONSTRUCT_MP5LD:
-! ----------------------------
-!
-! Subroutine reconstructs the interface states using the fifth order
-! low dissipation Monotonicity Preserving (MP) method.
+! Subroutine reconstructs the interface states using the 5th order
+! Monotonicity Preserving (MP) method with adjustable dissipation.
!
! Arguments are described in subroutine reconstruct().
!
@@ -4261,7 +3877,7 @@ module interpolations
!
!===============================================================================
!
- subroutine reconstruct_mp5ld(h, fc, fl, fr)
+ subroutine reconstruct_mp5(h, fc, fl, fr)
implicit none
@@ -4281,7 +3897,7 @@ module interpolations
!! === left-side interpolation ===
!!
do i = 3, n - 3
- u(i) = sum(ce5ld(:) * fc(i-2:i+3))
+ u(i) = sum(ce5a(:) * fc(i-2:i+3))
end do
u( 1) = sum(ce2(:) * fc( 1:2))
@@ -4299,7 +3915,7 @@ module interpolations
fi(1:n) = fc(n:1:-1)
do i = 3, n - 3
- u(i) = sum(ce5ld(:) * fi(i-2:i+3))
+ u(i) = sum(ce5a(:) * fi(i-2:i+3))
end do
u( 1) = sum(ce2(:) * fi( 1:2))
@@ -4320,25 +3936,25 @@ module interpolations
!-------------------------------------------------------------------------------
!
- end subroutine reconstruct_mp5ld
+ end subroutine reconstruct_mp5
!
!===============================================================================
!
-! subroutine RECONSTRUCT_MP7LD:
-! ----------------------------
+! subroutine RECONSTRUCT_MP7:
+! --------------------------
!
-! Subroutine reconstructs the interface states using the seventh order
-! low dissipation Monotonicity Preserving (MP) method.
+! Subroutine reconstructs the interface states using the 7th order
+! Monotonicity Preserving (MP) method with adjustable dissipation.
!
! Arguments are described in subroutine reconstruct().
!
! References:
!
-! - See RECONSTRUCT_MP5LD
+! - See RECONSTRUCT_MP5
!
!===============================================================================
!
- subroutine reconstruct_mp7ld(h, fc, fl, fr)
+ subroutine reconstruct_mp7(h, fc, fl, fr)
implicit none
@@ -4358,7 +3974,7 @@ module interpolations
!! === left-side interpolation ===
!!
do i = 4, n - 4
- u(i) = sum(ce7ld(:) * fc(i-3:i+4))
+ u(i) = sum(ce7a(:) * fc(i-3:i+4))
end do
u( 1) = sum(ce2(:) * fc( 1:2))
@@ -4378,7 +3994,7 @@ module interpolations
fi(1:n) = fc(n:1:-1)
do i = 4, n - 4
- u(i) = sum(ce7ld(:) * fi(i-3:i+4))
+ u(i) = sum(ce7a(:) * fi(i-3:i+4))
end do
u( 1) = sum(ce2(:) * fi( 1:2))
@@ -4401,15 +4017,15 @@ module interpolations
!-------------------------------------------------------------------------------
!
- end subroutine reconstruct_mp7ld
+ end subroutine reconstruct_mp7
!
!===============================================================================
!
-! subroutine RECONSTRUCT_MP9LD:
-! ----------------------------
+! subroutine RECONSTRUCT_MP9:
+! --------------------------
!
-! Subroutine reconstructs the interface states using the ninth order
-! low dissipation Monotonicity Preserving (MP) method.
+! Subroutine reconstructs the interface states using the 9th order
+! Monotonicity Preserving (MP) method with adjustable dissipation.
!
! Arguments are described in subroutine reconstruct().
!
@@ -4419,7 +4035,7 @@ module interpolations
!
!===============================================================================
!
- subroutine reconstruct_mp9ld(h, fc, fl, fr)
+ subroutine reconstruct_mp9(h, fc, fl, fr)
implicit none
@@ -4439,7 +4055,7 @@ module interpolations
!! === left-side interpolation ===
!!
do i = 5, n - 5
- u(i) = sum(ce9ld(:) * fc(i-4:i+5))
+ u(i) = sum(ce9a(:) * fc(i-4:i+5))
end do
u( 1) = sum(ce2(:) * fc( 1:2))
@@ -4461,7 +4077,7 @@ module interpolations
fi(1:n) = fc(n:1:-1)
do i = 5, n - 5
- u(i) = sum(ce9ld(:) * fi(i-4:i+5))
+ u(i) = sum(ce9a(:) * fi(i-4:i+5))
end do
u( 1) = sum(ce2(:) * fi( 1:2))
@@ -4486,7 +4102,7 @@ module interpolations
!-------------------------------------------------------------------------------
!
- end subroutine reconstruct_mp9ld
+ end subroutine reconstruct_mp9
!
!===============================================================================
!
From 71290e573123e053d37355229bd50c5f23b695b4 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 22:28:55 -0300
Subject: [PATCH 11/12] INTERPOLATIONS: Remove regular compact MP methods.
They are identical to low-dissipation methods with central_weight = 0.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 505 ++++---------------------------------
1 file changed, 50 insertions(+), 455 deletions(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 2172982..384a948 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -134,10 +134,10 @@ module interpolations
real(kind=8), dimension(10), save :: ce9a
real(kind=8), dimension( 8), save :: ce7a
real(kind=8), dimension( 6), save :: ce5a
- real(kind=8), dimension( 3), save :: di5ld, di7ld
- real(kind=8), dimension( 5), save :: di9ld
- real(kind=8), dimension( 4), save :: ci5ld
- real(kind=8), dimension( 6), save :: ci7ld, ci9ld
+ real(kind=8), dimension( 3), save :: di5a, di7a
+ real(kind=8), dimension( 5), save :: di9a
+ real(kind=8), dimension( 4), save :: ci5a
+ real(kind=8), dimension( 6), save :: ci7a, ci9a
real(kind=8), dimension(9), parameter :: &
ce9 = [ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03, &
@@ -375,24 +375,6 @@ module interpolations
reconstruct_states => reconstruct_crmp9
order = 9
nghosts = max(nghosts, 6)
- case ("crmp5l", "crmp5ld", "CRMP5L", "CRMP5LD")
- name_rec = "5th order Low-Dissipation Compact Monotonicity Preserving"
- interfaces => interfaces_dir
- reconstruct_states => reconstruct_crmp5ld
- order = 5
- nghosts = max(nghosts, 4)
- case ("crmp7l", "crmp7ld", "CRMP7L", "CRMP7LD")
- name_rec = "7th order Low-Dissipation Compact Monotonicity Preserving"
- interfaces => interfaces_dir
- reconstruct_states => reconstruct_crmp7ld
- order = 7
- nghosts = max(nghosts, 4)
- case ("crmp9l", "crmp9ld", "CRMP9L", "CRMP9LD")
- name_rec = "9th order Low-Dissipation Compact Monotonicity Preserving"
- interfaces => interfaces_dir
- reconstruct_states => reconstruct_crmp9ld
- order = 9
- nghosts = max(nghosts, 6)
case ("ocmp5", "OCMP5")
name_rec = "5th order Optimized Compact Monotonicity Preserving"
interfaces => interfaces_dir
@@ -568,23 +550,23 @@ module interpolations
-1.68d+02 * cweight -3.050d+02, 7.20d+01 * cweight +5.500d+01, &
-1.80d+01 * cweight -5.000d+00, 2.00d+00 * cweight ] / 2.52d+03
- di5ld = [ -(cweight -3.0d+00) / 6.0d+00, 1.0d+00, &
- (cweight +1.0d+00) / 6.0d+00 ]
- ci5ld = [ - cweight +2.0d+00, -9.0d+00 * cweight +3.8d+01, &
- 9.0d+00 * cweight +2.0d+01, cweight ] / 3.6d+01
+ di5a = [ -(cweight -3.0d+00) / 6.0d+00, 1.0d+00, &
+ (cweight +1.0d+00) / 6.0d+00 ]
+ ci5a = [ - cweight +2.0d+00, -9.0d+00 * cweight +3.8d+01, &
+ 9.0d+00 * cweight +2.0d+01, cweight ] / 3.6d+01
- di7ld = [ -(cweight -4.0d+00) / 8.0d+00, 1.0d+00, &
- (cweight +2.0d+00) / 8.0d+00 ]
- ci7ld = [ cweight -2.0d+00, -1.50d+01 * cweight +3.8d+01, &
- -8.0d+01 * cweight +4.78d+02, 8.0d+01 * cweight +3.18d+02, &
- 1.5d+01 * cweight +8.00d+00, -cweight ] / 4.8d+02
+ di7a = [ -(cweight -4.0d+00) / 8.0d+00, 1.0d+00, &
+ (cweight +2.0d+00) / 8.0d+00 ]
+ ci7a = [ cweight -2.0d+00, -1.50d+01 * cweight +3.8d+01, &
+ -8.0d+01 * cweight +4.78d+02, 8.0d+01 * cweight +3.18d+02, &
+ 1.5d+01 * cweight +8.00d+00, -cweight ] / 4.8d+02
- di9ld = [ -2.0d+00 * cweight +5.0d+00, -1.0d+01 * cweight +4.0d+01, &
- 6.0d+01, 1.0d+01 * cweight +2.0d+01, &
- 2.0d+00 * cweight +1.0d+00 ] / 6.0d+01
- ci9ld = [ -3.00d+00 * cweight +6.000d+00,-1.75d+02 * cweight +4.810d+02, &
- -3.00d+02 * cweight +1.881d+03, 3.00d+02 * cweight +1.281d+03, &
- 1.75d+02 * cweight +1.310d+02, 3.00d+00 * cweight ] / 1.8d+03
+ di9a = [ -2.0d+00 * cweight +5.0d+00, -1.0d+01 * cweight +4.0d+01, &
+ 6.0d+01, 1.0d+01 * cweight +2.0d+01, &
+ 2.0d+00 * cweight +1.0d+00 ] / 6.0d+01
+ ci9a = [ -3.00d+00 * cweight +6.000d+00,-1.75d+02 * cweight +4.810d+02, &
+ -3.00d+02 * cweight +1.881d+03, 3.00d+02 * cweight +1.281d+03, &
+ 1.75d+02 * cweight +1.310d+02, 3.00d+00 * cweight ] / 1.8d+03
#ifdef PROFILE
! stop accounting time for module initialization/finalization
@@ -4109,8 +4091,9 @@ module interpolations
! subroutine RECONSTRUCT_CRMP5:
! ----------------------------
!
-! Subroutine reconstructs the interface states using the fifth order
-! Compact Reconstruction Monotonicity Preserving (CRMP) method.
+! Subroutine reconstructs the interface states using the 5th order
+! Compact Reconstruction Monotonicity Preserving (CRMP) method with
+! adjustable dissipation.
!
! Arguments are described in subroutine reconstruct().
!
@@ -4122,12 +4105,12 @@ module interpolations
! 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
+! [2] Jian Fang, Yufeng Yao, Zhaorui Li, Lipeng Lu,
+! "Investigation of low-dissipation monotonicity-preserving scheme
+! for direct numerical simulation of compressible turbulent flows",
+! Computers & Fluids,
+! 2014, vol. 104, pp. 55-72,
+! http://dx.doi.org/10.1016/j.compfluid.2014.07.024
!
!===============================================================================
!
@@ -4147,11 +4130,6 @@ module interpolations
real(kind=8), dimension(size(fc)) :: a, b, c
real(kind=8), dimension(size(fc)) :: r
real(kind=8), dimension(size(fc)) :: u
-
- real(kind=8), dimension(3), parameter :: &
- di5 = [ 3.0d+00, 6.0d+00, 1.0d+00 ] / 6.0d+00
- real(kind=8), dimension(3), parameter :: &
- ci5 = [ 1.0d+00, 1.9d+01, 1.0d+01 ] / 1.8d+01
!
!-------------------------------------------------------------------------------
!
@@ -4163,9 +4141,9 @@ module interpolations
c(i) = 0.0d+00
end do
do i = ng + 1, n - ng - 1
- a(i) = di5(1)
- b(i) = di5(2)
- c(i) = di5(3)
+ a(i) = di5a(1)
+ b(i) = di5a(2)
+ c(i) = di5a(3)
end do
do i = n - ng, n
a(i) = 0.0d+00
@@ -4176,7 +4154,7 @@ module interpolations
!! === left-side interpolation ===
!!
do i = ng, n - ng + 1
- r(i) = sum(ci5(:) * fc(i-1:i+1))
+ r(i) = sum(ci5a(:) * fc(i-1:i+2))
end do
r( 1) = sum(ce2(:) * fc( 1: 2))
@@ -4201,7 +4179,7 @@ module interpolations
fi(1:n) = fc(n:1:-1)
do i = ng, n - ng + 1
- r(i) = sum(ci5(:) * fi(i-1:i+1))
+ r(i) = sum(ci5a(:) * fi(i-1:i+2))
end do ! i = ng, n - ng + 1
r( 1) = sum(ce2(:) * fi( 1: 2))
@@ -4236,8 +4214,9 @@ module interpolations
! subroutine RECONSTRUCT_CRMP7:
! ----------------------------
!
-! Subroutine reconstructs the interface states using the seventh order
-! Compact Reconstruction Monotonicity Preserving (CRMP) method.
+! Subroutine reconstructs the interface states using the 7th order
+! Compact Reconstruction Monotonicity Preserving (CRMP) method with
+! adjustable dissipation.
!
! Arguments are described in subroutine reconstruct().
!
@@ -4263,13 +4242,6 @@ module interpolations
real(kind=8), dimension(size(fc)) :: a, b, c
real(kind=8), dimension(size(fc)) :: r
real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
- real(kind=8), dimension(3), parameter :: &
- di7 = [ 5.0d-01, 1.0d+00, 2.5d-01 ]
- real(kind=8), dimension(5), parameter :: &
- ci7 = [-1.0d+00, 1.9d+01, 2.39d+02, 1.59d+02, 4.0d+00 ] / 2.4d+02
!
!-------------------------------------------------------------------------------
!
@@ -4281,9 +4253,9 @@ module interpolations
c(i) = 0.0d+00
end do
do i = ng + 1, n - ng - 1
- a(i) = di7(1)
- b(i) = di7(2)
- c(i) = di7(3)
+ a(i) = di7a(1)
+ b(i) = di7a(2)
+ c(i) = di7a(3)
end do
do i = n - ng, n
a(i) = 0.0d+00
@@ -4294,7 +4266,7 @@ module interpolations
!! === left-side interpolation ===
!!
do i = ng, n - ng + 1
- r(i) = sum(ci7(:) * fc(i-2:i+2))
+ r(i) = sum(ci7a(:) * fc(i-2:i+3))
end do
r( 1) = sum(ce2(:) * fc( 1: 2))
@@ -4321,7 +4293,7 @@ module interpolations
fi(1:n) = fc(n:1:-1)
do i = ng, n - ng + 1
- r(i) = sum(ci7(:) * fi(i-2:i+2))
+ r(i) = sum(ci7a(:) * fi(i-2:i+3))
end do ! i = ng, n - ng + 1
r( 1) = sum(ce2(:) * fi( 1: 2))
@@ -4359,7 +4331,8 @@ module interpolations
! ----------------------------
!
! Subroutine reconstructs the interface states using the 9th order
-! Compact Reconstruction Monotonicity Preserving (CRMP) method.
+! Compact Reconstruction Monotonicity Preserving (CRMP) method with
+! adjustable dissipation.
!
! Arguments are described in subroutine reconstruct().
!
@@ -4385,13 +4358,6 @@ module interpolations
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
-
-! local parameters
-!
- real(kind=8), dimension(5), parameter :: &
- di9 = [ 5.0d+00, 4.0d+01, 6.0d+01, 2.0d+01, 1.0d+00 ] / 6.0d+01
- real(kind=8), dimension(5), parameter :: &
- ci9 = [ 6.0d+00, 4.81d+02, 1.881d+03, 1.281d+03, 1.31d+02 ] / 1.8d+03
!
!-------------------------------------------------------------------------------
!
@@ -4405,11 +4371,11 @@ module interpolations
b(i) = 0.0d+00
end do
do i = ng + 1, n - ng - 1
- e(i) = di9(1)
- c(i) = di9(2)
- d(i) = di9(3)
- a(i) = di9(4)
- b(i) = di9(5)
+ e(i) = di9a(1)
+ c(i) = di9a(2)
+ d(i) = di9a(3)
+ a(i) = di9a(4)
+ b(i) = di9a(5)
end do
do i = n - ng, n
e(i) = 0.0d+00
@@ -4422,7 +4388,7 @@ module interpolations
!! === left-side interpolation ===
!!
do i = ng, n - ng + 1
- r(i) = sum(ci9(:) * fc(i-2:i+2))
+ r(i) = sum(ci9a(:) * fc(i-2:i+3))
end do
r( 1) = sum(ce2(:) * fc( 1: 2))
@@ -4451,7 +4417,7 @@ module interpolations
fi(1:n) = fc(n:1:-1)
do i = ng, n - ng + 1
- r(i) = sum(ci9(:) * fi(i-2:i+2))
+ r(i) = sum(ci9a(:) * fi(i-2:i+3))
end do ! i = ng, n - ng + 1
r( 1) = sum(ce2(:) * fi( 1: 2))
@@ -4487,377 +4453,6 @@ module interpolations
!
!===============================================================================
!
-! subroutine RECONSTRUCT_CRMP5LD:
-! ------------------------------
-!
-! Subroutine reconstructs the interface states using the fifth order
-! Low-Dissipation Compact Reconstruction Monotonicity Preserving (CRMP)
-! 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
-! [3] Ghosh, D. & Baeder, J.,
-! "Compact Reconstruction Schemes With Weighted ENO Limiting For
-! Hyperbolic Conservation Laws",
-! SIAM Journal on Scientific Computing,
-! 2012, vol. 34, no. 3, pp. A1678-A1705,
-! http://dx.doi.org/10.1137/110857659
-!
-!===============================================================================
-!
- subroutine reconstruct_crmp5ld(h, fc, fl, fr)
-
- use algebra, only : tridiag
-
- 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)) :: a, b, c
- real(kind=8), dimension(size(fc)) :: r
- real(kind=8), dimension(size(fc)) :: u
-!
-!-------------------------------------------------------------------------------
-!
- n = size(fc)
-
- 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) = di5ld(1)
- b(i) = di5ld(2)
- c(i) = di5ld(3)
- end do
- do i = n - ng, n
- a(i) = 0.0d+00
- b(i) = 1.0d+00
- c(i) = 0.0d+00
- end do
-
-!! === left-side interpolation ===
-!!
- do i = ng, n - ng + 1
- r(i) = sum(ci5ld(:) * fc(i-1:i+2))
- end do
-
- r( 1) = sum(ce2(:) * fc( 1: 2))
- r( 2) = sum(ce3(:) * fc( 1: 3))
- do i = 3, ng
- r(i) = sum(ce5(:) * fc(i-2:i+2))
- end do
- do i = n - ng, n - 2
- r(i) = sum(ce5(:) * fc(i-2:i+2))
- end do
- r(n-1) = sum(ce3(:) * fc(n-2: n))
- r(n ) = fc(n )
-
- call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
-
- 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(ci5ld(:) * fi(i-1:i+2))
- end do ! i = ng, n - ng + 1
-
- r( 1) = sum(ce2(:) * fi( 1: 2))
- r( 2) = sum(ce3(:) * fi( 1: 3))
- do i = 3, ng
- r(i) = sum(ce5(:) * fi(i-2:i+2))
- end do
- do i = n - ng, n - 2
- r(i) = sum(ce5(:) * fi(i-2:i+2))
- end do
- r(n-1) = sum(ce3(:) * fi(n-2: n))
- r(n ) = fi(n )
-
- call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
-
- 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_crmp5ld
-!
-!===============================================================================
-!
-! subroutine RECONSTRUCT_CRMP7LD:
-! ------------------------------
-!
-! Subroutine reconstructs the interface states using the seventh order
-! low dissipation Compact Reconstruction Monotonicity Preserving (CRMP)
-! method.
-!
-! Arguments are described in subroutine reconstruct().
-!
-! References:
-!
-! - See RECONSTRUCT_CRMP5LD
-!
-!===============================================================================
-!
- subroutine reconstruct_crmp7ld(h, fc, fl, fr)
-
- use algebra, only : tridiag
-
- 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)) :: a, b, c
- real(kind=8), dimension(size(fc)) :: r
- real(kind=8), dimension(size(fc)) :: u
-!
-!-------------------------------------------------------------------------------
-!
- n = size(fc)
-
- 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) = di7ld(1)
- b(i) = di7ld(2)
- c(i) = di7ld(3)
- end do
- do i = n - ng, n
- a(i) = 0.0d+00
- b(i) = 1.0d+00
- c(i) = 0.0d+00
- end do
-
-!! === left-side interpolation ===
-!!
- do i = ng, n - ng + 1
- r(i) = sum(ci7ld(:) * fc(i-2:i+3))
- 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 tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
-
- 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(ci7ld(:) * fi(i-2:i+3))
- 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 tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
-
- 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_crmp7ld
-!
-!===============================================================================
-!
-! subroutine RECONSTRUCT_CRMP9LD:
-! ------------------------------
-!
-! Subroutine reconstructs the interface states using the 9th order
-! low dissipation Compact Reconstruction Monotonicity Preserving (CRMP)
-! method.
-!
-! Arguments are described in subroutine reconstruct().
-!
-! References:
-!
-! - See RECONSTRUCT_CRMP5LD
-!
-!===============================================================================
-!
- subroutine reconstruct_crmp9ld(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
-!
-!-------------------------------------------------------------------------------
-!
- n = size(fc)
-
- 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) = di9ld(1)
- c(i) = di9ld(2)
- d(i) = di9ld(3)
- a(i) = di9ld(4)
- b(i) = di9ld(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(ci9ld(:) * fc(i-2:i+3))
- end do
-
- r( 1) = sum(ce2(:) * fc( 1: 2))
- r( 2) = sum(ce3(:) * fc( 1: 3))
- r( 3) = sum(ce5(:) * fc( 1: 5))
- r( 4) = sum(ce7(:) * fc( 1: 7))
- do i = 5, ng
- r(i) = sum(ce9(:) * fc(i-4:i+4))
- end do
- do i = n - ng, n - 4
- r(i) = sum(ce9(:) * fc(i-4:i+4))
- end do
- r(n-3) = sum(ce7(:) * fc(n-6: n))
- 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(ci9ld(:) * fi(i-2:i+3))
- 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))
- r( 4) = sum(ce7(:) * fi( 1: 7))
- do i = 5, ng
- r(i) = sum(ce9(:) * fi(i-4:i+4))
- end do
- do i = n - ng, n - 4
- r(i) = sum(ce9(:) * fi(i-4:i+4))
- end do
- r(n-3) = sum(ce7(:) * fi(n-6: n))
- 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_crmp9ld
-!
-!===============================================================================
-!
! subroutine RECONSTRUCT_OCMP5:
! -----------------------------
!
From 7bf0635c6ed7dde4c4bc79aa1b1447fa2ebcd48c Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 7 Oct 2020 22:30:50 -0300
Subject: [PATCH 12/12] INTERPOLATIONS: Set central_weight to zero by default.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/interpolations.F90 | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 384a948..9a8d6e4 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -116,7 +116,7 @@ module interpolations
! weight toward central scheme for low-dissipation methods
!
- real(kind=8), save :: cweight = 7.0d-01
+ real(kind=8), save :: cweight = 0.0d+00
! Gaussian process reconstruction coefficients vector
!