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 ! !===============================================================================