gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2024-10-07 10:42:24 -03:00
commit d140e34612
3 changed files with 924 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

672
sources/boundaries.F90
View File

@ -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
!
@ -110,6 +144,7 @@ module boundaries
use coordinates , only : periodic
use helpers , only : print_message
use interpolations, only : prolongation_method
#ifdef MPI
use mpitools , only : npmax
#endif /* MPI */
@ -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
!
!===============================================================================
!

12
sources/interpolations.F90
View File

@ -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,6 +819,8 @@ 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 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" )

264
sources/mesh.F90
View File

@ -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.
@ -97,6 +111,7 @@ module mesh
subroutine initialize_mesh(verbose, status)
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
!
!===============================================================================
!