Turn on castro.use_retry by default

CHANGES.md

@@ -1,5 +1,13 @@
 # 20.02
 
+   * We now default to use_retry = 1, instructing Castro to retry a
+     step with a smaller dt if there is a CFL violation, burning
+     failure, or negative timestep.  For the burning failure, we have
+     Castro set the Microphysics parameter abort_on_failure to .false.
+     at a high priority (so it overrides the Microphysics default).
+     We also check to make sure the combination of parameters makes
+     sense at runtime. (#724)
+
    * The parameter castro.hard_cfl_limit has been removed. (#723)
 
    * Some unnecessary clean_state calls were removed (#721)

Microphysics/networks/_parameters

@@ -1,4 +1,6 @@
 # cutoff for species mass fractions
 small_x                              real               1.d-30
 
-
+# override the Microphysics default -- this enables use_retry to take
+# action for burn failures
+abort_on_failure                     logical            .false.             100

Source/driver/Castro_F.H

@@ -40,6 +40,8 @@ extern "C"
 
   void ca_get_ngdnv(int* ngdnv);
 
+  void ca_get_abort_on_failure(int * abort_on_failure);
+
   void ca_amrinfo_init();
   void ca_amrinfo_finalize();
 

Source/driver/Castro_nd.F90

@@ -287,6 +287,26 @@ subroutine ca_get_ngdnv(ngdnv_in) bind(C, name="ca_get_ngdnv")
 
 end subroutine ca_get_ngdnv
 
+
+#ifdef REACTIONS
+subroutine ca_get_abort_on_failure(abort_on_failure_in) bind(C, name="ca_get_abort_on_failure")
+
+  use extern_probin_module, only : abort_on_failure
+
+  implicit none
+
+  integer, intent(inout) :: abort_on_failure_in
+
+  if (abort_on_failure) then
+     abort_on_failure_in = 1
+  else
+     abort_on_failure_in = 0
+  endif
+
+end subroutine ca_get_abort_on_failure
+#endif
+
+
 ! :::
 ! ::: ----------------------------------------------------------------
 ! :::

Source/driver/Castro_setup.cpp

@@ -153,16 +153,29 @@ Castro::variableSetUp ()
 
   BL_ASSERT(desc_lst.size() == 0);
 
-  // Get options, set phys_bc
+  // read the C++ parameters that are set in inputs and do other
+  // initializations (e.g., set phys_bc)
   read_params();
 
   // Initialize the runtime parameters for any of the external
-  // microphysics
+  // microphysics (these are the parameters that are in the &extern
+  // block of the probin file)
   extern_init();
 
   // Initialize the network
   network_init();
 
+
+  // some consistency checks on the parameters
+#ifdef REACTIONS
+  int abort_on_failure;
+  ca_get_abort_on_failure(&abort_on_failure);
+
+  if (!use_retry && !abort_on_failure) {
+    amrex::Error("use_retry = 0 and abort_on_failure = F is dangerous and not supported");
+  }
+#endif
+
 #ifdef REACTIONS
   // Initialize the burner
   burner_init();

Source/driver/_cpp_parameters

@@ -313,11 +313,12 @@ plot_per_is_exact            int           0
 # enforce that the AMR small plot interval must be hit exactly
 small_plot_per_is_exact      int           0
 
-# Retry a timestep if it violated the timestep-limiting criteria over
-# the course of an advance. The criteria will suggest a new timestep
-# that satisfies the criteria, and we will do subcycled timesteps
-# on the same level until we reach the original target time.
-use_retry                    int           0
+# Retry a timestep if it violated the timestep-limiting criteria or
+# other checks (negative density, burn failure) over the course of an
+# advance. The criteria will suggest a new timestep that satisfies the
+# criteria, and we will do subcycled timesteps on the same level until
+# we reach the original target time.
+use_retry                    int           1
 
 # Tolerance to use when evaluating whether to do a retry.
 # The timestep suggested by the retry will be multiplied by

Source/driver/meth_params.F90

@@ -543,16 +543,6 @@ subroutine ca_set_castro_method_params() bind(C, name="ca_set_castro_method_para
     point_mass = 0.0_rt;
     allocate(point_mass_fix_solution)
     point_mass_fix_solution = 0;
-#endif
-#ifdef DIFFUSION
-    allocate(diffuse_temp)
-    diffuse_temp = 0;
-    allocate(diffuse_cutoff_density)
-    diffuse_cutoff_density = -1.e200_rt;
-    allocate(diffuse_cutoff_density_hi)
-    diffuse_cutoff_density_hi = -1.e200_rt;
-    allocate(diffuse_cond_scale_fac)
-    diffuse_cond_scale_fac = 1.0_rt;
 #endif
     allocate(difmag)
     difmag = 0.1_rt;
@@ -704,6 +694,16 @@ subroutine ca_set_castro_method_params() bind(C, name="ca_set_castro_method_para
     grown_factor = 1;
     allocate(track_grid_losses)
     track_grid_losses = 0;
+#ifdef DIFFUSION
+    allocate(diffuse_temp)
+    diffuse_temp = 0;
+    allocate(diffuse_cutoff_density)
+    diffuse_cutoff_density = -1.e200_rt;
+    allocate(diffuse_cutoff_density_hi)
+    diffuse_cutoff_density_hi = -1.e200_rt;
+    allocate(diffuse_cond_scale_fac)
+    diffuse_cond_scale_fac = 1.0_rt;
+#endif
 
     call amrex_parmparse_build(pp, "castro")
 #ifdef ROTATION
@@ -721,12 +721,6 @@ subroutine ca_set_castro_method_params() bind(C, name="ca_set_castro_method_para
     call pp%query("use_point_mass", use_point_mass)
     call pp%query("point_mass", point_mass)
     call pp%query("point_mass_fix_solution", point_mass_fix_solution)
-#endif
-#ifdef DIFFUSION
-    call pp%query("diffuse_temp", diffuse_temp)
-    call pp%query("diffuse_cutoff_density", diffuse_cutoff_density)
-    call pp%query("diffuse_cutoff_density_hi", diffuse_cutoff_density_hi)
-    call pp%query("diffuse_cond_scale_fac", diffuse_cond_scale_fac)
 #endif
     call pp%query("difmag", difmag)
     call pp%query("small_dens", small_dens)
@@ -803,6 +797,26 @@ subroutine ca_set_castro_method_params() bind(C, name="ca_set_castro_method_para
     call pp%query("do_acc", do_acc)
     call pp%query("grown_factor", grown_factor)
     call pp%query("track_grid_losses", track_grid_losses)
+#ifdef DIFFUSION
+    call pp%query("diffuse_temp", diffuse_temp)
+    call pp%query("diffuse_cutoff_density", diffuse_cutoff_density)
+    call pp%query("diffuse_cutoff_density_hi", diffuse_cutoff_density_hi)
+    call pp%query("diffuse_cond_scale_fac", diffuse_cond_scale_fac)
+#endif
+    call amrex_parmparse_destroy(pp)
+
+
+    allocate(character(len=1)::gravity_type)
+    gravity_type = "fillme";
+    allocate(const_grav)
+    const_grav = 0.0_rt;
+    allocate(get_g_from_phi)
+    get_g_from_phi = 0;
+
+    call amrex_parmparse_build(pp, "gravity")
+    call pp%query("gravity_type", gravity_type)
+    call pp%query("const_grav", const_grav)
+    call pp%query("get_g_from_phi", get_g_from_phi)
     call amrex_parmparse_destroy(pp)
 
 
@@ -826,20 +840,6 @@ subroutine ca_set_castro_method_params() bind(C, name="ca_set_castro_method_para
     call amrex_parmparse_destroy(pp)
 
 
-    allocate(character(len=1)::gravity_type)
-    gravity_type = "fillme";
-    allocate(const_grav)
-    const_grav = 0.0_rt;
-    allocate(get_g_from_phi)
-    get_g_from_phi = 0;
-
-    call amrex_parmparse_build(pp, "gravity")
-    call pp%query("gravity_type", gravity_type)
-    call pp%query("const_grav", const_grav)
-    call pp%query("get_g_from_phi", get_g_from_phi)
-    call amrex_parmparse_destroy(pp)
-
-
 
     !$acc update &
     !$acc device(difmag, small_dens, small_temp) &

Source/driver/param_includes/castro_defaults.H

@@ -14,11 +14,6 @@ int         Castro::rot_source_type = 4;
 int         Castro::implicit_rotation_update = 1;
 int         Castro::rot_axis = 3;
 #endif
-#ifdef GRAVITY
-int         Castro::use_point_mass = 0;
-amrex::Real Castro::point_mass = 0.0;
-int         Castro::point_mass_fix_solution = 0;
-#endif
 #ifdef AMREX_PARTICLES
 int         Castro::do_tracer_particles = 0;
 #endif
@@ -28,6 +23,11 @@ amrex::Real Castro::diffuse_cutoff_density = -1.e200;
 amrex::Real Castro::diffuse_cutoff_density_hi = -1.e200;
 amrex::Real Castro::diffuse_cond_scale_fac = 1.0;
 #endif
+#ifdef GRAVITY
+int         Castro::use_point_mass = 0;
+amrex::Real Castro::point_mass = 0.0;
+int         Castro::point_mass_fix_solution = 0;
+#endif
 int         Castro::state_interp_order = 1;
 int         Castro::lin_limit_state_interp = 0;
 int         Castro::state_nghost = 0;
@@ -105,7 +105,7 @@ amrex::Real Castro::init_shrink = 1.0;
 amrex::Real Castro::change_max = 1.1;
 int         Castro::plot_per_is_exact = 0;
 int         Castro::small_plot_per_is_exact = 0;
-int         Castro::use_retry = 0;
+int         Castro::use_retry = 1;
 amrex::Real Castro::retry_tolerance = 0.02;
 amrex::Real Castro::retry_neg_dens_factor = 1.e-1;
 amrex::Real Castro::retry_subcycle_factor = 0.5;

Source/driver/param_includes/castro_job_info_tests.H

@@ -9,11 +9,6 @@ jobInfoFile << (Castro::rot_source_type == 4 ? "    " : "[*] ") << "castro.rot_s
 jobInfoFile << (Castro::implicit_rotation_update == 1 ? "    " : "[*] ") << "castro.implicit_rotation_update = " << Castro::implicit_rotation_update << std::endl;
 jobInfoFile << (Castro::rot_axis == 3 ? "    " : "[*] ") << "castro.rot_axis = " << Castro::rot_axis << std::endl;
 #endif
-#ifdef GRAVITY
-jobInfoFile << (Castro::use_point_mass == 0 ? "    " : "[*] ") << "castro.use_point_mass = " << Castro::use_point_mass << std::endl;
-jobInfoFile << (Castro::point_mass == 0.0 ? "    " : "[*] ") << "castro.point_mass = " << Castro::point_mass << std::endl;
-jobInfoFile << (Castro::point_mass_fix_solution == 0 ? "    " : "[*] ") << "castro.point_mass_fix_solution = " << Castro::point_mass_fix_solution << std::endl;
-#endif
 #ifdef AMREX_PARTICLES
 jobInfoFile << (Castro::do_tracer_particles == 0 ? "    " : "[*] ") << "castro.do_tracer_particles = " << Castro::do_tracer_particles << std::endl;
 #endif
@@ -23,6 +18,11 @@ jobInfoFile << (Castro::diffuse_cutoff_density == -1.e200 ? "    " : "[*] ") <<
 jobInfoFile << (Castro::diffuse_cutoff_density_hi == -1.e200 ? "    " : "[*] ") << "castro.diffuse_cutoff_density_hi = " << Castro::diffuse_cutoff_density_hi << std::endl;
 jobInfoFile << (Castro::diffuse_cond_scale_fac == 1.0 ? "    " : "[*] ") << "castro.diffuse_cond_scale_fac = " << Castro::diffuse_cond_scale_fac << std::endl;
 #endif
+#ifdef GRAVITY
+jobInfoFile << (Castro::use_point_mass == 0 ? "    " : "[*] ") << "castro.use_point_mass = " << Castro::use_point_mass << std::endl;
+jobInfoFile << (Castro::point_mass == 0.0 ? "    " : "[*] ") << "castro.point_mass = " << Castro::point_mass << std::endl;
+jobInfoFile << (Castro::point_mass_fix_solution == 0 ? "    " : "[*] ") << "castro.point_mass_fix_solution = " << Castro::point_mass_fix_solution << std::endl;
+#endif
 jobInfoFile << (Castro::state_interp_order == 1 ? "    " : "[*] ") << "castro.state_interp_order = " << Castro::state_interp_order << std::endl;
 jobInfoFile << (Castro::lin_limit_state_interp == 0 ? "    " : "[*] ") << "castro.lin_limit_state_interp = " << Castro::lin_limit_state_interp << std::endl;
 jobInfoFile << (Castro::state_nghost == 0 ? "    " : "[*] ") << "castro.state_nghost = " << Castro::state_nghost << std::endl;
@@ -100,7 +100,7 @@ jobInfoFile << (Castro::init_shrink == 1.0 ? "    " : "[*] ") << "castro.init_sh
 jobInfoFile << (Castro::change_max == 1.1 ? "    " : "[*] ") << "castro.change_max = " << Castro::change_max << std::endl;
 jobInfoFile << (Castro::plot_per_is_exact == 0 ? "    " : "[*] ") << "castro.plot_per_is_exact = " << Castro::plot_per_is_exact << std::endl;
 jobInfoFile << (Castro::small_plot_per_is_exact == 0 ? "    " : "[*] ") << "castro.small_plot_per_is_exact = " << Castro::small_plot_per_is_exact << std::endl;
-jobInfoFile << (Castro::use_retry == 0 ? "    " : "[*] ") << "castro.use_retry = " << Castro::use_retry << std::endl;
+jobInfoFile << (Castro::use_retry == 1 ? "    " : "[*] ") << "castro.use_retry = " << Castro::use_retry << std::endl;
 jobInfoFile << (Castro::retry_tolerance == 0.02 ? "    " : "[*] ") << "castro.retry_tolerance = " << Castro::retry_tolerance << std::endl;
 jobInfoFile << (Castro::retry_neg_dens_factor == 1.e-1 ? "    " : "[*] ") << "castro.retry_neg_dens_factor = " << Castro::retry_neg_dens_factor << std::endl;
 jobInfoFile << (Castro::retry_subcycle_factor == 0.5 ? "    " : "[*] ") << "castro.retry_subcycle_factor = " << Castro::retry_subcycle_factor << std::endl;

Source/driver/param_includes/castro_params.H

@@ -14,11 +14,6 @@ static int rot_source_type;
 static int implicit_rotation_update;
 static int rot_axis;
 #endif
-#ifdef GRAVITY
-static int use_point_mass;
-static amrex::Real point_mass;
-static int point_mass_fix_solution;
-#endif
 #ifdef AMREX_PARTICLES
 static int do_tracer_particles;
 #endif
@@ -28,6 +23,11 @@ static amrex::Real diffuse_cutoff_density;
 static amrex::Real diffuse_cutoff_density_hi;
 static amrex::Real diffuse_cond_scale_fac;
 #endif
+#ifdef GRAVITY
+static int use_point_mass;
+static amrex::Real point_mass;
+static int point_mass_fix_solution;
+#endif
 static int state_interp_order;
 static int lin_limit_state_interp;
 static int state_nghost;

Source/driver/param_includes/castro_queries.H

@@ -14,11 +14,6 @@ pp.query("rot_source_type", rot_source_type);
 pp.query("implicit_rotation_update", implicit_rotation_update);
 pp.query("rot_axis", rot_axis);
 #endif
-#ifdef GRAVITY
-pp.query("use_point_mass", use_point_mass);
-pp.query("point_mass", point_mass);
-pp.query("point_mass_fix_solution", point_mass_fix_solution);
-#endif
 #ifdef AMREX_PARTICLES
 pp.query("do_tracer_particles", do_tracer_particles);
 #endif
@@ -28,6 +23,11 @@ pp.query("diffuse_cutoff_density", diffuse_cutoff_density);
 pp.query("diffuse_cutoff_density_hi", diffuse_cutoff_density_hi);
 pp.query("diffuse_cond_scale_fac", diffuse_cond_scale_fac);
 #endif
+#ifdef GRAVITY
+pp.query("use_point_mass", use_point_mass);
+pp.query("point_mass", point_mass);
+pp.query("point_mass_fix_solution", point_mass_fix_solution);
+#endif
 pp.query("state_interp_order", state_interp_order);
 pp.query("lin_limit_state_interp", lin_limit_state_interp);
 pp.query("state_nghost", state_nghost);

sphinx_docs/source/inputs.rst

@@ -404,10 +404,6 @@ The following parameters affect the timestep choice:
   * ``castro.dt_cutoff``: time step below which calculation
     will abort (Real :math:`> 0`; default: 0.0)
 
-  * ``castro.hard_cfl_limit``: whether or not to abort the
-    simulation if the hydrodynamics update creates velocities that
-    violate the CFL criterion (Integer; default: 1)
-
 As an example, consider::
 
     castro.cfl = 0.9

sphinx_docs/source/retry.rst

@@ -6,8 +6,11 @@ time step on a level and restart with a smaller timestep, subcycling
 within the level to make up the full time step needed for that level.
 It is enabled by setting::
 
-   castro.use_retry      = 1
-   castro.hard_cfl_limit = 0
+   castro.use_retry = 1
+
+.. note::
+
+   The Castro retry mechanism is enabled by default.
 
 The number of subcycles to try in the level is controlled via the
 ``castro.max_subcycles`` parameter.  It is not really suggested to go
@@ -27,7 +30,13 @@ A retry can be triggered by a number of conditions:
       retry_burn = F
       abort_on_failure = F
 
-
+    This instructs the integration routine in Microphysics to not
+    abort when the integration fails, but instead to tell the calling
+    Castro routine that the integration failed so Castro can handle
+    the retry itself.
+
+    The combination of ``use_retry = 0`` and ``abort_on_failure = F``
+    is unsafe and not supported.