AMReX-Astro · zingale · Sep 23, 2023 · Sep 18, 2023 · Sep 18, 2023 · Sep 18, 2023
diff --git a/Docs/source/sdc.rst b/Docs/source/sdc.rst
@@ -42,23 +42,13 @@ The options that affect the nonlinear solve are:
   * 3 : use VODE for the first iteration and then Newton for the
     subsequent iterations.
 
-  In all cases, the type of Jacobian (analytic or numerical) is determined by
-  ``integrator.jacobian``.
-
-* ``sdc_solver_tol_dens`` : the relative error on the density in solving the nonlinear system.
-
-* ``sdc_solver_tol_spec`` : the relative error on the partial densities, :math:`(\rho X_k)`
-  in the nonlinear solve.
-
-* ``sdc_solver_tol_ener`` : the relative error of the energy in the nonlinear solve.
+  The tolerances for both the Newton and VODE solver are controlled by
+  the usual Microphysics parameters: ``integrator.rtol_spec``,
+  ``integrator.atol_spec``, ``integrator.rtol_enuc``,
+  ``integrator.atol_enuc``.
 
-* ``sdc_solver_atol`` : the absolute error in the mass fractions during the nonlinear solve.
-
-* ``sdc_solver_relax_factor`` : the factor by which to relax the
-  tolerances (i.e. increase them) for earlier iterations.  We reach
-  the desired tolerances on the final iteration.
-
-* ``sdc_solve_for_rhoe`` : whether we solve the system in terms of :math:`(\rho e)` or :math:`(\rho E)`.
+In all cases, the type of Jacobian (analytic or numerical) is determined by
+  ``integrator.jacobian``.
 
 
 

diff --git a/Exec/reacting_tests/bubble_convergence/inputs_2d.128 b/Exec/reacting_tests/bubble_convergence/inputs_2d.128
@@ -30,10 +30,6 @@ castro.time_integration_method=2
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_ener=1.e-5
-castro.sdc_solver_atol=1.e-10
 castro.sdc_solver=1
 
 castro.grav_source_type = 2
@@ -100,3 +96,8 @@ integrator.atol_spec = 1.e-10
 network.small_x = 1.e-10
 
 integrator.rtol_enuc = 1.e-10
+
+integrator.rtol_spec = 1.e-10
+integrator.rtol_enuc = 1.e-5
+integrator.atol_spec = 1.e-10
+integrator.atol_enuc = 1.e-10
diff --git a/Exec/reacting_tests/bubble_convergence/inputs_2d.256 b/Exec/reacting_tests/bubble_convergence/inputs_2d.256
@@ -30,10 +30,6 @@ castro.time_integration_method=2
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_ener=1.e-5
-castro.sdc_solver_atol=1.e-10
 castro.sdc_solver=1
 
 castro.grav_source_type = 2
@@ -99,4 +95,7 @@ integrator.atol_spec = 1.e-10
 
 network.small_x = 1.e-10
 
-integrator.rtol_enuc = 1.e-10
+integrator.rtol_spec = 1.e-10
+integrator.rtol_enuc = 1.e-5
+integrator.atol_spec = 1.e-10
+integrator.atol_enuc = 1.e-10
diff --git a/Exec/reacting_tests/bubble_convergence/inputs_2d.32 b/Exec/reacting_tests/bubble_convergence/inputs_2d.32
@@ -30,10 +30,6 @@ castro.time_integration_method=2
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_ener=1.e-5
-castro.sdc_solver_atol=1.e-10
 castro.sdc_solver=1
 
 castro.grav_source_type = 2
@@ -92,7 +88,7 @@ amr.refine.tempgrad.field_name = Temp
 
 integrator.rtol_spec = 1.e-10
 integrator.atol_spec = 1.e-10
+integrator.rtol_enuc = 1.e-5
+integrator.atol_enuc = 1.e-10
 
 network.small_x = 1.e-10
