From 6f23ae962ef132f1cd361fb0669debe6ac6c9276 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Feb 2022 19:25:18 -0300 Subject: [PATCH] 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 --- sources/boundaries.F90 | 213 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 187 insertions(+), 26 deletions(-) diff --git a/sources/boundaries.F90 b/sources/boundaries.F90 index e43ed03..5b6d668 100644 --- a/sources/boundaries.F90 +++ b/sources/boundaries.F90 @@ -3611,6 +3611,8 @@ module boundaries integer :: m #endif /* NDIMS == 2 */ + integer, dimension(NDIMS) :: s + real(kind=8), dimension(nn) :: x, y #if NDIMS == 3 real(kind=8), dimension(nn) :: z @@ -3652,13 +3654,15 @@ module boundaries m = 3 - n do j = 1, nsides + s(2) = j do i = 1, nsides - - if (.not. associated(pmeta%edges(i,j,m)%ptr)) & - call block_boundary_specific(n, [ i, j, k ], & - t, dt, x, y, z, & - pmeta%data%u, status) - + s(1) = i + if (.not. associated(pmeta%edges(i,j,m)%ptr)) then + call block_primitive_variables(n, pmeta%data, status) + call block_boundary_specific(n, [ i, j, k ], t, dt, & + x, y, z, pmeta%data%q, status) + call block_conservative_variables(n, s(n), pmeta%data) + end if end do end do @@ -3672,14 +3676,17 @@ module boundaries if (.not. periodic(n)) then do k = 1, nsides + s(3) = k do j = 1, nsides + s(2) = j do i = 1, nsides - - if (.not. associated(pmeta%faces(i,j,k,n)%ptr)) & - call block_boundary_specific(n, [ i, j, k ], & - t, dt, x, y, z, & - pmeta%data%u, status) - + s(1) = i + if (.not. associated(pmeta%faces(i,j,k,n)%ptr)) then + call block_primitive_variables(n, pmeta%data, status) + call block_boundary_specific(n, [ i, j, k ], t, dt, & + x, y, z, pmeta%data%q, status) + call block_conservative_variables(n, s(n), pmeta%data) + end if 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 : nb, ne, nbl, neu 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 - use equations , only : imz, ibz + use equations , only : ivz, ibz #endif /* NDIMS == 3 */ use equations , only : csnd2 use gravity , only : gravitational_acceleration @@ -4582,12 +4589,12 @@ module boundaries if (side(1) == 1) then do i = nbl, 1, -1 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 else do i = neu, nn 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 if @@ -4599,7 +4606,7 @@ module boundaries is = nbl + i 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 qn(ibx,it,jl:ju,kl:ku) = - qn(ibx,is,jl:ju,kl:ku) end if @@ -4610,7 +4617,7 @@ module boundaries is = neu - i 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 qn(ibx,it,jl:ju,kl:ku) = - qn(ibx,is,jl:ju,kl:ku) end if @@ -4741,12 +4748,12 @@ module boundaries if (side(2) == 1) then do j = nbl, 1, -1 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 else do j = neu, nn 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 if @@ -4758,7 +4765,7 @@ module boundaries js = nbl + j 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 qn(iby,il:iu,jt,kl:ku) = - qn(iby,il:iu,js,kl:ku) end if @@ -4769,7 +4776,7 @@ module boundaries js = neu - j 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 qn(iby,il:iu,jt,kl:ku) = - qn(iby,il:iu,js,kl:ku) end if @@ -4896,12 +4903,12 @@ module boundaries if (side(3) == 1) then do k = nbl, 1, -1 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 else do k = neu, nn 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 if @@ -4913,7 +4920,7 @@ module boundaries ks = nbl + k 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 qn(ibz,il:iu,jl:ju,kt) = - qn(ibz,il:iu,jl:ju,ks) end if @@ -4924,7 +4931,7 @@ module boundaries ks = neu - k 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 qn(ibz,il:iu,jl:ju,kt) = - qn(ibz,il:iu,jl:ju,ks) end if @@ -5038,6 +5045,160 @@ module boundaries ! 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 ! !===============================================================================