From 687659caa9fd06f4bff1be9b92e9b3cde665b70c Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Tue, 1 Apr 2025 09:30:04 -0400
Subject: [PATCH 01/23] Removing stop commands

---
 src/simulation/m_bubbles_EL.fpp         | 33 ++++++++++-------------
 src/simulation/m_bubbles_EL_kernels.fpp | 36 +++++--------------------
 src/simulation/m_checker.fpp            |  1 +
 src/simulation/m_start_up.fpp           | 10 ++++++-
 4 files changed, 30 insertions(+), 50 deletions(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index b48f941a39..ba0b1b164a 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -251,7 +251,7 @@ contains
                 end do
                 close (94)
             else
-                stop "if you include lagrange bubbles, you have to initialize them in input/lag_bubbles.dat"
+                stop "Lagrange bubbles: you have to initialize them in input/lag_bubbles.dat"
             end if
         else
             if (proc_rank == 0) print *, 'Restarting lagrange bubbles at save_count: ', save_count
@@ -338,6 +338,16 @@ contains
         cell = -buff_size
         call s_locate_cell(mtn_pos(bub_id, 1:3, 1), cell, mtn_s(bub_id, 1:3, 1))
 
+        ! Check if the bubble is located in the ghost cell of a symmetric boundary
+        if (bc_x%beg == -2 .and. cell(1) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_x%beg)."
+        if (bc_x%end == -2 .and. cell(1) > m) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_x%end)."
+        if (bc_y%beg == -2 .and. cell(2) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_y%beg)."
+        if (bc_y%end == -2 .and. cell(2) > n) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_y%end)."
+        if (p > 0) then
+            if (bc_z%beg == -2 .and. cell(3) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_z%beg)."
+            if (bc_z%end == -2 .and. cell(3) > p) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_z%end)."
+        end if
+
         ! If particle is in the ghost cells, find the closest non-ghost cell
         cell(1) = min(max(cell(1), 0), m)
         cell(2) = min(max(cell(2), 0), n)
@@ -366,7 +376,7 @@ contains
         gas_mv(bub_id, 1) = pv*volparticle*(1._wp/(R_v*Tw))*(massflag) ! vapermass
         gas_mg(bub_id) = (gas_p(bub_id, 1) - pv*(massflag))*volparticle*(1._wp/(R_n*Tw)) ! gasmass
         if (gas_mg(bub_id) <= 0._wp) then
-            stop 'the initial mass of gas inside the bubble is negative. Check your initial conditions'
+            stop "The initial mass of gas inside the bubble is negative. Check your initial conditions."
         end if
         totalmass = gas_mg(bub_id) + gas_mv(bub_id, 1) ! totalmass
 
@@ -374,9 +384,7 @@ contains
         concvap = gas_mv(bub_id, 1)/(gas_mv(bub_id, 1) + gas_mg(bub_id))
         omegaN = (3._wp*(gas_p(bub_id, 1) - pv*(massflag)) + 4._wp*(1._wp/Web)/bub_R0(bub_id))/rhol
         if (pv*(massflag) > gas_p(bub_id, 1)) then
-            print *, 'Not allowed: bubble initially located in a region with pressure below the vapor pressure'
-            print *, 'location:', mtn_pos(bub_id, 1:3, 1)
-            stop
+            stop "Lagrange bubble initially located in a region with pressure below the vapor pressure."
         end if
         omegaN = sqrt(omegaN/bub_R0(bub_id)**2._wp)
 
@@ -906,11 +914,7 @@ contains
             ! Bubble dynamic closure from Maeda and Colonius (2018)
 
             ! Range of cells included in Omega
-            if (lag_params%smooth_type == 1) then
-                mapCells_pinf = mapCells
-            else
-                stop "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
-            end if
+            mapCells_pinf = mapCells
 
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
             smearGrid = mapCells_pinf - (-mapCells_pinf) + 1
@@ -987,9 +991,6 @@ contains
             f_pinfl = charpres2/charvol2
             vol = charvol
             dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
-        else
-
-            stop "Check cluterflag. Exiting."
 
         end if
 
@@ -1038,7 +1039,6 @@ contains
                 mtn_vel(k, 1:3, 1) = mtn_vel(k, 1:3, 1) + dt*mtn_dveldt(k, 1:3, 1)
                 gas_p(k, 1) = gas_p(k, 1) + dt*gas_dpdt(k, 1)
                 gas_mv(k, 1) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1)
-                if (intfc_rad(k, 1) <= 0._wp) stop "Negative bubble radius encountered, please reduce dt"
             end do
 
             call s_transfer_data_to_tmp()
@@ -1061,7 +1061,6 @@ contains
                     mtn_vel(k, 1:3, 2) = mtn_vel(k, 1:3, 1) + dt*mtn_dveldt(k, 1:3, 1)
                     gas_p(k, 2) = gas_p(k, 1) + dt*gas_dpdt(k, 1)
                     gas_mv(k, 2) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1)
-                    if (intfc_rad(k, 2) <= 0._wp) stop "Negative bubble radius encountered, please reduce dt"
                 end do
 
             elseif (stage == 2) then
@@ -1074,7 +1073,6 @@ contains
                     mtn_vel(k, 1:3, 1) = mtn_vel(k, 1:3, 1) + dt*(mtn_dveldt(k, 1:3, 1) + mtn_dveldt(k, 1:3, 2))/2._wp
                     gas_p(k, 1) = gas_p(k, 1) + dt*(gas_dpdt(k, 1) + gas_dpdt(k, 2))/2._wp
                     gas_mv(k, 1) = gas_mv(k, 1) + dt*(gas_dmvdt(k, 1) + gas_dmvdt(k, 2))/2._wp
-                    if (intfc_rad(k, 1) <= 0._wp) stop "Negative bubble radius encountered, please reduce dt"
                 end do
 
                 call s_transfer_data_to_tmp()
@@ -1099,7 +1097,6 @@ contains
                     mtn_vel(k, 1:3, 2) = mtn_vel(k, 1:3, 1) + dt*mtn_dveldt(k, 1:3, 1)
                     gas_p(k, 2) = gas_p(k, 1) + dt*gas_dpdt(k, 1)
                     gas_mv(k, 2) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1)
-                    if (intfc_rad(k, 2) <= 0._wp) stop "Negative bubble radius encountered, please reduce dt"
                 end do
 
             elseif (stage == 2) then
@@ -1112,7 +1109,6 @@ contains
                     mtn_vel(k, 1:3, 2) = mtn_vel(k, 1:3, 1) + dt*(mtn_dveldt(k, 1:3, 1) + mtn_dveldt(k, 1:3, 2))/4._wp
                     gas_p(k, 2) = gas_p(k, 1) + dt*(gas_dpdt(k, 1) + gas_dpdt(k, 2))/4._wp
                     gas_mv(k, 2) = gas_mv(k, 1) + dt*(gas_dmvdt(k, 1) + gas_dmvdt(k, 2))/4._wp
-                    if (intfc_rad(k, 2) <= 0._wp) stop "Negative bubble radius encountered, please reduce dt"
                 end do
             elseif (stage == 3) then
                 !$acc parallel loop gang vector default(present) private(k)
@@ -1124,7 +1120,6 @@ contains
                     mtn_vel(k, 1:3, 1) = mtn_vel(k, 1:3, 1) + (2._wp/3._wp)*dt*(mtn_dveldt(k, 1:3, 1)/4._wp + mtn_dveldt(k, 1:3, 2)/4._wp + mtn_dveldt(k, 1:3, 3))
                     gas_p(k, 1) = gas_p(k, 1) + (2._wp/3._wp)*dt*(gas_dpdt(k, 1)/4._wp + gas_dpdt(k, 2)/4._wp + gas_dpdt(k, 3))
                     gas_mv(k, 1) = gas_mv(k, 1) + (2._wp/3._wp)*dt*(gas_dmvdt(k, 1)/4._wp + gas_dmvdt(k, 2)/4._wp + gas_dmvdt(k, 3))
-                    if (intfc_rad(k, 1) <= 0._wp) stop "Negative bubble radius encountered, please reduce dt"
                 end do
 
                 call s_transfer_data_to_tmp()
diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index 43246aad55..b45832cae9 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -317,51 +317,27 @@ contains
 
         ! x-dir
         if (bc_x%beg == -2 .and. (cell(1) <= mapCells - 1)) then
-            if (cell(1) >= 0) then
-                cellaux(1) = abs(cellaux(1)) - 1
-            else
-                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%beg)."
-            end if
+            cellaux(1) = abs(cellaux(1)) - 1
         end if
         if (bc_x%end == -2 .and. (cell(1) >= m + 1 - mapCells)) then
-            if (cell(1) <= m) then
-                cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
-            else
-                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%end)."
-            end if
+            cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
         end if
 
         !y-dir
         if (bc_y%beg == -2 .and. (cell(2) <= mapCells - 1)) then
-            if (cell(2) >= 0) then
-                cellaux(2) = abs(cellaux(2)) - 1
-            else
-                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%beg)."
-            end if
+            cellaux(2) = abs(cellaux(2)) - 1
         end if
         if (bc_y%end == -2 .and. (cell(2) >= n + 1 - mapCells)) then
-            if (cell(2) <= n) then
-                cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
-            else
-                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%end)."
-            end if
+            cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
         end if
 
         if (p > 0) then
             !z-dir
             if (bc_z%beg == -2 .and. (cell(3) <= mapCells - 1)) then
-                if (cell(3) >= 0) then
-                    cellaux(3) = abs(cellaux(3)) - 1
-                else
-                    stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%beg)."
-                end if
+                cellaux(3) = abs(cellaux(3)) - 1
             end if
             if (bc_z%end == -2 .and. (cell(3) >= p + 1 - mapCells)) then
-                if (cell(3) <= p) then
-                    cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
-                else
-                    stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%end)."
-                end if
+                cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
             end if
         end if
 
diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
index a7f6665df4..67a2d0c49d 100644
--- a/src/simulation/m_checker.fpp
+++ b/src/simulation/m_checker.fpp
@@ -330,6 +330,7 @@ contains
         @:PROHIBIT(bubbles_lagrange .and. file_per_process, "file_per_process must be false for bubbles_lagrange")
         @:PROHIBIT(bubbles_lagrange .and. n==0, "bubbles_lagrange accepts 2D and 3D simulations only")
         @:PROHIBIT(bubbles_lagrange .and. model_eqns==3, "The 6-equation flow model does not support bubbles_lagrange")
+        @:PROHIBIT(lag_params%cluster_type>=2 .and. lag_params%smooth_type/=1, "cluster_type=2 requires smooth_type=1")
     end subroutine s_check_inputs_bubbles_lagrange
 
     !> Checks miscellaneous constraints,
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index 39d3271c3c..9abf9ec205 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -1395,9 +1395,17 @@ contains
         end if
 
         if (bubbles_lagrange) then
+            !$acc update host(intfc_rad)
+            do i = 1, nBubs
+                if (ieee_is_nan(intfc_rad(i, 1)) .or. intfc_rad(i, 1) <= 0._wp) then
+                    print *, "Bubble radius is negative or NaN", proc_rank, t_step, i, intfc_rad(i, 1)
+                    error stop "Bubble radius is negative or NaN, please reduce dt"
+                end if
+            end do
+
             !$acc update host(q_beta%vf(1)%sf)
             call s_write_data_files(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, save_count, q_beta%vf(1))
