Addressing comments 3 and fixing Frontier

Diego Vaca · Diego Vaca · commit 81121a7282b4 · 2025-05-16T23:15:51.000-05:00
diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp
@@ -457,6 +457,7 @@ contains
         !!  @param fbeta_c Mass transfer coefficient (EL)
         !!  @param fbeta_t Heat transfer coefficient (EL)
         !!  @param fCson Speed of sound (EL)
+        !!  @param adap_dt_stop Fail-safe exit if max iteration count reached
     subroutine s_advance_step(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
                               fntait, fBtait, f_bub_adv_src, f_divu, &
                               bub_id, fmass_v, fmass_n, fbeta_c, &
@@ -471,17 +472,18 @@ contains
         real(wp), intent(in) :: fntait, fBtait, f_bub_adv_src, f_divu
         integer, intent(in) :: bub_id
         real(wp), intent(in) :: fmass_n, fbeta_c, fbeta_t, fCson
-        integer, intent(out) :: adap_dt_stop
+        integer, intent(inout) :: adap_dt_stop
 
         real(wp), dimension(5) :: err !< Error estimates for adaptive time stepping
         real(wp) :: t_new !< Updated time step size
         real(wp) :: h !< Time step size
         real(wp), dimension(4) :: myR_tmp1, myV_tmp1, myR_tmp2, myV_tmp2 !< Bubble radius, radial velocity, and radial acceleration for the inner loop
         real(wp), dimension(4) :: myPb_tmp1, myMv_tmp1, myPb_tmp2, myMv_tmp2 !< Gas pressure and vapor mass for the inner loop (EL)
+        real(wp) :: fR2, fV2, fpb2, fmass_v2
         integer :: iter_count
 
-        h = f_initial_substep_h(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
-                                fntait, fBtait, f_bub_adv_src, f_divu, fCson)
+        call s_initial_substep_h(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
+                                 fntait, fBtait, f_bub_adv_src, f_divu, fCson, h)
 
         ! Advancing one step
         t_new = 0._wp
@@ -500,27 +502,30 @@ contains
                 iter_count = iter_count + 1
 
                 ! Advance one sub-step
-                err(1) = f_advance_substep( &
-                         fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
-                         fntait, fBtait, f_bub_adv_src, f_divu, &
-                         bub_id, fmass_v, fmass_n, fbeta_c, &
-                         fbeta_t, fCson, h, &
-                         myR_tmp1, myV_tmp1, myPb_tmp1, myMv_tmp1)
+                call s_advance_substep(err(1), &
+                                       fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
+                                       fntait, fBtait, f_bub_adv_src, f_divu, &
+                                       bub_id, fmass_v, fmass_n, fbeta_c, &
+                                       fbeta_t, fCson, h, &
+                                       myR_tmp1, myV_tmp1, myPb_tmp1, myMv_tmp1)
 
                 ! Advance one sub-step by advancing two half steps
-                err(2) = f_advance_substep( &
-                         fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
-                         fntait, fBtait, f_bub_adv_src, f_divu, &
-                         bub_id, fmass_v, fmass_n, fbeta_c, &
-                         fbeta_t, fCson, 0.5_wp*h, &
-                         myR_tmp2, myV_tmp2, myPb_tmp2, myMv_tmp2)
-
-                err(3) = f_advance_substep( &
-                         fRho, fP, myR_tmp2(4), myV_tmp2(4), fR0, myPb_tmp2(4), fpbdot, alf, &
-                         fntait, fBtait, f_bub_adv_src, f_divu, &
-                         bub_id, myMv_tmp2(4), fmass_n, fbeta_c, &
-                         fbeta_t, fCson, 0.5_wp*h, &
-                         myR_tmp2, myV_tmp2, myPb_tmp2, myMv_tmp2)
+                call s_advance_substep(err(2), &
+                                       fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
+                                       fntait, fBtait, f_bub_adv_src, f_divu, &
+                                       bub_id, fmass_v, fmass_n, fbeta_c, &
+                                       fbeta_t, fCson, 0.5_wp*h, &
+                                       myR_tmp2, myV_tmp2, myPb_tmp2, myMv_tmp2)
+
+                fR2 = myR_tmp2(4); fV2 = myV_tmp2(4)
+                fpb2 = myPb_tmp2(4); fmass_v2 = myMv_tmp2(4)
+
+                call s_advance_substep(err(3), &
+                                       fRho, fP, fR2, fV2, fR0, fpb2, fpbdot, alf, &
+                                       fntait, fBtait, f_bub_adv_src, f_divu, &
+                                       bub_id, fmass_v2, fmass_n, fbeta_c, &
+                                       fbeta_t, fCson, 0.5_wp*h, &
+                                       myR_tmp2, myV_tmp2, myPb_tmp2, myMv_tmp2)
 
                 err(4) = abs((myR_tmp1(4) - myR_tmp2(4))/myR_tmp1(4))
                 err(5) = abs((myV_tmp1(4) - myV_tmp2(4))/myV_tmp1(4))
@@ -587,15 +592,20 @@ contains
         !!  @param f_bub_adv_src Source for bubble volume fraction
         !!  @param f_divu Divergence of velocity
         !!  @param fCson Speed of sound (EL)
-    function f_initial_substep_h(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
-                                 fntait, fBtait, f_bub_adv_src, f_divu, &
-                                 fCson)
+        !!  @param h Time step size
+    subroutine s_initial_substep_h(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
+                                   fntait, fBtait, f_bub_adv_src, f_divu, &
+                                   fCson, h)
+#ifdef _CRAYFTN
+        !DIR$ INLINEALWAYS s_initial_substep_h
+#else
         !$acc routine seq
+#endif
         real(wp), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
         real(wp), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu
         real(wp), intent(IN) :: fCson
+        real(wp), intent(OUT) :: h
 
-        real(wp) :: f_initial_substep_h
         real(wp), dimension(2) :: h_size !< Time step size (h0, h1)
         real(wp), dimension(3) :: d_norms !< norms (d_0, d_1, d_2)
         real(wp), dimension(2) :: myR_tmp, myV_tmp, myA_tmp !< Bubble radius, radial velocity, and radial acceleration
@@ -637,12 +647,13 @@ contains
             h_size(2) = (scale_guess/max(d_norms(2), d_norms(3)))**(1._wp/3._wp)
         end if
 
-        f_initial_substep_h = min(h_size(1)/scale_guess, h_size(2))
+        h = min(h_size(1)/scale_guess, h_size(2))
 
-    end function f_initial_substep_h
+    end subroutine s_initial_substep_h
 
     !>  Integrate bubble variables over the given time step size, h, using a
         !!      third-order accurate embedded Runge–Kutta scheme.
+        !!  @param err Estimated error
         !!  @param fRho Current density
         !!  @param fP Current driving pressure
         !!  @param fR Current bubble radius
@@ -666,20 +677,24 @@ contains
         !!  @param myV_tmp Bubble radial velocity at each stage
         !!  @param myPb_tmp Internal bubble pressure at each stage (EL)
         !!  @param myMv_tmp Mass of vapor in the bubble at each stage (EL)
-    function f_advance_substep(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
-                               fntait, fBtait, f_bub_adv_src, f_divu, &
-                               bub_id, fmass_v, fmass_n, fbeta_c, &
-                               fbeta_t, fCson, h, &
-                               myR_tmp, myV_tmp, myPb_tmp, myMv_tmp)
+    subroutine s_advance_substep(err, fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
+                                 fntait, fBtait, f_bub_adv_src, f_divu, &
+                                 bub_id, fmass_v, fmass_n, fbeta_c, &
+                                 fbeta_t, fCson, h, &
+                                 myR_tmp, myV_tmp, myPb_tmp, myMv_tmp)
+#ifdef _CRAYFTN
+        !DIR$ INLINEALWAYS s_advance_substep
+#else
         !$acc routine seq
+#endif
+        real(wp), intent(OUT) :: err
         real(wp), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
         real(wp), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu, h
         integer, intent(IN) :: bub_id
         real(wp), intent(IN) :: fmass_v, fmass_n, fbeta_c, fbeta_t, fCson
         real(wp), dimension(4), intent(OUT) :: myR_tmp, myV_tmp, myPb_tmp, myMv_tmp
 
         real(wp), dimension(4) :: myA_tmp, mydPbdt_tmp, mydMvdt_tmp
-        real(wp) :: f_advance_substep
         real(wp) :: err_R, err_V
 
         myPb_tmp(1:4) = fpb
@@ -751,9 +766,9 @@ contains
             f_approx_equal(myA_tmp(3), 0._wp) .and. f_approx_equal(myA_tmp(4), 0._wp)) then
             err_V = 0._wp
         end if