-
-integrator.rtol_enuc = 1.e-10
diff --git a/Exec/reacting_tests/bubble_convergence/inputs_2d.512 b/Exec/reacting_tests/bubble_convergence/inputs_2d.512
@@ -30,10 +30,6 @@ castro.time_integration_method=2
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_ener=1.e-5
-castro.sdc_solver_atol=1.e-10
 castro.sdc_solver=1
 
 castro.grav_source_type = 2
@@ -97,7 +93,7 @@ amr.refine.tempgrad.field_name = Temp
 
 integrator.rtol_spec = 1.e-10
 integrator.atol_spec = 1.e-10
+integrator.rtol_enuc = 1.e-5
+integrator.atol_enuc = 1.e-10
 
 network.small_x = 1.e-10
-
-integrator.rtol_enuc = 1.e-10
diff --git a/Exec/reacting_tests/bubble_convergence/inputs_2d.64 b/Exec/reacting_tests/bubble_convergence/inputs_2d.64
@@ -30,10 +30,6 @@ castro.time_integration_method=2
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_ener=1.e-5
-castro.sdc_solver_atol=1.e-10
 castro.sdc_solver=1
 
 castro.grav_source_type = 2
@@ -96,7 +92,7 @@ amr.refine.tempgrad.field_name = Temp
 
 integrator.rtol_spec = 1.e-10
 integrator.atol_spec = 1.e-10
+integrator.rtol_enuc = 1.e-5
+integrator.atol_enuc = 1.e-10
 
 network.small_x = 1.e-10
-
-integrator.rtol_enuc = 1.e-10
diff --git a/Exec/reacting_tests/bubble_convergence/job_scripts/cori.bubble_convergence.slurm b/Exec/reacting_tests/bubble_convergence/job_scripts/cori.bubble_convergence.slurm
@@ -22,10 +22,6 @@ castro.time_integration_method=2
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_ener=1.e-5
-castro.sdc_solver_atol=1.e-10
 castro.sdc_solver=1"
 
 srun ${CASTRO_EXEC} inputs_2d.64 ${RUNPARAMS} amr.plot_file=bubble_64_sdc4_plt &> 64.out

diff --git a/Exec/reacting_tests/reacting_convergence/convergence_sdc.sh b/Exec/reacting_tests/reacting_convergence/convergence_sdc.sh
@@ -11,9 +11,6 @@ castro.sdc_order=2
 castro.time_integration_method=2
 castro.ppm_type=0
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_ener=1.e-10
 castro.sdc_solver=1
 castro.use_retry=0"
 

diff --git a/Exec/reacting_tests/reacting_convergence/convergence_sdc4.sh b/Exec/reacting_tests/reacting_convergence/convergence_sdc4.sh
@@ -12,9 +12,6 @@ castro.time_integration_method=2
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_ener=1.e-10
 castro.sdc_solver=1
 castro.use_retry=0"
 

diff --git a/Exec/reacting_tests/reacting_convergence/inputs.128.testsuite b/Exec/reacting_tests/reacting_convergence/inputs.128.testsuite
@@ -77,5 +77,8 @@ amr.refine.denerr.field_name = density
 
 network.small_x = 1.e-10
 
-integrator.atol_spec = 1.e-12
+integrator.rtol_spec = 1.e-6
+integrator.atol_spec = 1.e-10
+integrator.rtol_enuc = 1.e-6
+integrator.atol_enuc = 1.e-10
 
diff --git a/Exec/reacting_tests/reacting_convergence/job_scripts/cori.reacting_convergence.slurm b/Exec/reacting_tests/reacting_convergence/job_scripts/cori.reacting_convergence.slurm
@@ -27,10 +27,6 @@ castro.time_integration_method=2
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_spec=1.e-10
-castro.sdc_solver_tol_dens=1.e-10
-castro.sdc_solver_tol_ener=1.e-5
-castro.sdc_solver_atol=1.e-10
 castro.sdc_solver=1"
 
 srun ${CASTRO_EXEC} inputs.64 amr.plot_file=react_converge_64_sdc4_plt ${RUNPARAMS}

