Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

switch true-SDC to completely use burn_t internally in the Newton solve #2605

Merged
merged 5 commits into from
Oct 9, 2023
Merged

switch true-SDC to completely use burn_t internally in the Newton solve #2605

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
a2994b1
switch true-SDC to completely use burn_t internally in the Newton solve
zingale Oct 7, 2023
2938549
move initial guess
zingale Oct 7, 2023
060404f
Merge branch 'development' into switch_newton_to_burn_t
zingale Oct 7, 2023
bc73dd3
Merge branch 'development' into switch_newton_to_burn_t
zingale Oct 8, 2023
3c10d2b
Merge branch 'development' into switch_newton_to_burn_t
zingale Oct 9, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
88 changes: 27 additions & 61 deletions Source/sdc/sdc_newton_solve.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,9 @@ constexpr int CONVERGENCE_FAILURE = -2;
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void
f_sdc_jac(const Real dt_m,
Array1D<Real, 1, NumSpec+1> const& U,
burn_t& burn_state,
Array1D<Real, 1, NumSpec+1>& f,
JacNetArray2D& Jac,
Array1D<Real, 1, NumSpec+1>& f_source,
const Real rho,
const Real T_old) {
JacNetArray2D& Jac) {

// This is used with the Newton solve and returns f and the Jacobian

Expand All @@ -29,30 +26,10 @@ 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;

// 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] = std::max(network_rp::small_x, std::min(1.0_rt, burn_state.y[SFS+n] / burn_state.y[SRHO]));
}
Expand All @@ -71,9 +48,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.

Expand Down Expand Up @@ -126,8 +103,6 @@ sdc_newton_solve(const Real dt_m,
// 1:NumSpec : species
// NumSpec+1 : (rho e)

Array1D<Real, 1, NumSpec+1> U_react;
Array1D<Real, 1, NumSpec+1> f_source;
Array1D<Real, 1, NumSpec+1> f;

const int MAX_ITER = 100;
Expand All @@ -142,37 +117,22 @@ 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];

// 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

Expand All @@ -183,8 +143,15 @@ 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

// initial guess
burn_state.T = U_old[UTEMP];

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
Expand All @@ -206,9 +173,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
Expand All @@ -219,10 +187,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);
Expand All @@ -242,19 +210,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];

}
Expand Down