From 2e42cc7231255d512a944c297f30030cd0072a31 Mon Sep 17 00:00:00 2001
From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com>
Date: Sat, 8 Mar 2025 00:14:40 -0500
Subject: [PATCH] boundary conditions refactor and MPI fix for elliptic
 smoothing

---
 .../m_boundary_conditions.fpp                 |  93 ++++----
 src/common/m_mpi_common.fpp                   |  13 +-
 src/pre_process/m_perturbation.fpp            | 209 +-----------------
 src/simulation/m_cbc.fpp                      |   4 +-
 src/simulation/m_global_parameters.fpp        |  15 +-
 src/simulation/m_ibm.fpp                      |   4 +-
 src/simulation/m_qbmm.fpp                     |  10 +-
 src/simulation/m_rhs.fpp                      |  18 +-
 src/simulation/m_riemann_solvers.fpp          |   8 +-
 src/simulation/m_surface_tension.fpp          |   4 +-
 src/simulation/m_viscous.fpp                  |   6 +-
 src/simulation/m_weno.fpp                     |   8 +-
 12 files changed, 89 insertions(+), 303 deletions(-)
 rename src/{simulation => common}/m_boundary_conditions.fpp (97%)

diff --git a/src/simulation/m_boundary_conditions.fpp b/src/common/m_boundary_conditions.fpp
similarity index 97%
rename from src/simulation/m_boundary_conditions.fpp
rename to src/common/m_boundary_conditions.fpp
index 073c62d07..252f68875 100644
--- a/src/simulation/m_boundary_conditions.fpp
+++ b/src/common/m_boundary_conditions.fpp
@@ -16,9 +16,11 @@ module m_boundary_conditions
 
     implicit none
 
-    private; 
-    public :: s_populate_variables_buffers, &
-              s_populate_capillary_buffers
+#ifdef MFC_SIMULATION
+    private; public :: s_populate_variables_buffers, s_populate_capillary_buffers
+#else
+    private; public :: s_populate_variables_buffers
+#endif
 
 contains
 
@@ -28,7 +30,8 @@ contains
     subroutine s_populate_variables_buffers(q_prim_vf, pb, mv)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+
+        real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
 
         integer :: bc_loc, bc_dir
 
@@ -66,6 +69,7 @@ contains
                 q_prim_vf, pb, mv, 1, 1)
         end select
 
+#ifdef MFC_SIMULATION
         if (qbmm .and. .not. polytropic) then
             select case (bc_x%beg)
             case (-13:-3) ! Ghost-cell extrap. BC at beginning
@@ -85,8 +89,7 @@ contains
                 call s_qbmm_extrapolation(pb, mv, 1, 1)
             end select
         end if
-
-        ! END: Population of Buffers in x-direction
+#endif
 
         ! Population of Buffers in y-direction
 
@@ -126,8 +129,8 @@ contains
                 q_prim_vf, pb, mv, 2, 1)
         end select
 
+#ifdef MFC_SIMULATION
         if (qbmm .and. .not. polytropic) then
-
             select case (bc_y%beg)
             case (-13:-3) ! Ghost-cell extrap. BC at beginning
                 call s_qbmm_extrapolation(pb, mv, 2, -1)
@@ -145,10 +148,8 @@ contains
             case (-16)    ! No-slip wall BC at end
                 call s_qbmm_extrapolation(pb, mv, 2, 1)
             end select
-
         end if
-
-        ! END: Population of Buffers in y-direction
+#endif
 
         ! Population of Buffers in z-direction
 
@@ -186,8 +187,8 @@ contains
                 q_prim_vf, pb, mv, 3, 1)
         end select
 
+#ifdef MFC_SIMULATION
         if (qbmm .and. .not. polytropic) then
-
             select case (bc_z%beg)
             case (-13:-3) ! Ghost-cell extrap. BC at beginning
                 call s_qbmm_extrapolation(pb, mv, 3, -1)
@@ -205,9 +206,8 @@ contains
             case (-16)    ! No-slip wall BC at end
                 call s_qbmm_extrapolation(pb, mv, 3, 1)
             end select
-
         end if
-
+#endif
         ! END: Population of Buffers in z-direction
 
     end subroutine s_populate_variables_buffers
