Addressing comments 2

Author: Diego Vaca

src/common/m_constants.fpp

@@ -53,7 +53,14 @@ module m_constants
     ! 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
-
+    ! Constants of the algorithm described by Heirer, E. Hairer S.P.Nørsett G. Wanner, Solving Ordinary Differential Equations I, Chapter II.4
+    ! to choose the initial time step size for the adaptive time stepping routine
+    real(wp), parameter :: threshold_first_guess = 1e-5_wp
+    real(wp), parameter :: threshold_second_guess = 1e-15_wp
+    real(wp), parameter :: scale_first_guess = 1e-3_wp
+    real(wp), parameter :: scale_guess = 1e-2_wp
+    real(wp), parameter :: small_guess = 1e-6_wp
+    
     ! Relativity
     integer, parameter :: relativity_cons_to_prim_max_iter = 100
 

src/simulation/m_bubbles.fpp

@@ -612,10 +612,10 @@ contains
         ! Compute d_0 = ||y0|| and d_1 = ||f(x0,y0)||
         d_norms(1) = sqrt((myR_tmp(1)**2._wp + myV_tmp(1)**2._wp)/2._wp)
         d_norms(2) = sqrt((myV_tmp(1)**2._wp + myA_tmp(1)**2._wp)/2._wp)
-        if (d_norms(1) < 1e-5_wp .or. d_norms(2) < 1e-5_wp) then
-            h_size(1) = 1e-6_wp
+        if (d_norms(1) < threshold_first_guess .or. d_norms(2) < threshold_first_guess) then
+            h_size(1) = small_guess
         else
-            h_size(1) = 1e-2_wp*(d_norms(1)/d_norms(2))
+            h_size(1) = scale_guess*(d_norms(1)/d_norms(2))
         end if
 
         ! Evaluate f(x0+h0,y0+h0*f(x0,y0))
@@ -631,13 +631,13 @@ contains
 
         ! Set h1 = (0.01/max(d_1,d_2))^{1/(p+1)}
         !      if max(d_1,d_2) < 1e-15_wp, h_size(2) = max(1e-6_wp, h0*1e-3_wp)
-        if (max(d_norms(2), d_norms(3)) < 1e-15_wp) then
-            h_size(2) = max(1e-6_wp, h_size(1)*1e-3_wp)
+        if (max(d_norms(2), d_norms(3)) < threshold_second_guess) then
+            h_size(2) = max(small_guess, h_size(1)*scale_first_guess)
         else
-            h_size(2) = (1e-2_wp/max(d_norms(2), d_norms(3)))**(1._wp/3._wp)
+            h_size(2) = (scale_guess/max(d_norms(2), d_norms(3)))**(1._wp/3._wp)
         end if
 
-        f_initial_substep_h = min(100._wp*h_size(1), h_size(2))
+        f_initial_substep_h = min(h_size(1)/scale_guess, h_size(2))
 
     end function f_initial_substep_h
 

src/simulation/m_bubbles_EE.fpp

@@ -282,14 +282,16 @@ contains
                         end if
 
                         ! Adaptive time stepping
+                        adap_dt_stop = 0
+
                         if (adap_dt) then
 
                             call s_advance_step(myRho, myP, myR, myV, R0(q), &
                                                 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, 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
 
@@ -302,6 +304,8 @@ contains
                             bub_r_src(j, k, l, q) = q_cons_vf(vs(q))%sf(j, k, l)
                         end if
 
+                        adap_dt_stop_max = max(adap_dt_stop_max, adap_dt_stop)
+
                         if (alf < 1.e-11_wp) then
                             bub_adv_src(j, k, l) = 0._wp
                             bub_r_src(j, k, l, q) = 0._wp

src/simulation/m_bubbles_EL.fpp

@@ -572,13 +572,14 @@ contains
             call s_compute_cson_from_pinf(k, q_prim_vf, myPinf, cell, myRho, gamma, pi_inf, myCson)
 
             ! Adaptive time stepping
+            adap_dt_stop = 0
+            
             if (adap_dt) then
 
                 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, adap_dt_stop)
-                adap_dt_stop_max = max(adap_dt_stop_max, adap_dt_stop)
 
                 ! Update bubble state
                 intfc_rad(k, 1) = myR
@@ -599,6 +600,8 @@ contains
 
             end if
 
+            adap_dt_stop_max = max(adap_dt_stop_max, adap_dt_stop)
+
         end do
 
         if (adap_dt .and. adap_dt_stop_max > 0) stop "Adaptive time stepping failed to converge"

src/simulation/m_checker.fpp

@@ -115,7 +115,7 @@ contains
         if (.not. cfl_dt) then
             @:PROHIBIT(dt <= 0)
         end if
-        @:PROHIBIT(time_stepper < 1 .or. time_stepper > 4)
+        @:PROHIBIT(time_stepper < 1 .or. time_stepper > 3)
     end subroutine s_check_inputs_time_stepping
 
     !> Checks constraints on parameters related to 6-equation model