diff --git a/Source/FieldSolver/ImplicitSolvers/DarwinLinearFieldOperator.H b/Source/FieldSolver/ImplicitSolvers/DarwinLinearFieldOperator.H index 3976e0cebb1..b2bc3a1a748 100644 --- a/Source/FieldSolver/ImplicitSolvers/DarwinLinearFieldOperator.H +++ b/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& a_ridx_g, diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H index 62a4329d616..8300a2eb5ed 100644 --- a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H +++ b/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: diff --git a/Source/NonlinearSolvers/CurlCurlMLMGPC.H b/Source/NonlinearSolvers/CurlCurlMLMGPC.H index 514d97a1b7a..8ba6262426d 100644 --- a/Source/NonlinearSolvers/CurlCurlMLMGPC.H +++ b/Source/NonlinearSolvers/CurlCurlMLMGPC.H @@ -81,7 +81,7 @@ class CurlCurlMLMGPC : public Preconditioner /** * \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::Define ( const T& a_U, } template -void CurlCurlMLMGPC::Update (const T& a_U) +void CurlCurlMLMGPC::Update () { BL_PROFILE("CurlCurlMLMGPC::Update()"); using namespace amrex; @@ -256,9 +256,6 @@ void CurlCurlMLMGPC::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.) { diff --git a/Source/NonlinearSolvers/JacobiPC.H b/Source/NonlinearSolvers/JacobiPC.H index c669c295dcf..b934ea8cda9 100644 --- a/Source/NonlinearSolvers/JacobiPC.H +++ b/Source/NonlinearSolvers/JacobiPC.H @@ -72,7 +72,7 @@ class JacobiPC : public Preconditioner 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 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> m_grids; + amrex::Vector> m_dmap; + const amrex::Vector>* m_bcoefs = nullptr; bool m_has_offdiag = false; @@ -269,13 +276,25 @@ void JacobiPC::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 -void JacobiPC::Update (const T& a_U) +void JacobiPC::Update () { BL_PROFILE("JacobiPC::Update()"); using namespace amrex; @@ -285,7 +304,6 @@ void JacobiPC::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::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); } } } diff --git a/Source/NonlinearSolvers/JacobianFunctionMF.H b/Source/NonlinearSolvers/JacobianFunctionMF.H index af9d3b2e999..5a34257ea27 100644 --- a/Source/NonlinearSolvers/JacobianFunctionMF.H +++ b/Source/NonlinearSolvers/JacobianFunctionMF.H @@ -52,9 +52,9 @@ class JacobianFunctionMF : public LinearFunction } 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 diff --git a/Source/NonlinearSolvers/LinearFunction.H b/Source/NonlinearSolvers/LinearFunction.H index 41cd71ab308..eaaf44b9d6f 100644 --- a/Source/NonlinearSolvers/LinearFunction.H +++ b/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&, diff --git a/Source/NonlinearSolvers/MatrixPC.H b/Source/NonlinearSolvers/MatrixPC.H index ad9d431d442..573a4a32368 100644 --- a/Source/NonlinearSolvers/MatrixPC.H +++ b/Source/NonlinearSolvers/MatrixPC.H @@ -104,7 +104,7 @@ class MatrixPC : public Preconditioner /** * \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 * 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 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::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::Define ( const T& a_U, } template -void MatrixPC::Update (const T& a_U) +void MatrixPC::Update () { BL_PROFILE("MatrixPC::Update()"); using namespace amrex; @@ -293,7 +291,7 @@ void MatrixPC::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::Update (const T& a_U) } template -int MatrixPC::Assemble (const T& a_U) +int MatrixPC::Assemble () { // Assemble the sparse matrix representation of the preconditioner // A = curl (alpha * curl []) + M @@ -347,11 +345,8 @@ int MatrixPC::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(); diff --git a/Source/NonlinearSolvers/NewtonSolver.H b/Source/NonlinearSolvers/NewtonSolver.H index 4debb5c9429..1ed7b40afa2 100644 --- a/Source/NonlinearSolvers/NewtonSolver.H +++ b/Source/NonlinearSolvers/NewtonSolver.H @@ -388,7 +388,7 @@ void NewtonSolver::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(); diff --git a/Source/NonlinearSolvers/Preconditioner.H b/Source/NonlinearSolvers/Preconditioner.H index 1a2ec025bcb..0822744be05 100644 --- a/Source/NonlinearSolvers/Preconditioner.H +++ b/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 diff --git a/Source/NonlinearSolvers/WarpX_PETSc.cpp b/Source/NonlinearSolvers/WarpX_PETSc.cpp index 124649bf6c0..1fb3b08a610 100644 --- a/Source/NonlinearSolvers/WarpX_PETSc.cpp +++ b/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*>(context->m_linop.get())->updatePreCondMat(context->m_U); + dynamic_cast*>(context->m_linop.get())->updatePreCondMat(); } PetscFunctionReturn(PETSC_SUCCESS); }