Robin BC: Abort if solver is not safe for reuse (#3788)

Because the coefficients stored inside LinOp objects are irreversibly
modified for Robin BC, it's unsafe to reuse the solver when there is
Robin BC, unless the scalars and the coefficients are reset.

---------

Co-authored-by: cgilet <cgilet@gmail.com>

Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.H

@@ -195,6 +195,9 @@ public:
     Vector<Vector<MF> > m_a_coeffs;
     Vector<Vector<Array<MF,AMREX_SPACEDIM> > > m_b_coeffs;
 
+    bool m_scalars_set = false;
+    bool m_acoef_set = false;
+
 protected:
 
     bool m_needs_update = true;
@@ -308,7 +311,9 @@ MLABecLaplacianT<MF>::setScalars (T1 a, T2 b) noexcept
         for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
             m_a_coeffs[amrlev][0].setVal(RT(0.0));
         }
+        m_acoef_set = true;
     }
+    m_scalars_set = true;
 }
 
 template <typename MF>
@@ -323,6 +328,7 @@ MLABecLaplacianT<MF>::setACoeffs (int amrlev, const AMF& alpha)
                               "MLABecLaplacian::setACoeffs: alpha is supposed to be single component.");
     m_a_coeffs[amrlev][0].LocalCopy(alpha, 0, 0, 1, IntVect(0));
     m_needs_update = true;
+    m_acoef_set = true;
 }
 
 template <typename MF>
@@ -333,6 +339,7 @@ MLABecLaplacianT<MF>::setACoeffs (int amrlev, T alpha)
 {
     m_a_coeffs[amrlev][0].setVal(RT(alpha));
     m_needs_update = true;
+    m_acoef_set = true;
 }
 
 
@@ -402,6 +409,8 @@ MLABecLaplacianT<MF>::update ()
     applyMetricTermsCoeffs();
 #endif
 
+    applyRobinBCTermsCoeffs();
+
     averageDownCoeffs();
 
     update_singular_flags();
@@ -488,6 +497,14 @@ void applyRobinBCTermsCoeffs (LP& linop)
     }
     const RT bovera = linop.m_b_scalar/linop.m_a_scalar;
 
+    if (!reset_alpha) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(linop.m_scalars_set && linop.m_acoef_set,
+                                         "To reuse solver With Robin BC, one must re-call setScalars (and setACoeffs if the scalar is not zero)");
+    }
+
+    linop.m_scalars_set = false;
+    linop.m_acoef_set = false;
+
     for (int amrlev = 0; amrlev < linop.NAMRLevels(); ++amrlev) {
         const int mglev = 0;
         const Box& domain = linop.Geom(amrlev,mglev).Domain();

Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.H

@@ -129,6 +129,9 @@ public:
     Vector<Vector<MultiFab> > m_a_coeffs;
     Vector<Vector<Array<MultiFab,AMREX_SPACEDIM> > > m_b_coeffs;
 
+    bool m_scalars_set = false;
+    bool m_acoef_set = false;
+
 protected:
 
     int m_ncomp = 1;

Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp

@@ -116,21 +116,25 @@ MLEBABecLap::setScalars (Real a, Real b)
         {
             m_a_coeffs[amrlev][0].setVal(0.0);
         }
+        m_acoef_set = true;
     }
+    m_scalars_set = true;
 }
 
 void
 MLEBABecLap::setACoeffs (int amrlev, const MultiFab& alpha)
 {
     MultiFab::Copy(m_a_coeffs[amrlev][0], alpha, 0, 0, 1, 0);
     m_needs_update = true;
+    m_acoef_set = true;
 }
 
 void
 MLEBABecLap::setACoeffs (int amrlev, Real alpha)
 {
     m_a_coeffs[amrlev][0].setVal(alpha);
     m_needs_update = true;
+    m_acoef_set = true;
 }
 
 void
@@ -1185,6 +1189,8 @@ MLEBABecLap::update ()
 {
     if (MLCellABecLap::needsUpdate()) { MLCellABecLap::update(); }
 
+    applyRobinBCTermsCoeffs();
+
     averageDownCoeffs();
 
     m_is_singular.clear();