Cleanup: remove unneeded argument for Preconditioner::Update

Source/FieldSolver/ImplicitSolvers/DarwinLinearFieldOperator.H

@@ -69,10 +69,7 @@ public:
     }
 
     inline
-    void updatePreCondMat ( const WarpXSolverVec& a_X ) override
-    {
-        amrex::ignore_unused(a_X);
-    }
+    void updatePreCondMat () override { }
 
     inline
     void getPCMatrix ( amrex::Gpu::DeviceVector<int>& a_ridx_g,

Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H

@@ -314,7 +314,7 @@ public:
     [[nodiscard]] auto numAMRLevels () const { return m_num_amr_levels; }
 
     // return DOFs object pointer
-    [[nodiscard]] inline const auto& getDOFsObject () const { return m_dofs; }
+    [[nodiscard]] inline static const auto& getDOFsObject () { return m_dofs; }
 
 private:
 

Source/NonlinearSolvers/CurlCurlMLMGPC.H

@@ -81,7 +81,7 @@ class CurlCurlMLMGPC : public Preconditioner<T,Ops>
         /**
          * \brief Update the preconditioner
          */
-        void Update (const T& a_U) override;
+        void Update () override;
 
         /**
          * \brief Apply (solve) the preconditioner given a RHS
@@ -247,7 +247,7 @@ void CurlCurlMLMGPC<T,Ops>::Define ( const T& a_U,
 }
 
 template <class T, class Ops>
-void CurlCurlMLMGPC<T,Ops>::Update (const T& a_U)
+void CurlCurlMLMGPC<T,Ops>::Update ()
 {
     BL_PROFILE("CurlCurlMLMGPC::Update()");
     using namespace amrex;
@@ -256,9 +256,6 @@ void CurlCurlMLMGPC<T,Ops>::Update (const T& a_U)
         IsDefined(),
         "CurlCurlMLMGPC::Update() called on undefined object" );
 
-    // a_U is not needed for a linear operator
-    amrex::ignore_unused(a_U);
-
     // set the alpha coefficient for the curl-curl op
     const RT thetaDt = m_ops->GetThetaForPC()*this->m_dt;
     if (thetaDt==0.) {

Source/NonlinearSolvers/JacobiPC.H

@@ -72,7 +72,7 @@ class JacobiPC : public Preconditioner<T,Ops>
 
         void Define (const T&, Ops*) override;
 
-        void Update (const T& a_U) override;
+        void Update () override;
 
         /**
          * \brief Solve (I + M) x = b via damped Jacobi iteration
@@ -110,6 +110,13 @@ class JacobiPC : public Preconditioner<T,Ops>
 
         int m_num_amr_levels = 0;
 
+        /**
+         * \brief Grid layout of the solver vector, cached in Define() so that
+         *  Update() can allocate its scratch without being handed a vector.
+         */
+        amrex::Vector<amrex::Array<amrex::BoxArray,3>> m_grids;
+        amrex::Vector<amrex::Array<amrex::DistributionMapping,3>> m_dmap;
+
         const amrex::Vector<amrex::Array<amrex::MultiFab*,3>>* m_bcoefs = nullptr;
 
         bool m_has_offdiag = false;
@@ -269,13 +276,25 @@ void JacobiPC<T,Ops>::Define ( const T& a_U,
     m_num_amr_levels = m_ops->numAMRLevels();
     m_bcoefs = m_ops->GetMassMatricesCoeff();
 
+    // Save the grid layout of the solver vector. This is then used in
+    // `Update` to allocate temporary MultiFabs
+    const auto& u_mfarrvec = a_U.getArrayVec();
+    m_grids.resize(m_num_amr_levels);
+    m_dmap.resize(m_num_amr_levels);
+    for (int n = 0; n < m_num_amr_levels; n++) {
+        for (int dim = 0; dim < 3; dim++) {
+            m_grids[n][dim] = u_mfarrvec[n][dim]->boxArray();
+            m_dmap[n][dim]  = u_mfarrvec[n][dim]->DistributionMap();
+        }
+    }
+
     readParameters();
 
     m_is_defined = true;
 }
 
 template <class T, class Ops>
-void JacobiPC<T,Ops>::Update (const T& a_U)
+void JacobiPC<T,Ops>::Update ()
 {
     BL_PROFILE("JacobiPC::Update()");
     using namespace amrex;
@@ -285,7 +304,6 @@ void JacobiPC<T,Ops>::Update (const T& a_U)
         "JacobiPC::Update() called on undefined object" );
 
     if (m_bcoefs != nullptr && !m_work_defined) {
-        auto& u_mfarrvec = a_U.getArrayVec();
         m_work.resize(m_num_amr_levels);
         m_x_ghost.resize(m_num_amr_levels);
 
@@ -307,15 +325,10 @@ void JacobiPC<T,Ops>::Update (const T& a_U)
         const amrex::IntVect nghost(m_stencil_width);
         for (int n = 0; n < m_num_amr_levels; n++) {
             for (int dim = 0; dim < 3; dim++) {
-                m_work[n][dim].define(
-                    u_mfarrvec[n][dim]->boxArray(),
-                    u_mfarrvec[n][dim]->DistributionMap(),
-                    1, 0);
+                m_work[n][dim].define(m_grids[n][dim], m_dmap[n][dim], 1, 0);
                 if (m_has_offdiag) {
-                    m_x_ghost[n][dim].define(
-                        u_mfarrvec[n][dim]->boxArray(),
-                        u_mfarrvec[n][dim]->DistributionMap(),
-                        1, nghost);
+                    m_x_ghost[n][dim].define(m_grids[n][dim], m_dmap[n][dim],
+                                             1, nghost);
                 }
             }
         }

Source/NonlinearSolvers/JacobianFunctionMF.H

@@ -52,9 +52,9 @@ class JacobianFunctionMF : public LinearFunction<T,Ops>
     }
 
     inline
-    void updatePreCondMat ( const T&  a_X ) override
+    void updatePreCondMat () override
     {
-        if (m_usePreCond) { m_preCond->Update(a_X); }
+        if (m_usePreCond) { m_preCond->Update(); }
     }
 
     inline

Source/NonlinearSolvers/LinearFunction.H

@@ -58,7 +58,7 @@ class LinearFunction
         virtual void precond ( T& a_U, const T& a_X ) = 0;
 
         //! update preconditioner
-        virtual void updatePreCondMat ( const T&  a_X ) = 0;
+        virtual void updatePreCondMat () = 0;
 
         //! get sparse matrix representation of preconditioner
         virtual void getPCMatrix( amrex::Gpu::DeviceVector<int>&,

Source/NonlinearSolvers/MatrixPC.H

@@ -104,7 +104,7 @@ class MatrixPC : public Preconditioner<T,Ops>
         /**
          * \brief Update the preconditioner
          */
-        void Update (const T& a_U) override;
+        void Update () override;
 
         /**
          * \brief Assemble the matrix
@@ -114,7 +114,7 @@ class MatrixPC : public Preconditioner<T,Ops>
          * non-zero elements (return value is difference in current number
          * of nonzero elements and the desired number).
          */
-        int Assemble (const T& a_U);
+        int Assemble ();
 
         /**
          * \brief Apply (solve) the preconditioner given a RHS
@@ -185,6 +185,7 @@ class MatrixPC : public Preconditioner<T,Ops>
 
         int m_ndofs_l = 0;
         int m_ndofs_g = 0;
+
         bool m_pc_diag_only = false;
         int m_pc_mat_nnz = 1;
         bool m_include_mass_matrices = false;
@@ -247,9 +248,6 @@ void MatrixPC<T,Ops>::Define ( const T& a_U,
     // read preconditioner parameters
     readParameters();
 
-    // a_U is not needed
-    amrex::ignore_unused(a_U);
-
     // Set number of AMR levels and create geometry, grids, and
     // distribution mapping vectors.
     m_num_amr_levels = m_ops->numAMRLevels();
@@ -282,7 +280,7 @@ void MatrixPC<T,Ops>::Define ( const T& a_U,
 }
 
 template <class T, class Ops>
-void MatrixPC<T,Ops>::Update (const T& a_U)
+void MatrixPC<T,Ops>::Update ()
 {
     BL_PROFILE("MatrixPC::Update()");
     using namespace amrex;
@@ -293,7 +291,7 @@ void MatrixPC<T,Ops>::Update (const T& a_U)
 
     while(true) {
 
-        auto nnz_diff = Assemble(a_U);
+        auto nnz_diff = Assemble();
         AMREX_ALWAYS_ASSERT(nnz_diff >= 0);
         if (nnz_diff) {
 
@@ -315,7 +313,7 @@ void MatrixPC<T,Ops>::Update (const T& a_U)
 }
 
 template <class T, class Ops>
-int MatrixPC<T,Ops>::Assemble (const T& a_U)
+int MatrixPC<T,Ops>::Assemble ()
 {
     // Assemble the sparse matrix representation of the preconditioner
     //      A = curl (alpha * curl []) + M
@@ -347,11 +345,8 @@ int MatrixPC<T,Ops>::Assemble (const T& a_U)
                        << "alpha = " << alpha << "\n";
     }
 
-    // Get DOF object from a_U
-    const auto& dofs_obj = a_U.getDOFsObject();
-    const auto& dofs_mfarrvec = dofs_obj->m_array;
-    AMREX_ALWAYS_ASSERT(m_ndofs_l == dofs_obj->m_nDoFs_l);
-    AMREX_ALWAYS_ASSERT(m_ndofs_g == dofs_obj->m_nDoFs_g);
+    // The DOF object is shared by every solver vector.
+    const auto& dofs_mfarrvec = T::getDOFsObject()->m_array;
 
     m_r_indices_g.clear();
     m_num_nz.clear();

Source/NonlinearSolvers/NewtonSolver.H

@@ -388,7 +388,7 @@ void NewtonSolver<Vec,Ops>::Solve (Vec&         a_U,
         m_ops->PreLinearSolve();
         m_linear_function->setBaseSolution(a_U);
         m_linear_function->setBaseRHS(m_R);
-        m_linear_function->updatePreCondMat(a_U);
+        m_linear_function->updatePreCondMat();
 
         // Solve linear system for Newton step [Jac]*dU = F
         m_dU.zero();

Source/NonlinearSolvers/Preconditioner.H

@@ -58,7 +58,7 @@ class Preconditioner
         /**
          * \brief Update the preconditioner
          */
-        virtual void Update ( const T& a_U ) = 0;
+        virtual void Update () = 0;
 
         /**
          * \brief Apply (solve) the preconditioner given a RHS

Source/NonlinearSolvers/WarpX_PETSc.cpp

@@ -106,7 +106,7 @@ PetscErrorCode RHSFunction( SNES a_solver, Vec a_U, Vec a_F, void* ctxt)
     VecAXPBY(a_F, 1.0, -1.0, a_U);
 
     if (!context->m_fd_jac_comput) {
-        dynamic_cast<JacobianFunctionMF<VecType,TIType>*>(context->m_linop.get())->updatePreCondMat(context->m_U);
+        dynamic_cast<JacobianFunctionMF<VecType,TIType>*>(context->m_linop.get())->updatePreCondMat();
     }
     PetscFunctionReturn(PETSC_SUCCESS);
 }