-            !$acc update host(Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_rad, intfc_vel)
+            !$acc update host(Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_vel)
             call s_write_restart_lag_bubbles(save_count) !parallel
             if (lag_params%write_bubbles_stats) call s_write_lag_bubble_stats()
         else

From 102bee30b5a917d4514e81b41b7fb22cadc026b6 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Tue, 1 Apr 2025 10:00:45 -0400
Subject: [PATCH 02/23] Remove stop commands from GPU regions.

---
 src/simulation/m_checker.fpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
index 67a2d0c49d..e5e00ca934 100644
--- a/src/simulation/m_checker.fpp
+++ b/src/simulation/m_checker.fpp
@@ -330,7 +330,8 @@ contains
         @:PROHIBIT(bubbles_lagrange .and. file_per_process, "file_per_process must be false for bubbles_lagrange")
         @:PROHIBIT(bubbles_lagrange .and. n==0, "bubbles_lagrange accepts 2D and 3D simulations only")
         @:PROHIBIT(bubbles_lagrange .and. model_eqns==3, "The 6-equation flow model does not support bubbles_lagrange")
-        @:PROHIBIT(lag_params%cluster_type>=2 .and. lag_params%smooth_type/=1, "cluster_type=2 requires smooth_type=1")
+        @:PROHIBIT(bubbles_lagrange .and. lag_params%cluster_type<=0, "cluster_type must be specified")
+        @:PROHIBIT(bubbles_lagrange .and. lag_params%cluster_type>=2 .and. lag_params%smooth_type/=1, "cluster_type=2 requires smooth_type=1")
     end subroutine s_check_inputs_bubbles_lagrange
 
     !> Checks miscellaneous constraints,

From fd5eb3b9d17ca585a6dc16795617075c1b87cca5 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Thu, 3 Apr 2025 10:13:44 -0400
Subject: [PATCH 03/23] Fix frontier tests: removing changes from m_checker.fpp

---
 src/simulation/m_checker.fpp | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
index e5e00ca934..a7f6665df4 100644
--- a/src/simulation/m_checker.fpp
+++ b/src/simulation/m_checker.fpp
@@ -330,8 +330,6 @@ contains
         @:PROHIBIT(bubbles_lagrange .and. file_per_process, "file_per_process must be false for bubbles_lagrange")
         @:PROHIBIT(bubbles_lagrange .and. n==0, "bubbles_lagrange accepts 2D and 3D simulations only")
         @:PROHIBIT(bubbles_lagrange .and. model_eqns==3, "The 6-equation flow model does not support bubbles_lagrange")
-        @:PROHIBIT(bubbles_lagrange .and. lag_params%cluster_type<=0, "cluster_type must be specified")
-        @:PROHIBIT(bubbles_lagrange .and. lag_params%cluster_type>=2 .and. lag_params%smooth_type/=1, "cluster_type=2 requires smooth_type=1")
     end subroutine s_check_inputs_bubbles_lagrange
 
     !> Checks miscellaneous constraints,

From 4a398c1ea21a1b17e4620d42482c7db11263299f Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Thu, 3 Apr 2025 11:17:32 -0400
Subject: [PATCH 04/23] Fix frontier tests: removing changes from
 m_start_up.fpp

---
 src/simulation/m_checker.fpp  |  2 ++
 src/simulation/m_start_up.fpp | 10 +---------
 2 files changed, 3 insertions(+), 9 deletions(-)

diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
index a7f6665df4..e5e00ca934 100644
--- a/src/simulation/m_checker.fpp
+++ b/src/simulation/m_checker.fpp
@@ -330,6 +330,8 @@ contains
         @:PROHIBIT(bubbles_lagrange .and. file_per_process, "file_per_process must be false for bubbles_lagrange")
         @:PROHIBIT(bubbles_lagrange .and. n==0, "bubbles_lagrange accepts 2D and 3D simulations only")
         @:PROHIBIT(bubbles_lagrange .and. model_eqns==3, "The 6-equation flow model does not support bubbles_lagrange")
+        @:PROHIBIT(bubbles_lagrange .and. lag_params%cluster_type<=0, "cluster_type must be specified")
+        @:PROHIBIT(bubbles_lagrange .and. lag_params%cluster_type>=2 .and. lag_params%smooth_type/=1, "cluster_type=2 requires smooth_type=1")
     end subroutine s_check_inputs_bubbles_lagrange
 
     !> Checks miscellaneous constraints,
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index 9abf9ec205..39d3271c3c 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -1395,17 +1395,9 @@ contains
         end if
 
         if (bubbles_lagrange) then
-            !$acc update host(intfc_rad)
-            do i = 1, nBubs
-                if (ieee_is_nan(intfc_rad(i, 1)) .or. intfc_rad(i, 1) <= 0._wp) then
-                    print *, "Bubble radius is negative or NaN", proc_rank, t_step, i, intfc_rad(i, 1)
-                    error stop "Bubble radius is negative or NaN, please reduce dt"
-                end if
-            end do
-
             !$acc update host(q_beta%vf(1)%sf)
             call s_write_data_files(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, save_count, q_beta%vf(1))
-            !$acc update host(Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_vel)
+            !$acc update host(Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_rad, intfc_vel)
             call s_write_restart_lag_bubbles(save_count) !parallel
             if (lag_params%write_bubbles_stats) call s_write_lag_bubble_stats()
         else

From 9cc0626a39bf8cd7ab61fbc4b424defb874bb74f Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Thu, 3 Apr 2025 12:37:15 -0400
Subject: [PATCH 05/23] Fix frontier tests: removing changes from
 m_bubbles_EL_kernels.fpp

---
 src/simulation/m_bubbles_EL_kernels.fpp | 36 ++++++++++++++++++++-----
 src/simulation/m_start_up.fpp           | 10 ++++++-
 2 files changed, 39 insertions(+), 7 deletions(-)

diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index b45832cae9..43246aad55 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -317,27 +317,51 @@ contains
 
         ! x-dir
         if (bc_x%beg == -2 .and. (cell(1) <= mapCells - 1)) then
-            cellaux(1) = abs(cellaux(1)) - 1
+            if (cell(1) >= 0) then
+                cellaux(1) = abs(cellaux(1)) - 1
+            else
+                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%beg)."
+            end if
         end if
         if (bc_x%end == -2 .and. (cell(1) >= m + 1 - mapCells)) then
-            cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
+            if (cell(1) <= m) then
+                cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
+            else
+                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%end)."
+            end if
         end if
 
         !y-dir
         if (bc_y%beg == -2 .and. (cell(2) <= mapCells - 1)) then
-            cellaux(2) = abs(cellaux(2)) - 1
+            if (cell(2) >= 0) then
+                cellaux(2) = abs(cellaux(2)) - 1
+            else
+                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%beg)."
+            end if
         end if
         if (bc_y%end == -2 .and. (cell(2) >= n + 1 - mapCells)) then
-            cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
+            if (cell(2) <= n) then
+                cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
+            else
+                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%end)."
+            end if
         end if
 
         if (p > 0) then
             !z-dir
             if (bc_z%beg == -2 .and. (cell(3) <= mapCells - 1)) then
-                cellaux(3) = abs(cellaux(3)) - 1
+                if (cell(3) >= 0) then
+                    cellaux(3) = abs(cellaux(3)) - 1
+                else
+                    stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%beg)."
+                end if
             end if
             if (bc_z%end == -2 .and. (cell(3) >= p + 1 - mapCells)) then
-                cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
+                if (cell(3) <= p) then
+                    cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
+                else
+                    stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%end)."
+                end if
             end if
         end if
 
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index 39d3271c3c..9abf9ec205 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -1395,9 +1395,17 @@ contains
         end if
 
         if (bubbles_lagrange) then
+            !$acc update host(intfc_rad)
+            do i = 1, nBubs
+                if (ieee_is_nan(intfc_rad(i, 1)) .or. intfc_rad(i, 1) <= 0._wp) then
+                    print *, "Bubble radius is negative or NaN", proc_rank, t_step, i, intfc_rad(i, 1)
+                    error stop "Bubble radius is negative or NaN, please reduce dt"
+                end if
+            end do
+
             !$acc update host(q_beta%vf(1)%sf)
             call s_write_data_files(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, save_count, q_beta%vf(1))
-            !$acc update host(Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_rad, intfc_vel)
+            !$acc update host(Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_vel)
             call s_write_restart_lag_bubbles(save_count) !parallel
             if (lag_params%write_bubbles_stats) call s_write_lag_bubble_stats()
         else

From 00fc3e90d1899c4b08f2e90c443962404364d47b Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-03.turing.wpi.edu>
Date: Thu, 3 Apr 2025 15:42:09 -0400
Subject: [PATCH 06/23] Fix frontier tests: removing changes in GPU regions
 from EL.fpp and EL_kernels.fpp

---
 src/simulation/m_bubbles_EL.fpp | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index ba0b1b164a..151eb9f40d 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -914,7 +914,11 @@ contains
             ! Bubble dynamic closure from Maeda and Colonius (2018)
 
             ! Range of cells included in Omega
-            mapCells_pinf = mapCells
+            if (lag_params%smooth_type == 1) then
+                mapCells_pinf = mapCells
+            else
+                stop "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
+            end if
 
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
             smearGrid = mapCells_pinf - (-mapCells_pinf) + 1
@@ -992,6 +996,10 @@ contains
             vol = charvol
             dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
 
+        else
+
+            stop "Check cluterflag. Exiting."
+
         end if
 
         if (lag_params%pressure_corrector) then

From 4330e06069df3263e86a2271ffa087fd71749519 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-03.turing.wpi.edu>
Date: Thu, 3 Apr 2025 17:25:28 -0400
Subject: [PATCH 07/23] Fixed frontier attempt 1

---
 src/simulation/m_bubbles_EL.fpp         | 12 ++-------
 src/simulation/m_bubbles_EL_kernels.fpp | 36 +++++--------------------
 2 files changed, 8 insertions(+), 40 deletions(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 151eb9f40d..d1dd3eaa3e 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -910,15 +910,11 @@ contains
             !R_Omega
             dc = (3._wp*vol/(4._wp*pi))**(1._wp/3._wp)
 
-        else if (lag_params%cluster_type >= 2) then
+        else
             ! Bubble dynamic closure from Maeda and Colonius (2018)
 
             ! Range of cells included in Omega
-            if (lag_params%smooth_type == 1) then
-                mapCells_pinf = mapCells
-            else
-                stop "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
-            end if
+            mapCells_pinf = mapCells
 
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
             smearGrid = mapCells_pinf - (-mapCells_pinf) + 1
@@ -996,10 +992,6 @@ contains
             vol = charvol
             dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
 
-        else
-
-            stop "Check cluterflag. Exiting."
-
         end if
 
         if (lag_params%pressure_corrector) then
diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index 43246aad55..b45832cae9 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -317,51 +317,27 @@ contains
 
         ! x-dir
         if (bc_x%beg == -2 .and. (cell(1) <= mapCells - 1)) then
-            if (cell(1) >= 0) then
-                cellaux(1) = abs(cellaux(1)) - 1
-            else
-                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%beg)."
-            end if
+            cellaux(1) = abs(cellaux(1)) - 1
         end if
         if (bc_x%end == -2 .and. (cell(1) >= m + 1 - mapCells)) then