diff --git a/Exec/science/Detonation/inputs-det-x.sdc.test b/Exec/science/Detonation/inputs-det-x.sdc.test
@@ -28,10 +28,6 @@ castro.time_integration_method = 2
 castro.sdc_order = 4
 
 castro.sdc_solver = 2
-castro.sdc_solver_tol_dens = 1.e-5
-castro.sdc_solver_tol_spec = 1.e-5
-castro.sdc_solver_tol_ener = 1.e-5
-castro.sdc_solver_atol = 1.e-10
 
 castro.ppm_type = 0
 
@@ -101,4 +97,10 @@ amr.refine.tempgrad.field_name = Temp
 # Microphysics
 
 integrator.call_eos_in_rhs = 1
+
+integrator.rtol_spec = 1.e-5
+integrator.rtol_enuc = 1.e-5
+integrator.atol_spec = 1.e-10
+integrator.atol_enuc = 1.e-10
+
 integrator.jacobian = 1
diff --git a/Exec/science/Detonation/nse_runs/inputs.template.true_sdc b/Exec/science/Detonation/nse_runs/inputs.template.true_sdc
@@ -39,10 +39,6 @@ castro.sdc_quadrature = @@SDC_QUADRATURE@@
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
-castro.sdc_solver_tol_dens=1.e-8
-castro.sdc_solver_tol_spec=1.e-8
-castro.sdc_solver_tol_ener=1.e-8
-castro.sdc_solver_atol = 1.e-8
 castro.sdc_solver=2
 
 
@@ -111,9 +107,10 @@ amr.refine.tempgrad.field_name = Temp
 
 integrator.call_eos_in_rhs = 1
 
-integrator.rtol_spec = 1.e-6
-integrator.atol_spec = 1.e-6
-integrator.rtol_enuc = 1.e-6
+integrator.rtol_spec = 1.e-8
+integrator.atol_spec = 1.e-8
+integrator.rtol_enuc = 1.e-8
+integrator.atol_enuc = 1.e-8
 
 integrator.jacobian = 2
 integrator.ode_max_steps = 1500000
diff --git a/Exec/science/flame/inputs.1d.sdc b/Exec/science/flame/inputs.1d.sdc
@@ -26,10 +26,6 @@ castro.time_integration_method = 2
 castro.sdc_order = 4
 
 castro.sdc_solve_for_rhoe = 1
-castro.sdc_solver_tol_dens = 1.e-10
-castro.sdc_solver_tol_spec = 1.e-10
-castro.sdc_solver_tol_ener = 1.e-6
-castro.sdc_solver_atol = 1.e-10
 castro.sdc_solver = 1
 castro.sdc_solver_relax_factor = 1
 
@@ -107,8 +103,10 @@ amr.refine.dengrad.field_name = density
 
 # Microphysics
 
-integrator.rtol_spec = 1.e-8
-integrator.atol_spec = 1.e-8
+integrator.rtol_spec = 1.e-10
+integrator.atol_spec = 1.e-10
+integrator.rtol_enuc = 1.e-6
+integrator.atol_enuc = 1.e-10
 
 network.small_x = 1.e-10
 

diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters
@@ -267,22 +267,6 @@ sdc_extra                    int           0
 # which SDC nonlinear solver to use?  1 = Newton, 2 = VODE, 3 = VODE for first iter
 sdc_solver                   int           1
 
-# relative tolerance for the nonlinear solve on rho X_k with SDC
-sdc_solver_tol_spec           Real          1.e-6
-
-# relative tolerance for the nonlinear solve on (rho e) or (rho E) with SDC
-sdc_solver_tol_ener           Real          1.e-6
-
-# absolute tolerance for species with SDC (this will be multiplied by
-# the current rho in the zone to define the absolute tolerance on (rho
-# X)).
-sdc_solver_atol               Real          1.e-10
-
-# factor by which to reduce the SDC solver tol for each iteration before the last
-# (e.g., for iteration k out of kmax iterations, the tol is
-# :math:`\epsilon/f^{(k_\mathrm{max} - k)}`
-sdc_solver_relax_factor      Real          1.0
-
 # do we solve for (rho e) or (rho E) in the SDC nonlinear solve?
 sdc_solve_for_rhoe           int           1
 

