INTERPOLATIONS: Implement tridiagonal version of OCMP9 scheme.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2023-02-22 14:31:49 -03:00
parent 7a99094427
commit 229759f1c0

View File

@ -354,6 +354,13 @@ module interpolations
reconstruct_states => reconstruct_ocmp7t
order = 7
nghosts = max(nghosts, 6)
case ("ocmp9t", "OCMP9t")
name_rec = "9th order Tridiagonal Optimized " // &
"Compact Monotonicity Preserving"
interfaces => interfaces_dir
reconstruct_states => reconstruct_ocmp9t
order = 7
nghosts = max(nghosts, 6)
case ("ocmp7", "OCMP7", "ocmp7p", "OCMP7P")
name_rec = "7th order Pentadiagonal Optimized " // &
"Compact Monotonicity Preserving"
@ -4528,6 +4535,141 @@ module interpolations
!
!===============================================================================
!
! subroutine RECONSTRUCT_OCMP9T:
! -----------------------------
!
! Subroutine reconstructs the interface states using the 9th order Optimized
! Compact Reconstruction Monotonicity Preserving (CRMP) method.
!
! The method uses an implicit tri-diagonal solver and its coefficients are
! optimized for the maximum value of imaginary part of the scaled modified
! wave number (Im(κ)) to be 2e-6 and the integrated error of compact
! discretization E to be 1.32e-11.
!
! Arguments are described in subroutine reconstruct().
!
!===============================================================================
!
subroutine reconstruct_ocmp9t(p, h, fc, fl, fr)
use algebra, only : tridiag
implicit none
logical , intent(in) :: p
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
real(kind=8), parameter :: a1 = 5.0087697436175480833840357d-01
real(kind=8), parameter :: a2 = 3.4091832657841234002695365d-01
real(kind=8), dimension(3), parameter :: di9 = [ a1, 1.0d+00, a2 ]
real(kind=8), dimension(9), parameter :: &
ci9 = [ - 5.000d+00 * a1 - 5.000d+00 * a2 + 4.000d+00, &
5.500d+01 * a1 + 4.900d+01 * a2 - 4.100d+01, &
- 3.050d+02 * a1 - 2.210d+02 * a2 + 1.990d+02, &
1.375d+03 * a1 + 6.190d+02 * a2 - 6.410d+02, &
1.879d+03 * a1 - 1.271d+03 * a2 + 1.879d+03, &
- 6.410d+02 * a1 + 2.509d+03 * a2 + 1.375d+03, &
1.990d+02 * a1 + 9.550d+02 * a2 - 3.050d+02, &
- 4.100d+01 * a1 - 1.250d+02 * a2 + 5.500d+01, &
4.000d+00 * a1 + 1.000d+01 * a2 - 5.000d+00 ] / 2.52d+03
!-------------------------------------------------------------------------------
!
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) = di9(1)
b(i) = di9(2)
c(i) = di9(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(ci9(:) * fc(i-4:i+4))
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 tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
call mp_limiting(p, n, 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-4:i+4))
end do
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 tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
call mp_limiting(p, n, fi(:), u(:))
fr(1:n-1) = u(n-1:1:-1)
i = n - 1
fl(1) = 5.0d-01 * (fc(1) + fc(2))
fr(i) = 5.0d-01 * (fc(i) + fc(n))
fl(n) = fc(n)
fr(n) = fc(n)
!-------------------------------------------------------------------------------
!
end subroutine reconstruct_ocmp9t
!
!===============================================================================
!
! subroutine RECONSTRUCT_OCMP7P:
! -----------------------------
!