-            if (cell(1) <= m) then
-                cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
-            else
-                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%end)."
-            end if
+            cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
         end if
 
         !y-dir
         if (bc_y%beg == -2 .and. (cell(2) <= mapCells - 1)) then
-            if (cell(2) >= 0) then
-                cellaux(2) = abs(cellaux(2)) - 1
-            else
-                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%beg)."
-            end if
+            cellaux(2) = abs(cellaux(2)) - 1
         end if
         if (bc_y%end == -2 .and. (cell(2) >= n + 1 - mapCells)) then
-            if (cell(2) <= n) then
-                cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
-            else
-                stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%end)."
-            end if
+            cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
         end if
 
         if (p > 0) then
             !z-dir
             if (bc_z%beg == -2 .and. (cell(3) <= mapCells - 1)) then
-                if (cell(3) >= 0) then
-                    cellaux(3) = abs(cellaux(3)) - 1
-                else
-                    stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%beg)."
-                end if
+                cellaux(3) = abs(cellaux(3)) - 1
             end if
             if (bc_z%end == -2 .and. (cell(3) >= p + 1 - mapCells)) then
-                if (cell(3) <= p) then
-                    cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
-                else
-                    stop "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%end)."
-                end if
+                cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
             end if
         end if
 

From 1874569e97ea4bc73d784fae97883912a12aed51 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Fri, 4 Apr 2025 10:15:42 -0400
Subject: [PATCH 08/23] Fixed frontier attempt 2

---
 src/simulation/m_bubbles_EL.fpp | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index d1dd3eaa3e..c719f0ca74 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -914,7 +914,9 @@ contains
             ! Bubble dynamic closure from Maeda and Colonius (2018)
 
             ! Range of cells included in Omega
-            mapCells_pinf = mapCells
+            if (lag_params%smooth_type == 1) then
+                mapCells_pinf = mapCells
+            end if
 
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
             smearGrid = mapCells_pinf - (-mapCells_pinf) + 1

From ce43b8241494983ffaa84a7219105c10df74b38c Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Fri, 4 Apr 2025 11:40:09 -0400
Subject: [PATCH 09/23] Fixed frontier attempt 3

---
 src/simulation/m_bubbles_EL.fpp | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index c719f0ca74..635994b7d6 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -910,12 +910,14 @@ contains
             !R_Omega
             dc = (3._wp*vol/(4._wp*pi))**(1._wp/3._wp)
 
-        else
+        else if (lag_params%cluster_type >= 2) then
             ! Bubble dynamic closure from Maeda and Colonius (2018)
 
             ! Range of cells included in Omega
             if (lag_params%smooth_type == 1) then
                 mapCells_pinf = mapCells
+            else
+                print *, "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
             end if
 
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
@@ -993,6 +995,9 @@ contains
             f_pinfl = charpres2/charvol2
             vol = charvol
             dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
+        else
+
+            print *, "Check cluterflag. Exiting."
 
         end if
 

From 35f71e2d1cec6e96699877490c27661378816c35 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Fri, 4 Apr 2025 13:06:40 -0400
Subject: [PATCH 10/23] Fixed frontier attempt 4

---
 src/simulation/m_bubbles_EL_kernels.fpp | 36 ++++++++++++++++++++-----
 1 file changed, 30 insertions(+), 6 deletions(-)

diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index b45832cae9..97d1ebefb1 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -317,27 +317,51 @@ contains
 
         ! x-dir
         if (bc_x%beg == -2 .and. (cell(1) <= mapCells - 1)) then
-            cellaux(1) = abs(cellaux(1)) - 1
+            if (cell(1) >= 0) then
+                cellaux(1) = abs(cellaux(1)) - 1
+            else
+                print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%beg)."
+            end if
         end if
         if (bc_x%end == -2 .and. (cell(1) >= m + 1 - mapCells)) then
-            cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
+            if (cell(1) <= m) then
+                cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
+            else
+                print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%end)."
+            end if
         end if
 
         !y-dir
         if (bc_y%beg == -2 .and. (cell(2) <= mapCells - 1)) then
-            cellaux(2) = abs(cellaux(2)) - 1
+            if (cell(2) >= 0) then
+                cellaux(2) = abs(cellaux(2)) - 1
+            else
+                print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%beg)."
+            end if
         end if
         if (bc_y%end == -2 .and. (cell(2) >= n + 1 - mapCells)) then
-            cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
+            if (cell(2) <= n) then
+                cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
+            else
+                print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%end)."
+            end if
         end if
 
         if (p > 0) then
             !z-dir
             if (bc_z%beg == -2 .and. (cell(3) <= mapCells - 1)) then
-                cellaux(3) = abs(cellaux(3)) - 1
+                if (cell(3) >= 0) then
+                    cellaux(3) = abs(cellaux(3)) - 1
+                else
+                    print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%beg)."
+                end if
             end if
             if (bc_z%end == -2 .and. (cell(3) >= p + 1 - mapCells)) then
-                cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
+                if (cell(3) <= p) then
+                    cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
+                else
+                    print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%end)."
+                end if
             end if
         end if
 

From 8fbecffb59baf8f9497114ff048197d9344b4a97 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Fri, 4 Apr 2025 15:18:31 -0400
Subject: [PATCH 11/23] Isolating bug 1

---
 src/simulation/m_bubbles_EL.fpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 635994b7d6..41e12e3a43 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -917,7 +917,7 @@ contains
             if (lag_params%smooth_type == 1) then
                 mapCells_pinf = mapCells
             else
-                print *, "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
+                stop "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
             end if
 
             ! Include the cell that contains the bubble (mapCells+1+mapCells)

From 0771c04d2b37019363f3c365414a79c088e62389 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-02.cm.cluster>
Date: Fri, 4 Apr 2025 16:15:35 -0400
Subject: [PATCH 12/23] Isolating bug 2

---
 src/simulation/m_bubbles_EL.fpp | 11 ++---------
 1 file changed, 2 insertions(+), 9 deletions(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 41e12e3a43..7ac64aaed1 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -825,7 +825,7 @@ contains
         real(wp) :: charvol, charpres, charvol2, charpres2
         integer, dimension(3) :: cellaux
         integer :: i, j, k
-        integer :: mapCells_pinf, smearGrid, smearGridz
+        integer :: smearGrid, smearGridz
         logical :: celloutside
 
         scoord = mtn_s(bub_id, 1:3, 2)
@@ -913,15 +913,8 @@ contains
         else if (lag_params%cluster_type >= 2) then
             ! Bubble dynamic closure from Maeda and Colonius (2018)
 
-            ! Range of cells included in Omega
-            if (lag_params%smooth_type == 1) then
-                mapCells_pinf = mapCells
-            else
-                stop "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
-            end if
-
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
-            smearGrid = mapCells_pinf - (-mapCells_pinf) + 1
+            smearGrid = mapCells - (-mapCells) + 1
             smearGridz = smearGrid
             if (p == 0) smearGridz = 1
 

From 20809e4610699bf5f03a59614a44513aa5bd2d7f Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-01.cm.cluster>
Date: Fri, 4 Apr 2025 17:21:53 -0400
Subject: [PATCH 13/23] Isolating bug 3

---
 src/simulation/m_bubbles_EL.fpp | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 7ac64aaed1..e0ca2f7460 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -825,7 +825,7 @@ contains
         real(wp) :: charvol, charpres, charvol2, charpres2
         integer, dimension(3) :: cellaux
         integer :: i, j, k
-        integer :: smearGrid, smearGridz
+        integer :: mapCells_pinf, smearGrid, smearGridz
         logical :: celloutside
 
         scoord = mtn_s(bub_id, 1:3, 2)
@@ -912,9 +912,14 @@ contains
 
         else if (lag_params%cluster_type >= 2) then
             ! Bubble dynamic closure from Maeda and Colonius (2018)
+            if (lag_params%smooth_type == 1) then
+                mapCells_pinf = mapCells
+            else
+                stop "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
+            end if
 
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
-            smearGrid = mapCells - (-mapCells) + 1
+            smearGrid = mapCells_pinf - (-mapCells_pinf) + 1
             smearGridz = smearGrid
             if (p == 0) smearGridz = 1
 
@@ -988,9 +993,6 @@ contains
             f_pinfl = charpres2/charvol2
             vol = charvol
             dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
-        else
-
-            print *, "Check cluterflag. Exiting."
 
         end if
 

From e9a47c36da4f19188c816803ce3dd9d891f0935c Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-03.turing.wpi.edu>
Date: Sat, 5 Apr 2025 09:52:29 -0400
Subject: [PATCH 14/23] Isolating bug 4

---
 src/simulation/m_bubbles_EL_kernels.fpp | 36 +++++--------------------
 1 file changed, 6 insertions(+), 30 deletions(-)

diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index 97d1ebefb1..b45832cae9 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -317,51 +317,27 @@ contains
 
         ! x-dir
         if (bc_x%beg == -2 .and. (cell(1) <= mapCells - 1)) then
-            if (cell(1) >= 0) then
-                cellaux(1) = abs(cellaux(1)) - 1
-            else
-                print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%beg)."
-            end if
+            cellaux(1) = abs(cellaux(1)) - 1
         end if
         if (bc_x%end == -2 .and. (cell(1) >= m + 1 - mapCells)) then
-            if (cell(1) <= m) then
-                cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
-            else
-                print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_x%end)."
-            end if
+            cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
         end if
 
         !y-dir
         if (bc_y%beg == -2 .and. (cell(2) <= mapCells - 1)) then
-            if (cell(2) >= 0) then
-                cellaux(2) = abs(cellaux(2)) - 1
-            else
-                print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%beg)."
-            end if
+            cellaux(2) = abs(cellaux(2)) - 1
         end if
         if (bc_y%end == -2 .and. (cell(2) >= n + 1 - mapCells)) then
-            if (cell(2) <= n) then
-                cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
-            else
-                print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_y%end)."
-            end if
+            cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
         end if
 
         if (p > 0) then
             !z-dir
             if (bc_z%beg == -2 .and. (cell(3) <= mapCells - 1)) then
-                if (cell(3) >= 0) then
-                    cellaux(3) = abs(cellaux(3)) - 1
-                else
-                    print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%beg)."
-                end if
+                cellaux(3) = abs(cellaux(3)) - 1
             end if
             if (bc_z%end == -2 .and. (cell(3) >= p + 1 - mapCells)) then
-                if (cell(3) <= p) then
-                    cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
-                else
-                    print *, "Lagrangian bubbles must not be located in the ghost cells of a symmetric boundary (bc_z%end)."
-                end if
+                cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
             end if
         end if
 

From b02693667a1347980e47bb0949194c197168231e Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-03.turing.wpi.edu>
Date: Sat, 5 Apr 2025 11:00:56 -0400
Subject: [PATCH 15/23] Fixing isolated bug 1

---
 src/simulation/m_bubbles_EL.fpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index e0ca2f7460..69ff604b33 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -915,7 +915,7 @@ contains
             if (lag_params%smooth_type == 1) then
                 mapCells_pinf = mapCells
             else
-                stop "lag_params%cluster_type: 2 requires lag_params%smooth_type: 1."
+                mapCells_pinf = 0
             end if
 
             ! Include the cell that contains the bubble (mapCells+1+mapCells)