@@ -215,7 +215,7 @@ contains
     subroutine s_ghost_cell_extrapolation(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
         integer, intent(in) :: bc_dir, bc_loc
         integer :: j, k, l, q, i
 
@@ -325,7 +325,7 @@ contains
     subroutine s_symmetry(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
         integer, intent(in) :: bc_dir, bc_loc
 
         integer :: j, k, l, q, i
@@ -362,7 +362,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -380,6 +380,7 @@ contains
                         end do
                     end do
                 end if
+#endif
 
             else !< bc_x%end
 
@@ -411,7 +412,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -429,7 +430,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             end if
 
             !< y-direction
@@ -463,7 +464,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -481,7 +482,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             else !< bc_y%end
 
                 !$acc parallel loop collapse(3) gang vector default(present)
@@ -510,7 +511,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -528,7 +529,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             end if
 
             !< z-direction
@@ -562,7 +563,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -580,7 +581,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             else !< bc_z%end
 
                 !$acc parallel loop collapse(3) gang vector default(present)
@@ -609,7 +610,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -627,7 +628,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             end if
 
         end if
@@ -637,7 +638,7 @@ contains
     subroutine s_periodic(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
         integer, intent(in) :: bc_dir, bc_loc
 
         integer :: j, k, l, q, i
@@ -658,7 +659,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -676,7 +677,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             else !< bc_x%end
 
                 !$acc parallel loop collapse(4) gang vector default(present)
@@ -690,7 +691,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -708,7 +709,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             end if
 
             !< y-direction
@@ -727,7 +728,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(4) gang vector default(present)
                     do i = 1, nb
@@ -745,7 +746,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             else !< bc_y%end
 
                 !$acc parallel loop collapse(4) gang vector default(present)
@@ -759,7 +760,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -777,7 +778,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             end if
 
             !< z-direction
@@ -796,7 +797,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -814,7 +815,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             else !< bc_z%end
 
                 !$acc parallel loop collapse(4) gang vector default(present)
@@ -828,7 +829,7 @@ contains
                         end do
                     end do
                 end do
-
+#ifdef MFC_SIMULATION
                 if (qbmm .and. .not. polytropic) then
                     !$acc parallel loop collapse(5) gang vector default(present)
                     do i = 1, nb
@@ -846,7 +847,7 @@ contains
                         end do
                     end do
                 end if
-
+#endif
             end if
 
         end if
@@ -856,7 +857,7 @@ contains
     subroutine s_axis(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
         integer, intent(in) :: bc_dir, bc_loc
 
         integer :: j, k, l, q, i
@@ -905,7 +906,7 @@ contains
                 end do
             end do
         end do
-
+#ifdef MFC_SIMULATION
         if (qbmm .and. .not. polytropic) then
             !$acc parallel loop collapse(5) gang vector default(present)
             do i = 1, nb
@@ -923,13 +924,13 @@ contains
                 end do
             end do
         end if
-
+#endif
     end subroutine s_axis
 
     subroutine s_slip_wall(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
         integer, intent(in) :: bc_dir, bc_loc
 
         integer :: j, k, l, q, i
@@ -1070,7 +1071,7 @@ contains
     subroutine s_no_slip_wall(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
         integer, intent(in) :: bc_dir, bc_loc
 
         integer :: j, k, l, q, i
@@ -1246,7 +1247,7 @@ contains
 
     subroutine s_qbmm_extrapolation(pb, mv, bc_dir, bc_loc)
 
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
         integer, intent(in) :: bc_dir, bc_loc
 
         integer :: j, k, l, q, i
@@ -1378,6 +1379,7 @@ contains
 
     end subroutine s_qbmm_extrapolation
 
+#ifdef MFC_SIMULATION
     subroutine s_populate_capillary_buffers(c_divs)
 
         type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs
@@ -1659,5 +1661,6 @@ contains
         end if
 
     end subroutine s_populate_capillary_buffers
+#endif
 
 end module m_boundary_conditions
diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp
index 29d55ba8c..4eb4373b7 100644
--- a/src/common/m_mpi_common.fpp
+++ b/src/common/m_mpi_common.fpp
@@ -613,7 +613,6 @@ contains
 
     end subroutine s_mpi_finalize
 
-
     !>  The goal of this procedure is to populate the buffers of
         !!      the cell-average conservative variables by communicating
         !!      with the neighboring processors.
@@ -621,16 +620,13 @@ contains
         !!  @param mpi_dir MPI communication coordinate direction
         !!  @param pbc_loc Processor boundary condition (PBC) location
     subroutine s_mpi_sendrecv_variables_buffers(q_cons_vf, &
-#ifdef MFC_SIMULATION
                                                 pb, mv, &
-#endif
                                                 mpi_dir, &
                                                 pbc_loc)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
-#ifdef MFC_SIMULATION
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
-#endif
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
+
         integer, intent(in) :: mpi_dir, pbc_loc
 
         integer :: i, j, k, l, r, q !< Generic loop iterators
@@ -850,12 +846,13 @@ contains
         #:endfor
         call nvtxEndRange ! Packbuf
 
+        p_send => q_prims_buff_send(0)
+        p_recv => q_prims_buff_recv(0)
+
 #ifdef MFC_SIMULATION
         ! Send/Recv
         #:for rdma_mpi in [False, True]
             if (rdma_mpi .eqv. ${'.true.' if rdma_mpi else '.false.'}$) then
-                p_send => q_prims_buff_send(0)
-                p_recv => q_prims_buff_recv(0)
                 #:if rdma_mpi
                     !$acc data attach(p_send, p_recv)
                     !$acc host_data use_device(p_send, p_recv)
diff --git a/src/pre_process/m_perturbation.fpp b/src/pre_process/m_perturbation.fpp
index ff3926add..1dd71295c 100644
--- a/src/pre_process/m_perturbation.fpp
+++ b/src/pre_process/m_perturbation.fpp
@@ -13,6 +13,8 @@ module m_perturbation
     use m_mpi_proxy              !< Message passing interface (MPI) module proxy
 
     use m_eigen_solver          ! Subroutines to solve eigenvalue problem for
+
+    use m_boundary_conditions   ! Boundary conditions module
     ! complex general matrix
 
     use ieee_arithmetic
@@ -623,211 +625,8 @@ contains
 
         do q = 1, elliptic_smoothing_iters
 
-            if (bcxb >= 0) then
-                call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, -1)
-            else if (bcxb == -1) then
-                do l = 0, p
-                    do k = 0, n
-                        do j = 1, buff_size
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(-j, k, l) = q_prim_vf(i)%sf(m - j + 1, k, l)
-                            end do
-                        end do
-                    end do
-                end do
-            else if (bcxb == -2) then
-                do l = 0, p
-                    do k = 0, n
-                        do j = 1, buff_size
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(-j, k, l) = q_prim_vf(i)%sf(j - 1, k, l)
-                            end do
-                        end do
-                    end do
-                end do
-            else
-                do l = 0, p
-                    do k = 0, n
-                        do j = 1, buff_size
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(-j, k, l) = q_prim_vf(i)%sf(0, k, l)
-                            end do
-                        end do
-                    end do
-                end do
-            end if
-
-            if (bcxe >= 0) then
-                call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, 1)
-            else if (bcxe == -1) then
-                do l = 0, p
-                    do k = 0, n
-                        do j = 1, buff_size
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(m + j, k, l) = q_prim_vf(i)%sf(j - 1, k, l)
-                            end do
-                        end do
-                    end do
-                end do
-            else if (bcxe == -2) then
-                do l = 0, p
-                    do k = 0, n
-                        do j = 1, buff_size
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(m + j, k, l) = q_prim_vf(i)%sf(m - (j - 1), k, l)
-                            end do
-                        end do
-                    end do
-                end do
-            else
-                do l = 0, p
-                    do k = 0, n
-                        do j = 1, buff_size
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(m + j, k, l) = q_prim_vf(i)%sf(m, k, l)
-                            end do
-                        end do
-                    end do
-                end do
-            end if
-
-            if (bcyb >= 0) then
-                call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, -1)
-            else if (bcyb == -1) then
-                do l = 0, p
-                    do k = 1, buff_size
-                        do j = 0, m
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(j, -k, l) = q_prim_vf(i)%sf(j, n - k + 1, l)
-                            end do
-                        end do
-                    end do
-                end do
-            else if (bcyb == -2) then
-                do l = 0, p
-                    do k = 1, buff_size
-                        do j = 0, m
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(j, -k, l) = q_prim_vf(i)%sf(j, k - 1, l)
-                            end do
-                        end do
-                    end do
-                end do
-            else
-                do l = 0, p
-                    do k = 1, buff_size
-                        do j = 0, m
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(j, -k, l) = q_prim_vf(i)%sf(j, 0, l)
-                            end do
-                        end do
-                    end do
-                end do
-            end if
-
-            if (bcye >= 0) then
-                call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, 1)
-            else if (bcye == -1) then
-                do l = 0, p
-                    do k = 1, buff_size
-                        do j = 0, m
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(j, n + k, l) = q_prim_vf(i)%sf(j, k - 1, l)
-                            end do
-                        end do
-                    end do
-                end do
-            else if (bcye == -2) then
-                do l = 0, p
-                    do k = 1, buff_size
-                        do j = 0, m
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(j, n + k, l) = q_prim_vf(i)%sf(j, n - (k - 1), l)
-                            end do
-                        end do
-                    end do
-                end do
-            else
-                do l = 0, p
-                    do k = 1, buff_size
-                        do j = 0, m
-                            do i = 1, sys_size
-                                q_prim_vf(i)%sf(j, n + k, l) = q_prim_vf(i)%sf(j, n, l)
-                            end do
-                        end do
-                    end do
-                end do
-            end if
-
-            if (p > 0) then
-                if (bczb >= 0) then
-                    call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, -1)
-                else if (bczb == -1) then
-                    do l = 1, buff_size
-                        do k = 0, n
-                            do j = 0, m
-                                do i = 1, sys_size
-                                    q_prim_vf(i)%sf(j, k, -l) = q_prim_vf(i)%sf(j, k, p - l + 1)
-                                end do
-                            end do
-                        end do
-                    end do
-                else if (bczb == -2) then
-                    do l = 1, buff_size
-                        do k = 0, n
-                            do j = 0, m
-                                do i = 1, sys_size
-                                    q_prim_vf(i)%sf(j, k, -l) = q_prim_vf(i)%sf(j, k, l - 1)
-                                end do
-                            end do
-                        end do
-                    end do
-                else
-                    do l = 1, buff_size
-                        do k = 0, n
-                            do j = 0, m
-                                do i = 1, sys_size
-                                    q_prim_vf(i)%sf(j, k, -l) = q_prim_vf(i)%sf(j, k, 0)
-                                end do
-                            end do
-                        end do
-                    end do
-                end if
-
-                if (bcze >= 0) then
-                    call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, 1)
-                else if (bcze == -1) then
-                    do l = 1, buff_size
-                        do k = 0, n
-                            do j = 0, m
-                                do i = 1, sys_size
-                                    q_prim_vf(i)%sf(j, k, p + l) = q_prim_vf(i)%sf(j, k, l - 1)
-                                end do
-                            end do
-                        end do
-                    end do
-                else if (bcze == -2) then
-                    do l = 1, buff_size
-                        do k = 0, n
-                            do j = 0, m
-                                do i = 1, sys_size
-                                    q_prim_vf(i)%sf(j, k, p + l) = q_prim_vf(i)%sf(j, k, p - (l - 1))
-                                end do
-                            end do
-                        end do
-                    end do
-                else
-                    do l = 1, buff_size
-                        do k = 0, n
-                            do j = 0, m
-                                do i = 1, sys_size
-                                    q_prim_vf(i)%sf(j, k, p + l) = q_prim_vf(i)%sf(j, k, p)
-                                end do
-                            end do
-                        end do
-                    end do
-                end if
-            end if
+            ! Communication of buffer regions and apply boundary conditions
+            call s_populate_variables_buffers(q_prim_vf)
 
             ! Perform smoothing and store in temp array
             if (p == 0) then
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
index 8ab1a97eb..1632899cd 100644
--- a/src/simulation/m_cbc.fpp
+++ b/src/simulation/m_cbc.fpp
@@ -681,7 +681,7 @@ contains
                     call s_convert_primitive_to_flux_variables(q_prim_rs${XYZ}$_vf, &
                                                                F_rs${XYZ}$_vf, &
                                                                F_src_rs${XYZ}$_vf, &
-                                                               is1, is2, is3, starty, startz)
+                                                               is1, is2, is3, idwbuff(2)%beg, idwbuff(3)%beg)
 
                     !$acc parallel loop collapse(3) gang vector default(present)
                     do i = 1, advxe
