addressing comments

Author: Diego Vaca

src/common/m_constants.fpp

@@ -49,6 +49,10 @@ module m_constants
     integer, parameter :: mapCells = 3 !< Number of cells around the bubble where the smoothening function will have effect
     real(wp), parameter :: R_uni = 8314._wp ! Universal gas constant - J/kmol/K
 
+    ! Strang Splitting constants
+    real(wp), parameter :: dflt_adap_dt_tol = 1e-4_wp !< Default tolerance for adaptive step size
+    integer, parameter :: adap_dt_max_iters = 100 !< Maximum number of iterations
+
     ! Relativity
     integer, parameter :: relativity_cons_to_prim_max_iter = 100
 

src/simulation/m_bubbles.fpp

@@ -65,6 +65,8 @@ contains
         @:ALLOCATE(bub_p_src(0:m, 0:n, 0:p, 1:nb))
         @:ALLOCATE(bub_m_src(0:m, 0:n, 0:p, 1:nb))
 
+        if (adap_dt .and. f_is_default(adap_dt_tol)) adap_dt_tol = dflt_adap_dt_tol
+
     end subroutine s_initialize_bubbles_EE_module
 
     ! Compute the bubble volume fraction alpha from the bubble number density n
@@ -166,6 +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 :: dmBub_id !< Dummy variables for unified subgrid bubble subroutines
         real(wp) :: dmMass_v, dmMass_n, dmBeta_c, dmBeta_t, dmCson
 
@@ -186,7 +189,9 @@ contains
             end do
         end do
 
-        !$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp, myalpha_rho, myalpha)
+        adap_dt_stop_max = 0
+        !$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp, myalpha_rho, myalpha)  &
+        !$acc reduction(MAX:adap_dt_stop_max) copy(adap_dt_stop_max)
         do l = 0, p
             do k = 0, n
                 do j = 0, m
@@ -283,8 +288,8 @@ contains
                                                 pb, pbdot, alf, n_tait, B_tait, &
                                                 bub_adv_src(j, k, l), divu%sf(j, k, l), &
                                                 dmBub_id, dmMass_v, dmMass_n, dmBeta_c, &
-                                                dmBeta_t, dmCson)
-
+                                                dmBeta_t, dmCson, adap_dt_stop)
+                            adap_dt_stop_max = max(adap_dt_stop_max, adap_dt_stop)
                             q_cons_vf(rs(q))%sf(j, k, l) = nbub*myR
                             q_cons_vf(vs(q))%sf(j, k, l) = nbub*myV
 
@@ -311,6 +316,8 @@ contains
             end do
         end do
 
+        if (adap_dt .and. adap_dt_stop_max > 0) stop "Adaptive time stepping failed to converge"
+
         if (.not. adap_dt) then
             !$acc parallel loop collapse(3) gang vector default(present)
             do l = 0, p

src/simulation/m_bubbles_EE.fpp

@@ -135,6 +135,8 @@ contains
         @:ALLOCATE(mtn_dposdt(1:nBubs_glb, 1:3, 1:lag_num_ts))
         @:ALLOCATE(mtn_dveldt(1:nBubs_glb, 1:3, 1:lag_num_ts))
 
+        if (adap_dt .and. f_is_default(adap_dt_tol)) adap_dt_tol = dflt_adap_dt_tol
+
         ! Starting bubbles
         call s_start_lagrange_inputs()
         call s_read_input_bubbles(q_cons_vf)
@@ -504,6 +506,8 @@ contains
         real(wp), dimension(contxe) :: myalpha_rho, myalpha
         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
         real(wp) :: dmalf, dmntait, dmBtait, dm_bub_adv_src, dm_divu !< Dummy variables for unified subgrid bubble subroutines
 
         integer :: i, k, l
@@ -533,7 +537,9 @@ contains
         end if
 
         ! Radial motion model
-        !$acc parallel loop gang vector default(present) private(k, myalpha_rho, myalpha, Re, cell) copyin(stage)
+        adap_dt_stop_max = 0
+        !$acc parallel loop gang vector default(present) private(k, myalpha_rho, myalpha, Re, cell) &
+        !$acc reduction(MAX:adap_dt_stop_max) copy(adap_dt_stop_max) copyin(stage)
         do k = 1, nBubs
             ! Keller-Miksis model
 
@@ -571,7 +577,8 @@ contains
                 call s_advance_step(myRho, myPinf, myR, myV, myR0, myPb, myPbdot, dmalf, &
                                     dmntait, dmBtait, dm_bub_adv_src, dm_divu, &
                                     k, myMass_v, myMass_n, myBeta_c, &
-                                    myBeta_t, myCson)
+                                    myBeta_t, myCson, adap_dt_stop)
+                adap_dt_stop_max = max(adap_dt_stop_max, adap_dt_stop)
 
                 ! Update bubble state
                 intfc_rad(k, 1) = myR