From 654bc3cfc96ded1395c16b327a1ffb1a1d28672e Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-01.cm.cluster>
Date: Sat, 5 Apr 2025 12:29:21 -0400
Subject: [PATCH 16/23] Fixing isolated bug 2

---
 src/simulation/m_bubbles_EL.fpp | 16 +++++-----------
 1 file changed, 5 insertions(+), 11 deletions(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 69ff604b33..0e07f2ead7 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -825,7 +825,7 @@ contains
         real(wp) :: charvol, charpres, charvol2, charpres2
         integer, dimension(3) :: cellaux
         integer :: i, j, k
-        integer :: mapCells_pinf, smearGrid, smearGridz
+        integer :: smearGrid, smearGrid_zDir
         logical :: celloutside
 
         scoord = mtn_s(bub_id, 1:3, 2)
@@ -912,16 +912,10 @@ contains
 
         else if (lag_params%cluster_type >= 2) then
             ! Bubble dynamic closure from Maeda and Colonius (2018)
-            if (lag_params%smooth_type == 1) then
-                mapCells_pinf = mapCells
-            else
-                mapCells_pinf = 0
-            end if
-
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
-            smearGrid = mapCells_pinf - (-mapCells_pinf) + 1
-            smearGridz = smearGrid
-            if (p == 0) smearGridz = 1
+            smearGrid = 2*mapCells + 1
+            smearGrid_zDir = smearGrid
+            if (p == 0) smearGrid_zDir = 1
 
             charvol = 0._wp
             charpres = 0._wp
@@ -934,7 +928,7 @@ contains
                 !$acc loop seq
                 do j = 1, smearGrid
                     !$acc loop seq
-                    do k = 1, smearGridz
+                    do k = 1, smearGrid_zDir
                         cellaux(1) = cell(1) + i - (mapCells + 1)
                         cellaux(2) = cell(2) + j - (mapCells + 1)
                         cellaux(3) = cell(3) + k - (mapCells + 1)

From e2c7f7c43186a6beeb3adddf5d6fe2bf61e2ff1d Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-03.turing.wpi.edu>
Date: Sun, 6 Apr 2025 10:41:25 -0400
Subject: [PATCH 17/23] p_inf from subroutine to function

---
 src/simulation/m_bubbles_EL.fpp | 65 +++++++++++++++------------------
 1 file changed, 30 insertions(+), 35 deletions(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 0e07f2ead7..573e76feef 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -518,10 +518,10 @@ contains
         integer, intent(in) :: stage
 
         real(wp) :: vaporflux, pliqint
-        real(wp) :: preterm1, term2, paux, pint, Romega, term1_fac, Rb
+        real(wp) :: preterm1, term2, pint, Romega, term1_fac, Rb
         real(wp) :: conc_v, R_m, gamma_m, fpb, fmass_n, fmass_v
         real(wp) :: fR, fV, fbeta_c, fbeta_t, fR0, fpbdt
-        real(wp) :: pinf, aux1, aux2, velint, cson, rhol
+        real(wp) :: pinf, velint, cson, rhol
         real(wp) :: gamma, pi_inf, qv
         real(wp), dimension(contxe) :: myalpha_rho, myalpha
         real(wp), dimension(2) :: Re
@@ -538,7 +538,7 @@ contains
             ! Calculate velocity potentials (valid for one bubble per cell)
             !$acc parallel loop gang vector default(present) private(k, cell)
             do k = 1, nBubs
-                call s_get_pinf(k, q_prim_vf, 2, paux, cell, preterm1, term2, Romega)
+                pinf = f_pinf(k, q_prim_vf, 2, cell, preterm1, term2, Romega)
                 fR0 = bub_R0(k)
                 fR = intfc_rad(k, 2)
                 fV = intfc_vel(k, 2)
@@ -546,11 +546,10 @@ contains
                 pint = f_cpbw_KM(fR0, fR, fV, fpb)
                 pint = pint + 0.5_wp*fV**2._wp
                 if (lag_params%cluster_type == 2) then
-                    bub_dphidt(k) = (paux - pint) + term2
+                    bub_dphidt(k) = (pinf - pint) + term2
                     ! Accounting for the potential induced by the bubble averaged over the control volume
                     ! Note that this is based on the incompressible flow assumption near the bubble.
-                    Rb = intfc_rad(k, 2)
-                    term1_fac = 3._wp/2._wp*(Rb*(Romega**2._wp - Rb**2._wp))/(Romega**3._wp - Rb**3._wp)
+                    term1_fac = 3._wp/2._wp*(fR*(Romega**2._wp - fR**2._wp))/(Romega**3._wp - fR**3._wp)
                     bub_dphidt(k) = bub_dphidt(k)/(1._wp - term1_fac)
                 end if
             end do
@@ -587,7 +586,7 @@ contains
                 pliqint = f_cpbw_KM(fR0, fR, fV, fpb)
 
                 ! Obtaining driving pressure
-                call s_get_pinf(k, q_prim_vf, 1, pinf, cell, aux1, aux2)
+                pinf = f_pinf(k, q_prim_vf, 1, cell, preterm1, term2, Romega)
 
                 ! Obtain liquid density and computing speed of sound from pinf
                 !$acc loop seq
@@ -804,20 +803,16 @@ contains
         !! @param bub_id Particle identifier
         !! @param q_prim_vf  Primitive variables
         !! @param ptype 1: p at infinity, 2: averaged P at the bubble location
-        !! @param f_pinfl Driving pressure
         !! @param cell Bubble cell
         !! @param Romega Control volume radius
-    subroutine s_get_pinf(bub_id, q_prim_vf, ptype, f_pinfl, cell, preterm1, term2, Romega)
-#ifdef _CRAYFTN
-        !DIR$ INLINEALWAYS s_get_pinf
-#else
+    function f_pinf(bub_id, q_prim_vf, ptype, cell, preterm1, term2, Romega)
         !$acc routine seq
-#endif
         integer, intent(in) :: bub_id, ptype
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
-        real(wp), intent(out) :: f_pinfl
         integer, dimension(3), intent(out) :: cell
-        real(wp), intent(out), optional :: preterm1, term2, Romega
+        real(wp), intent(out) :: preterm1, term2, Romega
+
+        real(wp) :: f_pinf
 
         real(wp), dimension(3) :: scoord, psi
         real(wp) :: dc, vol, aux
@@ -825,11 +820,11 @@ contains
         real(wp) :: charvol, charpres, charvol2, charpres2
         integer, dimension(3) :: cellaux
         integer :: i, j, k
-        integer :: smearGrid, smearGrid_zDir
+        integer :: smearGrid, smearGridz
         logical :: celloutside
 
         scoord = mtn_s(bub_id, 1:3, 2)
-        f_pinfl = 0._wp
+        f_pinf = 0._wp
 
         !< Find current bubble cell
         cell(:) = int(scoord(:))
@@ -892,19 +887,19 @@ contains
 
             !< Perform bilinear interpolation
             if (p == 0) then  !2D
-                f_pinfl = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)
+                f_pinf = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)
             else              !3D
-                f_pinfl = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))*(1._wp - psi(3))
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))*(1._wp - psi(3))
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)*(1._wp - psi(3))
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)*(1._wp - psi(3))
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1)*(1._wp - psi(1))*(1._wp - psi(2))*psi(3)
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3) + 1)*psi(1)*(1._wp - psi(2))*psi(3)
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3) + 1)*psi(1)*psi(2)*psi(3)
-                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3) + 1)*(1._wp - psi(1))*psi(2)*psi(3)
+                f_pinf = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))*(1._wp - psi(3))
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))*(1._wp - psi(3))
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)*(1._wp - psi(3))
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)*(1._wp - psi(3))
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1)*(1._wp - psi(1))*(1._wp - psi(2))*psi(3)
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3) + 1)*psi(1)*(1._wp - psi(2))*psi(3)
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3) + 1)*psi(1)*psi(2)*psi(3)
+                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3) + 1)*(1._wp - psi(1))*psi(2)*psi(3)
             end if
 
             !R_Omega
@@ -914,8 +909,8 @@ contains
             ! Bubble dynamic closure from Maeda and Colonius (2018)
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
             smearGrid = 2*mapCells + 1
-            smearGrid_zDir = smearGrid
-            if (p == 0) smearGrid_zDir = 1
+            smearGridz = smearGrid
+            if (p == 0) smearGridz = 1
 
             charvol = 0._wp
             charpres = 0._wp
@@ -928,7 +923,7 @@ contains
                 !$acc loop seq
                 do j = 1, smearGrid
                     !$acc loop seq
-                    do k = 1, smearGrid_zDir
+                    do k = 1, smearGridz
                         cellaux(1) = cell(1) + i - (mapCells + 1)
                         cellaux(2) = cell(2) + j - (mapCells + 1)
                         cellaux(3) = cell(3) + k - (mapCells + 1)
@@ -984,7 +979,7 @@ contains
                 end do
             end do
 
-            f_pinfl = charpres2/charvol2
+            f_pinf = charpres2/charvol2
             vol = charvol
             dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
 
@@ -1009,12 +1004,12 @@ contains
 
             ! Getting p_inf
             if (ptype == 1) then
-                f_pinfl = f_pinfl + preterm1*term1 + term2
+                f_pinf = f_pinf + preterm1*term1 + term2
             end if
 
         end if
 
-    end subroutine s_get_pinf
+    end function f_pinf
 
     !>  This subroutine updates the Lagrange variables using the tvd RK time steppers.
         !!      The time derivative of the bubble variables must be stored at every stage to avoid precision errors.

From 665321bd387db0d116c227c4ac40fdf2c854fe7f Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-01.cm.cluster>
Date: Mon, 7 Apr 2025 10:28:15 -0400
Subject: [PATCH 18/23] printing statements in rkck