@@ -712,7 +712,7 @@ contains
                     call s_convert_primitive_to_flux_variables(q_prim_rs${XYZ}$_vf, &
                                                                F_rs${XYZ}$_vf, &
                                                                F_src_rs${XYZ}$_vf, &
-                                                               is1, is2, is3, starty, startz)
+                                                               is1, is2, is3, idwbuff(2)%beg, idwbuff(3)%beg)
 
                     !$acc parallel loop collapse(4) gang vector default(present)
                     do i = 1, advxe
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
index 6963d6dfe..fd44d0f4a 100644
--- a/src/simulation/m_global_parameters.fpp
+++ b/src/simulation/m_global_parameters.fpp
@@ -287,9 +287,7 @@ module m_global_parameters
     !! conditions data to march the solution in the physical computational domain
     !! to the next time-step.
 
-    integer :: startx, starty, startz
-
-    !$acc declare create(sys_size, buff_size, startx, starty, startz, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, b_size, tensor_size, xi_idx, species_idx)
+    !$acc declare create(sys_size, buff_size, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, b_size, tensor_size, xi_idx, species_idx)
 
     ! END: Simulation Algorithm Parameters
 
@@ -1140,18 +1138,7 @@ contains
                 & idwbuff(3)%beg:idwbuff(3)%end))
         end if
 
