From dee6ca80f492dbbe73991b69aff72310dde2d73c Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 30 Apr 2025 07:32:56 -0400 Subject: [PATCH] fix(operator): changed MULT (residuals) to deal with negated M properly --- src/poly/utils/private/operator.cpp | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp index 6dda6c2..d37e8df 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/operator.cpp @@ -26,9 +26,6 @@ #include "../../../../utils/debugUtils/MFEMAnalysisUtils/MFEMAnalysis-cpp/src/include/mfem_smout.h" -static int s_PolyOperatorMultCalledCount = 0; -static int s_PolyOperatorGradCalledCount = 0; - void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr& invMat, const std::string& name="matrix") { // PERF: This likely can be made much more efficient and will probably be called in tight loops, a good // PERF: place for some easy optimization might be here. @@ -99,7 +96,8 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) { // m_negM_op = std::make_unique(m_Mmat.get(), -1.0); // m_negQ_op = std::make_unique(m_Qmat.get(), -1.0); - m_Mmat->operator*=(-1.0); + // TODO Replace this with a scaled operator when I am done writing these out to disk + m_Mmat->operator*=(-1.0); // Sparse matrix only defines an inplace operator*= not an operator* so we need to do this m_Qmat->operator*=(-1.0); // Set up the constant parts of the jacobian now @@ -165,7 +163,7 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult"); m_f->Mult(x_theta, f_term); // fixme: this may be the wrong way to assemble m_f? - m_M->Mult(x_phi, Mphi_term); + m_M->Mult(x_phi, Mphi_term); // Does the order of operations work out between this and the next subtract to make Mφ negative properly? I think so but double check m_D->Mult(x_phi, Dphi_term); m_Q->Mult(x_theta, Qtheta_term); @@ -175,16 +173,16 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { // -- Apply essential boundary conditions -- for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) { if (int idx = m_theta_ess_tdofs.first[i]; idx >= 0 && idx < y_R0.Size()) { - // const double &targetValue = m_theta_ess_tdofs.second[i]; - // y_block.GetBlock(0)[idx] = targetValue - x_theta(idx); // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising - y_block.GetBlock(0)[idx] = 0; // Zero out the essential theta dofs in the bi-linear form // TODO Check if this is double zeroing (i.e if they were already removed maybe I am removing something that should not be removed here) + const double &targetValue = m_theta_ess_tdofs.second[i]; + y_block.GetBlock(0)[idx] = x_theta(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising } } std::cout << "θ Norm: " << y_block.GetBlock(0).Norml2() << std::endl; for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) { if (int idx = m_phi_ess_tdofs.first[i]; idx >= 0 && idx < y_R1.Size()) { - y_block.GetBlock(1)[idx] = 0; // Zero out the essential phi dofs in the bi-linear form // TODO Check if this is double zeroing (i.e if they were already removed maybe I am removing something that should not be removed here) + const double &targetValue = m_phi_ess_tdofs.second[i]; + y_block.GetBlock(1)[idx] = x_phi(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising } } @@ -205,14 +203,15 @@ void PolytropeOperator::updateInverseSchurCompliment() const { if (m_invNonlinearJacobian == nullptr) { MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian"); } - mfem::SparseMatrix* schurCompliment(&m_D->SpMat()); // This is now a copy of D + mfem::SparseMatrix* schurCompliment(&m_D->SpMat()); // Represents S in the preconditioner, starts as a copy of D before being modified in place // Calculate S = D - Q df^{-1} M - mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(m_Q->SpMat(), *m_invNonlinearJacobian); // Q * df^{-1} - const mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, m_M->SpMat()); // Q * df^{-1} * M + mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(*m_Qmat, *m_invNonlinearJacobian); // Q * df^{-1} + const mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, *m_Mmat); // Q * df^{-1} * M // PERF: Could potentially add some caching in here to only update the preconditioner when some condition has been met - schurCompliment->Add(-1, *temp_QxdfInvxM); // D - Q * df^{-1} * M + schurCompliment->Add(1, *temp_QxdfInvxM); // D - Q * df^{-1} * M + saveSparseMatrixBinary(*schurCompliment, "schurCompliment.bin"); approxJacobiInvert(*schurCompliment, m_invSchurCompliment, "Schur Compliment"); if (m_schurPreconditioner == nullptr) {