---
 src/simulation/m_bubbles_EL.fpp    | 20 ++++++++++++++------
 src/simulation/m_time_steppers.fpp | 19 +++++++++++++++++++
 2 files changed, 33 insertions(+), 6 deletions(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 573e76feef..79b911d08d 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -533,12 +533,14 @@ contains
 
         !< BUBBLE DYNAMICS
 
+        print *, 'bub_dphidt kernel loop'
         ! Subgrid p_inf model based on Maeda and Colonius (2018).
         if (lag_params%pressure_corrector) then
             ! Calculate velocity potentials (valid for one bubble per cell)
             !$acc parallel loop gang vector default(present) private(k, cell)
             do k = 1, nBubs
                 pinf = f_pinf(k, q_prim_vf, 2, cell, preterm1, term2, Romega)
+                print *, 'f_pinf'
                 fR0 = bub_R0(k)
                 fR = intfc_rad(k, 2)
                 fV = intfc_vel(k, 2)
@@ -555,6 +557,7 @@ contains
             end do
         end if
 
+        print *, 'bub dynamics kernel loop'
         ! Radial motion model
         if (bubble_model == 2) then
             !$acc parallel loop gang vector default(present) private(k, myalpha_rho, myalpha, Re, cell) copyin(stage)
@@ -580,14 +583,15 @@ contains
                 ! Vapor and heat fluxes
                 vaporflux = f_vflux(fR, fV, fmass_v, k, fmass_n, fbeta_c, conc_v)
                 fpbdt = f_bpres_dot(vaporflux, fR, fV, fpb, fmass_v, k, fbeta_t, R_m, gamma_m, conc_v)
+                print *, 'f_vflux and f_bpres_dot'
                 gas_dmvdt(k, stage) = 4._wp*pi*fR**2._wp*vaporflux
-
+                print *, 'updated gas_dmvdt', stage
                 ! Pressure at the bubble wall
                 pliqint = f_cpbw_KM(fR0, fR, fV, fpb)
-
+                print *, 'f_cpbw_KM'
                 ! Obtaining driving pressure
                 pinf = f_pinf(k, q_prim_vf, 1, cell, preterm1, term2, Romega)
-
+                print *, 'f_pinf'
                 ! Obtain liquid density and computing speed of sound from pinf
                 !$acc loop seq
                 do i = 1, contxe
@@ -598,14 +602,15 @@ contains
                 call s_convert_species_to_mixture_variables_acc(rhol, gamma, pi_inf, qv, myalpha, &
                                                                 myalpha_rho, Re, cell(1), cell(2), cell(3))
                 call s_compute_cson_from_pinf(k, q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
-
+                print *, 'cson and rhol'
                 ! Velocity correction due to massflux
                 velint = fV - gas_dmvdt(k, stage)/(4._wp*pi*fR**2._wp*rhol)
-
+                print *, 'velint'
                 ! Interphase acceleration and update vars
                 intfc_dveldt(k, stage) = f_rddot_KM(fpbdt, pinf, pliqint, rhol, fR, velint, fR0, cson)
                 gas_dpdt(k, stage) = fpbdt
                 intfc_draddt(k, stage) = fV
+                print *, 'update intfc_dveldt, gas_dpdt, intfc_draddt'
 
             end do
         else
@@ -618,6 +623,7 @@ contains
             end do
         end if
 
+        print *, 'bub position kernel loop'
         ! Bubbles remain in a fixed position
         !$acc parallel loop collapse(2) gang vector default(present) private(k) copyin(stage)
         do k = 1, nBubs
@@ -630,7 +636,9 @@ contains
         call nvtxEndRange
 
         !< EULER-LAGRANGE COUPLING
+        print *, 's_smear_voidfraction'
         call s_smear_voidfraction()
+        print *, 's_add_rhs_sources'
         if (lag_params%solver_approach == 2) call s_add_rhs_sources(q_cons_vf, q_prim_vf, rhs_vf)
 
     end subroutine s_compute_EL_coupled_solver
@@ -1000,7 +1008,7 @@ contains
             preterm1 = 3._wp/2._wp*Rbeq*(dc**2._wp - Rbeq**2._wp)/(aux*denom)
 
             !Control volume radius
-            if (ptype == 2) Romega = dc
+            Romega = dc
 
             ! Getting p_inf
             if (ptype == 1) then
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index 2ffbc8e58b..2e9a2f2a6c 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -1075,11 +1075,16 @@ contains
 #ifdef DEBUG
             if (proc_rank == 0) print *, 'RKCK 1st time-stage at', rkck_time_tmp
 #endif
+            print *, 's_compute_rhs', RKstep
             call s_compute_rhs(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, rhs_ts_rkck(1)%vf, pb_ts(1)%sf, rhs_pb, mv_ts(1)%sf, rhs_mv, t_step, time_avg)
+            print *, 's_compute_EL_coupled_solver', RKstep
             call s_compute_EL_coupled_solver(q_cons_ts(1)%vf, q_prim_vf, rhs_ts_rkck(1)%vf, RKstep)
+            print *, 's_update_tmp_rkck', RKstep
             call s_update_tmp_rkck(RKstep, q_cons_ts, rhs_ts_rkck, lag_largestep)
+            print *, 's_compute_rkck_dt', RKstep, lag_largestep
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
+            print *, 'done with stage 1'
 
             ! SECOND TIME-STAGE
             RKstep = 2
@@ -1089,11 +1094,16 @@ contains
 #ifdef DEBUG
             if (proc_rank == 0) print *, 'RKCK 2nd time-stage at', rkck_time_tmp
 #endif
+            print *, 's_compute_rhs', RKstep
             call s_compute_rhs(q_cons_ts(2)%vf, q_T_sf, q_prim_vf, rhs_ts_rkck(2)%vf, pb_ts(1)%sf, rhs_pb, mv_ts(1)%sf, rhs_mv, t_step, time_avg)
+            print *, 's_compute_EL_coupled_solver', RKstep
             call s_compute_EL_coupled_solver(q_cons_ts(2)%vf, q_prim_vf, rhs_ts_rkck(2)%vf, RKstep)
+            print *, 's_update_tmp_rkck', RKstep
             call s_update_tmp_rkck(RKstep, q_cons_ts, rhs_ts_rkck, lag_largestep)
+            print *, 's_compute_rkck_dt', RKstep, lag_largestep
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
+            print *, 'done with stage 2'
 
             ! THIRD TIME-STAGE
             RKstep = 3
@@ -1108,6 +1118,7 @@ contains
             call s_update_tmp_rkck(RKstep, q_cons_ts, rhs_ts_rkck, lag_largestep)
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
+            print *, 'done with stage 3'
 
             ! FOURTH TIME-STAGE
             RKstep = 4
@@ -1122,6 +1133,7 @@ contains
             call s_update_tmp_rkck(RKstep, q_cons_ts, rhs_ts_rkck, lag_largestep)
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
+            print *, 'done with stage 4'
 
             ! FIFTH TIME-STAGE
             RKstep = 5
@@ -1136,6 +1148,7 @@ contains
             call s_update_tmp_rkck(5, q_cons_ts, rhs_ts_rkck, lag_largestep)
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
+            print *, 'done with stage 5'
 
             ! SIXTH TIME-STAGE
             RKstep = 6
@@ -1150,6 +1163,7 @@ contains
             call s_update_tmp_rkck(6, q_cons_ts, rhs_ts_rkck, lag_largestep)
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
+            print *, 'done with stage 6'
 
             dt_did = dt
 
@@ -1158,9 +1172,12 @@ contains
 #ifdef DEBUG
                 if (proc_rank == 0) print *, 'Computing truncation error (4th/5th RKCK)'
 #endif
+                print *, 's_calculate_rkck_truncation_error'
                 call s_calculate_rkck_truncation_error(rkck_errmax)
+                print *, 's_compute_rkck_dt'
                 call s_compute_rkck_dt(lag_largestep, restart_rkck_step, rkck_errmax)
                 if (restart_rkck_step) cycle
+                print *, 'done with truncation error', dt
             end if
 
         end do
@@ -1168,6 +1185,7 @@ contains
         !> Update values
         mytime = mytime + dt_did
         call s_update_rkck(q_cons_ts)
+        print *, 'done with s_update_rkck'
 
         call s_write_void_evol(mytime)
         if (lag_params%write_bubbles_stats) call s_calculate_lag_bubble_stats()
@@ -1181,6 +1199,7 @@ contains
             call s_write_run_time_information(q_prim_vf, t_step)
         end if
 
+        print *, 'done with s_4th_5th_order_rkck'
     end subroutine s_4th_5th_order_rkck
 
     !> Module deallocation and/or disassociation procedures

From d1a90501489ba1ef68c6cb50033c31fb5a3bd2e9 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-01.cm.cluster>
Date: Mon, 7 Apr 2025 11:44:14 -0400
Subject: [PATCH 19/23] printing statements in smear_voidfraction

---
 src/simulation/m_bubbles_EL.fpp         |  3 +++
 src/simulation/m_bubbles_EL_kernels.fpp | 11 ++++++++++-
 2 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 79b911d08d..0ccc2ca69e 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -685,6 +685,7 @@ contains
 
         call nvtxStartRange("BUBBLES-LAGRANGE-KERNELS")
 
+        print *, 'Zeroing vars'
         !$acc parallel loop collapse(4) gang vector default(present)
         do i = 1, q_beta_idx
             do l = idwbuff(3)%beg, idwbuff(3)%end
@@ -696,9 +697,11 @@ contains
             end do
         end do
 
+        print *, 's_smoothfunction'
         call s_smoothfunction(nBubs, intfc_rad, intfc_vel, &
                               mtn_s, mtn_pos, q_beta)
 
+        print *, '1-beta'
         !Store 1-beta
         !$acc parallel loop collapse(3) gang vector default(present)
         do l = idwbuff(3)%beg, idwbuff(3)%end
diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index b45832cae9..5e36f90874 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -120,6 +120,8 @@ contains
         smearGridz = smearGrid
         if (p == 0) smearGridz = 1
 
+        print *, 's_gaussian'
+
         !$acc parallel loop gang vector default(present) private(nodecoord, l, s_coord, cell, center) copyin(smearGrid, smearGridz)
         do l = 1, nBubs
             nodecoord(1:3) = 0
@@ -128,9 +130,10 @@ contains
             s_coord(1:3) = lbk_s(l, 1:3, 2)
             center(1:2) = lbk_pos(l, 1:2, 2)
             if (p > 0) center(3) = lbk_pos(l, 3, 2)
+            print *, 'reading initial state'
             call s_get_cell(s_coord, cell)
             call s_compute_stddsv(cell, volpart, stddsv)
-
+            print *, 's_compute_stddsv'
             strength_vol = volpart
             strength_vel = 4._wp*pi*lbk_rad(l, 2)**2._wp*lbk_vel(l, 2)
 
@@ -145,6 +148,7 @@ contains
 
                         !Check if the cells intended to smear the bubbles in are in the computational domain
                         !and redefine the cells for symmetric boundary
+                        print *, 's_check_celloutside'
                         call s_check_celloutside(cellaux, celloutside)
 
                         if (.not. celloutside) then
@@ -152,6 +156,7 @@ contains
                             nodecoord(1) = x_cc(cellaux(1))
                             nodecoord(2) = y_cc(cellaux(2))
                             if (p > 0) nodecoord(3) = z_cc(cellaux(3))
+                            print *, 's_applygaussian'
                             call s_applygaussian(center, cellaux, nodecoord, stddsv, 0._wp, func)
                             if (lag_params%cluster_type >= 4) call s_applygaussian(center, cellaux, nodecoord, stddsv, 1._wp, func2)
 
@@ -159,6 +164,7 @@ contains
                             if (bc_x%beg == -2 .or. bc_x%end == -2 .or. bc_y%beg == -2 .or. bc_y%end == -2 &
                                 .or. bc_z%beg == -2 .or. bc_z%end == -2) then
                                 call s_shift_cell_symmetric_bc(cellaux, cell)
+                                print *, 's_shift_cell_symmetric_bc'
                             end if
                         else
                             func = 0._wp
@@ -169,6 +175,7 @@ contains
                             if (p == 0) cellaux(3) = 0
                         end if
 
+                        print *, 'Update 1'
                         !Update void fraction field
                         addFun1 = func*strength_vol
                         !$acc atomic update
@@ -176,6 +183,7 @@ contains
                             updatedvar%vf(1)%sf(cellaux(1), cellaux(2), cellaux(3)) &
                             + addFun1
 
+                        print *, 'Update 2'
                         !Update time derivative of void fraction
                         addFun2 = func*strength_vel
                         !$acc atomic update
@@ -192,6 +200,7 @@ contains
                                 updatedvar%vf(5)%sf(cellaux(1), cellaux(2), cellaux(3)) &
                                 + addFun3
                         end if
+                        print *, 'Update 3'
                     end do
                 end do
             end do

From 6736d6c2883eb2a335d52daf55ab93e03279d4c1 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-01.cm.cluster>
Date: Mon, 7 Apr 2025 13:04:47 -0400
Subject: [PATCH 20/23] re-defining s_gaussian kernel

---
 src/simulation/m_bubbles_EL_kernels.fpp | 62 ++++++++++++-------------
 1 file changed, 31 insertions(+), 31 deletions(-)

diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index 5e36f90874..6ae6e6188f 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -116,55 +116,55 @@ contains
         logical :: celloutside
         integer :: smearGrid, smearGridz
 
-        smearGrid = mapCells - (-mapCells) + 1 ! Include the cell that contains the bubble (3+1+3)
+        smearGrid = mapCells - (-mapCells) ! Include the cell that contains the bubble (3+1+3)
         smearGridz = smearGrid
-        if (p == 0) smearGridz = 1
+        if (p == 0) smearGridz = 0
 
         print *, 's_gaussian'
 
-        !$acc parallel loop gang vector default(present) private(nodecoord, l, s_coord, cell, center) copyin(smearGrid, smearGridz)
+        !$acc parallel loop collapse(4) gang vector default(present) private(l, s_coord, cell, center, cellaux, nodecoord) &
+        !$acc copyin(smearGrid, smearGridz)
         do l = 1, nBubs
-            nodecoord(1:3) = 0
-            center(1:3) = 0._wp
-            volpart = 4._wp/3._wp*pi*lbk_rad(l, 2)**3._wp
-            s_coord(1:3) = lbk_s(l, 1:3, 2)
-            center(1:2) = lbk_pos(l, 1:2, 2)
-            if (p > 0) center(3) = lbk_pos(l, 3, 2)
-            print *, 'reading initial state'
-            call s_get_cell(s_coord, cell)
-            call s_compute_stddsv(cell, volpart, stddsv)
-            print *, 's_compute_stddsv'
-            strength_vol = volpart
-            strength_vel = 4._wp*pi*lbk_rad(l, 2)**2._wp*lbk_vel(l, 2)
-
-            !$acc loop collapse(3) private(cellaux, nodecoord)
-            do i = 1, smearGrid
-                do j = 1, smearGrid
-                    do k = 1, smearGridz
-                        cellaux(1) = cell(1) + i - (mapCells + 1)
-                        cellaux(2) = cell(2) + j - (mapCells + 1)
-                        cellaux(3) = cell(3) + k - (mapCells + 1)
+            do i = 0, smearGrid
+                do j = 0, smearGrid
+                    do k = 0, smearGridz
+
+                        nodecoord(1:3) = 0
+                        center(1:3) = 0._wp
+                        volpart = 4._wp/3._wp*pi*lbk_rad(l, 2)**3._wp
+                        s_coord(1:3) = lbk_s(l, 1:3, 2)
+                        center(1:2) = lbk_pos(l, 1:2, 2)
+                        print *, 'reading initial state'
+                        if (p > 0) center(3) = lbk_pos(l, 3, 2)
+                        call s_get_cell(s_coord, cell)
+                        call s_compute_stddsv(cell, volpart, stddsv)
+                        strength_vol = volpart
+                        strength_vel = 4._wp*pi*lbk_rad(l, 2)**2._wp*lbk_vel(l, 2)
+                        print *, 's_compute_stddsv'
+
+                        cellaux(1) = cell(1) + i - mapCells
+                        cellaux(2) = cell(2) + j - mapCells
+                        cellaux(3) = cell(3) + k - mapCells
                         if (p == 0) cellaux(3) = 0
 
                         !Check if the cells intended to smear the bubbles in are in the computational domain
                         !and redefine the cells for symmetric boundary
-                        print *, 's_check_celloutside'
                         call s_check_celloutside(cellaux, celloutside)
+                        print *, 's_check_celloutside'
 
                         if (.not. celloutside) then
 
                             nodecoord(1) = x_cc(cellaux(1))
                             nodecoord(2) = y_cc(cellaux(2))
                             if (p > 0) nodecoord(3) = z_cc(cellaux(3))
-                            print *, 's_applygaussian'
                             call s_applygaussian(center, cellaux, nodecoord, stddsv, 0._wp, func)
-                            if (lag_params%cluster_type >= 4) call s_applygaussian(center, cellaux, nodecoord, stddsv, 1._wp, func2)
+                            if (p == 0) call s_applygaussian(center, cellaux, nodecoord, stddsv, 1._wp, func2)
+                            print *, 's_applygaussian'
 
                             ! Relocate cells for bubbles intersecting symmetric boundaries
                             if (bc_x%beg == -2 .or. bc_x%end == -2 .or. bc_y%beg == -2 .or. bc_y%end == -2 &
                                 .or. bc_z%beg == -2 .or. bc_z%end == -2) then
                                 call s_shift_cell_symmetric_bc(cellaux, cell)
-                                print *, 's_shift_cell_symmetric_bc'
                             end if
                         else
                             func = 0._wp
@@ -175,7 +175,6 @@ contains
                             if (p == 0) cellaux(3) = 0
                         end if
 
-                        print *, 'Update 1'
                         !Update void fraction field
                         addFun1 = func*strength_vol
                         !$acc atomic update
@@ -183,7 +182,6 @@ contains
                             updatedvar%vf(1)%sf(cellaux(1), cellaux(2), cellaux(3)) &
                             + addFun1
 
-                        print *, 'Update 2'
                         !Update time derivative of void fraction
                         addFun2 = func*strength_vel
                         !$acc atomic update
@@ -193,14 +191,16 @@ contains
 
                         !Product of two smeared functions
                         !Update void fraction * time derivative of void fraction
-                        if (lag_params%cluster_type >= 4) then
+                        if (p == 0) then
                             addFun3 = func2*strength_vol*strength_vel
                             !$acc atomic update
                             updatedvar%vf(5)%sf(cellaux(1), cellaux(2), cellaux(3)) = &
                                 updatedvar%vf(5)%sf(cellaux(1), cellaux(2), cellaux(3)) &
                                 + addFun3
                         end if
-                        print *, 'Update 3'
+
+                        print *, 'Update 1-2-3'
+
                     end do
                 end do
             end do

From 3773d10fdff3dd4edb86e2c0030dab82b2668b26 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-01.cm.cluster>
Date: Mon, 7 Apr 2025 15:08:32 -0400
Subject: [PATCH 21/23] removing bc_xyz from gaussian

---
 src/simulation/m_bubbles_EL_kernels.fpp | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index 6ae6e6188f..849ae7a57e 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -159,13 +159,12 @@ contains
                             if (p > 0) nodecoord(3) = z_cc(cellaux(3))
                             call s_applygaussian(center, cellaux, nodecoord, stddsv, 0._wp, func)
                             if (p == 0) call s_applygaussian(center, cellaux, nodecoord, stddsv, 1._wp, func2)
-                            print *, 's_applygaussian'
 
                             ! Relocate cells for bubbles intersecting symmetric boundaries
-                            if (bc_x%beg == -2 .or. bc_x%end == -2 .or. bc_y%beg == -2 .or. bc_y%end == -2 &
-                                .or. bc_z%beg == -2 .or. bc_z%end == -2) then
-                                call s_shift_cell_symmetric_bc(cellaux, cell)
-                            end if
+                            ! if (bc_x%beg == -2 .or. bc_x%end == -2 .or. bc_y%beg == -2 .or. bc_y%end == -2 &
+                            !     .or. bc_z%beg == -2 .or. bc_z%end == -2) then
+                            !     call s_shift_cell_symmetric_bc(cellaux, cell)
+                            ! end if
                         else
                             func = 0._wp
                             func2 = 0._wp

From d88de9f2ffed9cc0678d24c03d122436dc148138 Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-01.cm.cluster>
Date: Thu, 10 Apr 2025 11:50:11 -0400
Subject: [PATCH 22/23] printing in s_gaussian

---
 src/simulation/m_bubbles_EL_kernels.fpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index 849ae7a57e..09446e01ae 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -131,6 +131,7 @@ contains
 
                         nodecoord(1:3) = 0
                         center(1:3) = 0._wp
+                        print *, 'Starting loop'
                         volpart = 4._wp/3._wp*pi*lbk_rad(l, 2)**3._wp
                         s_coord(1:3) = lbk_s(l, 1:3, 2)
                         center(1:2) = lbk_pos(l, 1:2, 2)

From a7c09308f7aa24530b1b18a7c1e796dbedd6e2ad Mon Sep 17 00:00:00 2001
From: Diego Vaca <dgvacarevelo@login-03.turing.wpi.edu>
Date: Fri, 11 Apr 2025 12:31:33 -0400
Subject: [PATCH 23/23] remove all stops in kernels, do not relocate

---
 src/simulation/m_bubbles_EL.fpp         | 98 +++++++++++--------------
 src/simulation/m_bubbles_EL_kernels.fpp | 89 +++++++++++-----------
 src/simulation/m_checker.fpp            |  2 -
 src/simulation/m_start_up.fpp           | 10 +--
 src/simulation/m_time_steppers.fpp      | 19 -----
 5 files changed, 90 insertions(+), 128 deletions(-)

diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 0ccc2ca69e..03b00d27c9 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -338,16 +338,6 @@ contains
         cell = -buff_size
         call s_locate_cell(mtn_pos(bub_id, 1:3, 1), cell, mtn_s(bub_id, 1:3, 1))
 
-        ! Check if the bubble is located in the ghost cell of a symmetric boundary
-        if (bc_x%beg == -2 .and. cell(1) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_x%beg)."
-        if (bc_x%end == -2 .and. cell(1) > m) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_x%end)."
-        if (bc_y%beg == -2 .and. cell(2) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_y%beg)."
-        if (bc_y%end == -2 .and. cell(2) > n) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_y%end)."
-        if (p > 0) then
-            if (bc_z%beg == -2 .and. cell(3) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_z%beg)."
-            if (bc_z%end == -2 .and. cell(3) > p) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_z%end)."
-        end if
-
         ! If particle is in the ghost cells, find the closest non-ghost cell
         cell(1) = min(max(cell(1), 0), m)
         cell(2) = min(max(cell(2), 0), n)
