Skip to content

Commit

Permalink
some changes to the SDC Newton solver to make it more robust (#2586)
Browse files Browse the repository at this point in the history 
we now reject a solve if the mass fractions are too far out of [0, 1]
also do the total energy update conservatively
  • Loading branch information
@zingale
zingale authored Oct 12, 2023
1 parent d2ace11 commit f8446cd
Show file tree
Hide file tree
Showing 4 changed files with 58 additions and 29 deletions.
18 changes: 9 additions & 9 deletions Exec/reacting_tests/reacting_convergence/ci-benchmarks/react_converge_128_true_sdc.out
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,16 +21,16 @@
rho_e_sdc_source -6.1984065959e+22 6.3082607681e+21
Temp_sdc_source 0 0
rho_He4_sdc_source -543231.29383 78510.973919
rho_C12_sdc_source -1694.4890575 10.073635142
rho_C12_sdc_source -1694.4890575 10.073635144
rho_O16_sdc_source -0.0098093776861 7.1506350196e-06
rho_Fe56_sdc_source -5.4492579275e-05 7.8513247889e-06
pressure 1.4155320805e+22 4.2608130858e+22
kineng 5.228354476e-13 1.6647276226e+18
kineng 3.6057617076e-14 1.6647276226e+18
soundspeed 212069503.63 257404342.21
Gamma_1 1.5601126452 1.5885713572
MachNumber 6.8192135042e-18 0.0086982771596
MachNumber 1.7908132003e-18 0.0086982771596
entropy 348439780.9 349209883.57
magvort 6.3108872418e-30 0.00018541051835
magvort 0 0.00018541051835
divu -0.1163387912 0.55816306007
eint_E 4.5262357143e+16 6.9937847678e+16
eint_e 4.5262357143e+16 6.9937843728e+16
Expand All @@ -49,11 +49,11 @@
z_velocity 0 0
t_sound_t_enuc 0.00023710316663 0.0195732693
enuc 2.9131490847e+15 4.5102586513e+17
magvel 1.446147223e-09 2067446.1363
radvel -0.00067837194793 2067446.1363
magvel 3.7977686648e-10 2067446.1363
radvel -0.00067839201193 2067446.1363
circvel 0 11.820144682
magmom 0.00072307361144 1.6422547233e+12
angular_momentum_x -0 -0
magmom 0.00018988843323 1.6422547233e+12
angular_momentum_x 0 0
angular_momentum_y 0 0
angular_momentum_z -1.2410862734e+14 1.2410862747e+14
angular_momentum_z -1.241086281e+14 1.2410862748e+14

16 changes: 14 additions & 2 deletions Source/sdc/Castro_sdc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,8 +208,6 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt)
{
normalize_species_sdc(i, j, k, U_center_arr);
});
// ca_normalize_species(AMREX_INT_ANYD(bx1.loVect()), AMREX_INT_ANYD(bx1.hiVect()),
// BL_TO_FORTRAN_ANYD(U_center));

// convert the C source to cell-centers
C_center.resize(bx1, NUM_STATE);
Expand All @@ -234,6 +232,13 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt)
sdc_update_centers_o4(i, j, k, U_center_arr, U_new_center_arr, C_center_arr, dt_m, sdc_iteration);
});

// enforce that the species sum to one after the reaction solve
amrex::ParallelFor(bx1,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
normalize_species_sdc(i, j, k, U_new_center_arr);
});

// compute R_i and in 1 ghost cell and then convert to <R> in
// place (only for the interior)
R_new.resize(bx1, NUM_STATE);
Expand Down Expand Up @@ -350,6 +355,13 @@ Castro::construct_old_react_source(MultiFab& U_state,

make_cell_center(obx, U_state.array(mfi), U_center_arr, domain_lo, domain_hi);

// sometimes the Laplacian can make the species go negative near discontinuities
amrex::ParallelFor(obx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
normalize_species_sdc(i, j, k, U_center_arr);
});

