Merge branch 'master' into flux-tubes
This commit is contained in:
commit
22ae707b99
@ -138,15 +138,6 @@ module interpolations
|
||||
real(kind=8), dimension(2), parameter :: &
|
||||
ce2 = [ 5.0d-01, 5.0d-01 ]
|
||||
|
||||
! coefficients for the optimized compact schemes
|
||||
!
|
||||
logical :: oc_scheme = .false.
|
||||
logical :: oc_stable = .false.
|
||||
real(kind=8), dimension(2) :: oc_a = 0.0d+00
|
||||
real(kind=8), dimension(3) :: di5
|
||||
real(kind=8), dimension(5) :: di7, di9, ci5, ci7
|
||||
real(kind=8), dimension(7) :: ci9
|
||||
|
||||
! by default everything is private
|
||||
!
|
||||
private
|
||||
@ -231,15 +222,6 @@ module interpolations
|
||||
call get_parameter("central_weight" , cweight )
|
||||
call get_parameter("cfl" , cfl )
|
||||
|
||||
stmp = 'minimum'
|
||||
call get_parameter("ocmp_scheme_mode", stmp)
|
||||
select case(trim(stmp))
|
||||
case ("stable", "STABLE", "stability", "STABILITY")
|
||||
oc_stable = .true.
|
||||
case default
|
||||
oc_stable = .false.
|
||||
end select
|
||||
|
||||
! calculate κ = (1 - ν) / ν
|
||||
!
|
||||
kappa = min(kappa, (1.0d+00 - cfl) / cfl)
|
||||
@ -365,60 +347,34 @@ module interpolations
|
||||
reconstruct_states => reconstruct_ocmp5
|
||||
order = 5
|
||||
nghosts = max(nghosts, 4)
|
||||
oc_scheme = .true.
|
||||
if (oc_stable) then
|
||||
oc_a = [ 4.98486328125000044409d-01, 2.50756835937499977796d-01 ]
|
||||
else
|
||||
oc_a = [ 5.01779599466318559919d-01, 2.54377898581479855444d-01 ]
|
||||
end if
|
||||
di5 = [ oc_a(1), 1.0d+00, oc_a(2) ]
|
||||
ci5 = [- 3.0d+00 * oc_a(1) - 3.0d+00 * oc_a(2) + 2.0d+00, &
|
||||
2.7d+01 * oc_a(1) + 1.7d+01 * oc_a(2) - 1.3d+01, &
|
||||
4.7d+01 * oc_a(1) - 4.3d+01 * oc_a(2) + 4.7d+01, &
|
||||
- 1.3d+01 * oc_a(1) + 7.7d+01 * oc_a(2) + 2.7d+01, &
|
||||
2.0d+00 * oc_a(1) + 1.2d+01 * oc_a(2) - 3.0d+00 ] / 6.0d+01
|
||||
case ("ocmp7", "OCMP7")
|
||||
name_rec = "7th order Optimized Compact Monotonicity Preserving"
|
||||
case ("ocmp7t", "OCMP7t")
|
||||
name_rec = "7th order Tridiagonal Optimized " // &
|
||||
"Compact Monotonicity Preserving"
|
||||
interfaces => interfaces_dir
|
||||
reconstruct_states => reconstruct_ocmp7
|
||||
reconstruct_states => reconstruct_ocmp7t
|
||||
order = 7
|
||||
nghosts = max(nghosts, 6)
|
||||
oc_scheme = .true.
|
||||
if (oc_stable) then
|
||||
oc_a = [ 6.64390055338541563046d-01, 3.34471638997395903647d-01 ]
|
||||
else
|
||||
oc_a = [ 6.71157762607103580699d-01, 3.40895359923885754583d-01 ]
|
||||
end if
|
||||
di7 = [ (3.0d+00 * oc_a(1) + 2.0d+00 * oc_a(2) - 2.0d+00) / 8.0d+00, &
|
||||
oc_a(1), 1.0d+00, oc_a(2), &
|
||||
( oc_a(1) + 6.0d+00 * oc_a(2) - 2.0d+00) / 4.0d+01 ]
|
||||
ci7 = [ 1.80d+01 * oc_a(1) + 1.80d+01 * oc_a(2) - 1.60d+01, &
|
||||
5.43d+02 * oc_a(1) + 2.68d+02 * oc_a(2) - 2.91d+02, &
|
||||
3.43d+02 * oc_a(1) - 3.32d+02 * oc_a(2) + 5.09d+02, &
|
||||
-1.07d+02 * oc_a(1) + 5.68d+02 * oc_a(2) + 3.09d+02, &
|
||||
4.30d+01 * oc_a(1) + 3.18d+02 * oc_a(2) - 9.10d+01 ] / 6.0d+02
|
||||
case ("ocmp9", "OCMP9")
|
||||
name_rec = "9th order Optimized Compact Monotonicity Preserving"
|
||||
case ("ocmp9t", "OCMP9t")
|
||||
name_rec = "9th order Tridiagonal Optimized " // &
|
||||
"Compact Monotonicity Preserving"
|
||||
interfaces => interfaces_dir
|
||||
reconstruct_states => reconstruct_ocmp9
|
||||
reconstruct_states => reconstruct_ocmp9t
|
||||
order = 7
|
||||
nghosts = max(nghosts, 6)
|
||||
case ("ocmp7", "OCMP7", "ocmp7p", "OCMP7P")
|
||||
name_rec = "7th order Pentadiagonal Optimized " // &
|
||||
"Compact Monotonicity Preserving"
|
||||
interfaces => interfaces_dir
|
||||
reconstruct_states => reconstruct_ocmp7p
|
||||
order = 7
|
||||
nghosts = max(nghosts, 6)
|
||||
case ("ocmp9", "OCMP9", "ocmp9p", "OCMP9P")
|
||||
name_rec = "9th order Pentadiagonal Optimized " // &
|
||||
"Compact Monotonicity Preserving"
|
||||
interfaces => interfaces_dir
|
||||
reconstruct_states => reconstruct_ocmp9p
|
||||
order = 9
|
||||
nghosts = max(nghosts, 6)
|
||||
oc_scheme = .true.
|
||||
if (oc_stable) then
|
||||
oc_a = [ 6.64322884877522779057d-01, 4.01406269073486310361d-01 ]
|
||||
else
|
||||
oc_a = [ 6.70794426721399439373d-01, 4.09731006175920509094d-01 ]
|
||||
end if
|
||||
di9 = [ (9.0d+00 * oc_a(1) + 5.0d+00 * oc_a(2) - 6.0d+00) / 2.0d+01, &
|
||||
oc_a(1), 1.0d+00, oc_a(2), &
|
||||
( oc_a(1) + 5.0d+00 * oc_a(2) - 2.0d+00) / 2.0d+01 ]
|
||||
ci9 = [-1.000d+01 * oc_a(1) -1.000d+01 * oc_a(2) +1.000d+01, &
|
||||
2.420d+02 * oc_a(1) +2.000d+02 * oc_a(2) -2.140d+02, &
|
||||
4.085d+03 * oc_a(1) +1.635d+03 * oc_a(2) -2.146d+03, &
|
||||
2.335d+03 * oc_a(1) -1.865d+03 * oc_a(2) +3.454d+03, &
|
||||
-8.150d+02 * oc_a(1) +3.385d+03 * oc_a(2) +2.404d+03, &
|
||||
4.450d+02 * oc_a(1) +2.895d+03 * oc_a(2) -9.560d+02, &
|
||||
1.800d+01 * oc_a(1) +6.000d+01 * oc_a(2) -3.200d+01 ] / 4.2d+03
|
||||
case ("gp", "GP")
|
||||
write(stmp, '(f16.1)') sgp
|
||||
write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp &
|
||||
@ -709,13 +665,6 @@ module interpolations
|
||||
|
||||
call print_section(verbose, "Interpolations")
|
||||
call print_parameter(verbose, "reconstruction" , name_rec )
|
||||
if (oc_scheme) then
|
||||
if (oc_stable) then
|
||||
call print_parameter(verbose, "scheme optimized for", "stability")
|
||||
else
|
||||
call print_parameter(verbose, "scheme optimized for", "minimum error")
|
||||
end if
|
||||
end if
|
||||
call print_parameter(verbose, "TVD limiter" , name_tlim)
|
||||
call print_parameter(verbose, "prolongation limiter", name_plim)
|
||||
if (mlp) then
|
||||
@ -4050,7 +3999,7 @@ module interpolations
|
||||
|
||||
do i = ng, n - ng + 1
|
||||
r(i) = sum(ci5a(:) * fi(i-1:i+2))
|
||||
end do ! i = ng, n - ng + 1
|
||||
end do
|
||||
|
||||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||||
@ -4165,7 +4114,7 @@ module interpolations
|
||||
|
||||
do i = ng, n - ng + 1
|
||||
r(i) = sum(ci7a(:) * fi(i-2:i+3))
|
||||
end do ! i = ng, n - ng + 1
|
||||
end do
|
||||
|
||||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||||
@ -4290,7 +4239,7 @@ module interpolations
|
||||
|
||||
do i = ng, n - ng + 1
|
||||
r(i) = sum(ci9a(:) * fi(i-2:i+3))
|
||||
end do ! i = ng, n - ng + 1
|
||||
end do
|
||||
|
||||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||||
@ -4326,11 +4275,16 @@ module interpolations
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RECONSTRUCT_OCMP5:
|
||||
! -----------------------------
|
||||
! ----------------------------
|
||||
!
|
||||
! Subroutine reconstructs the interface states using the 5th 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 6.4e-12.
|
||||
!
|
||||
! Arguments are described in subroutine reconstruct().
|
||||
!
|
||||
! References:
|
||||
@ -4362,6 +4316,16 @@ module interpolations
|
||||
real(kind=8), dimension(size(fc)) :: r
|
||||
real(kind=8), dimension(size(fc)) :: u
|
||||
|
||||
real(kind=8), parameter :: a1 = 5.0148583220385203994203557d-01
|
||||
real(kind=8), parameter :: a2 = 2.5452055302783787784904423d-01
|
||||
real(kind=8), dimension(3), parameter :: di5 = [ a1, 1.0d+00, a2 ]
|
||||
real(kind=8), dimension(5), parameter :: &
|
||||
ci5 = [- 3.0d+00 * a1 - 3.0d+00 * a2 + 2.0d+00, &
|
||||
2.7d+01 * a1 + 1.7d+01 * a2 - 1.3d+01, &
|
||||
4.7d+01 * a1 - 4.3d+01 * a2 + 4.7d+01, &
|
||||
- 1.3d+01 * a1 + 7.7d+01 * a2 + 2.7d+01, &
|
||||
2.0d+00 * a1 + 1.2d+01 * a2 - 3.0d+00 ] / 6.0d+01
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
n = size(fc)
|
||||
@ -4411,7 +4375,7 @@ module interpolations
|
||||
|
||||
do i = ng, n - ng + 1
|
||||
r(i) = sum(ci5(:) * fi(i-2:i+2))
|
||||
end do ! i = ng, n - ng + 1
|
||||
end do
|
||||
|
||||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||||
@ -4431,8 +4395,8 @@ module interpolations
|
||||
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(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)
|
||||
|
||||
@ -4442,26 +4406,286 @@ module interpolations
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RECONSTRUCT_OCMP7:
|
||||
! subroutine RECONSTRUCT_OCMP7T:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine reconstructs the interface states using the 7th 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.0e-11.
|
||||
!
|
||||
! Arguments are described in subroutine reconstruct().
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Myeong-Hwan Ahn, Duck-Joo Lee,
|
||||
! "Modified Monotonicity Preserving Constraints for High-Resolution
|
||||
! Optimized Compact Scheme",
|
||||
! Journal of Scientific Computing,
|
||||
! 2020, vol. 83, p. 34
|
||||
! https://doi.org/10.1007/s10915-020-01221-0
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine reconstruct_ocmp7(p, h, fc, fl, fr)
|
||||
subroutine reconstruct_ocmp7t(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.0119544102147968726471782d-01
|
||||
real(kind=8), parameter :: a2 = 3.0649545768561641522180678d-01
|
||||
real(kind=8), dimension(3), parameter :: di7 = [ a1, 1.0d+00, a2 ]
|
||||
real(kind=8), dimension(7), parameter :: &
|
||||
ci7 = [ 4.00d+00 * a1 + 4.00d+00 * a2 - 3.00d+00, &
|
||||
- 3.80d+01 * a1 - 3.10d+01 * a2 + 2.50d+01, &
|
||||
2.14d+02 * a1 + 1.09d+02 * a2 - 1.01d+02, &
|
||||
3.19d+02 * a1 - 2.41d+02 * a2 + 3.19d+02, &
|
||||
- 1.01d+02 * a1 + 4.59d+02 * a2 + 2.14d+02, &
|
||||
2.50d+01 * a1 + 1.30d+02 * a2 - 3.80d+01, &
|
||||
- 3.00d+00 * a1 - 1.00d+01 * a2 + 4.00d+00 ] / 4.2d+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-3: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(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(ci7(:) * fi(i-3:i+3))
|
||||
end do
|
||||
|
||||
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(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_ocmp7t
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! 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:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine reconstructs the interface states using the 7th order Optimized
|
||||
! Compact Reconstruction Monotonicity Preserving (CRMP) method.
|
||||
!
|
||||
! The method uses an implicit penta-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.38e-11.
|
||||
!
|
||||
! Arguments are described in subroutine reconstruct().
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine reconstruct_ocmp7p(p, h, fc, fl, fr)
|
||||
|
||||
use algebra, only : pentadiag
|
||||
|
||||
@ -4479,12 +4703,23 @@ module interpolations
|
||||
real(kind=8), dimension(size(fc)) :: r
|
||||
real(kind=8), dimension(size(fc)) :: u
|
||||
|
||||
real(kind=8), parameter :: a1 = 6.7100477315625665203995935d-01
|
||||
real(kind=8), parameter :: a2 = 3.4097048071278782662422942d-01
|
||||
real(kind=8), dimension(5), parameter :: &
|
||||
di7 = [ (3.0d+00 * a1 + 2.0d+00 * a2 - 2.0d+00) / 8.0d+00, &
|
||||
a1, 1.0d+00, a2, &
|
||||
( a1 + 6.0d+00 * a2 - 2.0d+00) / 4.0d+01 ]
|
||||
real(kind=8), dimension(5), parameter :: &
|
||||
ci7 = [ 1.80d+01 * a1 + 1.80d+01 * a2 - 1.60d+01, &
|
||||
5.43d+02 * a1 + 2.68d+02 * a2 - 2.91d+02, &
|
||||
3.43d+02 * a1 - 3.32d+02 * a2 + 5.09d+02, &
|
||||
-1.07d+02 * a1 + 5.68d+02 * a2 + 3.09d+02, &
|
||||
4.30d+01 * a1 + 3.18d+02 * a2 - 9.10d+01 ] / 6.0d+02
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
n = size(fc)
|
||||
|
||||
! prepare the diagonals of the tridiagonal matrix
|
||||
!
|
||||
do i = 1, ng
|
||||
e(i) = 0.0d+00
|
||||
c(i) = 0.0d+00
|
||||
@ -4538,7 +4773,7 @@ module interpolations
|
||||
|
||||
do i = ng, n - ng + 1
|
||||
r(i) = sum(ci7(:) * fi(i-2:i+2))
|
||||
end do ! i = ng, n - ng + 1
|
||||
end do
|
||||
|
||||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||||
@ -4560,37 +4795,33 @@ module interpolations
|
||||
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(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_ocmp7
|
||||
end subroutine reconstruct_ocmp7p
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RECONSTRUCT_OCMP9:
|
||||
! subroutine RECONSTRUCT_OCMP9P:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine reconstructs the interface states using the 9th order Optimized
|
||||
! Compact Reconstruction Monotonicity Preserving (CRMP) method.
|
||||
!
|
||||
! The method uses an implicit penta-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.54e-11.
|
||||
!
|
||||
! Arguments are described in subroutine reconstruct().
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Myeong-Hwan Ahn, Duck-Joo Lee,
|
||||
! "Modified Monotonicity Preserving Constraints for High-Resolution
|
||||
! Optimized Compact Scheme",
|
||||
! Journal of Scientific Computing,
|
||||
! 2020, vol. 83, p. 34
|
||||
! https://doi.org/10.1007/s10915-020-01221-0
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine reconstruct_ocmp9(p, h, fc, fl, fr)
|
||||
subroutine reconstruct_ocmp9p(p, h, fc, fl, fr)
|
||||
|
||||
use algebra, only : pentadiag
|
||||
|
||||
@ -4608,12 +4839,25 @@ module interpolations
|
||||
real(kind=8), dimension(size(fc)) :: r
|
||||
real(kind=8), dimension(size(fc)) :: u
|
||||
|
||||
real(kind=8), parameter :: a1 = 6.7028117788580347100541103d-01
|
||||
real(kind=8), parameter :: a2 = 4.1002518667894543684869155d-01
|
||||
real(kind=8), dimension(5), parameter :: &
|
||||
di9 = [ (9.0d+00 * a1 + 5.0d+00 * a2 - 6.0d+00) / 2.0d+01, &
|
||||
a1, 1.0d+00, a2, &
|
||||
( a1 + 5.0d+00 * a2 - 2.0d+00) / 2.0d+01 ]
|
||||
real(kind=8), dimension(7), parameter :: &
|
||||
ci9 = [-1.000d+01 * a1 -1.000d+01 * a2 +1.000d+01, &
|
||||
2.420d+02 * a1 +2.000d+02 * a2 -2.140d+02, &
|
||||
4.085d+03 * a1 +1.635d+03 * a2 -2.146d+03, &
|
||||
2.335d+03 * a1 -1.865d+03 * a2 +3.454d+03, &
|
||||
-8.150d+02 * a1 +3.385d+03 * a2 +2.404d+03, &
|
||||
4.450d+02 * a1 +2.895d+03 * a2 -9.560d+02, &
|
||||
1.800d+01 * a1 +6.000d+01 * a2 -3.200d+01 ] / 4.2d+03
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
n = size(fc)
|
||||
|
||||
! prepare the diagonals of the tridiagonal matrix
|
||||
!
|
||||
do i = 1, ng
|
||||
e(i) = 0.0d+00
|
||||
c(i) = 0.0d+00
|
||||
@ -4669,7 +4913,7 @@ module interpolations
|
||||
|
||||
do i = ng, n - ng + 1
|
||||
r(i) = sum(ci9(:) * fi(i-3:i+3))
|
||||
end do ! i = ng, n - ng + 1
|
||||
end do
|
||||
|
||||
r( 1) = sum(ce2(:) * fi( 1: 2))
|
||||
r( 2) = sum(ce3(:) * fi( 1: 3))
|
||||
@ -4693,14 +4937,14 @@ module interpolations
|
||||
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(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_ocmp9
|
||||
end subroutine reconstruct_ocmp9p
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user