@@ -594,6 +601,8 @@ contains
 
         end do
 
+        if (adap_dt .and. adap_dt_stop_max > 0) stop "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)
         do k = 1, nBubs

src/simulation/m_bubbles_EL.fpp

@@ -385,6 +385,7 @@ module m_global_parameters
     logical :: polydisperse !< Polydisperse bubbles
     logical :: adv_n        !< Solve the number density equation and compute alpha from number density
     logical :: adap_dt      !< Adaptive step size control
+    real(wp) :: adap_dt_tol !< Tolerance to control adaptive step size
 
     integer :: bubble_model !< Gilmore or Keller--Miksis bubble model
     integer :: thermal      !< Thermal behavior. 1 = adiabatic, 2 = isotherm, 3 = transfer
@@ -406,7 +407,7 @@ module m_global_parameters
         !$acc declare create(nb)
     #:endif
 
-    !$acc declare create(R0ref, Ca, Web, Re_inv, bubbles_euler, polytropic, polydisperse, qbmm, nmomsp, nmomtot, R0_type, bubble_model, thermal, poly_sigma, adv_n, adap_dt, pi_fac)
+    !$acc declare create(R0ref, Ca, Web, Re_inv, bubbles_euler, polytropic, polydisperse, qbmm, nmomsp, nmomtot, R0_type, bubble_model, thermal, poly_sigma, adv_n, adap_dt, adap_dt_tol, pi_fac)
 
     type(scalar_field), allocatable, dimension(:) :: mom_sp
     type(scalar_field), allocatable, dimension(:, :, :) :: mom_3d
@@ -637,6 +638,7 @@ contains
 
         adv_n = .false.
         adap_dt = .false.
+        adap_dt_tol = dflt_real
 
         pi_fac = 1._wp
 

src/simulation/m_global_parameters.fpp

@@ -175,7 +175,7 @@ contains
             & 'x_domain%beg', 'x_domain%end', 'y_domain%beg', 'y_domain%end',    &
             & 'z_domain%beg', 'z_domain%end', 'x_a', 'x_b', 'y_a', 'y_b', 'z_a', &
             & 'z_b', 't_stop', 't_save', 'cfl_target', 'Bx0',  &
-            & 'tau_star', 'cont_damage_s', 'alpha_bar' ]
+            & 'tau_star', 'cont_damage_s', 'alpha_bar', 'adap_dt_tol' ]
             call MPI_BCAST(${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 

src/simulation/m_mpi_proxy.fpp

@@ -168,7 +168,8 @@ contains
             relax, relax_model, &
             palpha_eps, ptgalpha_eps, &
             R0_type, file_per_process, sigma, &
-            pi_fac, adv_n, adap_dt, bf_x, bf_y, bf_z, &
+            pi_fac, adv_n, adap_dt, adap_dt_tol, &
+            bf_x, bf_y, bf_z, &
             k_x, k_y, k_z, w_x, w_y, w_z, p_x, p_y, p_z, &
             g_x, g_y, g_z, n_start, t_save, t_stop, &
             cfl_adap_dt, cfl_const_dt, cfl_target, &
@@ -1621,7 +1622,7 @@ contains
         if (chemistry) then
             !$acc update device(q_T_sf%sf)
         end if
-        !$acc update device(nb, R0ref, Ca, Web, Re_inv, weight, R0, V0, bubbles_euler, polytropic, polydisperse, qbmm, R0_type, ptil, bubble_model, thermal, poly_sigma, adv_n, adap_dt, n_idx, pi_fac, low_Mach)
+        !$acc update device(nb, R0ref, Ca, Web, Re_inv, weight, R0, V0, bubbles_euler, polytropic, polydisperse, qbmm, R0_type, ptil, bubble_model, thermal, poly_sigma, adv_n, adap_dt, adap_dt_tol, n_idx, pi_fac, low_Mach)
         !$acc update device(R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v, k_n, k_v, pb0, mass_n0, mass_v0, Pe_T, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN , mul0, ss, gamma_v, mu_v, gamma_m, gamma_n, mu_n, gam)
 
         !$acc update device(acoustic_source, num_source)

src/simulation/m_start_up.fpp

@@ -244,6 +244,7 @@ def analytic(self):
     'ptgalpha_eps': ParamType.REAL,
     'pi_fac': ParamType.REAL,
     'adap_dt': ParamType.LOG,
+    'adap_dt_tol': ParamType.REAL,
     'ib': ParamType.LOG,
     'num_ibs': ParamType.INT,
     'n_start': ParamType.INT,