-        f_advance_substep = sqrt((err_R**2._wp + err_V**2._wp)/2._wp)
+        err = sqrt((err_R**2._wp + err_V**2._wp)/2._wp)
 
-    end function f_advance_substep
+    end subroutine s_advance_substep
 
     !>  Changes of pressure and vapor mass in the lagrange bubbles.
         !!  @param bub_id Bubble identifier
@@ -766,8 +781,8 @@ contains
         !!  @param fMv_tmp Mass of vapor in the bubble
         !!  @param fdPbdt_tmp Rate of change of the internal bubble pressure
         !!  @param fdMvdt_tmp Rate of change of the mass of vapor in the bubble
-    function f_advance_EL(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, fmass_n, fbeta_c, fbeta_t, &
-                          fdPbdt_tmp)
+    function f_advance_EL(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, &
+                          fmass_n, fbeta_c, fbeta_t, fdPbdt_tmp)
         !$acc routine seq
         real(wp), intent(IN) :: fR_tmp, fV_tmp, fPb_tmp, fMv_tmp
         real(wp), intent(IN) :: fmass_n, fbeta_c, fbeta_t
diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp
@@ -168,7 +168,7 @@ contains
 
         integer :: i, j, k, l, q, ii !< Loop variables
 
-        integer :: adap_dt_stop_max, adap_dt_stop !< fail-safe exit if max iteration count reached
+        integer :: adap_dt_stop_max, adap_dt_stop !< Fail-safe exit if max iteration count reached
         integer :: dmBub_id !< Dummy variables for unified subgrid bubble subroutines
         real(wp) :: dmMass_v, dmMass_n, dmBeta_c, dmBeta_t, dmCson
 