-        startx = -buff_size
-        starty = 0
-        startz = 0
-        if (n > 0) then
-            starty = -buff_size
-        end if
-        if (p > 0) then
-            startz = -buff_size
-        end if
-
         !$acc update device(fd_order,fd_number)
-        !$acc update device(startx, starty, startz)
 
         if (cyl_coord .neqv. .true.) then ! Cartesian grid
             grid_geometry = 1
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
index 21f6bb27f..b6967b913 100644
--- a/src/simulation/m_ibm.fpp
+++ b/src/simulation/m_ibm.fpp
@@ -123,7 +123,7 @@ contains
             dimension(sys_size), &
             intent(INOUT) :: q_prim_vf !< Primitive Variables
 
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), optional, intent(INOUT) :: pb, mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), optional, intent(INOUT) :: pb, mv
 
         integer :: i, j, k, l, q, r!< Iterator variables
         integer :: patch_id !< Patch ID of ghost point
@@ -730,7 +730,7 @@ contains
         type(scalar_field), &
             dimension(sys_size), &
             intent(IN) :: q_prim_vf !< Primitive Variables
-        real(wp), optional, dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
+        real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(INOUT) :: pb, mv
 
         type(ghost_point), intent(IN) :: gp
         real(wp), intent(INOUT) :: pres_IP
diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp
index 11298e4f7..c93a4be2f 100644
--- a/src/simulation/m_qbmm.fpp
+++ b/src/simulation/m_qbmm.fpp
@@ -417,8 +417,8 @@ contains
         type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf, q_prim_vf
         type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
         type(scalar_field), dimension(sys_size), intent(in) :: flux_n_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, rhs_pb
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: mv, rhs_mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, rhs_pb
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: mv, rhs_mv
 
         integer :: i, j, k, l, q
 
@@ -820,10 +820,10 @@ contains
         type(scalar_field), dimension(:), intent(inout) :: q_cons_vf, q_prim_vf
         type(scalar_field), dimension(:), intent(inout) :: momsp
         type(scalar_field), dimension(0:, 0:, :), intent(inout) :: moms3d
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, rhs_pb
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: mv, rhs_mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, rhs_pb
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: mv, rhs_mv
         type(int_bounds_info), intent(in) :: ix, iy, iz
-        real(wp), dimension(startx:, starty:, startz:), intent(inout) :: nbub_sc
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:), intent(inout) :: nbub_sc
 
         real(wp), dimension(nmom) :: moms, msum
         real(wp), dimension(nnode, nb) :: wght, abscX, abscY, wght_pb, wght_mv, wght_ht, ht
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
index d9617f930..56ee36a99 100644
--- a/src/simulation/m_rhs.fpp
+++ b/src/simulation/m_rhs.fpp
@@ -585,9 +585,9 @@ contains
         !$acc parallel loop collapse(4) gang vector default(present)
         do id = 1, num_dims
             do i = 1, sys_size
-                do l = startz, p - startz
-                    do k = starty, n - starty
-                        do j = startx, m - startx
+                do l = idwbuff(3)%beg, idwbuff(3)%end
+                    do k = idwbuff(2)%beg, idwbuff(2)%end
+                        do j = idwbuff(1)%beg, idwbuff(1)%end
                             flux_gsrc_n(id)%vf(i)%sf(j, k, l) = 0._wp
                         end do
                     end do
@@ -607,8 +607,8 @@ contains
         type(scalar_field), intent(inout) :: q_T_sf
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
         type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, rhs_pb
-        real(wp), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: mv, rhs_mv
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, rhs_pb
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: mv, rhs_mv
         integer, intent(in) :: t_step
         real(wp), intent(inout) :: time_avg
 
@@ -1971,8 +1971,8 @@ contains
                                                   norm_dir)
 
         type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_x, vR_y, vR_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: vL_x, vL_y, vL_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: vR_x, vR_y, vR_z
         integer, intent(in) :: norm_dir
 
         integer :: weno_dir !< Coordinate direction of the WENO reconstruction