diff --git a/Source/sdc/sdc_newton_solve.H b/Source/sdc/sdc_newton_solve.H
@@ -100,7 +100,6 @@ f_sdc_jac(const Real dt_m,
 
     // get the Jacobian.
 
-    /// First we get the dR/dU expression.
     // Instead of the decomposition into dw/dU and dR/dw
     // written out in the original paper Appendix A, we instead use the
     // form from the simplified-SDC paper (Zingale et al. 2022).  Note:
@@ -157,12 +156,6 @@ sdc_newton_solve(const Real dt_m,
 
     ierr = NEWTON_SUCCESS;
 
-    // the tolerance we are solving to may depend on the
-    // iteration
-    Real relax_fac = std::pow(sdc_solver_relax_factor, sdc_order - sdc_iteration - 1);
-    Real tol_spec = sdc_solver_tol_spec * relax_fac;
-    Real tol_ener = sdc_solver_tol_ener * relax_fac;
-
     // update the density and momenta for this zone -- they don't react
 
     U_new[URHO] = U_old[URHO] + dt_m * C[URHO];
@@ -264,10 +257,11 @@ sdc_newton_solve(const Real dt_m,
         // multiply by density to get a partial density limit
 
         for (int n = 0; n < NumSpec; ++n) {
-            eps_tot(1 + n) = tol_spec * std::abs(U_react(1 + n)) +
-                sdc_solver_atol * std::abs(U_new[URHO]);
+            eps_tot(1 + n) = integrator_rp::rtol_spec * std::abs(U_react(1 + n)) +
+                integrator_rp::atol_spec * std::abs(U_new[URHO]);
         }
-        eps_tot(NumSpec+1) = tol_ener * std::abs(U_react(NumSpec+1)) + sdc_solver_atol;
+        eps_tot(NumSpec+1) = integrator_rp::rtol_enuc * std::abs(U_react(NumSpec+1)) +
+            integrator_rp::atol_enuc;
 
         // compute the norm of the weighted error, where the
         // weights are 1/eps_tot

diff --git a/Source/sdc/vode_rhs_true_sdc.H b/Source/sdc/vode_rhs_true_sdc.H
@@ -128,12 +128,6 @@ sdc_vode_solve(const Real dt_m,
 
 #if (INTEGRATOR == 0)
 
-    // The tolerance we are solving to may depend on the
-    // iteration
-    auto relax_fac = std::pow(sdc_solver_relax_factor, sdc_order - sdc_iteration - 1);
-    auto tol_spec = sdc_solver_tol_spec * relax_fac;
-    auto tol_ener = sdc_solver_tol_ener * relax_fac;
-
     // Update the momenta for this zone -- they don't react
 
     for (int n = 0; n < 3; ++n) {
@@ -207,12 +201,12 @@ sdc_vode_solve(const Real dt_m,
     dvode_state.jacobian_type = integrator_rp::jacobian;
 
     // relative tolerances
-    dvode_state.rtol_spec = tol_spec;
-    dvode_state.rtol_enuc = tol_ener;
+    dvode_state.rtol_spec = integrator_rp::rtol_spec;
+    dvode_state.rtol_enuc = integrator_rp::rtol_enuc;
 
     // absolute tolerances
-    dvode_state.atol_spec = sdc_solver_atol * U_old[URHO]; // this way, atol is the minimum x
-    dvode_state.atol_enuc = sdc_solver_atol * U_old[UEINT];
+    dvode_state.atol_spec = integrator_rp::atol_spec * U_old[URHO]; // this way, atol is the minimum x
+    dvode_state.atol_enuc = integrator_rp::atol_enuc * U_old[UEINT];
 
     int istate = dvode(burn_state, dvode_state);