@@ -320,7 +320,7 @@ contains
             end do
         end do
 
-        if (adap_dt .and. adap_dt_stop_max > 0) stop "Adaptive time stepping failed to converge"
+        if (adap_dt .and. adap_dt_stop_max > 0) call s_mpi_abort("Adaptive time stepping failed to converge.")
 
         if (.not. adap_dt) then
             !$acc parallel loop collapse(3) gang vector default(present)
diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
@@ -224,7 +224,7 @@ contains
                     indomain = particle_in_domain(inputBubble(1:3))
                     id = id + 1
                     if (id > lag_params%nBubs_glb .and. proc_rank == 0) then
-                        call s_mpi_abort('Current number of bubbles is larger than nBubs_glb')
+                        call s_mpi_abort("Current number of bubbles is larger than nBubs_glb")
                     end if
                     if (indomain) then
                         bub_id = bub_id + 1
@@ -236,7 +236,7 @@ contains
                 end do
                 close (94)
             else
-                stop "Lagrange bubbles: you have to initialize them in input/lag_bubbles.dat"
+                call s_mpi_abort("Initialize the lagrange bubbles in input/lag_bubbles.dat")
             end if
         else
             if (proc_rank == 0) print *, 'Restarting lagrange bubbles at save_count: ', save_count