@@ -518,10 +508,10 @@ contains
         integer, intent(in) :: stage
 
         real(wp) :: vaporflux, pliqint
-        real(wp) :: preterm1, term2, pint, Romega, term1_fac, Rb
+        real(wp) :: preterm1, term2, paux, pint, Romega, term1_fac, Rb
         real(wp) :: conc_v, R_m, gamma_m, fpb, fmass_n, fmass_v
         real(wp) :: fR, fV, fbeta_c, fbeta_t, fR0, fpbdt
-        real(wp) :: pinf, velint, cson, rhol
+        real(wp) :: pinf, aux1, aux2, velint, cson, rhol
         real(wp) :: gamma, pi_inf, qv
         real(wp), dimension(contxe) :: myalpha_rho, myalpha
         real(wp), dimension(2) :: Re
@@ -533,14 +523,12 @@ contains
 
         !< BUBBLE DYNAMICS
 
-        print *, 'bub_dphidt kernel loop'
         ! Subgrid p_inf model based on Maeda and Colonius (2018).
         if (lag_params%pressure_corrector) then
             ! Calculate velocity potentials (valid for one bubble per cell)
             !$acc parallel loop gang vector default(present) private(k, cell)
             do k = 1, nBubs
-                pinf = f_pinf(k, q_prim_vf, 2, cell, preterm1, term2, Romega)
-                print *, 'f_pinf'
+                call s_get_pinf(k, q_prim_vf, 2, paux, cell, preterm1, term2, Romega)
                 fR0 = bub_R0(k)
                 fR = intfc_rad(k, 2)
                 fV = intfc_vel(k, 2)
