From f5dea85db1bb47e66d811279dc54765838700b62 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 22 May 2025 11:30:24 -0400 Subject: [PATCH] fix(PolytropeOperator): seperated boundary aware and unaware operators for M, Q, and D residual calculation needs to be done when boundary degrees of freedom have not been removed (since their removal takes place in the Mult step in order to introduce the proper restoring force). Whereas jacobian calculation needs to always work from the boundary aware operators (these are sparse matrices) The other thing to note here is that this seems contrary to a matrix free design. While true it is common practive to assemble linear terms since they are cheap. We still never assemble the non linear matrix form. --- src/poly/utils/private/operator.cpp | 48 +++++++++++++++++------------ src/poly/utils/public/operator.h | 46 ++++++++++++++------------- 2 files changed, 53 insertions(+), 41 deletions(-) diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp index 16c659b..f4f23c2 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/operator.cpp @@ -54,24 +54,34 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) { return; } - m_M -> EliminateTestEssentialBC(m_theta_ess_tdofs.first); - m_M -> EliminateTrialEssentialBC(m_phi_ess_tdofs.first); + m_Mmat = std::make_unique(m_M->SpMat()); + m_Qmat = std::make_unique(m_Q->SpMat()); + m_Dmat = std::make_unique(m_D->SpMat()); - m_Q -> EliminateTestEssentialBC(m_phi_ess_tdofs.first); - m_Q -> EliminateTrialEssentialBC(m_theta_ess_tdofs.first); + // Remove the essential dofs from the constant matrices + for (const auto& dof: m_theta_ess_tdofs.first) { + m_Mmat->EliminateRow(dof); + m_Qmat->EliminateCol(dof); + } - m_D -> EliminateEssentialBC(m_phi_ess_tdofs.first); + for (const auto& dof: m_phi_ess_tdofs.first) { + m_Mmat->EliminateCol(dof); + m_Qmat->EliminateRow(dof); + m_Dmat->EliminateRowCol(dof); + } - m_negM_op = std::make_unique(m_M.get(), -1.0); - m_negQ_op = std::make_unique(m_Q.get(), -1.0); + m_negM_mat = std::make_unique(m_Mmat.get(), -1.0); + m_negQ_mat = std::make_unique(m_Qmat.get(), -1.0); - m_schurCompliment = std::make_unique(*m_Q, *m_D, *m_M); + m_schurCompliment = std::make_unique(*m_Qmat, *m_Dmat, *m_Mmat); // Set up the constant parts of the jacobian now + + // We use the assembled matrices with their boundary conditions already removed for the jacobian m_jacobian = std::make_unique(m_blockOffsets); - m_jacobian->SetBlock(0, 1, m_negM_op.get()); //<- -M (constant) - m_jacobian->SetBlock(1, 0, m_negQ_op.get()); //<- -Q (constant) - m_jacobian->SetBlock(1, 1, m_D.get()); //<- D (constant) + m_jacobian->SetBlock(0, 1, m_negM_mat.get()); //<- -M (constant) + m_jacobian->SetBlock(1, 0, m_negQ_mat.get()); //<- -Q (constant) + m_jacobian->SetBlock(1, 1, m_Dmat.get()); //<- D (constant) m_invSchurCompliment = std::make_unique(*m_schurCompliment); @@ -240,9 +250,9 @@ void GMRESInverter::Mult(const mfem::Vector &x, mfem::Vector &y) const { SchurCompliment::SchurCompliment( - const mfem::MixedBilinearForm &QOp, - const mfem::BilinearForm &DOp, - const mfem::MixedBilinearForm &MOp, + const mfem::SparseMatrix &QOp, + const mfem::SparseMatrix &DOp, + const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp) : mfem::Operator(DOp.Height(), DOp.Width()) { @@ -252,9 +262,9 @@ mfem::Operator(DOp.Height(), DOp.Width()) } SchurCompliment::SchurCompliment( - const mfem::MixedBilinearForm &QOp, - const mfem::BilinearForm &DOp, - const mfem::MixedBilinearForm &MOp) : + const mfem::SparseMatrix &QOp, + const mfem::SparseMatrix &DOp, + const mfem::SparseMatrix &MOp) : mfem::Operator(DOp.Height(), DOp.Width()) { updateConstantTerms(QOp, DOp, MOp); @@ -262,7 +272,7 @@ mfem::Operator(DOp.Height(), DOp.Width()) m_nTheta = m_MOp->Height(); } -void SchurCompliment::SetOperator(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp, const mfem::Solver &GradInvOp) { +void SchurCompliment::SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp) { updateConstantTerms(QOp, DOp, MOp); updateInverseNonlinearJacobian(GradInvOp); } @@ -271,7 +281,7 @@ void SchurCompliment::updateInverseNonlinearJacobian(const mfem::Solver &gradInv m_GradInvOp = &gradInv; } -void SchurCompliment::updateConstantTerms(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp) { +void SchurCompliment::updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp) { m_QOp = &QOp; m_DOp = &DOp; m_MOp = &MOp; diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h index 1cffe05..051845f 100644 --- a/src/poly/utils/public/operator.h +++ b/src/poly/utils/public/operator.h @@ -28,22 +28,23 @@ class SchurCompliment final : public mfem::Operator { public: - SchurCompliment(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp, const mfem::Solver &GradInvOp); - SchurCompliment(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp); + SchurCompliment(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp); + SchurCompliment(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp); ~SchurCompliment() override = default; void Mult(const mfem::Vector &x, mfem::Vector &y) const override; - void SetOperator(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp, const mfem::Solver &GradInvOp); + void SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver & + GradInvOp); void updateInverseNonlinearJacobian(const mfem::Solver &gradInv); private: - void updateConstantTerms(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp); + void updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp); private: // Note that these are not owned by this class - const mfem::MixedBilinearForm* m_QOp = nullptr; - const mfem::BilinearForm* m_DOp = nullptr; - const mfem::MixedBilinearForm* m_MOp = nullptr; + const mfem::SparseMatrix* m_QOp = nullptr; + const mfem::SparseMatrix* m_DOp = nullptr; + const mfem::SparseMatrix* m_MOp = nullptr; const mfem::Solver* m_GradInvOp = nullptr; int m_nPhi = 0; int m_nTheta = 0; @@ -101,32 +102,33 @@ private: std::unique_ptr m_D; std::unique_ptr m_f; + // These are used to store the matrix representations + // for the bi-linear forms. This might seem counterintuitive + // for a matrix-free approach. However, these will be computed + // regardless due to MFEM's implementation of these operators. + // Further since these do not change it is not a performance issue. + + // We need to store these separately because the jacobian and preconditioner + // must be computed from the boundary aware operators (which will be these + // matrices) whereas the residuals must be computed from the raw, physical, + // operators + std::unique_ptr m_Mmat; + std::unique_ptr m_Qmat; + std::unique_ptr m_Dmat; + const mfem::Array m_blockOffsets; SSE::MFEMArrayPair m_theta_ess_tdofs; SSE::MFEMArrayPair m_phi_ess_tdofs; - // std::unique_ptr m_Mmat, m_Qmat, m_Dmat; - std::unique_ptr m_negM_op; - std::unique_ptr m_negQ_op; + std::unique_ptr m_negM_mat; + std::unique_ptr m_negQ_mat; mutable std::unique_ptr m_jacobian; - // mutable std::unique_ptr m_invSchurCompliment; - // mutable std::unique_ptr m_invNonlinearJacobian; mutable std::unique_ptr m_schurCompliment; mutable std::unique_ptr m_invSchurCompliment; mutable std::unique_ptr m_invNonlinearJacobian; - /* - * The schur preconditioner has the form - * - * ⎡ḟ(θ)^-1 0 ⎤ - * ⎣ 0 S^-1 ⎦ - * - * Where S is the Schur compliment of the system - * - */ - mutable std::unique_ptr m_schurPreconditioner; bool m_isFinalized = false;