From a2994b1a7364e021dc75ec39b71538000d33f6ce Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 7 Oct 2023 10:54:43 -0400 Subject: [PATCH 1/2] switch true-SDC to completely use burn_t internally in the Newton solve --- Source/sdc/sdc_newton_solve.H | 88 +++++++++++------------------------ 1 file changed, 27 insertions(+), 61 deletions(-) diff --git a/Source/sdc/sdc_newton_solve.H b/Source/sdc/sdc_newton_solve.H index d7a1b96ab2..c4b811af1f 100644 --- a/Source/sdc/sdc_newton_solve.H +++ b/Source/sdc/sdc_newton_solve.H @@ -15,12 +15,9 @@ constexpr int CONVERGENCE_FAILURE = -2; AMREX_GPU_HOST_DEVICE AMREX_INLINE void f_sdc_jac(const Real dt_m, - Array1D const& U, + burn_t& burn_state, Array1D& f, - JacNetArray2D& Jac, - Array1D& f_source, - const Real rho, - const Real T_old) { + JacNetArray2D& Jac) { // This is used with the Newton solve and returns f and the Jacobian @@ -29,24 +26,6 @@ f_sdc_jac(const Real dt_m, // create a burn_t -- this will hold the full state, reconstructed // from the current solution U - burn_t burn_state; - - burn_state.y[SRHO] = rho; - - for (int n = 0; n < NumSpec; ++n) { - burn_state.y[SFS+n] = U(n+1); - } - - burn_state.y[SEINT] = U(NumSpec+1); - - // we are not solving or accessing the momentum or total energy, - // so we'll just set those to zero here - - for (int n = 0; n < 3; ++n) { - burn_state.y[SMX+n] = 0.0_rt; - } - - burn_state.y[SEDEN] = 0.0_rt; // normalize the species auto sum_rhoX = 0.0_rt; @@ -63,8 +42,6 @@ f_sdc_jac(const Real dt_m, // compute the temperature -- this can be removed shortly, since // we already do this in single_zone_react_source - burn_state.rho = burn_state.y[SRHO]; - burn_state.T = T_old; // initial guess for (int n = 0; n < NumSpec; ++n) { burn_state.xn[n] = burn_state.y[SFS+n] / burn_state.y[SRHO]; } @@ -83,9 +60,9 @@ f_sdc_jac(const Real dt_m, // note: we store -f for (int n = 1; n <= NumSpec; ++n) { - f(n) = -U(n) + dt_m * R_full[UFS-1+n] + f_source(n); + f(n) = -burn_state.y[SFS-1+n] + dt_m * R_full[UFS-1+n] + burn_state.ydot_a[SFS-1+n]; } - f(NumSpec+1) = -U(NumSpec+1) + dt_m * R_full[UEINT] + f_source(NumSpec+1); + f(NumSpec+1) = -burn_state.y[SEINT] + dt_m * R_full[UEINT] + burn_state.ydot_a[SEINT]; // get the Jacobian. @@ -138,8 +115,6 @@ sdc_newton_solve(const Real dt_m, // 1:NumSpec : species // NumSpec+1 : (rho e) - Array1D U_react; - Array1D f_source; Array1D f; const int MAX_ITER = 100; @@ -154,37 +129,25 @@ sdc_newton_solve(const Real dt_m, U_new[UMX+n] = U_old[UMX+n] + dt_m * C[UMX+n]; } - // now only save the subset that participates in the - // nonlinear solve -- note: we include the old state in - // f_source + burn_t burn_state; + + copy_cons_to_burn_type(U_new, burn_state); + burn_state.rho = U_new[URHO]; + + // initial guess + burn_state.T = U_old[UTEMP]; // for the Jacobian solve, we are solving // f(U) = U - dt R(U) - U_old - dt C = 0 // we define f_source = U_old + dt C so we are solving // f(U) = U - dt R(U) - f_source = 0 + // + // we'll store this in the burn_t ydot_a[] for (int n = 0; n < NumSpec; ++n) { - f_source(1 + n) = U_old[UFS + n] + dt_m * C[UFS + n]; - } - f_source(NumSpec+1) = U_old[UEINT] + dt_m * C[UEINT]; - - // temperature will be used as an initial guess in the EOS - - Real T_old = U_old[UTEMP]; - - // we need the density too - - Real rho_new = U_new[URHO]; - - // store the subset for the nonlinear solve - // We use an initial guess if possible - - for (int n = 0; n < NumSpec; ++n) { - U_react(1+n) = U_new[UFS+n]; + burn_state.ydot_a[SFS + n] = U_old[UFS + n] + dt_m * C[UFS + n]; } - U_react(NumSpec+1) = U_new[UEINT]; - -#if (INTEGRATOR == 0) + burn_state.ydot_a[SEINT] = U_old[UEINT] + dt_m * C[UEINT]; // do a simple Newton solve @@ -195,8 +158,12 @@ sdc_newton_solve(const Real dt_m, bool converged = false; while (!converged && iter < MAX_ITER) { + + // burn_state.y[] will always contain the current guess for the solution + // only the species and internal energy are updated though + int info = 0; - f_sdc_jac(dt_m, U_react, f, Jac, f_source, rho_new, T_old); + f_sdc_jac(dt_m, burn_state, f, Jac); // solve the linear system: Jac dU = -f #ifdef NEW_NETWORK_IMPLEMENTATION @@ -218,9 +185,10 @@ sdc_newton_solve(const Real dt_m, #endif // on output, f is the solution (dU) - for (int n = 1; n <= NumSpec+1; ++n) { - U_react(n) += f(n); + for (int n = 0; n < NumSpec; ++n) { + burn_state.y[SFS+n] += f(n + 1); } + burn_state.y[SEINT] += f(NumSpec+1); // compute the norm of the weighted error, where the // weights are 1/eps_tot @@ -231,10 +199,10 @@ sdc_newton_solve(const Real dt_m, if (n < NumSpec+1) { // for species, atol is the mass fraction limit, so we // multiply by density to get a partial density limit - eps = integrator_rp::rtol_spec * std::abs(U_react(n)) + + eps = integrator_rp::rtol_spec * std::abs(burn_state.y[SFS-1+n]) + integrator_rp::atol_spec * std::abs(U_new[URHO]); } else { - eps = integrator_rp::rtol_enuc * std::abs(U_react(NumSpec+1)) + + eps = integrator_rp::rtol_enuc * std::abs(burn_state.y[SEINT]) + integrator_rp::atol_enuc; } err_sum += f(n) * f(n) / (eps * eps); @@ -254,19 +222,17 @@ sdc_newton_solve(const Real dt_m, return; } -#endif - // update the full U_new for (int n = 0; n < NumSpec; ++n) { - U_new[UFS+n] = U_react(1+n); + U_new[UFS+n] = burn_state.y[SFS+n]; } auto v2 = 0.0_rt; for (int m = 0; m < 3; ++m) { v2 += U_new[UMX+m] * U_new[UMX+m]; } - U_new[UEINT] = U_react(NumSpec+1); + U_new[UEINT] = burn_state.y[SEINT]; U_new[UEDEN] = U_new[UEINT] + 0.5_rt * v2 / U_new[URHO]; } From 2938549b4edf6430e38b3a6b98d540b3c76ba4dd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 7 Oct 2023 11:23:26 -0400 Subject: [PATCH 2/2] move initial guess --- Source/sdc/sdc_newton_solve.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/sdc/sdc_newton_solve.H b/Source/sdc/sdc_newton_solve.H index c4b811af1f..03847ae105 100644 --- a/Source/sdc/sdc_newton_solve.H +++ b/Source/sdc/sdc_newton_solve.H @@ -134,9 +134,6 @@ sdc_newton_solve(const Real dt_m, copy_cons_to_burn_type(U_new, burn_state); burn_state.rho = U_new[URHO]; - // initial guess - burn_state.T = U_old[UTEMP]; - // for the Jacobian solve, we are solving // f(U) = U - dt R(U) - U_old - dt C = 0 // we define f_source = U_old + dt C so we are solving @@ -162,6 +159,9 @@ sdc_newton_solve(const Real dt_m, // burn_state.y[] will always contain the current guess for the solution // only the species and internal energy are updated though + // initial guess + burn_state.T = U_old[UTEMP]; + int info = 0; f_sdc_jac(dt_m, burn_state, f, Jac);