@@ -548,16 +536,16 @@ contains
                 pint = f_cpbw_KM(fR0, fR, fV, fpb)
                 pint = pint + 0.5_wp*fV**2._wp
                 if (lag_params%cluster_type == 2) then
-                    bub_dphidt(k) = (pinf - pint) + term2
+                    bub_dphidt(k) = (paux - pint) + term2
                     ! Accounting for the potential induced by the bubble averaged over the control volume
                     ! Note that this is based on the incompressible flow assumption near the bubble.
-                    term1_fac = 3._wp/2._wp*(fR*(Romega**2._wp - fR**2._wp))/(Romega**3._wp - fR**3._wp)
+                    Rb = intfc_rad(k, 2)
+                    term1_fac = 3._wp/2._wp*(Rb*(Romega**2._wp - Rb**2._wp))/(Romega**3._wp - Rb**3._wp)
                     bub_dphidt(k) = bub_dphidt(k)/(1._wp - term1_fac)
                 end if
             end do
         end if
 
-        print *, 'bub dynamics kernel loop'
         ! Radial motion model
         if (bubble_model == 2) then
             !$acc parallel loop gang vector default(present) private(k, myalpha_rho, myalpha, Re, cell) copyin(stage)
@@ -583,15 +571,14 @@ contains
                 ! Vapor and heat fluxes
                 vaporflux = f_vflux(fR, fV, fmass_v, k, fmass_n, fbeta_c, conc_v)
                 fpbdt = f_bpres_dot(vaporflux, fR, fV, fpb, fmass_v, k, fbeta_t, R_m, gamma_m, conc_v)
-                print *, 'f_vflux and f_bpres_dot'
                 gas_dmvdt(k, stage) = 4._wp*pi*fR**2._wp*vaporflux
-                print *, 'updated gas_dmvdt', stage
+
                 ! Pressure at the bubble wall
                 pliqint = f_cpbw_KM(fR0, fR, fV, fpb)
-                print *, 'f_cpbw_KM'
+
                 ! Obtaining driving pressure
-                pinf = f_pinf(k, q_prim_vf, 1, cell, preterm1, term2, Romega)
-                print *, 'f_pinf'
+                call s_get_pinf(k, q_prim_vf, 1, pinf, cell, aux1, aux2)
+
                 ! Obtain liquid density and computing speed of sound from pinf
                 !$acc loop seq
                 do i = 1, contxe
@@ -602,15 +589,14 @@ contains
                 call s_convert_species_to_mixture_variables_acc(rhol, gamma, pi_inf, qv, myalpha, &
                                                                 myalpha_rho, Re, cell(1), cell(2), cell(3))
                 call s_compute_cson_from_pinf(k, q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
-                print *, 'cson and rhol'
+
                 ! Velocity correction due to massflux
                 velint = fV - gas_dmvdt(k, stage)/(4._wp*pi*fR**2._wp*rhol)
-                print *, 'velint'
+
                 ! Interphase acceleration and update vars
                 intfc_dveldt(k, stage) = f_rddot_KM(fpbdt, pinf, pliqint, rhol, fR, velint, fR0, cson)
                 gas_dpdt(k, stage) = fpbdt
                 intfc_draddt(k, stage) = fV
-                print *, 'update intfc_dveldt, gas_dpdt, intfc_draddt'
 
             end do
         else
@@ -623,7 +609,6 @@ contains
             end do
         end if
 
-        print *, 'bub position kernel loop'
         ! Bubbles remain in a fixed position
         !$acc parallel loop collapse(2) gang vector default(present) private(k) copyin(stage)
         do k = 1, nBubs
@@ -636,9 +621,7 @@ contains
         call nvtxEndRange
 
         !< EULER-LAGRANGE COUPLING
-        print *, 's_smear_voidfraction'
         call s_smear_voidfraction()
-        print *, 's_add_rhs_sources'
         if (lag_params%solver_approach == 2) call s_add_rhs_sources(q_cons_vf, q_prim_vf, rhs_vf)
 
     end subroutine s_compute_EL_coupled_solver
@@ -685,7 +668,6 @@ contains
 
         call nvtxStartRange("BUBBLES-LAGRANGE-KERNELS")
 
-        print *, 'Zeroing vars'
         !$acc parallel loop collapse(4) gang vector default(present)
         do i = 1, q_beta_idx
             do l = idwbuff(3)%beg, idwbuff(3)%end
@@ -697,11 +679,9 @@ contains
             end do
         end do
 
-        print *, 's_smoothfunction'
         call s_smoothfunction(nBubs, intfc_rad, intfc_vel, &
                               mtn_s, mtn_pos, q_beta)
 
-        print *, '1-beta'
         !Store 1-beta
         !$acc parallel loop collapse(3) gang vector default(present)
         do l = idwbuff(3)%beg, idwbuff(3)%end
@@ -814,16 +794,20 @@ contains
         !! @param bub_id Particle identifier
         !! @param q_prim_vf  Primitive variables
         !! @param ptype 1: p at infinity, 2: averaged P at the bubble location
+        !! @param f_pinfl Driving pressure
         !! @param cell Bubble cell
         !! @param Romega Control volume radius
-    function f_pinf(bub_id, q_prim_vf, ptype, cell, preterm1, term2, Romega)
+    subroutine s_get_pinf(bub_id, q_prim_vf, ptype, f_pinfl, cell, preterm1, term2, Romega)
+#ifdef _CRAYFTN
+        !DIR$ INLINEALWAYS s_get_pinf
+#else
         !$acc routine seq
+#endif
         integer, intent(in) :: bub_id, ptype
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+        real(wp), intent(out) :: f_pinfl
         integer, dimension(3), intent(out) :: cell
-        real(wp), intent(out) :: preterm1, term2, Romega
-
-        real(wp) :: f_pinf
+        real(wp), intent(out), optional :: preterm1, term2, Romega
 
         real(wp), dimension(3) :: scoord, psi
         real(wp) :: dc, vol, aux
@@ -831,11 +815,11 @@ contains
         real(wp) :: charvol, charpres, charvol2, charpres2
         integer, dimension(3) :: cellaux
         integer :: i, j, k
-        integer :: smearGrid, smearGridz
+        integer :: mapCells_pinf, smearGrid, smearGridz
         logical :: celloutside
 
         scoord = mtn_s(bub_id, 1:3, 2)
-        f_pinf = 0._wp
+        f_pinfl = 0._wp
 
         !< Find current bubble cell
         cell(:) = int(scoord(:))
@@ -898,19 +882,19 @@ contains
 
             !< Perform bilinear interpolation
             if (p == 0) then  !2D
-                f_pinf = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)
+                f_pinfl = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)
             else              !3D
-                f_pinf = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))*(1._wp - psi(3))
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))*(1._wp - psi(3))
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)*(1._wp - psi(3))
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)*(1._wp - psi(3))
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1)*(1._wp - psi(1))*(1._wp - psi(2))*psi(3)
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3) + 1)*psi(1)*(1._wp - psi(2))*psi(3)
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3) + 1)*psi(1)*psi(2)*psi(3)
-                f_pinf = f_pinf + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3) + 1)*(1._wp - psi(1))*psi(2)*psi(3)
+                f_pinfl = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))*(1._wp - psi(3))
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))*(1._wp - psi(3))
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)*(1._wp - psi(3))
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)*(1._wp - psi(3))
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1)*(1._wp - psi(1))*(1._wp - psi(2))*psi(3)
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3) + 1)*psi(1)*(1._wp - psi(2))*psi(3)
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3) + 1)*psi(1)*psi(2)*psi(3)
+                f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3) + 1)*(1._wp - psi(1))*psi(2)*psi(3)
             end if
 
             !R_Omega
@@ -918,8 +902,12 @@ contains
 
         else if (lag_params%cluster_type >= 2) then
             ! Bubble dynamic closure from Maeda and Colonius (2018)
+            if (lag_params%smooth_type == 1) then
+                mapCells_pinf = mapCells
+            end if
+
             ! Include the cell that contains the bubble (mapCells+1+mapCells)
-            smearGrid = 2*mapCells + 1
+            smearGrid = mapCells_pinf - (-mapCells_pinf) + 1
             smearGridz = smearGrid
             if (p == 0) smearGridz = 1
 
@@ -990,7 +978,7 @@ contains
                 end do
             end do
 
-            f_pinf = charpres2/charvol2
+            f_pinfl = charpres2/charvol2
             vol = charvol
             dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
 
@@ -1011,16 +999,16 @@ contains
             preterm1 = 3._wp/2._wp*Rbeq*(dc**2._wp - Rbeq**2._wp)/(aux*denom)
 
             !Control volume radius
-            Romega = dc
+            if (ptype == 2) Romega = dc
 
             ! Getting p_inf
             if (ptype == 1) then
-                f_pinf = f_pinf + preterm1*term1 + term2
+                f_pinfl = f_pinfl + preterm1*term1 + term2
             end if
 
         end if
 
-    end function f_pinf
+    end subroutine s_get_pinf
 
     !>  This subroutine updates the Lagrange variables using the tvd RK time steppers.
         !!      The time derivative of the bubble variables must be stored at every stage to avoid precision errors.
diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
index 09446e01ae..7bab026e67 100644
--- a/src/simulation/m_bubbles_EL_kernels.fpp
+++ b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -116,42 +116,36 @@ contains
         logical :: celloutside
         integer :: smearGrid, smearGridz
 
-        smearGrid = mapCells - (-mapCells) ! Include the cell that contains the bubble (3+1+3)
+        smearGrid = mapCells - (-mapCells) + 1 ! Include the cell that contains the bubble (3+1+3)
         smearGridz = smearGrid
-        if (p == 0) smearGridz = 0
+        if (p == 0) smearGridz = 1
 
-        print *, 's_gaussian'
-
-        !$acc parallel loop collapse(4) gang vector default(present) private(l, s_coord, cell, center, cellaux, nodecoord) &
-        !$acc copyin(smearGrid, smearGridz)
+        !$acc parallel loop gang vector default(present) private(nodecoord, l, s_coord, cell, center) copyin(smearGrid, smearGridz)
         do l = 1, nBubs
