Merge branch 'master' into reconnection

This commit is contained in:
Grzegorz Kowal 2020-09-28 12:17:33 -03:00
commit b4bc20036c

View File

@ -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))