diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 5b77afc8f7..cca4a8e3b3 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -272,6 +272,15 @@ sdc_extra int 0 # which SDC nonlinear solver to use? 1 = Newton, 2 = VODE, 3 = VODE for first iter sdc_solver int 1 +# do we apply damping to the corrections in the SDC Newton solver? +sdc_newton_use_damping int 0 + +# for the SDC Newton solver with damping, what is the mass fraction +# threshold below which we ignore species in the damping factor +# calculation? +sdc_newton_damping_xcutoff Real 1.e-4 + + # Do we include geometry source terms due to local unit vectors in non-Cartesian Coord? # We currently support R-Z cylinderical 2D (Bernand-Champmartin) and R-THETA spherical 2D use_geom_source bool 1 diff --git a/Source/sdc/Castro_sdc_util.H b/Source/sdc/Castro_sdc_util.H index 5fd6c5a5ee..8b8f9a3131 100644 --- a/Source/sdc/Castro_sdc_util.H +++ b/Source/sdc/Castro_sdc_util.H @@ -1,6 +1,8 @@ #ifndef CASTRO_SDC_UTIL_H #define CASTRO_SDC_UTIL_H +#include + #include #ifdef REACTIONS #ifdef NEW_NETWORK_IMPLEMENTATION @@ -56,11 +58,17 @@ sdc_solve(const Real dt_m, // We are going to assume we already have a good guess // for the solve in U_new and just pass the solve onto // the main Newton solve - sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr); + std::vector err_history{}; + sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_history, ierr); // failing? if (ierr != newton::NEWTON_SUCCESS) { - amrex::Abort("Newton subcycling failed in sdc_solve"); + std::cout << std::format("Newton subcycling failed in sdc_solve with ierr = {}\n", + ierr); + for (auto e : err_history) { + std::cout << e << std::endl; + } + amrex::Abort("Newton failure"); } } else if (sdc_solver == VODE_SOLVE) { // Use VODE to do the solution @@ -75,7 +83,8 @@ sdc_solve(const Real dt_m, // Now U_new is the update that VODE predicts, so we // will use that as the initial guess to the Newton solve - sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr); + std::vector err_history{}; + sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_history, ierr); // Failing? if (ierr != newton::NEWTON_SUCCESS) { diff --git a/Source/sdc/sdc_newton_solve.H b/Source/sdc/sdc_newton_solve.H index bf5be1e2d0..cb41f9348c 100644 --- a/Source/sdc/sdc_newton_solve.H +++ b/Source/sdc/sdc_newton_solve.H @@ -11,13 +11,15 @@ namespace newton { constexpr int CONVERGENCE_FAILURE = -2; constexpr int BAD_MASS_FRACTIONS = -3; - constexpr Real species_failure_tolerance = 1.e-2_rt; + constexpr amrex::Real species_failure_tolerance = 1.e-2_rt; - constexpr int MAX_ITER = 100; - constexpr int MAX_NSUB = 64; + constexpr int MAX_ITER = 50; + constexpr int MAX_NSUB = 2048; }; + + #ifdef REACTIONS AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -170,10 +172,53 @@ sdc_newton_solve(const Real dt_m, #endif // on output, f is the correction to the solution (dU) + + // we want to limit the change by some amount k, such that + // + // 1/k < (U + dU) / U < k + // + // we'll write this as + // + // 1/k < 1 + φ < k + // + // where φ = dU/U + // + // we want to limit this, so we will introduce a term α, such that + // + // 1/k < 1 + α φ < k + // + // if φ > 0, then we need α < (k - 1) / φ + // + // if φ < 0, then we need α < |1/k - 1| / |φ| + + amrex::Real alpha{1.0_rt}; + constexpr amrex::Real k = 10.0; + + if (castro::sdc_newton_use_damping) { + for (int n = 0; n < NumSpec; ++n) { + const amrex::Real rhoX = burn_state.y[SFS+n]; + if (rhoX < burn_state.rho * castro::sdc_newton_damping_xcutoff) { + continue; + } + const amrex::Real drhoX = f(n + 1); + const amrex::Real phi = drhoX / rhoX; + + if (phi > 0.0_rt) { + alpha = std::min(alpha, (k - 1.0_rt) / phi); + } else if (phi < 0.0_rt) { + alpha = std::min(alpha, (1.0_rt / k - 1.0_rt) / phi); + } + } + + // make sure we apply at least some correction + alpha = std::max(alpha, 5.e-2_rt); + } + for (int n = 0; n < NumSpec; ++n) { - burn_state.y[SFS+n] += f(n + 1); + burn_state.y[SFS+n] += alpha * f(n + 1); + } - burn_state.y[SEINT] += f(NumSpec+1); + burn_state.y[SEINT] += alpha * f(NumSpec+1); // compute the norm of the weighted error, where the // weights are 1/eps_tot @@ -190,7 +235,7 @@ sdc_newton_solve(const Real dt_m, eps = integrator_rp::rtol_enuc * std::abs(burn_state.y[SEINT]) + integrator_rp::atol_enuc; } - err_sum += f(n) * f(n) / (eps * eps); + err_sum += (alpha * f(n)) * (alpha * f(n)) / (eps * eps); } err = std::sqrt(err_sum / static_cast(NumSpec+1)); @@ -231,7 +276,7 @@ sdc_newton_subdivide(const Real dt_m, GpuArray& U_new, GpuArray const& C, const int sdc_iteration, - Real& err_out, + std::vector& err_out, int& ierr) { // This is the driver for solving the nonlinear update for the @@ -282,7 +327,9 @@ sdc_newton_subdivide(const Real dt_m, U_begin[UFS + n] *= U_begin[URHO] / sum_rhoX; } - sdc_newton_solve(dt_sub, U_begin, U_new, C, sdc_iteration, err_out, ierr); + amrex::Real err; + sdc_newton_solve(dt_sub, U_begin, U_new, C, sdc_iteration, err, ierr); + err_out.push_back(err); // our solve may have resulted in mass fractions outside // of [0, 1] -- reject if this is the case