-            do i = 0, smearGrid
-                do j = 0, smearGrid
-                    do k = 0, smearGridz
-
-                        nodecoord(1:3) = 0
-                        center(1:3) = 0._wp
-                        print *, 'Starting loop'
-                        volpart = 4._wp/3._wp*pi*lbk_rad(l, 2)**3._wp
-                        s_coord(1:3) = lbk_s(l, 1:3, 2)
-                        center(1:2) = lbk_pos(l, 1:2, 2)
-                        print *, 'reading initial state'
-                        if (p > 0) center(3) = lbk_pos(l, 3, 2)
-                        call s_get_cell(s_coord, cell)
-                        call s_compute_stddsv(cell, volpart, stddsv)
-                        strength_vol = volpart
-                        strength_vel = 4._wp*pi*lbk_rad(l, 2)**2._wp*lbk_vel(l, 2)
-                        print *, 's_compute_stddsv'
-
-                        cellaux(1) = cell(1) + i - mapCells
-                        cellaux(2) = cell(2) + j - mapCells
-                        cellaux(3) = cell(3) + k - mapCells
+            nodecoord(1:3) = 0
+            center(1:3) = 0._wp
+            volpart = 4._wp/3._wp*pi*lbk_rad(l, 2)**3._wp
+            s_coord(1:3) = lbk_s(l, 1:3, 2)
+            center(1:2) = lbk_pos(l, 1:2, 2)
+            if (p > 0) center(3) = lbk_pos(l, 3, 2)
+            call s_get_cell(s_coord, cell)
+            call s_compute_stddsv(cell, volpart, stddsv)
+
+            strength_vol = volpart
+            strength_vel = 4._wp*pi*lbk_rad(l, 2)**2._wp*lbk_vel(l, 2)
+
+            !$acc loop collapse(3) private(cellaux, nodecoord)
+            do i = 1, smearGrid
+                do j = 1, smearGrid
+                    do k = 1, smearGridz
+                        cellaux(1) = cell(1) + i - (mapCells + 1)
+                        cellaux(2) = cell(2) + j - (mapCells + 1)
+                        cellaux(3) = cell(3) + k - (mapCells + 1)
                         if (p == 0) cellaux(3) = 0
 
                         !Check if the cells intended to smear the bubbles in are in the computational domain
                         !and redefine the cells for symmetric boundary
                         call s_check_celloutside(cellaux, celloutside)
-                        print *, 's_check_celloutside'
 
                         if (.not. celloutside) then
 
@@ -159,13 +153,13 @@ contains
                             nodecoord(2) = y_cc(cellaux(2))
                             if (p > 0) nodecoord(3) = z_cc(cellaux(3))
                             call s_applygaussian(center, cellaux, nodecoord, stddsv, 0._wp, func)
-                            if (p == 0) call s_applygaussian(center, cellaux, nodecoord, stddsv, 1._wp, func2)
+                            if (lag_params%cluster_type >= 4) call s_applygaussian(center, cellaux, nodecoord, stddsv, 1._wp, func2)
 
                             ! Relocate cells for bubbles intersecting symmetric boundaries
-                            ! if (bc_x%beg == -2 .or. bc_x%end == -2 .or. bc_y%beg == -2 .or. bc_y%end == -2 &
-                            !     .or. bc_z%beg == -2 .or. bc_z%end == -2) then
-                            !     call s_shift_cell_symmetric_bc(cellaux, cell)
-                            ! end if
+                            if (bc_x%beg == -2 .or. bc_x%end == -2 .or. bc_y%beg == -2 .or. bc_y%end == -2 &
+                                .or. bc_z%beg == -2 .or. bc_z%end == -2) then
+                                call s_shift_cell_symmetric_bc(cellaux, cell)
+                            end if
                         else
                             func = 0._wp
                             func2 = 0._wp
@@ -191,16 +185,13 @@ contains
 
                         !Product of two smeared functions
                         !Update void fraction * time derivative of void fraction
-                        if (p == 0) then
+                        if (lag_params%cluster_type >= 4) then
                             addFun3 = func2*strength_vol*strength_vel
                             !$acc atomic update
                             updatedvar%vf(5)%sf(cellaux(1), cellaux(2), cellaux(3)) = &
                                 updatedvar%vf(5)%sf(cellaux(1), cellaux(2), cellaux(3)) &
                                 + addFun3
                         end if
-
-                        print *, 'Update 1-2-3'
-
                     end do
                 end do
             end do
@@ -326,27 +317,39 @@ contains
 
         ! x-dir
         if (bc_x%beg == -2 .and. (cell(1) <= mapCells - 1)) then
-            cellaux(1) = abs(cellaux(1)) - 1
+            if (cell(1) >= 0) then
+                cellaux(1) = abs(cellaux(1)) - 1
+            end if
         end if
         if (bc_x%end == -2 .and. (cell(1) >= m + 1 - mapCells)) then
-            cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
+            if (cell(1) <= m) then
+                cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
+            end if
         end if
 
         !y-dir
         if (bc_y%beg == -2 .and. (cell(2) <= mapCells - 1)) then
-            cellaux(2) = abs(cellaux(2)) - 1
+            if (cell(2) >= 0) then
+                cellaux(2) = abs(cellaux(2)) - 1
+            end if
         end if
         if (bc_y%end == -2 .and. (cell(2) >= n + 1 - mapCells)) then
-            cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
+            if (cell(2) <= n) then
+                cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
+            end if
         end if
 
         if (p > 0) then
             !z-dir
             if (bc_z%beg == -2 .and. (cell(3) <= mapCells - 1)) then
-                cellaux(3) = abs(cellaux(3)) - 1
+                if (cell(3) >= 0) then
+                    cellaux(3) = abs(cellaux(3)) - 1
+                end if
             end if
             if (bc_z%end == -2 .and. (cell(3) >= p + 1 - mapCells)) then
-                cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
+                if (cell(3) <= p) then
+                    cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
+                end if
             end if
         end if
 
diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
index e5e00ca934..a7f6665df4 100644
--- a/src/simulation/m_checker.fpp
+++ b/src/simulation/m_checker.fpp
@@ -330,8 +330,6 @@ contains
         @:PROHIBIT(bubbles_lagrange .and. file_per_process, "file_per_process must be false for bubbles_lagrange")
         @:PROHIBIT(bubbles_lagrange .and. n==0, "bubbles_lagrange accepts 2D and 3D simulations only")
         @:PROHIBIT(bubbles_lagrange .and. model_eqns==3, "The 6-equation flow model does not support bubbles_lagrange")
-        @:PROHIBIT(bubbles_lagrange .and. lag_params%cluster_type<=0, "cluster_type must be specified")
-        @:PROHIBIT(bubbles_lagrange .and. lag_params%cluster_type>=2 .and. lag_params%smooth_type/=1, "cluster_type=2 requires smooth_type=1")
     end subroutine s_check_inputs_bubbles_lagrange
 
     !> Checks miscellaneous constraints,
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index 9abf9ec205..39d3271c3c 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -1395,17 +1395,9 @@ contains
         end if
 
         if (bubbles_lagrange) then
-            !$acc update host(intfc_rad)
-            do i = 1, nBubs
-                if (ieee_is_nan(intfc_rad(i, 1)) .or. intfc_rad(i, 1) <= 0._wp) then
-                    print *, "Bubble radius is negative or NaN", proc_rank, t_step, i, intfc_rad(i, 1)
-                    error stop "Bubble radius is negative or NaN, please reduce dt"
-                end if
-            end do
-
             !$acc update host(q_beta%vf(1)%sf)
             call s_write_data_files(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, save_count, q_beta%vf(1))
-            !$acc update host(Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_vel)
+            !$acc update host(Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_rad, intfc_vel)
             call s_write_restart_lag_bubbles(save_count) !parallel
             if (lag_params%write_bubbles_stats) call s_write_lag_bubble_stats()
         else
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index 2e9a2f2a6c..2ffbc8e58b 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -1075,16 +1075,11 @@ contains
 #ifdef DEBUG
             if (proc_rank == 0) print *, 'RKCK 1st time-stage at', rkck_time_tmp
 #endif
-            print *, 's_compute_rhs', RKstep
             call s_compute_rhs(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, rhs_ts_rkck(1)%vf, pb_ts(1)%sf, rhs_pb, mv_ts(1)%sf, rhs_mv, t_step, time_avg)
-            print *, 's_compute_EL_coupled_solver', RKstep
             call s_compute_EL_coupled_solver(q_cons_ts(1)%vf, q_prim_vf, rhs_ts_rkck(1)%vf, RKstep)
-            print *, 's_update_tmp_rkck', RKstep
             call s_update_tmp_rkck(RKstep, q_cons_ts, rhs_ts_rkck, lag_largestep)
-            print *, 's_compute_rkck_dt', RKstep, lag_largestep
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
-            print *, 'done with stage 1'
 
             ! SECOND TIME-STAGE
             RKstep = 2
@@ -1094,16 +1089,11 @@ contains
 #ifdef DEBUG
             if (proc_rank == 0) print *, 'RKCK 2nd time-stage at', rkck_time_tmp
 #endif
-            print *, 's_compute_rhs', RKstep
             call s_compute_rhs(q_cons_ts(2)%vf, q_T_sf, q_prim_vf, rhs_ts_rkck(2)%vf, pb_ts(1)%sf, rhs_pb, mv_ts(1)%sf, rhs_mv, t_step, time_avg)
-            print *, 's_compute_EL_coupled_solver', RKstep
             call s_compute_EL_coupled_solver(q_cons_ts(2)%vf, q_prim_vf, rhs_ts_rkck(2)%vf, RKstep)
-            print *, 's_update_tmp_rkck', RKstep
             call s_update_tmp_rkck(RKstep, q_cons_ts, rhs_ts_rkck, lag_largestep)
-            print *, 's_compute_rkck_dt', RKstep, lag_largestep
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
-            print *, 'done with stage 2'
 
             ! THIRD TIME-STAGE
             RKstep = 3
@@ -1118,7 +1108,6 @@ contains
             call s_update_tmp_rkck(RKstep, q_cons_ts, rhs_ts_rkck, lag_largestep)
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
-            print *, 'done with stage 3'
 
             ! FOURTH TIME-STAGE
             RKstep = 4
@@ -1133,7 +1122,6 @@ contains
             call s_update_tmp_rkck(RKstep, q_cons_ts, rhs_ts_rkck, lag_largestep)
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
-            print *, 'done with stage 4'
 
             ! FIFTH TIME-STAGE
             RKstep = 5
@@ -1148,7 +1136,6 @@ contains
             call s_update_tmp_rkck(5, q_cons_ts, rhs_ts_rkck, lag_largestep)
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
-            print *, 'done with stage 5'
 
             ! SIXTH TIME-STAGE
             RKstep = 6
@@ -1163,7 +1150,6 @@ contains
             call s_update_tmp_rkck(6, q_cons_ts, rhs_ts_rkck, lag_largestep)
             if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step)
             if (restart_rkck_step) cycle
-            print *, 'done with stage 6'
 
             dt_did = dt
 
@@ -1172,12 +1158,9 @@ contains
 #ifdef DEBUG
                 if (proc_rank == 0) print *, 'Computing truncation error (4th/5th RKCK)'
 #endif
-                print *, 's_calculate_rkck_truncation_error'
                 call s_calculate_rkck_truncation_error(rkck_errmax)
-                print *, 's_compute_rkck_dt'
                 call s_compute_rkck_dt(lag_largestep, restart_rkck_step, rkck_errmax)
                 if (restart_rkck_step) cycle
-                print *, 'done with truncation error', dt
             end if
 
         end do
@@ -1185,7 +1168,6 @@ contains
         !> Update values
         mytime = mytime + dt_did
         call s_update_rkck(q_cons_ts)
-        print *, 'done with s_update_rkck'
 
         call s_write_void_evol(mytime)
         if (lag_params%write_bubbles_stats) call s_calculate_lag_bubble_stats()
@@ -1199,7 +1181,6 @@ contains
             call s_write_run_time_information(q_prim_vf, t_step)
         end if
 
-        print *, 'done with s_4th_5th_order_rkck'
     end subroutine s_4th_5th_order_rkck
 
     !> Module deallocation and/or disassociation procedures