@@ -318,13 +318,18 @@ contains
         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 == BC_REFLECTIVE .and. cell(1) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_x%beg)."
-        if (bc_x%end == BC_REFLECTIVE .and. cell(1) > m) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_x%end)."
-        if (bc_y%beg == BC_REFLECTIVE .and. cell(2) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_y%beg)."
-        if (bc_y%end == BC_REFLECTIVE .and. cell(2) > n) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_y%end)."
+        if ((bc_x%beg == BC_REFLECTIVE .and. cell(1) < 0) .or. &
+            (bc_x%end == BC_REFLECTIVE .and. cell(1) > m) .or. &
+            (bc_y%beg == BC_REFLECTIVE .and. cell(2) < 0) .or. &
+            (bc_y%end == BC_REFLECTIVE .and. cell(2) > n)) then
+            call s_mpi_abort("Lagrange bubble is in the ghost cells of a symmetric boundary.")
+        end if
+
         if (p > 0) then
-            if (bc_z%beg == BC_REFLECTIVE .and. cell(3) < 0) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_z%beg)."
-            if (bc_z%end == BC_REFLECTIVE .and. cell(3) > p) stop "Lagrange bubble is in the ghost cells of a symmetric boundary (bc_z%end)."
+            if ((bc_z%beg == BC_REFLECTIVE .and. cell(3) < 0) .or. &
+                (bc_z%end == BC_REFLECTIVE .and. cell(3) > p)) then
+                call s_mpi_abort("Lagrange bubble is in the ghost cells of a symmetric boundary.")
+            end if
         end if
 
         ! If particle is in the ghost cells, find the closest non-ghost cell
@@ -355,15 +360,15 @@ 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'
+            call s_mpi_abort("The initial mass of gas inside the bubble is negative. Check the 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
-            stop "Lagrange bubble initially located in a region with pressure below the vapor pressure."
+            call s_mpi_abort("Lagrange bubble initially located in a region with pressure below the vapor pressure.")
         end if
         omegaN = sqrt(omegaN/bub_R0(bub_id)**2._wp)
 
@@ -379,7 +384,9 @@ contains
         call s_transcoeff(1._wp, PeG, Re_trans, Im_trans)
         gas_betaC(bub_id) = Re_trans*(massflag)*lag_params%diffcoefvap
 
-        if (gas_mg(bub_id) <= 0._wp) stop "Negative gas mass in the bubble, check if the bubble is in the domain."
+        if (gas_mg(bub_id) <= 0._wp) then
+            call s_mpi_abort("Negative gas mass in the bubble, check if the bubble is in the domain.")
+        end if
 
     end subroutine s_add_bubbles
 
@@ -507,7 +514,7 @@ contains
         real(wp), dimension(2) :: Re
         integer, dimension(3) :: cell
 
-        integer :: adap_dt_stop_max, adap_dt_stop !< fail-safe exit if max iteration count reached
+        integer :: adap_dt_stop_max, adap_dt_stop !< Fail-safe exit if max iteration count reached
         real(wp) :: dmalf, dmntait, dmBtait, dm_bub_adv_src, dm_divu !< Dummy variables for unified subgrid bubble subroutines
 
         integer :: i, k, l
@@ -604,7 +611,7 @@ contains
 
         end do
 
-        if (adap_dt .and. adap_dt_stop_max > 0) stop "Adaptive time stepping failed to converge"
+        if (adap_dt .and. adap_dt_stop_max > 0) call s_mpi_abort("Adaptive time stepping failed to converge.")
 
         ! Bubbles remain in a fixed position
         !$acc parallel loop collapse(2) gang vector default(present) private(k) copyin(stage)
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
@@ -1431,8 +1431,7 @@ contains
             !$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"
+                    call s_mpi_abort("Bubble radius is negative or NaN, please reduce dt.")
                 end if
             end do
 
