diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index fc2b2b4..94e0881 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -156,7 +156,7 @@ std::unique_ptr PolySolver::buildIndividualForms(const mfem::Arrayf->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex)); - return std::move(forms); + return forms; } void PolySolver::assembleAndFinalizeForm(auto &f) { @@ -264,11 +264,11 @@ void PolySolver::setInitialGuess() const { [this](const mfem::Vector &x) { const double r = x.Norml2(); const double radius = Probe::getMeshRadius(*m_mesh); - const double u = 1/radius; + // const double u = 1/radius; // 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); + return laneEmden::thetaSeriesExpansion(r, m_polytropicIndex, 10); } ); @@ -399,7 +399,8 @@ solverBundle PolySolver::setupNewtonSolver() const { solver.solver.SetMaxIter(gmresMaxIter); solver.solver.SetPrintLevel(gmresPrintLevel); - solver.solver.SetPreconditioner(m_polytropOperator->GetPreconditioner()); + // Preconditioner turned off because the polytrope operator seems *very* well conditioned without it + // solver.solver.SetPreconditioner(m_polytropOperator->GetPreconditioner()); // --- Set up the Newton solver --- solver.newton.SetRelTol(newtonRelTol); solver.newton.SetAbsTol(newtonAbsTol); diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build index 9b73678..8fc7d68 100644 --- a/src/poly/utils/meson.build +++ b/src/poly/utils/meson.build @@ -30,6 +30,7 @@ dependencies = [ quill_dep, config_dep, types_dep, + mfemanalysis_dep, ] libpolyutils = static_library('polyutils', diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp index ea363d5..6dda6c2 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/operator.cpp @@ -21,70 +21,13 @@ #include "operator.h" #include "4DSTARTypes.h" #include "mfem.hpp" +#include "mfem_smout.h" #include +#include "../../../../utils/debugUtils/MFEMAnalysisUtils/MFEMAnalysis-cpp/src/include/mfem_smout.h" + 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."); - } - - std::ofstream outfile(filename); - if (!outfile.is_open()) { - throw std::runtime_error("Failed to open file: " + filename); - } - - - int height = mat->Height(); - int width = mat->Width(); - - // Set precision for floating-point output - outfile << std::fixed << std::setprecision(precision); - - for (int i = 0; i < width; i++) { - outfile << i; - if (i < width - 1) { - outfile << ","; - } - else { - outfile << "\n"; - } - } - - // Iterate through rows - for (int i = 0; i < height; ++i) { - for (int j = 0; j < width; ++j) { - outfile << mat->Elem(i, j); - if (j < width - 1) { - outfile << ","; - } - } - outfile << std::endl; - } - - outfile.close(); -} - -/** - * @brief Writes the dense representation of an MFEM Operator (if it's a SparseMatrix) to a CSV file. - * - * @param op The MFEM Operator to write. - * @param filename The name of the output CSV file. - * @param precision Number of decimal places for floating-point values. - */ - void writeOperatorToCSV(const mfem::Operator &op, - const std::string &filename, - const int precision = 6) // Add precision argument -{ - // Attempt to cast the Operator to a SparseMatrix - const auto *sparse_mat = dynamic_cast(&op); - if (!sparse_mat) { - throw std::runtime_error("The operator is not a SparseMatrix."); - } - const mfem::DenseMatrix *mat = sparse_mat->ToDenseMatrix(); - writeDenseMatrixToCSV(filename, precision, mat); -} +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 @@ -138,14 +81,32 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) { return; // do nothing if already finalized } - m_negM_op = std::make_unique(m_M.get(), -1.0); - m_negQ_op = std::make_unique(m_Q.get(), -1.0); + 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 auto& dof: m_theta_ess_tdofs.first) { + m_Mmat->EliminateRow(dof); + m_Qmat->EliminateCol(dof); + } + + 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_Mmat.get(), -1.0); + // m_negQ_op = std::make_unique(m_Qmat.get(), -1.0); + + m_Mmat->operator*=(-1.0); + m_Qmat->operator*=(-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_D.get()); //<- D (constant) + m_jacobian->SetBlock(0, 1, m_Mmat.get()); //<- -M (constant) + m_jacobian->SetBlock(1, 0, m_Qmat.get()); //<- -Q (constant) + m_jacobian->SetBlock(1, 1, m_Dmat.get()); //<- D (constant) // Build the initial preconditioner based on some initial guess const auto &grad = m_f->GetGradient(initTheta); @@ -214,28 +175,20 @@ 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]; + // 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 + 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) + } + } + 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) } } - // 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; - + std::cout << "φ Norm: " << y_block.GetBlock(1).Norml2() << std::endl; } @@ -270,6 +223,7 @@ void PolytropeOperator::updateInverseSchurCompliment() const { // ⎣0 S^-1⎦ m_schurPreconditioner->SetDiagonalBlock(0, m_invNonlinearJacobian.get()); m_schurPreconditioner->SetDiagonalBlock(1, m_invSchurCompliment.get()); + } void PolytropeOperator::updatePreconditioner(const mfem::Operator &grad) const { @@ -287,11 +241,10 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const { const mfem::Vector& x_theta = x_block.GetBlock(0); auto &grad = m_f->GetGradient(x_theta); - updatePreconditioner(grad); + // auto *gradPtr = &grad; + // updatePreconditioner(grad); m_jacobian->SetBlock(0, 0, &grad); - // The other blocks are set up in finalize since they are constant. Only J00 depends on the current state of theta - return *m_jacobian; } @@ -305,26 +258,6 @@ 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 1c84db0..e046480 100644 --- a/src/poly/utils/public/operator.h +++ b/src/poly/utils/public/operator.h @@ -70,6 +70,7 @@ private: 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; mutable std::unique_ptr m_jacobian;