// burn, including one ghost cell
R_center.resize(obx, NUM_STATE);
Elixir elix_r_center = R_center.elixir();
Expand Down
11 changes: 2 additions & 9 deletions Source/sdc/Castro_sdc_util.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,21 +52,14 @@ sdc_solve(const Real dt_m,
int ierr;
Real err_out;

// for debugging
GpuArray<Real, NUM_STATE> U_orig;

for (int n = 0; n < NUM_STATE; ++n) {
U_orig[n] = U_old[n];
}

if (sdc_solver == NEWTON_SOLVE) {
// 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);

// failing?
if (ierr != NEWTON_SUCCESS) {
if (ierr != newton::NEWTON_SUCCESS) {
Abort("Newton subcycling failed in sdc_solve");
}
} else if (sdc_solver == VODE_SOLVE) {
Expand All @@ -85,7 +78,7 @@ sdc_solve(const Real dt_m,
sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr);

// Failing?
if (ierr != NEWTON_SUCCESS) {
if (ierr != newton::NEWTON_SUCCESS) {
Abort("Newton failure in sdc_solve");
}
}
Expand Down
42 changes: 33 additions & 9 deletions Source/sdc/sdc_newton_solve.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,14 @@
#include <linpack.H>

// error codes
constexpr int NEWTON_SUCCESS = 0;
constexpr int SINGULAR_MATRIX = -1;
constexpr int CONVERGENCE_FAILURE = -2;
namespace newton {
constexpr int NEWTON_SUCCESS = 0;
constexpr int SINGULAR_MATRIX = -1;
constexpr int CONVERGENCE_FAILURE = -2;
constexpr int BAD_MASS_FRACTIONS = -3;

constexpr Real species_failure_tolerance = 1.e-2_rt;
};

#ifdef REACTIONS

Expand Down Expand Up @@ -107,7 +111,7 @@ sdc_newton_solve(const Real dt_m,

const int MAX_ITER = 100;

ierr = NEWTON_SUCCESS;
ierr = newton::NEWTON_SUCCESS;

// update the density and momenta for this zone -- they don't react

Expand Down Expand Up @@ -162,7 +166,7 @@ sdc_newton_solve(const Real dt_m,
dgefa<NumSpec+1>(Jac, ipvt, info);
#endif
if (info != 0) {
ierr = SINGULAR_MATRIX;
ierr = newton::SINGULAR_MATRIX;
return;
}

Expand Down Expand Up @@ -206,7 +210,7 @@ sdc_newton_solve(const Real dt_m,
err_out = err;

if (!converged) {
ierr = CONVERGENCE_FAILURE;
ierr = newton::CONVERGENCE_FAILURE;
return;
}

Expand All @@ -221,7 +225,13 @@ sdc_newton_solve(const Real dt_m,
}

U_new[UEINT] = burn_state.y[SEINT];
U_new[UEDEN] = U_new[UEINT] + 0.5_rt * v2 / U_new[URHO];

// we want to do a conservative update for (rho E), so first figure out the
// energy generation rate

Real rho_Sdot = (U_new[UEINT] - U_old[UEINT]) / dt_m - C[UEINT];

U_new[UEDEN] = U_old[UEDEN] + dt_m * (C[UEDEN] + rho_Sdot);

}

Expand Down Expand Up @@ -251,13 +261,13 @@ sdc_newton_subdivide(const Real dt_m,
// use the old time solution.

int nsub = 1;
ierr = CONVERGENCE_FAILURE;
ierr = newton::CONVERGENCE_FAILURE;

for (int n = 0; n < NUM_STATE; ++n) {
U_begin[n] = U_old[n];
}

while (nsub < MAX_NSUB && ierr != NEWTON_SUCCESS) {
while (nsub < MAX_NSUB && ierr != newton::NEWTON_SUCCESS) {
if (nsub > 1) {
for (int n = 0; n < NUM_STATE; ++n) {
U_new[n] = U_old[n];
Expand All @@ -278,6 +288,20 @@ sdc_newton_subdivide(const Real dt_m,

sdc_newton_solve(dt_sub, U_begin, U_new, C, sdc_iteration, err_out, ierr);

// our solve may have resulted in mass fractions outside
// of [0, 1] -- reject if this is the case
for (int n = 0; n < NumSpec; ++n) {
if (U_new[UFS+n] < -newton::species_failure_tolerance * U_new[URHO] ||
U_new[UFS+n] > (1.0_rt + newton::species_failure_tolerance) * U_new[URHO]) {
ierr = newton::BAD_MASS_FRACTIONS;
}
}

if (ierr != newton::NEWTON_SUCCESS) {
// no point in continuing this subdivision if one stage failed
break;
}

for (int n = 0; n < NUM_STATE; ++n) {
U_begin[n] = U_new[n];
}
Expand Down

0 comments on commit f8446cd

Please sign in to comment.