From 2acc0371110659a5158c32f58956c1b3f51d781b Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 25 Apr 2025 10:32:06 -0400 Subject: [PATCH] fix(poly): phi boundary values now coorespond to theta flux through polytrope surface --- src/poly/coeff/private/polyCoeff.cpp | 10 ++ src/poly/coeff/public/polyCoeff.h | 32 ++++--- src/poly/solver/private/polySolver.cpp | 66 +++++++------ src/poly/solver/public/polySolver.h | 2 + src/poly/utils/private/integrators.cpp | 2 + src/poly/utils/private/operator.cpp | 124 ++++++++++++++----------- src/poly/utils/public/operator.h | 4 - 7 files changed, 146 insertions(+), 94 deletions(-) diff --git a/src/poly/coeff/private/polyCoeff.cpp b/src/poly/coeff/private/polyCoeff.cpp index c018fd8..d0d84c9 100644 --- a/src/poly/coeff/private/polyCoeff.cpp +++ b/src/poly/coeff/private/polyCoeff.cpp @@ -38,6 +38,8 @@ namespace polycoeff{ // v = -std::pow(r, 2); } + //PERF: One area of future optimization might be turning these into matrix operations since they are + //PERF: fundamentally just the dot product of [a b c ...] * [n^0 n^1 n^2 ...] double x1(const double n) { double r = 0; @@ -47,4 +49,12 @@ namespace polycoeff{ return r; } + double thetaSurfaceFlux(const double n) { + double surfaceFlux = 0; + for (int i = 0; i < dThetaInterpCoeff::numTerms; i++) { + surfaceFlux += dThetaInterpCoeff::coeff[i] * std::pow(n, dThetaInterpCoeff::power[i]); + } + return surfaceFlux; + } + } \ No newline at end of file diff --git a/src/poly/coeff/public/polyCoeff.h b/src/poly/coeff/public/polyCoeff.h index 3933a63..39393b6 100644 --- a/src/poly/coeff/public/polyCoeff.h +++ b/src/poly/coeff/public/polyCoeff.h @@ -42,15 +42,25 @@ namespace polycoeff void diffusionCoeff(const mfem::Vector &x, mfem::Vector &v); double x1(const double n); - /** - * @brief Coefficients for the interpolations of the surface location of a polytrope - * @param numTerms Number of terms in the polynomial interpolator - * @param power Array of the powers of the polynomial interpolator - * @param coeff Array of the coefficients of the polynomial interpolator - */ - struct x1InterpCoeff { - constexpr static int numTerms = 51; - constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50}; - constexpr static double coeff[51] = {2.449260049071003e+00, 4.340987403815444e-01, 2.592462347475860e+00, -2.283794512491552e+01, 1.135535208066129e+02, -3.460093673587010e+02, 6.948139716983651e+02, -9.564368645158952e+02, 9.184546798891624e+02, -6.130308987803503e+02, 2.747228735318193e+02, -7.416795903196430e+01, 7.318638538580859e+00, 1.749441306797260e+00, -4.240148582456829e-01, -5.818809544982156e-02, 1.514877217199105e-02, 3.228634707578998e-03, -2.862524323980516e-04, -1.622486968261819e-04, -1.253644717104076e-05, 4.334945141292894e-06, 1.296452565229763e-06, 8.634802209209870e-08, -3.337511676486084e-08, -1.094796628367775e-08, -1.228178760540410e-09, 1.416744125622751e-10, 8.513777351265677e-11, 1.624582561364811e-11, 8.377207519041114e-13, -4.363812865112836e-13, -1.535862757816461e-13, -2.485085045037669e-14, -6.566281276491033e-16, 8.405047965853478e-16, 2.673804441025638e-16, 4.176337890285142e-17, 8.150073570140493e-19, -1.531673805257016e-18, -4.746996933716653e-19, -6.976825828390195e-20, 8.807513368604331e-22, 3.307739310180278e-21, 8.593260940093030e-22, 7.385969061093440e-23, -2.211130577977291e-23, -8.557291048388455e-24, -4.169359901215994e-25, 4.609379358875657e-25, -3.590870335035984e-26}; - }; + + double thetaSurfaceFlux(const double n); + + /** + * @brief Coefficients for the interpolations of the surface location of a polytrope + * @param numTerms Number of terms in the polynomial interpolator + * @param power Array of the powers of the polynomial interpolator + * @param coeff Array of the coefficients of the polynomial interpolator + */ + const struct x1InterpCoeff { + constexpr static int numTerms = 51; + constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50}; + constexpr static double coeff[51] = {2.449260049071003e+00, 4.340987403815444e-01, 2.592462347475860e+00, -2.283794512491552e+01, 1.135535208066129e+02, -3.460093673587010e+02, 6.948139716983651e+02, -9.564368645158952e+02, 9.184546798891624e+02, -6.130308987803503e+02, 2.747228735318193e+02, -7.416795903196430e+01, 7.318638538580859e+00, 1.749441306797260e+00, -4.240148582456829e-01, -5.818809544982156e-02, 1.514877217199105e-02, 3.228634707578998e-03, -2.862524323980516e-04, -1.622486968261819e-04, -1.253644717104076e-05, 4.334945141292894e-06, 1.296452565229763e-06, 8.634802209209870e-08, -3.337511676486084e-08, -1.094796628367775e-08, -1.228178760540410e-09, 1.416744125622751e-10, 8.513777351265677e-11, 1.624582561364811e-11, 8.377207519041114e-13, -4.363812865112836e-13, -1.535862757816461e-13, -2.485085045037669e-14, -6.566281276491033e-16, 8.405047965853478e-16, 2.673804441025638e-16, 4.176337890285142e-17, 8.150073570140493e-19, -1.531673805257016e-18, -4.746996933716653e-19, -6.976825828390195e-20, 8.807513368604331e-22, 3.307739310180278e-21, 8.593260940093030e-22, 7.385969061093440e-23, -2.211130577977291e-23, -8.557291048388455e-24, -4.169359901215994e-25, 4.609379358875657e-25, -3.590870335035984e-26}; + }; + + const struct dThetaInterpCoeff { + constexpr static int numTerms = 51; + constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50}; + constexpr static double coeff[51] = {-8.164634353536671e-01, 8.651804355981230e-01, -6.391939109867470e-01, 4.882118355278944e-01, -3.898475553686603e-01, 3.102560043891329e-01, -2.365565113144737e-01, 1.639832755778444e-01, -9.645271947133248e-02, 4.450140636696046e-02, -1.460885194751961e-02, 2.816293786806504e-03, -8.178583546493117e-05, -8.584494556484958e-05, 8.476252127593868e-06, 2.923593421315422e-06, -2.206768214995963e-07, -1.203227957690371e-07, -3.381181730985542e-09, 4.022824706790907e-09, 7.041049107708875e-10, -3.562885681170365e-11, -3.281525407784209e-11, -5.031807464141896e-12, 2.034136401832885e-13, 2.361284283178230e-13, 4.602774507763180e-14, 1.809170850970874e-15, -1.333813332262995e-15, -4.045891156434286e-16, -5.197949512114809e-17, 2.222220713310119e-18, 2.625223897130583e-18, 6.226001466529447e-19, 6.571419077089260e-20, -6.672159423054950e-21, -4.476242224056620e-21, -9.790792477821165e-22, -9.222211318122281e-23, 1.427942034536028e-23, 7.759197090219954e-24, 1.546886518887300e-24, 9.585471984274525e-26, -4.005276449706623e-26, -1.459299762834743e-26, -1.870491354620814e-27, 2.271573838802745e-28, 1.360979028415734e-28, 1.102172718357361e-29, -6.919347646474293e-30, 4.875282352118995e-31}; + }; + } // namespace polyCoeff \ No newline at end of file diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 26aa293..fc2b2b4 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -139,7 +139,7 @@ std::unique_ptr PolySolver::buildIndividualForms(const mfem::Array M --- - mfem::Vector negOneVec(m_fePhi->GetVDim()); + mfem::Vector negOneVec(m_mesh->SpaceDimension()); negOneVec = -1.0; m_negationCoeff = std::make_unique(negOneVec); @@ -184,9 +184,6 @@ void PolySolver::solve() const { throw std::runtime_error("PolytropeOperator is not finalized. Cannot solve."); } - GetDofCoordinates(*m_feTheta, "dof2posTheta.csv"); - GetDofCoordinates(*m_fePhi, "dof2posPhi.csv"); - // It's safer to get the offsets directly from the operator after finalization const mfem::Array& block_offsets = m_polytropOperator->GetBlockOffsets(); // Assuming a getter exists or accessing member if public/friend mfem::BlockVector state_vector(block_offsets); @@ -194,23 +191,14 @@ void PolySolver::solve() const { state_vector.GetBlock(1) = phiVec; // NOLINT(*-slicing) - mfem::Vector zero_rhs(block_offsets.Last()); + mfem::Vector zero_rhs(block_offsets.Last()); // TODO: Is this right? zero_rhs = 0.0; const solverBundle sb = setupNewtonSolver(); - - // EMB 2025: Calling Mult gets the gradient of the operator for the NewtonSolver - // This then is set as the operator for the solver for NewtonSolver - // The solver (assuming it is an iterative solver) then sets the - // operator for its preconditioner to this. - // What this means is that there is no need to manually deal - // with updating the preconditioner at every newton step as the - // changes to the jacobian are automatically propagated through the - // solving chain. This is at least true with MFEM 4.8-rc0 - std::string windowTitle = "testWindow"; - std::string keyset = ""; sb.newton.Mult(zero_rhs, state_vector); + // Ax = b for x + // --- Save and view the solution --- saveAndViewSolution(state_vector); } @@ -219,28 +207,31 @@ SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const { mfem::Array theta_ess_tdof_list; mfem::Array phi_ess_tdof_list; - mfem::Array thetaCenterDofs, phiCenterDofs; - mfem::Array thetaCenterVals, phiCenterVals; - std::tie(thetaCenterDofs, phi_ess_tdof_list) = findCenterElement(); + mfem::Array thetaCenterDofs, phiCenterDofs; // phiCenterDofs are not used + mfem::Array thetaCenterVals; + std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement(); thetaCenterVals.SetSize(thetaCenterDofs.Size()); - phiCenterVals.SetSize(phi_ess_tdof_list.Size()); thetaCenterVals = 1.0; - phiCenterVals = 0.0; mfem::Array ess_brd(m_mesh->bdr_attributes.Max()); ess_brd = 1; - mfem::Array thetaSurfaceVals; + mfem::Array thetaSurfaceVals, phiSurfaceVals; m_feTheta->GetEssentialTrueDofs(ess_brd, theta_ess_tdof_list); + m_fePhi->GetEssentialTrueDofs(ess_brd, phi_ess_tdof_list); + thetaSurfaceVals.SetSize(theta_ess_tdof_list.Size()); thetaSurfaceVals = 0.0; + phiSurfaceVals.SetSize(phi_ess_tdof_list.Size()); + phiSurfaceVals = polycoeff::thetaSurfaceFlux(m_polytropicIndex); + // combine the essential dofs with the center dofs theta_ess_tdof_list.Append(thetaCenterDofs); thetaSurfaceVals.Append(thetaCenterVals); SSE::MFEMArrayPair thetaPair = std::make_pair(theta_ess_tdof_list, thetaSurfaceVals); - SSE::MFEMArrayPair phiPair = std::make_pair(phi_ess_tdof_list, phiCenterVals); + SSE::MFEMArrayPair phiPair = std::make_pair(phi_ess_tdof_list, phiSurfaceVals); SSE::MFEMArrayPairSet pairSet = std::make_pair(thetaPair, phiPair); return pairSet; @@ -275,16 +266,39 @@ void PolySolver::setInitialGuess() const { const double radius = Probe::getMeshRadius(*m_mesh); const double u = 1/radius; - return (-1.0/radius) * r + 1; + // return (-1.0/radius) * r + 1; // return -std::pow((u*r), 2)+1.0; // The series expansion is a better guess; however, this is cheaper and ensures that the value at the surface is very close to zero in a way that the series expansion does not + return laneEmden::thetaSeriesExpansion(r, m_polytropicIndex, 5); } ); + + mfem::VectorFunctionCoefficient phiSurfaceVectors (m_mesh->SpaceDimension(), + [this](const mfem::Vector &x, mfem::Vector &y) { + const double r = x.Norml2(); + mfem::Vector xh(x); + xh /= r; // Normalize the vector + y.SetSize(m_mesh->SpaceDimension()); + y = xh; + y *= polycoeff::thetaSurfaceFlux(m_polytropicIndex); + } + ); + // We want to apply specific boundary conditions to the surface + mfem::Array ess_brd(m_mesh->bdr_attributes.Max()); + ess_brd = 1; + + // θ = 0 at surface + mfem::ConstantCoefficient surfacePotential(0); + m_theta->ProjectCoefficient(thetaInitGuess); + m_theta->ProjectBdrCoefficient(surfacePotential, ess_brd); + mfem::GradientGridFunctionCoefficient phiInitGuess (m_theta.get()); m_phi->ProjectCoefficient(phiInitGuess); + m_phi->ProjectBdrCoefficientNormal(phiSurfaceVectors, ess_brd); + if (m_config.get("Poly:Solver:ViewInitialGuess", false)) { Probe::glVisView(*m_theta, *m_mesh, "θ init"); - Probe::glVisView(*m_phi, *m_mesh, "ɸ init"); + Probe::glVisView(*m_phi, *m_mesh, "φ init"); } } @@ -335,7 +349,7 @@ void PolySolver::LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); } -void PolySolver::GetDofCoordinates(mfem::FiniteElementSpace &fes, const std::string& filename) const { +void PolySolver::GetDofCoordinates(const mfem::FiniteElementSpace &fes, const std::string& filename) { mfem::Mesh *mesh = fes.GetMesh(); double r = Probe::getMeshRadius(*mesh); std::ofstream outputFile(filename, std::ios::out | std::ios::trunc); diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 952ffe8..b463973 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -111,6 +111,8 @@ private: // Private methods * trying to account for the true size of the system. VSize on the other hand refers to the total number of dofs. * * @return blockOffsets The offsets for the blocks in the operator + * + * @pre m_feTheta and m_fePhi must be valid, populated FiniteElementSpaces. */ mfem::Array computeBlockOffsets() const; diff --git a/src/poly/utils/private/integrators.cpp b/src/poly/utils/private/integrators.cpp index dab8c36..e433d62 100644 --- a/src/poly/utils/private/integrators.cpp +++ b/src/poly/utils/private/integrators.cpp @@ -97,6 +97,8 @@ namespace polyMFEMUtils { } // Calculate the Jacobian + + // TODO: Check for when theta ~ 0? const double u_safe = std::max(u_val, 0.0); const double d_u_nl = m_polytropicIndex * std::pow(u_safe, m_polytropicIndex - 1); const double x2_d_u_nl = d_u_nl; diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp index 5f16560..813524d 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/operator.cpp @@ -28,6 +28,8 @@ #include "linalg/tensor.hpp" +static int s_PolyOperatorMultCalledCount = 0; + void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfem::DenseMatrix *mat) { if (!mat) { throw std::runtime_error("The operator is not a SparseMatrix."); @@ -78,7 +80,7 @@ void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfe */ void writeOperatorToCSV(const mfem::Operator &op, const std::string &filename, - int precision = 6) // Add precision argument + const int precision = 6) // Add precision argument { // Attempt to cast the Operator to a SparseMatrix const auto *sparse_mat = dynamic_cast(&op); @@ -138,34 +140,19 @@ PolytropeOperator::PolytropeOperator( void PolytropeOperator::finalize(const mfem::Vector &initTheta) { if (m_isFinalized) { - return; + return; // do nothing if already finalized } - m_Mmat = std::make_unique(m_M->SpMat()); - m_Qmat = std::make_unique(m_Q->SpMat()); - m_Dmat = std::make_unique(m_D->SpMat()); - for (const int thetaDof : m_theta_ess_tdofs.first) { - m_Mmat->EliminateRow(thetaDof); - m_Qmat->EliminateCol(thetaDof); - } - // These are commented out because they theoretically are wrong (need to think more about how to apply essential dofs to a vector div field) - // for (const int phiDof : m_phi_ess_tdofs.first) { - // if (phiDof >=0 && phiDof < m_Dmat->Height()) { - // m_Dmat->EliminateRowCol(phiDof); - // m_Qmat->EliminateRow(phiDof); - // m_Mmat->EliminateCol(phiDof); - // } - // } - - m_negM_op = std::make_unique(m_Mmat.get(), -1.0); - m_negQ_op = std::make_unique(m_Qmat.get(), -1.0); + m_negM_op = std::make_unique(m_M.get(), -1.0); + m_negQ_op = std::make_unique(m_Q.get(), -1.0); // Set up the constant parts of the jacobian now 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_Dmat.get()); // D (constant) + 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) + // Build the initial preconditioner based on some initial guess const auto &grad = m_f->GetGradient(initTheta); updatePreconditioner(grad); @@ -173,15 +160,22 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) { } const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const { - if (m_jacobian == nullptr) { - MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call."); - } - return *m_jacobian; + if (m_jacobian == nullptr) { + MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call."); + } + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator not finalized prior to call to GetJacobianOperator()."); + } + + return *m_jacobian; } mfem::BlockDiagonalPreconditioner& PolytropeOperator::GetPreconditioner() const { + if (m_schurPreconditioner == nullptr) { + MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call."); + } if (!m_isFinalized) { - MFEM_ABORT("PolytropeOperator::GetPreconditioner called before finalize"); + MFEM_ABORT("PolytropeOperator not finalized prior to call to GetPreconditioner()."); } return *m_schurPreconditioner; } @@ -213,38 +207,39 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { // R1 = Dɸ - Qθ MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult"); - MFEM_ASSERT(m_Mmat.get() != nullptr, "SparseMatrix m_Mmat is null in PolytropeOperator::Mult"); - MFEM_ASSERT(m_Dmat.get() != nullptr, "SparseMatrix m_Dmat is null in PolytropeOperator::Mult"); - MFEM_ASSERT(m_Qmat.get() != nullptr, "SparseMatrix m_Qmat is null in PolytropeOperator::Mult"); - m_f->Mult(x_theta, f_term); - m_Mmat->Mult(x_phi, Mphi_term); - m_Dmat->Mult(x_phi, Dphi_term); - m_Qmat->Mult(x_theta, Qtheta_term); + 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_D->Mult(x_phi, Dphi_term); + m_Q->Mult(x_theta, Qtheta_term); subtract(f_term, Mphi_term, y_R0); subtract(Dphi_term, Qtheta_term, y_R1); - // -- Apply essential boundary conditions -- for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) { - int idx = m_theta_ess_tdofs.first[i]; - if (idx >= 0 && idx < y_R0.Size()) { + 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); // Zero out the essential theta dofs in the bi-linear form - // y_block.GetBlock(0)[idx] = 0; // Zero out the essential theta dofs in the bi-linear form + // 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 look into how the true dof -> vector component works - // for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) { - // int idx = m_phi_ess_tdofs.first[i]; - // if (idx >= 0 && idx < y_R1.Size()) { - // // const double &targetValue = m_phi_ess_tdofs.second[i]; - // // y_block.GetBlock(1)[idx] = targetValue - x_phi(idx); // Zero out the essential phi dofs in the bi-linear form - // y_block.GetBlock(1)[idx] = 0; // Zero out the essential phi dofs in the bi-linear form - // } - // } + // TODO: Are the residuals for φ being calculated correctly? + + std::ofstream outfileθ("PolyOperatorMultCallTheta_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc); + outfileθ << "dof,R_θ\n"; + for (int i = 0; i < y_R0.Size(); i++) { + outfileθ << i << "," << y_R0(i) << "\n"; + } + outfileθ.close(); + std::ofstream outfileφ("PolyOperatorMultCallPhi_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc); + outfileφ << "dof,R_ɸ\n"; + for (int i = 0; i < y_R1.Size(); i++) { + outfileφ << i << "," << y_R1(i) << "\n"; + } + outfileφ.close(); + ++s_PolyOperatorMultCalledCount; } @@ -253,7 +248,7 @@ void PolytropeOperator::updateInverseNonlinearJacobian(const mfem::Operator &gra if (const auto *sparse_mat = dynamic_cast(&grad); sparse_mat != nullptr) { approxJacobiInvert(*sparse_mat, m_invNonlinearJacobian, "Nonlinear Jacobian"); } else { - MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear jacobian where nonlinear jacobian is not trivially castable to a sparse matrix"); + MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear jacobian where nonlinear jacobian is not dynamically castable to a sparse matrix"); } } @@ -262,11 +257,11 @@ void PolytropeOperator::updateInverseSchurCompliment() const { if (m_invNonlinearJacobian == nullptr) { MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian"); } - mfem::SparseMatrix* schurCompliment(m_Dmat.get()); // This is now a copy of D + mfem::SparseMatrix* schurCompliment(&m_D->SpMat()); // This is now a copy of D // Calculate S = D - Q df^{-1} M - mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(*m_Qmat, *m_invNonlinearJacobian); // Q * df^{-1} - mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, *m_Mmat); // 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 // 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 @@ -275,6 +270,9 @@ void PolytropeOperator::updateInverseSchurCompliment() const { if (m_schurPreconditioner == nullptr) { m_schurPreconditioner = std::make_unique(m_blockOffsets); } + + // ⎡ḟ(θ)^-1 0⎤ + // ⎣0 S^-1⎦ m_schurPreconditioner->SetDiagonalBlock(0, m_invNonlinearJacobian.get()); m_schurPreconditioner->SetDiagonalBlock(1, m_invSchurCompliment.get()); } @@ -312,6 +310,26 @@ void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess } else { MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs"); } + + if (m_M) { + m_M->EliminateTestEssentialBC(theta_ess_tdofs.first); + m_M->EliminateTrialEssentialBC(phi_ess_tdofs.first); + } else { + MFEM_ABORT("m_M is null in PolytropeOperator::SetEssentialTrueDofs"); + } + + if (m_Q) { + m_Q->EliminateTrialEssentialBC(theta_ess_tdofs.first); + m_Q->EliminateTestEssentialBC(phi_ess_tdofs.first); + } else { + MFEM_ABORT("m_Q is null in PolytropeOperator::SetEssentialTrueDofs"); + } + + if (m_D) { + m_D->EliminateEssentialBC(phi_ess_tdofs.first); + } else { + MFEM_ABORT("m_D is null in PolytropeOperator::SetEssentialTrueDofs"); + } } void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) { diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h index d06b978..1c84db0 100644 --- a/src/poly/utils/public/operator.h +++ b/src/poly/utils/public/operator.h @@ -70,10 +70,6 @@ private: SSE::MFEMArrayPair m_theta_ess_tdofs; SSE::MFEMArrayPair m_phi_ess_tdofs; - std::unique_ptr m_Mmat; - std::unique_ptr m_Qmat; - std::unique_ptr m_Dmat; - std::unique_ptr m_negM_op; std::unique_ptr m_negQ_op; mutable std::unique_ptr m_jacobian;