@@ -2023,8 +2023,8 @@ contains
                                                               norm_dir)
 
         type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_x, vR_y, vR_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: vL_x, vL_y, vL_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: vR_x, vR_y, vR_z
         integer, intent(in) :: norm_dir
 
         integer :: recon_dir !< Coordinate direction of the WENO reconstruction
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
index 0122ba9ff..5314012a4 100644
--- a/src/simulation/m_riemann_solvers.fpp
+++ b/src/simulation/m_riemann_solvers.fpp
@@ -156,7 +156,7 @@ contains
                                 flux_gsrc_vf, &
                                 norm_dir, ix, iy, iz)
 
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(INOUT) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
         type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
 
         type(scalar_field), allocatable, dimension(:), intent(INOUT) :: qL_prim_vf, qR_prim_vf
@@ -277,7 +277,7 @@ contains
                                     flux_gsrc_vf, &
                                     norm_dir, ix, iy, iz)
 
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
 
         type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
@@ -967,7 +967,7 @@ contains
                                      flux_gsrc_vf, &
                                      norm_dir, ix, iy, iz)
 
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
         type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
 
@@ -2903,7 +2903,7 @@ contains
         qR_prim_vf, &
         norm_dir, ix, iy, iz)
 
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
 
         type(scalar_field), &
             allocatable, dimension(:), &
diff --git a/src/simulation/m_surface_tension.fpp b/src/simulation/m_surface_tension.fpp
index 06b230c60..86777e174 100644
--- a/src/simulation/m_surface_tension.fpp
+++ b/src/simulation/m_surface_tension.fpp
@@ -303,8 +303,8 @@ contains
 
         type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
 
-        real(wp), dimension(startx:, starty:, startz:, iv%beg:), intent(out) :: vL_x, vL_y, vL_z
-        real(wp), dimension(startx:, starty:, startz:, iv%beg:), intent(out) :: vR_x, vR_y, vR_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, iv%beg:), intent(out) :: vL_x, vL_y, vL_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, iv%beg:), intent(out) :: vR_x, vR_y, vR_z
         integer, intent(in) :: norm_dir
 
         integer :: recon_dir !< Coordinate direction of the WENO reconstruction
diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp
index 44d6bbf62..0abd1e978 100644
--- a/src/simulation/m_viscous.fpp
+++ b/src/simulation/m_viscous.fpp
@@ -522,7 +522,7 @@ contains
                              dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp, &
                              ix, iy, iz)
 
-        real(wp), dimension(startx:, starty:, startz:, 1:), &
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), &
             intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf, &
                              qL_prim_rsy_vf, qR_prim_rsy_vf, &
                              qL_prim_rsz_vf, qR_prim_rsz_vf
@@ -968,7 +968,7 @@ contains
         type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
         type(scalar_field), dimension(iv%beg:iv%end), intent(inout) :: vL_prim_vf, vR_prim_vf
 
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
         integer, intent(in) :: norm_dir
         type(int_bounds_info), intent(in) :: ix, iy, iz
 
@@ -1064,7 +1064,7 @@ contains
                                                              norm_dir, vL_prim_vf, vR_prim_vf, ix, iy, iz)
 
         type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
-        real(wp), dimension(startx:, starty:, startz:, iv%beg:), intent(inout) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, iv%beg:), intent(inout) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
         type(scalar_field), dimension(iv%beg:iv%end), intent(inout) :: vL_prim_vf, vR_prim_vf
         type(int_bounds_info), intent(in) :: ix, iy, iz
 
diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp
index 0036a5eee..3c787afd8 100644
--- a/src/simulation/m_weno.fpp
+++ b/src/simulation/m_weno.fpp
@@ -643,8 +643,8 @@ contains
                       is1_weno_d, is2_weno_d, is3_weno_d)
 
         type(scalar_field), dimension(1:), intent(in) :: v_vf
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z
         integer, intent(in) :: norm_dir
         integer, intent(in) :: weno_dir
         type(int_bounds_info), intent(in) :: is1_weno_d, is2_weno_d, is3_weno_d
@@ -1235,8 +1235,8 @@ contains
         !!  @param l Thire-coordinate cell index
     subroutine s_preserve_monotonicity(v_rs_ws, vL_rs_vf, vR_rs_vf)
 
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(IN) :: v_rs_ws
-        real(wp), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: vL_rs_vf, vR_rs_vf
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(IN) :: v_rs_ws
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(INOUT) :: vL_rs_vf, vR_rs_vf
 
         integer :: i, j, k, l