From 0380d0fcbbe05fa665b8389b1a36de4b76b87758 Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Fri, 5 Jul 2024 14:49:14 -0400 Subject: [PATCH 1/2] Avoid dividing by zero in constexpr linear algebra routines This fixes the crash in `Detonation-collision-retry-SDC` in the Castro test suite. --- integration/VODE/vode_dvjac.H | 3 +-- networks/rhs.H | 6 +++++- unit_test/test_linear_algebra/test_linear_algebra.H | 5 +++-- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/integration/VODE/vode_dvjac.H b/integration/VODE/vode_dvjac.H index 1e326aa9d5..82e016ebe9 100644 --- a/integration/VODE/vode_dvjac.H +++ b/integration/VODE/vode_dvjac.H @@ -169,8 +169,7 @@ void dvjac (int& IERPJ, BurnT& state, DvodeT& vstate) int IER{}; #ifdef NEW_NETWORK_IMPLEMENTATION - RHS::dgefa(vstate.jac); - IER = 0; + RHS::dgefa(vstate.jac, IER); #else if (integrator_rp::linalg_do_pivoting == 1) { constexpr bool allow_pivot{true}; diff --git a/networks/rhs.H b/networks/rhs.H index dddf396c61..9ca21cb1de 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -307,7 +307,7 @@ void dgesl (const RArray2D& a, RArray1D& b) } AMREX_GPU_HOST_DEVICE AMREX_INLINE -void dgefa (RArray2D& a) +void dgefa (RArray2D& a, int& info) { // LU factorization in-place without pivoting. @@ -317,6 +317,10 @@ void dgefa (RArray2D& a) [[maybe_unused]] constexpr int k = n1; // compute multipliers + if (a(k, k) == 0.0_rt) { + info = k; + return; // same as continue in a normal loop + } amrex::Real t = -1.0_rt / a(k,k); amrex::constexpr_for([&] (auto n2) diff --git a/unit_test/test_linear_algebra/test_linear_algebra.H b/unit_test/test_linear_algebra/test_linear_algebra.H index 0762163da3..485d65e141 100644 --- a/unit_test/test_linear_algebra/test_linear_algebra.H +++ b/unit_test/test_linear_algebra/test_linear_algebra.H @@ -101,7 +101,9 @@ void linear_algebra() { // solve the linear system with the templated solver - RHS::dgefa(A); + int info; + + RHS::dgefa(A, info); RHS::dgesl(A, b); std::cout << std::setprecision(16); @@ -122,7 +124,6 @@ void linear_algebra() { b = Ax(A, x); IArray1D pivot; - int info; constexpr bool allow_pivot{true}; From b779e957d4b300e724abbe2b9586476908b6dbed Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Sat, 6 Jul 2024 18:43:00 -0400 Subject: [PATCH 2/2] Have RHS::dgefa() return the info parameter directly --- integration/VODE/vode_dvjac.H | 2 +- networks/rhs.H | 6 +++++- unit_test/test_linear_algebra/test_linear_algebra.H | 5 ++--- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/integration/VODE/vode_dvjac.H b/integration/VODE/vode_dvjac.H index 82e016ebe9..ac66444490 100644 --- a/integration/VODE/vode_dvjac.H +++ b/integration/VODE/vode_dvjac.H @@ -169,7 +169,7 @@ void dvjac (int& IERPJ, BurnT& state, DvodeT& vstate) int IER{}; #ifdef NEW_NETWORK_IMPLEMENTATION - RHS::dgefa(vstate.jac, IER); + IER = RHS::dgefa(vstate.jac); #else if (integrator_rp::linalg_do_pivoting == 1) { constexpr bool allow_pivot{true}; diff --git a/networks/rhs.H b/networks/rhs.H index 9ca21cb1de..26e3c096e1 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -307,11 +307,13 @@ void dgesl (const RArray2D& a, RArray1D& b) } AMREX_GPU_HOST_DEVICE AMREX_INLINE -void dgefa (RArray2D& a, int& info) +int dgefa (RArray2D& a) { // LU factorization in-place without pivoting. + int info = 0; + amrex::constexpr_for<1, INT_NEQS>([&] (auto n1) { [[maybe_unused]] constexpr int k = n1; @@ -344,6 +346,8 @@ void dgefa (RArray2D& a, int& info) }); }); }); + + return info; } // Calculate the density dependence term for tabulated rates. The RHS has a term diff --git a/unit_test/test_linear_algebra/test_linear_algebra.H b/unit_test/test_linear_algebra/test_linear_algebra.H index 485d65e141..0762163da3 100644 --- a/unit_test/test_linear_algebra/test_linear_algebra.H +++ b/unit_test/test_linear_algebra/test_linear_algebra.H @@ -101,9 +101,7 @@ void linear_algebra() { // solve the linear system with the templated solver - int info; - - RHS::dgefa(A, info); + RHS::dgefa(A); RHS::dgesl(A, b); std::cout << std::setprecision(16); @@ -124,6 +122,7 @@ void linear_algebra() { b = Ax(A, x); IArray1D pivot; + int info; constexpr bool allow_pivot{true};