diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index 635741d..338e2a7 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -367,6 +367,12 @@ module interpolations reconstruct_states => reconstruct_crmp7ld order = 7 nghosts = max(nghosts, 4) + case ("ocmp5", "OCMP5") + name_rec = "5th order Optimized Compact Monotonicity Preserving" + interfaces => interfaces_dir + reconstruct_states => reconstruct_ocmp5 + order = 5 + nghosts = max(nghosts, 4) case ("gp", "GP") write(stmp, '(f16.1)') sgp write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp & @@ -5318,6 +5324,160 @@ module interpolations ! !=============================================================================== ! +! subroutine RECONSTRUCT_OCMP5: +! ----------------------------- +! +! Subroutine reconstructs the interface states using the 5th order Optimized +! Compact Reconstruction Monotonicity Preserving (CRMP) method. +! +! Arguments are described in subroutine reconstruct(). +! +! References: +! +! [1] Myeong-Hwan Ahn, Duck-Joo Lee, +! "Modified Monotonicity Preserving Constraints for High-Resolution +! Optimized Compact Scheme", +! Journal of Scientific Computing, +! 2020, vol. 83, p. 34 +! https://doi.org/10.1007/s10915-020-01221-0 +! +!=============================================================================== +! + subroutine reconstruct_ocmp5(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, iret + + 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), 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 ] + 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 ] +! +!------------------------------------------------------------------------------- +! + 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 right-hand side of the linear system +! + do i = ng, n - ng + 1 + r(i) = sum(ci5(:) * fc(i-2:i+2)) + end do + +! use explicit methods for ghost zones +! + r( 1) = sum(ce2(:) * fc( 1: 2)) + r( 2) = sum(ce3(:) * fc( 1: 3)) + do i = 3, ng + 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 ) + +! solve the tridiagonal system of equations +! + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + +! apply the monotonicity preserving limiter +! + call mp_limiting(fc(:), u(:)) + +! return the interpolated values of the left state +! + fl(1:n) = u(1:n) + +!! === right-side interpolation === +!! +! invert the cell-centered integrals +! + fi(1:n) = fc(n:1:-1) + +! prepare the right-hand side of the linear system +! + do i = ng, n - ng + 1 + r(i) = sum(ci5(:) * fi(i-2:i+2)) + end do ! i = ng, n - ng + 1 + +! use explicit methods for ghost zones +! + r( 1) = sum(ce2(:) * fi( 1: 2)) + r( 2) = sum(ce3(:) * fi( 1: 3)) + do i = 3, ng + 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 ) + +! solve the tridiagonal system of equations +! + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + +! apply the monotonicity preserving limiter +! + call mp_limiting(fi(:), u(:)) + +! return the interpolated values of the right state +! + fr(1:n-1) = u(n-1:1:-1) + +! update the extremum points +! + i = n - 1 + fl(1) = 0.5d+00 * (fc(1) + fc(2)) + fr(i) = 0.5d+00 * (fc(i) + fc(n)) + fl(n) = fc(n) + fr(n) = fc(n) + +!------------------------------------------------------------------------------- +! + end subroutine reconstruct_ocmp5 +! +!=============================================================================== +! ! subroutine PREPARE_GP: ! --------------------- ! @@ -6104,6 +6264,12 @@ module interpolations ! Science China Physics, Mechanics and Astronomy, ! Volume 54, Issue 3, pp. 511-522, ! http://dx.doi.org/10.1007/s11433-010-4220-x +! [3] 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 ! !=============================================================================== ! @@ -6120,8 +6286,9 @@ module interpolations ! local variables ! - integer :: n, i, im1, ip1, ip2 - real(kind=8) :: dq, ds, dc0, dc4, dm1, dp1, dml, dmr + logical :: test + integer :: n, i, im2, im1, ip1, ip2 + real(kind=8) :: dq, ds, dc0, dc4, dm1, dp1, dml, dmr, bt real(kind=8) :: qlc, qmd, qmp, qmn, qmx, qul ! local vectors @@ -6129,12 +6296,10 @@ module interpolations real(kind=8), dimension(0:size(qc)+2) :: dm ! !------------------------------------------------------------------------------- -! -! get the input vector size ! n = size(qc) -! calculate derivatives +! 1st order derivatives ! dm(0 ) = 0.0d+00 dm(1 ) = 0.0d+00 @@ -6142,7 +6307,7 @@ module interpolations dm(n+1) = 0.0d+00 dm(n+2) = 0.0d+00 -! check monotonicity condition for all elements and apply limiting if required +! check the monotonicity condition and apply limiting if necessary ! do i = 1, n - 1 @@ -6159,6 +6324,7 @@ module interpolations if (ds > eps) then + im2 = i - 2 im1 = i - 1 ip2 = i + 2 @@ -6174,6 +6340,17 @@ module interpolations qul = qc(i) + dq qlc = qc(i) + 0.5d+00 * dq + dml + test = qc(i) > max(qc(im1), qc(ip1)) .and. & + min(qc(im1), qc(ip1)) > max(qc(im2), qc(ip2)) + test = test .or. qc(i) < min(qc(im1), qc(ip1)) .and. & + max(qc(im1), qc(ip1)) < min(qc(im2), qc(ip2)) + if (test) then + if (qc(im2) <= qc(ip1) .and. qc(ip2) <= qc(im1)) then + bt = dc0 / (qc(ip2) + qc(im2) - 2.0d+00 * qc(i)) + if (bt >= 3.0d-01) qlc = qc(im2) + end if + end if + qmx = max(min(qc(i), qc(ip1), qmd), min(qc(i), qul, qlc)) qmn = min(max(qc(i), qc(ip1), qmd), max(qc(i), qul, qlc))