BOUNDARIES: Prepare primitive variables for specific boundaries.

The specific boundaries (open, outflow, reflective, custom, etc.) work
on the primitive variables. Since the boundaries are updated using the
conservative variables, we have to update the primitive variables in the
block to which the specific boundary is applied. Then, the conservative
variables are update at the ghost zone where the specific boundary
condition was applied.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-02-09 19:25:18 -03:00
parent 47ebb9703c
commit 6f23ae962e

View File

@ -3611,6 +3611,8 @@ module boundaries
integer :: m integer :: m
#endif /* NDIMS == 2 */ #endif /* NDIMS == 2 */
integer, dimension(NDIMS) :: s
real(kind=8), dimension(nn) :: x, y real(kind=8), dimension(nn) :: x, y
#if NDIMS == 3 #if NDIMS == 3
real(kind=8), dimension(nn) :: z real(kind=8), dimension(nn) :: z
@ -3652,13 +3654,15 @@ module boundaries
m = 3 - n m = 3 - n
do j = 1, nsides do j = 1, nsides
s(2) = j
do i = 1, nsides do i = 1, nsides
s(1) = i
if (.not. associated(pmeta%edges(i,j,m)%ptr)) & if (.not. associated(pmeta%edges(i,j,m)%ptr)) then
call block_boundary_specific(n, [ i, j, k ], & call block_primitive_variables(n, pmeta%data, status)
t, dt, x, y, z, & call block_boundary_specific(n, [ i, j, k ], t, dt, &
pmeta%data%u, status) x, y, z, pmeta%data%q, status)
call block_conservative_variables(n, s(n), pmeta%data)
end if
end do end do
end do end do
@ -3672,14 +3676,17 @@ module boundaries
if (.not. periodic(n)) then if (.not. periodic(n)) then
do k = 1, nsides do k = 1, nsides
s(3) = k
do j = 1, nsides do j = 1, nsides
s(2) = j
do i = 1, nsides do i = 1, nsides
s(1) = i
if (.not. associated(pmeta%faces(i,j,k,n)%ptr)) & if (.not. associated(pmeta%faces(i,j,k,n)%ptr)) then
call block_boundary_specific(n, [ i, j, k ], & call block_primitive_variables(n, pmeta%data, status)
t, dt, x, y, z, & call block_boundary_specific(n, [ i, j, k ], t, dt, &
pmeta%data%u, status) x, y, z, pmeta%data%q, status)
call block_conservative_variables(n, s(n), pmeta%data)
end if
end do end do
end do end do
end do end do
@ -4518,9 +4525,9 @@ module boundaries
use coordinates , only : nn => bcells, nh => bcells_half, ng => nghosts use coordinates , only : nn => bcells, nh => bcells_half, ng => nghosts
use coordinates , only : nb, ne, nbl, neu use coordinates , only : nb, ne, nbl, neu
use equations , only : magnetized use equations , only : magnetized
use equations , only : idn, ipr, imx, imy, ibx, iby use equations , only : idn, ipr, ivx, ivy, ibx, iby
#if NDIMS == 3 #if NDIMS == 3
use equations , only : imz, ibz use equations , only : ivz, ibz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
use equations , only : csnd2 use equations , only : csnd2
use gravity , only : gravitational_acceleration use gravity , only : gravitational_acceleration
@ -4582,12 +4589,12 @@ module boundaries
if (side(1) == 1) then if (side(1) == 1) then
do i = nbl, 1, -1 do i = nbl, 1, -1
qn( : ,i,jl:ju,kl:ku) = qn( : ,nb,jl:ju,kl:ku) qn( : ,i,jl:ju,kl:ku) = qn( : ,nb,jl:ju,kl:ku)
qn(imx,i,jl:ju,kl:ku) = min(0.0d+00, qn(imx,nb,jl:ju,kl:ku)) qn(ivx,i,jl:ju,kl:ku) = min(0.0d+00, qn(ivx,nb,jl:ju,kl:ku))
end do end do
else else
do i = neu, nn do i = neu, nn
qn( : ,i,jl:ju,kl:ku) = qn( : ,ne,jl:ju,kl:ku) qn( : ,i,jl:ju,kl:ku) = qn( : ,ne,jl:ju,kl:ku)
qn(imx,i,jl:ju,kl:ku) = max(0.0d+00, qn(imx,ne,jl:ju,kl:ku)) qn(ivx,i,jl:ju,kl:ku) = max(0.0d+00, qn(ivx,ne,jl:ju,kl:ku))
end do end do
end if end if
@ -4599,7 +4606,7 @@ module boundaries
is = nbl + i is = nbl + i
qn( : ,it,jl:ju,kl:ku) = qn( : ,is,jl:ju,kl:ku) qn( : ,it,jl:ju,kl:ku) = qn( : ,is,jl:ju,kl:ku)
qn(imx,it,jl:ju,kl:ku) = - qn(imx,is,jl:ju,kl:ku) qn(ivx,it,jl:ju,kl:ku) = - qn(ivx,is,jl:ju,kl:ku)
if (magnetized) then if (magnetized) then
qn(ibx,it,jl:ju,kl:ku) = - qn(ibx,is,jl:ju,kl:ku) qn(ibx,it,jl:ju,kl:ku) = - qn(ibx,is,jl:ju,kl:ku)
end if end if
@ -4610,7 +4617,7 @@ module boundaries
is = neu - i is = neu - i
qn( : ,it,jl:ju,kl:ku) = qn( : ,is,jl:ju,kl:ku) qn( : ,it,jl:ju,kl:ku) = qn( : ,is,jl:ju,kl:ku)
qn(imx,it,jl:ju,kl:ku) = - qn(imx,is,jl:ju,kl:ku) qn(ivx,it,jl:ju,kl:ku) = - qn(ivx,is,jl:ju,kl:ku)
if (magnetized) then if (magnetized) then
qn(ibx,it,jl:ju,kl:ku) = - qn(ibx,is,jl:ju,kl:ku) qn(ibx,it,jl:ju,kl:ku) = - qn(ibx,is,jl:ju,kl:ku)
end if end if
@ -4741,12 +4748,12 @@ module boundaries
if (side(2) == 1) then if (side(2) == 1) then
do j = nbl, 1, -1 do j = nbl, 1, -1
qn( : ,il:iu,j,kl:ku) = qn( : ,il:iu,nb,kl:ku) qn( : ,il:iu,j,kl:ku) = qn( : ,il:iu,nb,kl:ku)
qn(imy,il:iu,j,kl:ku) = min(0.0d+00, qn(imy,il:iu,nb,kl:ku)) qn(ivy,il:iu,j,kl:ku) = min(0.0d+00, qn(ivy,il:iu,nb,kl:ku))
end do end do
else else
do j = neu, nn do j = neu, nn
qn( : ,il:iu,j,kl:ku) = qn( : ,il:iu,ne,kl:ku) qn( : ,il:iu,j,kl:ku) = qn( : ,il:iu,ne,kl:ku)
qn(imy,il:iu,j,kl:ku) = max(0.0d+00, qn(imy,il:iu,ne,kl:ku)) qn(ivy,il:iu,j,kl:ku) = max(0.0d+00, qn(ivy,il:iu,ne,kl:ku))
end do end do
end if end if
@ -4758,7 +4765,7 @@ module boundaries
js = nbl + j js = nbl + j
qn( : ,il:iu,jt,kl:ku) = qn( : ,il:iu,js,kl:ku) qn( : ,il:iu,jt,kl:ku) = qn( : ,il:iu,js,kl:ku)
qn(imy,il:iu,jt,kl:ku) = - qn(imy,il:iu,js,kl:ku) qn(ivy,il:iu,jt,kl:ku) = - qn(ivy,il:iu,js,kl:ku)
if (magnetized) then if (magnetized) then
qn(iby,il:iu,jt,kl:ku) = - qn(iby,il:iu,js,kl:ku) qn(iby,il:iu,jt,kl:ku) = - qn(iby,il:iu,js,kl:ku)
end if end if
@ -4769,7 +4776,7 @@ module boundaries
js = neu - j js = neu - j
qn( : ,il:iu,jt,kl:ku) = qn( : ,il:iu,js,kl:ku) qn( : ,il:iu,jt,kl:ku) = qn( : ,il:iu,js,kl:ku)
qn(imy,il:iu,jt,kl:ku) = - qn(imy,il:iu,js,kl:ku) qn(ivy,il:iu,jt,kl:ku) = - qn(ivy,il:iu,js,kl:ku)
if (magnetized) then if (magnetized) then
qn(iby,il:iu,jt,kl:ku) = - qn(iby,il:iu,js,kl:ku) qn(iby,il:iu,jt,kl:ku) = - qn(iby,il:iu,js,kl:ku)
end if end if
@ -4896,12 +4903,12 @@ module boundaries
if (side(3) == 1) then if (side(3) == 1) then
do k = nbl, 1, -1 do k = nbl, 1, -1
qn( : ,il:iu,jl:ju,k) = qn( : ,il:iu,jl:ju,nb) qn( : ,il:iu,jl:ju,k) = qn( : ,il:iu,jl:ju,nb)
qn(imz,il:iu,jl:ju,k) = min(0.0d+00, qn(imz,il:iu,jl:ju,nb)) qn(ivz,il:iu,jl:ju,k) = min(0.0d+00, qn(ivz,il:iu,jl:ju,nb))
end do end do
else else
do k = neu, nn do k = neu, nn
qn( : ,il:iu,jl:ju,k) = qn( : ,il:iu,jl:ju,ne) qn( : ,il:iu,jl:ju,k) = qn( : ,il:iu,jl:ju,ne)
qn(imz,il:iu,jl:ju,k) = max(0.0d+00, qn(imz,il:iu,jl:ju,ne)) qn(ivz,il:iu,jl:ju,k) = max(0.0d+00, qn(ivz,il:iu,jl:ju,ne))
end do end do
end if end if
@ -4913,7 +4920,7 @@ module boundaries
ks = nbl + k ks = nbl + k
qn( : ,il:iu,jl:ju,kt) = qn( : ,il:iu,jl:ju,ks) qn( : ,il:iu,jl:ju,kt) = qn( : ,il:iu,jl:ju,ks)
qn(imz,il:iu,jl:ju,kt) = - qn(imz,il:iu,jl:ju,ks) qn(ivz,il:iu,jl:ju,kt) = - qn(ivz,il:iu,jl:ju,ks)
if (magnetized) then if (magnetized) then
qn(ibz,il:iu,jl:ju,kt) = - qn(ibz,il:iu,jl:ju,ks) qn(ibz,il:iu,jl:ju,kt) = - qn(ibz,il:iu,jl:ju,ks)
end if end if
@ -4924,7 +4931,7 @@ module boundaries
ks = neu - k ks = neu - k
qn( : ,il:iu,jl:ju,kt) = qn( : ,il:iu,jl:ju,ks) qn( : ,il:iu,jl:ju,kt) = qn( : ,il:iu,jl:ju,ks)
qn(imz,il:iu,jl:ju,kt) = - qn(imz,il:iu,jl:ju,ks) qn(ivz,il:iu,jl:ju,kt) = - qn(ivz,il:iu,jl:ju,ks)
if (magnetized) then if (magnetized) then
qn(ibz,il:iu,jl:ju,kt) = - qn(ibz,il:iu,jl:ju,ks) qn(ibz,il:iu,jl:ju,kt) = - qn(ibz,il:iu,jl:ju,ks)
end if end if
@ -5038,6 +5045,160 @@ module boundaries
! OTHER BOUNDARY SUBROUTINES ! OTHER BOUNDARY SUBROUTINES
! !
!=============================================================================== !===============================================================================
!
!===============================================================================
!
! subroutine BLOCK_PRIMITIVE_VARIABLES:
! ------------------------------------
!
! Subroutine updates primitive variables in the block interior excluding the
! ghost zones along the specified direction.
!
! Arguments:
!
! dir - the direction;
! side - the side of the boundary;
! pdata - the data block pointer;
! status - the call status;
!
!===============================================================================
!
subroutine block_primitive_variables(dir, pdata, status)
use blocks , only : block_data
use coordinates, only : nn => bcells, nb, ne
use equations , only : cons2prim
implicit none
integer , intent(in) :: dir
integer , intent(out) :: status
type(block_data), pointer, intent(inout) :: pdata
integer :: il, iu
integer :: jl, ju, j
integer :: kl, ku, k
!-------------------------------------------------------------------------------
!
if (dir == 1) then
il = nb
iu = ne
else
il = 1
iu = nn
end if
if (dir == 2) then
jl = nb
ju = ne
else
jl = 1
ju = nn
end if
#if NDIMS == 3
if (dir == 3) then
kl = nb
ku = ne
else
kl = 1
ku = nn
end if
#else /* NDIMS == 3 */
kl = 1
ku = 1
#endif /* NDIMS == 3 */
do k = kl, ku
do j = jl, ju
call cons2prim(pdata%u(:,il:iu,j,k), pdata%q(:,il:iu,j,k), status)
end do
end do
!-------------------------------------------------------------------------------
!
end subroutine block_primitive_variables
!
!===============================================================================
!
! subroutine BLOCK_CONSERVATIVE_VARIABLES:
! ---------------------------------------
!
! Subroutine updates conservative variables in the ghost region specified by
! the direction and side.
!
! Arguments:
!
! dir - the direction;
! side - the side of the boundary;
! pdata - the data block pointer;
!
!===============================================================================
!
subroutine block_conservative_variables(dir, side, pdata)
use blocks , only : block_data
use coordinates, only : nn => bcells, nb, ne
use equations , only : prim2cons
implicit none
integer , intent(in) :: dir, side
type(block_data), pointer, intent(inout) :: pdata
integer :: il, iu
integer :: jl, ju, j
integer :: kl, ku, k
logical , save :: first = .true.
integer, dimension(2,2), save :: nlims
!-------------------------------------------------------------------------------
!
if (first) then
nlims(1,1) = 1
nlims(1,2) = nb - 1
nlims(2,1) = ne + 1
nlims(2,2) = nn
first = .false.
end if
if (dir == 1) then
il = nlims(side,1)
iu = nlims(side,2)
else
il = 1
iu = nn
end if
if (dir == 2) then
jl = nlims(side,1)
ju = nlims(side,2)
else
jl = 1
ju = nn
end if
#if NDIMS == 3
if (dir == 3) then
kl = nlims(side,1)
ku = nlims(side,2)
else
kl = 1
ku = nn
end if
#else /* NDIMS == 3 */
kl = 1
ku = 1
#endif /* NDIMS == 3 */
do k = kl, ku
do j = jl, ju
call prim2cons(pdata%q(:,il:iu,j,k), pdata%u(:,il:iu,j,k), .true.)
end do
end do
!-------------------------------------------------------------------------------
!
end subroutine block_conservative_variables
#ifdef MPI #ifdef MPI
! !
!=============================================================================== !===============================================================================