From 3263a08ec04a39c37a67c605f77f24160e0e1677 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Mon, 7 Oct 2024 10:35:04 -0300
Subject: [PATCH] BOUNDARIES, MESH: Implement PPM block prolongation.
Implement Colella and Woodward's Piecewise Parabolic Method (PPM) for
block prolongation. This method is also used for the prolongation of
boundaries from blocks at lower levels to higher levels.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
sources/boundaries.F90 | 680 +++++++++++++++++++++++++++++++++++--
sources/interpolations.F90 | 14 +-
sources/mesh.F90 | 268 ++++++++++++++-
3 files changed, 924 insertions(+), 38 deletions(-)
diff --git a/sources/boundaries.F90 b/sources/boundaries.F90
index 649db61..82c7708 100644
--- a/sources/boundaries.F90
+++ b/sources/boundaries.F90
@@ -45,6 +45,32 @@ module boundaries
real(kind=8), dimension(:) , intent(in) :: x, y, z
real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
end subroutine
+
+#if NDIMS == 3
+ subroutine block_face_prolong_iface(dir, pos, qn, qb, status)
+ implicit none
+ integer , intent(in) :: dir
+ integer , dimension(3) , intent(in) :: pos
+ real(kind=8), dimension(:,:,:,:), intent(in) :: qn
+ real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+ integer , intent(out) :: status
+ end subroutine
+#endif /* NDIMS == 3 */
+ subroutine block_edge_prolong_iface(dir, pos, qn, qb, status)
+ implicit none
+ integer , intent(in) :: dir
+ integer , dimension(3) , intent(in) :: pos
+ real(kind=8), dimension(:,:,:,:), intent(in) :: qn
+ real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+ integer , intent(out) :: status
+ end subroutine
+ subroutine block_corner_prolong_iface(pos, qn, qb, status)
+ implicit none
+ integer , dimension(3) , intent(in) :: pos
+ real(kind=8), dimension(:,:,:,:), intent(in) :: qn
+ real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+ integer , intent(out) :: status
+ end subroutine
end interface
! pointers to custom boundary conditions, one per direction handling both sides
@@ -53,6 +79,14 @@ module boundaries
custom_boundary_x => null(), &
custom_boundary_y => null(), &
custom_boundary_z => null()
+#if NDIMS == 3
+ procedure(block_face_prolong_iface), pointer, save :: &
+ block_face_prolong => block_face_prolong_tvd
+#endif /* NDIMS == 3 */
+ procedure(block_edge_prolong_iface), pointer, save :: &
+ block_edge_prolong => block_edge_prolong_tvd
+ procedure(block_corner_prolong_iface), pointer, save :: &
+ block_corner_prolong => block_corner_prolong_tvd
! supported boundary types
!
@@ -108,12 +142,13 @@ module boundaries
!
subroutine initialize_boundaries(status)
- use coordinates, only : periodic
- use helpers , only : print_message
+ use coordinates , only : periodic
+ use helpers , only : print_message
+ use interpolations, only : prolongation_method
#ifdef MPI
- use mpitools , only : npmax
+ use mpitools , only : npmax
#endif /* MPI */
- use parameters , only : get_parameter
+ use parameters , only : get_parameter
implicit none
@@ -129,6 +164,21 @@ module boundaries
!
status = 0
+ select case(trim(prolongation_method))
+ case ("ppm", "PPM")
+#if NDIMS == 3
+ block_face_prolong => block_face_prolong_ppm
+#endif /* NDIMS == 3 */
+ block_edge_prolong => block_edge_prolong_ppm
+ block_corner_prolong => block_corner_prolong_ppm
+ case default
+#if NDIMS == 3
+ block_face_prolong => block_face_prolong_tvd
+#endif /* NDIMS == 3 */
+ block_edge_prolong => block_edge_prolong_tvd
+ block_corner_prolong => block_corner_prolong_tvd
+ end select
+
do n = 1, NDIMS
do s = 1, 2
@@ -3793,11 +3843,12 @@ module boundaries
!
!===============================================================================
!
-! subroutine BLOCK_FACE_PROLONG:
-! -----------------------------
+! subroutine BLOCK_FACE_PROLONG_TVD:
+! ---------------------------------
!
-! Subroutine prolongs the region of qn specified by the face position
-! and inserts it in the corresponding face boundary region of qb.
+! This subroutine uses Total Variation Diminishing Method (TVD) to prolong
+! a region of qn, specified by the face position, and inserts it into
+! the corresponding face boundary region of qb.
!
! Arguments:
!
@@ -3809,7 +3860,7 @@ module boundaries
!
!===============================================================================
!
- subroutine block_face_prolong(dir, pos, qn, qb, status)
+ subroutine block_face_prolong_tvd(dir, pos, qn, qb, status)
use coordinates , only : nb, ne, ncells_half, nghosts_half
use equations , only : nv, positive
@@ -3835,7 +3886,7 @@ module boundaries
logical , save :: first = .true.
integer, dimension(2,2), save :: nlims, tlims
- character(len=*), parameter :: loc = 'BOUNDARIES::block_face_prolong()'
+ character(len=*), parameter :: loc = 'BOUNDARIES::block_face_prolong_tvd()'
!-------------------------------------------------------------------------------
!
@@ -3943,7 +3994,190 @@ module boundaries
!-------------------------------------------------------------------------------
!
- end subroutine block_face_prolong
+ end subroutine block_face_prolong_tvd
+!
+!===============================================================================
+!
+! subroutine BLOCK_FACE_PROLONG_PPM:
+! ---------------------------------
+!
+! This subroutine uses Colella and Woodward's Piecewise Parabolic Method (PPM)
+! to prolong a region of qn, specified by the face position, and inserts it
+! into the corresponding face boundary region of qb.
+!
+! Arguments:
+!
+! dir - the face direction;
+! pos - the face position;
+! qn - the array to prolong from;
+! qb - the array of the face boundary region;
+! status - the call status;
+!
+! References:
+!
+! [1] Colella, P., Woodward, P. R.,
+! "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
+! Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
+! https://doi.org/10.1016/0021-9991(84)90143-8
+!
+!===============================================================================
+!
+ subroutine block_face_prolong_ppm(dir, pos, qn, qb, status)
+
+ use coordinates , only : nb, ne, ncells_half, nghosts_half
+ use equations , only : nv, positive
+ use helpers , only : print_message
+ use interpolations, only : limiter_prol
+
+ implicit none
+
+ integer , intent(in) :: dir
+ integer , dimension(3) , intent(in) :: pos
+ real(kind=8), dimension(:,:,:,:), intent(in) :: qn
+ real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+ integer , intent(out) :: status
+
+ integer :: p, n
+ integer :: i, il, iu, is, it, im, ip
+ integer :: j, jl, ju, js, jt, jm, jp
+ integer :: k, kl, ku, ks, kt, km, kp
+ real(kind=8) :: qc, dq1, dq2, dq3, dq4
+
+ real(kind=8), dimension(3) :: ql, qr, dq
+
+ logical , save :: first = .true.
+ integer, dimension(2,2), save :: nlims, tlims
+
+ character(len=*), parameter :: loc = 'BOUNDARIES::block_face_prolong_ppm()'
+
+!-------------------------------------------------------------------------------
+!
+ if (first) then
+ nlims(1,1) = ne - nghosts_half + 1
+ nlims(1,2) = ne
+ nlims(2,1) = nb
+ nlims(2,2) = nb + nghosts_half - 1
+
+ tlims(1,1) = nb
+ tlims(1,2) = nb + ncells_half + nghosts_half - 1
+ tlims(2,1) = ne - ncells_half - nghosts_half + 1
+ tlims(2,2) = ne
+
+ first = .false.
+ end if
+
+ status = 0
+
+ select case(dir)
+ case(1)
+ il = nlims(pos(1),1)
+ iu = nlims(pos(1),2)
+ jl = tlims(pos(2),1)
+ ju = tlims(pos(2),2)
+ kl = tlims(pos(3),1)
+ ku = tlims(pos(3),2)
+ case(2)
+ il = tlims(pos(1),1)
+ iu = tlims(pos(1),2)
+ jl = nlims(pos(2),1)
+ ju = nlims(pos(2),2)
+ kl = tlims(pos(3),1)
+ ku = tlims(pos(3),2)
+ case default
+ il = tlims(pos(1),1)
+ iu = tlims(pos(1),2)
+ jl = tlims(pos(2),1)
+ ju = tlims(pos(2),2)
+ kl = nlims(pos(3),1)
+ ku = nlims(pos(3),2)
+ end select
+
+ do k = kl, ku
+ km = k - 1
+ kp = k + 1
+ ks = 2 * (k - kl) + 1
+ kt = ks + 1
+ do j = jl, ju
+ jm = j - 1
+ jp = j + 1
+ js = 2 * (j - jl) + 1
+ jt = js + 1
+ do i = il, iu
+ im = i - 1
+ ip = i + 1
+ is = 2 * (i - il) + 1
+ it = is + 1
+
+ do p = 1, nv
+
+ qc = qn(p,i,j,k)
+
+ ql(1) = (7.0d+00 * (qn(p,i-1,j,k) + qn(p,i ,j,k)) &
+ - (qn(p,i-2,j,k) + qn(p,i+1,j,k))) / 1.2d+01
+ qr(1) = (7.0d+00 * (qn(p,i ,j,k) + qn(p,i+1,j,k)) &
+ - (qn(p,i-1,j,k) + qn(p,i+2,j,k))) / 1.2d+01
+ ql(2) = (7.0d+00 * (qn(p,i,j-1,k) + qn(p,i,j ,k)) &
+ - (qn(p,i,j-2,k) + qn(p,i,j+1,k))) / 1.2d+01
+ qr(2) = (7.0d+00 * (qn(p,i,j ,k) + qn(p,i,j+1,k)) &
+ - (qn(p,i,j-1,k) + qn(p,i,j+2,k))) / 1.2d+01
+#if NDIMS == 3
+ ql(3) = (7.0d+00 * (qn(p,i,j,k-1) + qn(p,i,j,k )) &
+ - (qn(p,i,j,k-2) + qn(p,i,j,k+1))) / 1.2d+01
+ qr(3) = (7.0d+00 * (qn(p,i,j,k ) + qn(p,i,j,k+1)) &
+ - (qn(p,i,j,k-1) + qn(p,i,j,k+2))) / 1.2d+01
+#endif /* NDIMS == 3 */
+
+ do n = 1, NDIMS
+ if ((qr(n) - qc) * (ql(n) - qc) > 0.0d+00) then
+ ql(n) = qc
+ qr(n) = qc
+ else
+ if (abs(qr(n) - qc) >= 2.0d+00 * abs(ql(n) - qc)) then
+ qr(n) = qc - 2.0d+00 * (ql(n) - qc)
+ end if
+ if (abs(ql(n) - qc) >= 2.0d+00 * abs(qr(n) - qc)) then
+ ql(n) = qc - 2.0d+00 * (qr(n) - qc)
+ end if
+ end if
+ end do
+
+ dq(:) = 2.5d-01 * (qr(:) - ql(:))
+
+ if (positive(p)) then
+ if (qn(p,i,j,k) > 0.0d+00) then
+ dq1 = sum(abs(dq(1:NDIMS)))
+ if (qn(p,i,j,k) <= dq1) &
+ dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:)
+ else
+ call print_message(loc, "Positive variable is not positive!")
+ status = 1
+ return
+ end if
+ end if
+
+ dq1 = dq(1) + dq(2) + dq(3)
+ dq2 = dq(1) - dq(2) - dq(3)
+ dq3 = dq(1) - dq(2) + dq(3)
+ dq4 = dq(1) + dq(2) - dq(3)
+
+ qb(p,is,js,ks) = qn(p,i,j,k) - dq1
+ qb(p,it,js,ks) = qn(p,i,j,k) + dq2
+ qb(p,is,jt,ks) = qn(p,i,j,k) - dq3
+ qb(p,it,jt,ks) = qn(p,i,j,k) + dq4
+ qb(p,is,js,kt) = qn(p,i,j,k) - dq4
+ qb(p,it,js,kt) = qn(p,i,j,k) + dq3
+ qb(p,is,jt,kt) = qn(p,i,j,k) - dq2
+ qb(p,it,jt,kt) = qn(p,i,j,k) + dq1
+
+ end do
+
+ end do
+ end do
+ end do
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine block_face_prolong_ppm
#endif /* NDIMS == 3 */
!
!===============================================================================
@@ -4052,11 +4286,12 @@ module boundaries
!
!===============================================================================
!
-! subroutine BLOCK_EDGE_PROLONG:
-! -----------------------------
+! subroutine BLOCK_EDGE_PROLONG_TVD:
+! ---------------------------------
!
-! Subroutine prolongss the region of qn specified by the edge position
-! and inserts it in the corresponding edge boundary region of qb.
+! This subroutine uses Total Variation Diminishing Method (TVD) to prolong
+! a region of qn, specified by the edge position, and inserts it into
+! the corresponding edge boundary region of qb.
!
! Arguments:
!
@@ -4068,7 +4303,7 @@ module boundaries
!
!===============================================================================
!
- subroutine block_edge_prolong(dir, pos, qn, qb, status)
+ subroutine block_edge_prolong_tvd(dir, pos, qn, qb, status)
use coordinates , only : nb, ne, ncells_half, nghosts_half
use equations , only : nv, positive
@@ -4101,7 +4336,7 @@ module boundaries
logical , save :: first = .true.
integer, dimension(2,2), save :: nlims, tlims
- character(len=*), parameter :: loc = 'BOUNDARIES::block_edge_prolong()'
+ character(len=*), parameter :: loc = 'BOUNDARIES::block_edge_prolong_tvd()'
!-------------------------------------------------------------------------------
!
@@ -4227,7 +4462,213 @@ module boundaries
!-------------------------------------------------------------------------------
!
- end subroutine block_edge_prolong
+ end subroutine block_edge_prolong_tvd
+!
+!===============================================================================
+!
+! subroutine BLOCK_EDGE_PROLONG_PPM:
+! ---------------------------------
+!
+! This subroutine uses Colella and Woodward's Piecewise Parabolic Method (PPM)
+! to prolong a region of qn, specified by the edge position, and inserts it
+! into the corresponding edge boundary region of qb.
+!
+! Arguments:
+!
+! dir - the edge direction;
+! pos - the edge position;
+! qn - the array to prolong from;
+! qb - the array of the edge boundary region;
+! status - the call status;
+!
+! References:
+!
+! [1] Colella, P., Woodward, P. R.,
+! "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
+! Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
+! https://doi.org/10.1016/0021-9991(84)90143-8
+!
+!===============================================================================
+!
+ subroutine block_edge_prolong_ppm(dir, pos, qn, qb, status)
+
+ use coordinates , only : nb, ne, ncells_half, nghosts_half
+ use equations , only : nv, positive
+ use helpers , only : print_message
+ use interpolations, only : limiter_prol
+
+ implicit none
+
+ integer , intent(in) :: dir
+ integer , dimension(3) , intent(in) :: pos
+ real(kind=8), dimension(:,:,:,:), intent(in) :: qn
+ real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+ integer , intent(out) :: status
+
+ integer :: p, n
+ integer :: i, il, iu, is, it, im, ip
+ integer :: j, jl, ju, js, jt, jm, jp
+#if NDIMS == 3
+ integer :: k, kl, ku, ks, kt, km, kp
+#else /* NDIMS == 3 */
+ integer :: k
+#endif /* NDIMS == 3 */
+ real(kind=8) :: qc, dq1, dq2
+#if NDIMS == 3
+ real(kind=8) :: dq3, dq4
+#endif /* NDIMS == 3 */
+
+ real(kind=8), dimension(NDIMS) :: ql, qr, dq
+
+ logical , save :: first = .true.
+ integer, dimension(2,2), save :: nlims, tlims
+
+ character(len=*), parameter :: loc = 'BOUNDARIES::block_edge_prolong_ppm()'
+
+!-------------------------------------------------------------------------------
+!
+ if (first) then
+ nlims(1,1) = ne - nghosts_half + 1
+ nlims(1,2) = ne
+ nlims(2,1) = nb
+ nlims(2,2) = nb + nghosts_half - 1
+
+ tlims(1,1) = nb
+ tlims(1,2) = nb + ncells_half + nghosts_half - 1
+ tlims(2,1) = ne - ncells_half - nghosts_half + 1
+ tlims(2,2) = ne
+
+ first = .false.
+ end if
+
+ status = 0
+
+ il = nlims(pos(1),1)
+ iu = nlims(pos(1),2)
+ jl = nlims(pos(2),1)
+ ju = nlims(pos(2),2)
+#if NDIMS == 3
+ kl = nlims(pos(3),1)
+ ku = nlims(pos(3),2)
+#endif /* NDIMS == 3 */
+
+ select case(dir)
+ case(1)
+ il = tlims(pos(1),1)
+ iu = tlims(pos(1),2)
+ case(2)
+ jl = tlims(pos(2),1)
+ ju = tlims(pos(2),2)
+#if NDIMS == 3
+ case(3)
+ kl = tlims(pos(3),1)
+ ku = tlims(pos(3),2)
+#endif /* NDIMS == 3 */
+ end select
+
+#if NDIMS == 3
+ do k = kl, ku
+ km = k - 1
+ kp = k + 1
+ ks = 2 * (k - kl) + 1
+ kt = ks + 1
+#else /* NDIMS == 3 */
+ k = 1
+#endif /* NDIMS == 3 */
+ do j = jl, ju
+ jm = j - 1
+ jp = j + 1
+ js = 2 * (j - jl) + 1
+ jt = js + 1
+ do i = il, iu
+ im = i - 1
+ ip = i + 1
+ is = 2 * (i - il) + 1
+ it = is + 1
+
+ do p = 1, nv
+
+ qc = qn(p,i,j,k)
+
+ ql(1) = (7.0d+00 * (qn(p,i-1,j,k) + qn(p,i ,j,k)) &
+ - (qn(p,i-2,j,k) + qn(p,i+1,j,k))) / 1.2d+01
+ qr(1) = (7.0d+00 * (qn(p,i ,j,k) + qn(p,i+1,j,k)) &
+ - (qn(p,i-1,j,k) + qn(p,i+2,j,k))) / 1.2d+01
+ ql(2) = (7.0d+00 * (qn(p,i,j-1,k) + qn(p,i,j ,k)) &
+ - (qn(p,i,j-2,k) + qn(p,i,j+1,k))) / 1.2d+01
+ qr(2) = (7.0d+00 * (qn(p,i,j ,k) + qn(p,i,j+1,k)) &
+ - (qn(p,i,j-1,k) + qn(p,i,j+2,k))) / 1.2d+01
+#if NDIMS == 3
+ ql(3) = (7.0d+00 * (qn(p,i,j,k-1) + qn(p,i,j,k )) &
+ - (qn(p,i,j,k-2) + qn(p,i,j,k+1))) / 1.2d+01
+ qr(3) = (7.0d+00 * (qn(p,i,j,k ) + qn(p,i,j,k+1)) &
+ - (qn(p,i,j,k-1) + qn(p,i,j,k+2))) / 1.2d+01
+#endif /* NDIMS == 3 */
+
+ do n = 1, NDIMS
+ if ((qr(n) - qc) * (ql(n) - qc) > 0.0d+00) then
+ ql(n) = qc
+ qr(n) = qc
+ else
+ if (abs(qr(n) - qc) >= 2.0d+00 * abs(ql(n) - qc)) then
+ qr(n) = qc - 2.0d+00 * (ql(n) - qc)
+ end if
+ if (abs(ql(n) - qc) >= 2.0d+00 * abs(qr(n) - qc)) then
+ ql(n) = qc - 2.0d+00 * (qr(n) - qc)
+ end if
+ end if
+ end do
+
+ dq(:) = 2.5d-01 * (qr(:) - ql(:))
+
+ if (positive(p)) then
+ if (qn(p,i,j,k) > 0.0d+00) then
+ dq1 = sum(abs(dq(1:NDIMS)))
+ if (qn(p,i,j,k) <= dq1) &
+ dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:)
+ else
+ call print_message(loc, "Positive variable is not positive!")
+ status = 1
+ return
+ end if
+ end if
+
+#if NDIMS == 2
+ dq1 = dq(1) + dq(2)
+ dq2 = dq(1) - dq(2)
+
+ qb(p,is,js,k ) = qn(p,i,j,k) - dq1
+ qb(p,it,js,k ) = qn(p,i,j,k) + dq2
+ qb(p,is,jt,k ) = qn(p,i,j,k) - dq2
+ qb(p,it,jt,k ) = qn(p,i,j,k) + dq1
+#endif /* NDIMS == 2 */
+#if NDIMS == 3
+ dq1 = dq(1) + dq(2) + dq(3)
+ dq2 = dq(1) - dq(2) - dq(3)
+ dq3 = dq(1) - dq(2) + dq(3)
+ dq4 = dq(1) + dq(2) - dq(3)
+
+ qb(p,is,js,ks) = qn(p,i,j,k) - dq1
+ qb(p,it,js,ks) = qn(p,i,j,k) + dq2
+ qb(p,is,jt,ks) = qn(p,i,j,k) - dq3
+ qb(p,it,jt,ks) = qn(p,i,j,k) + dq4
+ qb(p,is,js,kt) = qn(p,i,j,k) - dq4
+ qb(p,it,js,kt) = qn(p,i,j,k) + dq3
+ qb(p,is,jt,kt) = qn(p,i,j,k) - dq2
+ qb(p,it,jt,kt) = qn(p,i,j,k) + dq1
+#endif /* NDIMS == 3 */
+
+ end do
+
+ end do
+ end do
+#if NDIMS == 3
+ end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine block_edge_prolong_ppm
!
!===============================================================================
!
@@ -4322,11 +4763,12 @@ module boundaries
!
!===============================================================================
!
-! subroutine BLOCK_CORNER_PROLONG:
-! -------------------------------
+! subroutine BLOCK_CORNER_PROLONG_TVD:
+! -----------------------------------
!
-! Subroutine prolongs the region of qn specified by the corner position
-! and inserts it in the corresponding corner boundary region of qb.
+! This subroutine uses Total Variation Diminishing Method (TVD) to prolong
+! a region of qn, specified by the corner position, and inserts it into
+! the corresponding corner boundary region of qb.
!
! Remarks:
!
@@ -4342,7 +4784,7 @@ module boundaries
!
!===============================================================================
!
- subroutine block_corner_prolong(pos, qn, qb, status)
+ subroutine block_corner_prolong_tvd(pos, qn, qb, status)
use coordinates , only : nb, ne, nghosts_half
use equations , only : nv, positive
@@ -4374,7 +4816,7 @@ module boundaries
logical , save :: first = .true.
integer, dimension(2,2), save :: nlims
- character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong()'
+ character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong_tvd()'
!-------------------------------------------------------------------------------
!
@@ -4479,7 +4921,195 @@ module boundaries
!-------------------------------------------------------------------------------
!
- end subroutine block_corner_prolong
+ end subroutine block_corner_prolong_tvd
+!
+!===============================================================================
+!
+! subroutine BLOCK_CORNER_PROLONG_PPM:
+! -----------------------------------
+!
+! This subroutine uses Colella and Woodward's Piecewise Parabolic Method (PPM)
+! to prolong a region of qn, specified by the corner position, and inserts it
+! into the corresponding corner boundary region of qb.
+!
+! Remarks:
+!
+! This subroutine requires the ghost zone regions of the qn already
+! updated in order to perform a consistent interpolation.
+!
+! Arguments:
+!
+! pos - the corner position;
+! qn - the array to prolong from;
+! qb - the array of the corner boundary region;
+! status - the call status;
+!
+! References:
+!
+! [1] Colella, P., Woodward, P. R.,
+! "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
+! Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
+! https://doi.org/10.1016/0021-9991(84)90143-8
+!
+!===============================================================================
+!
+ subroutine block_corner_prolong_ppm(pos, qn, qb, status)
+
+ use coordinates , only : nb, ne, nghosts_half
+ use equations , only : nv, positive
+ use helpers , only : print_message
+ use interpolations, only : limiter_prol
+
+ implicit none
+
+ integer , dimension(3) , intent(in) :: pos
+ real(kind=8), dimension(:,:,:,:), intent(in) :: qn
+ real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+ integer , intent(out) :: status
+
+ integer :: p, n
+ integer :: i, il, iu, im, ip, is, it
+ integer :: j, jl, ju, jm, jp, js, jt
+#if NDIMS == 3
+ integer :: k, kl, ku, km, kp, ks, kt
+#else /* NDIMS == 3 */
+ integer :: k
+#endif /* NDIMS == 3 */
+ real(kind=8) :: qc, dq1, dq2
+#if NDIMS == 3
+ real(kind=8) :: dq3, dq4
+#endif /* NDIMS == 3 */
+
+ real(kind=8), dimension(NDIMS) :: ql, qr, dq
+
+ logical , save :: first = .true.
+ integer, dimension(2,2), save :: nlims
+
+ character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong_ppm()'
+
+!-------------------------------------------------------------------------------
+!
+ if (first) then
+ nlims(1,1) = ne - nghosts_half + 1
+ nlims(1,2) = ne
+ nlims(2,1) = nb
+ nlims(2,2) = nb + nghosts_half - 1
+
+ first = .false.
+ end if
+
+ status = 0
+
+ il = nlims(pos(1),1)
+ iu = nlims(pos(1),2)
+ jl = nlims(pos(2),1)
+ ju = nlims(pos(2),2)
+#if NDIMS == 3
+ kl = nlims(pos(3),1)
+ ku = nlims(pos(3),2)
+
+ do k = kl, ku
+ km = k - 1
+ kp = k + 1
+ ks = 2 * (k - kl) + 1
+ kt = ks + 1
+#else /* NDIMS == 3 */
+ k = 1
+#endif /* NDIMS == 3 */
+ do j = jl, ju
+ jm = j - 1
+ jp = j + 1
+ js = 2 * (j - jl) + 1
+ jt = js + 1
+ do i = il, iu
+ im = i - 1
+ ip = i + 1
+ is = 2 * (i - il) + 1
+ it = is + 1
+
+ do p = 1, nv
+
+ qc = qn(p,i,j,k)
+
+ ql(1) = (7.0d+00 * (qn(p,i-1,j,k) + qn(p,i ,j,k)) &
+ - (qn(p,i-2,j,k) + qn(p,i+1,j,k))) / 1.2d+01
+ qr(1) = (7.0d+00 * (qn(p,i ,j,k) + qn(p,i+1,j,k)) &
+ - (qn(p,i-1,j,k) + qn(p,i+2,j,k))) / 1.2d+01
+ ql(2) = (7.0d+00 * (qn(p,i,j-1,k) + qn(p,i,j ,k)) &
+ - (qn(p,i,j-2,k) + qn(p,i,j+1,k))) / 1.2d+01
+ qr(2) = (7.0d+00 * (qn(p,i,j ,k) + qn(p,i,j+1,k)) &
+ - (qn(p,i,j-1,k) + qn(p,i,j+2,k))) / 1.2d+01
+#if NDIMS == 3
+ ql(3) = (7.0d+00 * (qn(p,i,j,k-1) + qn(p,i,j,k )) &
+ - (qn(p,i,j,k-2) + qn(p,i,j,k+1))) / 1.2d+01
+ qr(3) = (7.0d+00 * (qn(p,i,j,k ) + qn(p,i,j,k+1)) &
+ - (qn(p,i,j,k-1) + qn(p,i,j,k+2))) / 1.2d+01
+#endif /* NDIMS == 3 */
+
+ do n = 1, NDIMS
+ if ((qr(n) - qc) * (ql(n) - qc) > 0.0d+00) then
+ ql(n) = qc
+ qr(n) = qc
+ else
+ if (abs(qr(n) - qc) >= 2.0d+00 * abs(ql(n) - qc)) then
+ qr(n) = qc - 2.0d+00 * (ql(n) - qc)
+ end if
+ if (abs(ql(n) - qc) >= 2.0d+00 * abs(qr(n) - qc)) then
+ ql(n) = qc - 2.0d+00 * (qr(n) - qc)
+ end if
+ end if
+ end do
+
+ dq(:) = 2.5d-01 * (qr(:) - ql(:))
+
+ if (positive(p)) then
+ if (qn(p,i,j,k) > 0.0d+00) then
+ dq1 = sum(abs(dq(1:NDIMS)))
+ if (qn(p,i,j,k) <= dq1) &
+ dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:)
+ else
+ call print_message(loc, "Positive variable is not positive!")
+ status = 1
+ return
+ end if
+ end if
+
+#if NDIMS == 2
+ dq1 = dq(1) + dq(2)
+ dq2 = dq(1) - dq(2)
+
+ qb(p,is,js,k ) = qn(p,i,j,k) - dq1
+ qb(p,it,js,k ) = qn(p,i,j,k) + dq2
+ qb(p,is,jt,k ) = qn(p,i,j,k) - dq2
+ qb(p,it,jt,k ) = qn(p,i,j,k) + dq1
+#endif /* NDIMS == 2 */
+#if NDIMS == 3
+ dq1 = dq(1) + dq(2) + dq(3)
+ dq2 = dq(1) - dq(2) - dq(3)
+ dq3 = dq(1) - dq(2) + dq(3)
+ dq4 = dq(1) + dq(2) - dq(3)
+
+ qb(p,is,js,ks) = qn(p,i,j,k) - dq1
+ qb(p,it,js,ks) = qn(p,i,j,k) + dq2
+ qb(p,is,jt,ks) = qn(p,i,j,k) - dq3
+ qb(p,it,jt,ks) = qn(p,i,j,k) + dq4
+ qb(p,is,js,kt) = qn(p,i,j,k) - dq4
+ qb(p,it,js,kt) = qn(p,i,j,k) + dq3
+ qb(p,is,jt,kt) = qn(p,i,j,k) - dq2
+ qb(p,it,jt,kt) = qn(p,i,j,k) + dq1
+#endif /* NDIMS == 3 */
+
+ end do
+
+ end do
+ end do
+#if NDIMS == 3
+ end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine block_corner_prolong_ppm
!
!===============================================================================
!
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 6b937bf..6f7b3a9 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -79,6 +79,7 @@ module interpolations
! method names
!
+ character(len=4) :: prolongation_method = "TVD"
character(len=255), save :: name_rec = ""
character(len=255), save :: name_tlim = ""
character(len=255), save :: name_plim = ""
@@ -169,7 +170,7 @@ module interpolations
public :: print_interpolations
public :: interfaces, reconstruct, limiter_prol
public :: fix_positivity
- public :: order
+ public :: order, prolongation_method
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
@@ -242,6 +243,13 @@ module interpolations
call get_parameter("ppm_const" , ppm_const )
call get_parameter("central_weight" , cweight )
call get_parameter("cfl" , cfl )
+ call get_parameter("prolongation_method" , prolongation_method)
+ select case(trim(prolongation_method))
+ case ("ppm", "PPM")
+ prolongation_method = 'PPM'
+ case default
+ prolongation_method = 'TVD'
+ end select
! calculate κ = (1 - ν) / ν
!
@@ -811,7 +819,9 @@ module interpolations
call print_parameter(verbose, "WENO sensitivity", weno_sensitivity)
end if
call print_parameter(verbose, "TVD limiter" , name_tlim)
- call print_parameter(verbose, "prolongation limiter", name_plim)
+ call print_parameter(verbose, "prolongation method" , prolongation_method)
+ if (prolongation_method == 'TVD') &
+ call print_parameter(verbose, "prolongation limiter", name_plim)
if (mlp) then
call print_parameter(verbose, "MLP limiting" , "on" )
else
diff --git a/sources/mesh.F90 b/sources/mesh.F90
index 2e3b143..9351a8c 100644
--- a/sources/mesh.F90
+++ b/sources/mesh.F90
@@ -46,6 +46,15 @@ module mesh
end subroutine
end interface
+ abstract interface
+ subroutine prolong_block_iface(pmeta, status)
+ use blocks, only : block_meta
+ implicit none
+ type(block_meta), pointer, intent(inout) :: pmeta
+ integer , intent(out) :: status
+ end subroutine
+ end interface
+
! pointer to the problem domain setup subroutine
!
procedure(setup_domain_iface) , pointer, save :: setup_domain => null()
@@ -54,6 +63,11 @@ module mesh
!
procedure(setup_problem_iface), pointer, save :: setup_problem => null()
+! pointer to the block prolongation subroutine
+!
+ procedure(prolong_block_iface), pointer, save :: &
+ prolong_block => prolong_block_tvd
+
! the flag indicating that the block structure or distribution has changed
!
logical, save :: changed = .true.
@@ -96,8 +110,9 @@ module mesh
!
subroutine initialize_mesh(verbose, status)
- use coordinates, only : ncells, nghosts
- use refinement , only : initialize_refinement
+ use coordinates , only : ncells, nghosts
+ use interpolations, only : prolongation_method
+ use refinement , only : initialize_refinement
implicit none
@@ -108,6 +123,13 @@ module mesh
!
status = 0
+ select case(trim(prolongation_method))
+ case ("ppm", "PPM")
+ prolong_block => prolong_block_ppm
+ case default
+ prolong_block => prolong_block_tvd
+ end select
+
setup_domain => setup_domain_default
pm(1:NDIMS) = 2 * (ncells + nghosts)
@@ -1774,13 +1796,11 @@ module mesh
!
!===============================================================================
!
-! subroutine PROLONG_BLOCK:
-! ------------------------
+! subroutine PROLONG_BLOCK_TVD:
+! ----------------------------
!
-! Subroutine prolongs variables from a data blocks linked to the input
-! meta block and copy the resulting array of variables to
-! the newly created children data block. The process of data restriction
-! conserves stored variables.
+! This subroutine prolongs conserved variables from a parent data block
+! to its children using the TVD method with a selected prolongation limiter.
!
! Arguments:
!
@@ -1789,7 +1809,7 @@ module mesh
!
!===============================================================================
!
- subroutine prolong_block(pmeta, status)
+ subroutine prolong_block_tvd(pmeta, status)
use blocks , only : block_meta, block_data, nchildren
use coordinates , only : ni => ncells, nn => bcells
@@ -1836,7 +1856,7 @@ module mesh
!$ integer :: omp_get_thread_num
!$omp threadprivate(first, nt, tmp)
- character(len=*), parameter :: loc = 'MESH::prolong_block()'
+ character(len=*), parameter :: loc = 'MESH::prolong_block_tvd()'
!-------------------------------------------------------------------------------
!
@@ -1972,7 +1992,233 @@ module mesh
!-------------------------------------------------------------------------------
!
- end subroutine prolong_block
+ end subroutine prolong_block_tvd
+!
+!===============================================================================
+!
+! subroutine PROLONG_BLOCK_PPM:
+! ----------------------------
+!
+! This subroutine prolongs conserved variables from a parent data block to
+! its children using Colella and Woodward's Piecewise Parabolic Method (PPM).
+!
+! Arguments:
+!
+! pmeta - the input meta block;
+! status - the subroutine call status: 0 for success, otherwise failure;
+!
+! References:
+!
+! [1] Colella, P., Woodward, P. R.,
+! "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
+! Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
+! https://doi.org/10.1016/0021-9991(84)90143-8
+!
+!===============================================================================
+!
+ subroutine prolong_block_ppm(pmeta, status)
+
+ use blocks , only : block_meta, block_data, nchildren
+ use coordinates , only : ni => ncells, nn => bcells
+ use coordinates , only : nh => nghosts_half
+ use coordinates , only : nb, ne
+ use equations , only : nv, positive
+ use helpers , only : print_message
+ use interpolations, only : limiter_prol
+ use workspace , only : resize_workspace, work, work_in_use
+
+ implicit none
+
+ type(block_meta), pointer, intent(inout) :: pmeta
+ integer , intent(out) :: status
+
+ type(block_meta), pointer :: pchild
+ type(block_data), pointer :: pdata
+
+ logical, save :: first = .true.
+
+ integer :: n, p, nl, nu
+ integer :: i, im, ip
+ integer :: j, jm, jp
+#if NDIMS == 3
+ integer :: k, km, kp
+#else /* NDIMS == 3 */
+ integer :: k
+#endif /* NDIMS == 3 */
+ real(kind=8) :: uc
+#if NDIMS == 3
+ real(kind=8) :: du1, du2, du3, du4
+#else /* NDIMS == 3 */
+ real(kind=8) :: du1, du2
+#endif /* NDIMS == 3 */
+
+ character(len=80) :: msg
+
+ integer , dimension(NDIMS) :: l, u
+ real(kind=8), dimension(NDIMS) :: ul, ur, du
+
+ real(kind=8), dimension(:,:,:), pointer, save :: tmp
+
+ integer, save :: nt
+!$ integer :: omp_get_thread_num
+!$omp threadprivate(first, nt, tmp)
+
+ character(len=*), parameter :: loc = 'MESH::prolong_block_ppm()'
+
+!-------------------------------------------------------------------------------
+!
+ nt = 0
+!$ nt = omp_get_thread_num()
+ status = 0
+
+ if (first) then
+ n = product(pm(:))
+
+ call resize_workspace(n, status)
+ if (status /= 0) then
+ call print_message(loc, "Could not resize the workspace!")
+ go to 100
+ end if
+
+ tmp(1:pm(1),1:pm(2),1:pm(3)) => work(1:n,nt)
+
+ first = .false.
+ end if
+
+ if (work_in_use(nt)) &
+ call print_message(loc, "Workspace is being used right now! " // &
+ "Corruptions can occur!")
+
+ work_in_use(nt) = .true.
+
+#if NDIMS == 2
+ k = 1
+#endif /* NDIMS == 2 */
+
+ nl = nb - nh
+ nu = ne + nh
+
+ pdata => pmeta%data
+
+ do n = 1, nv
+
+#if NDIMS == 3
+ do k = nl, nu
+ km = k - 1
+ kp = k + 1
+ l(3) = 2 * (k - nl) + 1
+ u(3) = l(3) + 1
+#endif /* NDIMS == 3 */
+ do j = nl, nu
+ jm = j - 1
+ jp = j + 1
+ l(2) = 2 * (j - nl) + 1
+ u(2) = l(2) + 1
+ do i = nl, nu
+ im = i - 1
+ ip = i + 1
+ l(1) = 2 * (i - nl) + 1
+ u(1) = l(1) + 1
+
+ uc = pdata%u(n,i,j,k)
+
+ ul(1) = (7.0d+00 * (pdata%u(n,i-1,j,k) + pdata%u(n,i ,j,k)) &
+ - (pdata%u(n,i-2,j,k) + pdata%u(n,i+1,j,k))) / 1.2d+01
+ ur(1) = (7.0d+00 * (pdata%u(n,i ,j,k) + pdata%u(n,i+1,j,k)) &
+ - (pdata%u(n,i-1,j,k) + pdata%u(n,i+2,j,k))) / 1.2d+01
+ ul(2) = (7.0d+00 * (pdata%u(n,i,j-1,k) + pdata%u(n,i,j ,k)) &
+ - (pdata%u(n,i,j-2,k) + pdata%u(n,i,j+1,k))) / 1.2d+01
+ ur(2) = (7.0d+00 * (pdata%u(n,i,j ,k) + pdata%u(n,i,j+1,k)) &
+ - (pdata%u(n,i,j-1,k) + pdata%u(n,i,j+2,k))) / 1.2d+01
+#if NDIMS == 3
+ ul(3) = (7.0d+00 * (pdata%u(n,i,j,k-1) + pdata%u(n,i,j,k )) &
+ - (pdata%u(n,i,j,k-2) + pdata%u(n,i,j,k+1))) / 1.2d+01
+ ur(3) = (7.0d+00 * (pdata%u(n,i,j,k ) + pdata%u(n,i,j,k+1)) &
+ - (pdata%u(n,i,j,k-1) + pdata%u(n,i,j,k+2))) / 1.2d+01
+#endif /* NDIMS == 3 */
+
+ do p = 1, NDIMS
+ if ((ur(p) - uc) * (ul(p) - uc) > 0.0d+00) then
+ ul(p) = uc
+ ur(p) = uc
+ else
+ if (abs(ur(p) - uc) >= 2.0d+00 * abs(ul(p) - uc)) then
+ ur(p) = uc - 2.0d+00 * (ul(p) - uc)
+ end if
+ if (abs(ul(p) - uc) >= 2.0d+00 * abs(ur(p) - uc)) then
+ ul(p) = uc - 2.0d+00 * (ur(p) - uc)
+ end if
+ end if
+ end do
+
+ du(:) = 2.5d-01 * (ur(:) - ul(:))
+
+ if (positive(n) .and. pdata%u(n,i,j,k) < sum(abs(du(1:NDIMS)))) then
+ if (pdata%u(n,i,j,k) > 0.0d+00) then
+ do while (pdata%u(n,i,j,k) <= sum(abs(du(1:NDIMS))))
+ du(:) = 0.5d+00 * du(:)
+ end do
+ else
+ write(msg,"(a,3i4,a)") &
+ "Positive variable is not positive at (", i, j, k, " )"
+ call print_message(loc, msg)
+ du(:) = 0.0d+00
+ end if
+ end if
+
+#if NDIMS == 2
+ du1 = du(1) + du(2)
+ du2 = du(1) - du(2)
+ tmp(l(1),l(2),k) = pdata%u(n,i,j,k) - du1
+ tmp(u(1),l(2),k) = pdata%u(n,i,j,k) + du2
+ tmp(l(1),u(2),k) = pdata%u(n,i,j,k) - du2
+ tmp(u(1),u(2),k) = pdata%u(n,i,j,k) + du1
+#endif /* NDIMS == 2 */
+
+#if NDIMS == 3
+ du1 = du(1) + du(2) + du(3)
+ du2 = du(1) - du(2) - du(3)
+ du3 = du(1) - du(2) + du(3)
+ du4 = du(1) + du(2) - du(3)
+ tmp(l(1),l(2),l(3)) = pdata%u(n,i,j,k) - du1
+ tmp(u(1),l(2),l(3)) = pdata%u(n,i,j,k) + du2
+ tmp(l(1),u(2),l(3)) = pdata%u(n,i,j,k) - du3
+ tmp(u(1),u(2),l(3)) = pdata%u(n,i,j,k) + du4
+ tmp(l(1),l(2),u(3)) = pdata%u(n,i,j,k) - du4
+ tmp(u(1),l(2),u(3)) = pdata%u(n,i,j,k) + du3
+ tmp(l(1),u(2),u(3)) = pdata%u(n,i,j,k) - du2
+ tmp(u(1),u(2),u(3)) = pdata%u(n,i,j,k) + du1
+#endif /* NDIMS == 3 */
+ end do
+ end do
+#if NDIMS == 3
+ end do
+#endif /* NDIMS == 3 */
+
+ do p = 1, nchildren
+
+ pchild => pmeta%child(p)%ptr
+
+ l(1:NDIMS) = pchild%pos(1:NDIMS) * ni + 1
+ u(1:NDIMS) = l(1:NDIMS) + nn - 1
+
+#if NDIMS == 2
+ pchild%data%u(n,:,:,1) = tmp(l(1):u(1),l(2):u(2),k)
+#endif /* NDIMS == 2 */
+#if NDIMS == 3
+ pchild%data%u(n,:,:,:) = tmp(l(1):u(1),l(2):u(2),l(3):u(3))
+#endif /* NDIMS == 3 */
+
+ end do ! nchildren
+ end do ! n = 1, nv
+
+ work_in_use(nt) = .false.
+
+ 100 continue
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine prolong_block_ppm
!
!===============================================================================
!