Fix #797

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,17 +376,16 @@ 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
 
         ! Bubble natural frequency
         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 +915,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 +992,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 +1040,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 +1062,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 +1074,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 +1098,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 +1110,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 +1121,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()

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
 

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(lag_params%cluster_type<=0, "cluster_type must be specified")
+        @: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,

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