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

some changes to the SDC Newton solver to make it more robust #2586

Merged
merged 16 commits into from
Oct 12, 2023
Merged

some changes to the SDC Newton solver to make it more robust #2586

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
da85e1a
some changes to the SDC Newton solver to make it more robust
zingale Sep 25, 2023
bbfffcd
fix check
zingale Sep 25, 2023
58efc8f
some more robustness
zingale Sep 28, 2023
7306214
fix
zingale Sep 28, 2023
4561f85
add another normalization
zingale Sep 29, 2023
1bba5fa
Merge branch 'development' into sdc_newton_robustness
zingale Oct 5, 2023
b1b30b7
Merge branch 'development' into sdc_newton_robustness
zingale Oct 6, 2023
4ad1016
Merge branch 'development' into sdc_newton_robustness
zingale Oct 7, 2023
aa90c9d
Merge branch 'development' into sdc_newton_robustness
zingale Oct 8, 2023
2b41b89
some reverts
zingale Oct 9, 2023
c2ab502
Merge branch 'development' into sdc_newton_robustness
zingale Oct 9, 2023
fc43f16
update benchmark
zingale Oct 9, 2023
b512a60
Merge branch 'development' into sdc_newton_robustness
zingale Oct 9, 2023
4bee43e
add a break
zingale Oct 12, 2023
c6dab80
Merge branch 'sdc_newton_robustness' of github.com:zingale/Castro int…
zingale Oct 12, 2023
1bbc944
fix location of check
zingale Oct 12, 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
18 changes: 9 additions & 9 deletions Exec/reacting_tests/reacting_convergence/ci-benchmarks/react_converge_128_true_sdc.out
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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 @@ -132,7 +136,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 @@ -195,7 +199,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 @@ -238,7 +242,7 @@ sdc_newton_solve(const Real dt_m,
err_out = err;

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

Expand All @@ -255,7 +259,13 @@ sdc_newton_solve(const Real dt_m,
}

U_new[UEINT] = U_react(NumSpec+1);
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 @@ -285,13 +295,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 @@ -312,6 +322,20 @@ sdc_newton_subdivide(const Real dt_m,

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

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

// 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;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does it makes sense to keep doing the sub-stepping if we hit this error case?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

oh, you mean I should do a break here?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Something like that, yes

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

break added

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well I'm not as worried about a break inside the species loop, I mean a break after the species loop before we finish up anything else we do below this

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah, okay

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just moved it to before the previous break. Thanks. I missed that.

}
}

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