From 2a91d57ad712d51a1919b3c73f9e8ee2b605aeaf Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 9 Jun 2025 10:19:18 -0400 Subject: [PATCH] fix(poly): fixed numerous bugs related to inconsistent system sizing with the reduced operator this has restored the symmetry which we relied on before. --- src/poly/solver/private/polySolver.cpp | 122 ++++- src/poly/utils/private/polytropeOperator.cpp | 546 ++++++++++--------- src/poly/utils/private/utilities.cpp | 9 +- src/poly/utils/public/polytropeOperator.h | 86 +-- src/poly/utils/public/utilities.h | 34 +- 5 files changed, 481 insertions(+), 316 deletions(-) diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 29c221c..50c6cc9 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -36,6 +36,7 @@ #include "probe.h" #include "resourceManager.h" #include "resourceManagerTypes.h" +#include "utilities.h" #include "quill/LogMacros.h" @@ -113,6 +114,19 @@ void PolySolver::assembleBlockSystem() { const std::unique_ptr forms = buildIndividualForms(blockOffsets); + // const double penalty_param = m_config.get("Poly::Solver::ZeroDerivativePenalty", 1e10); + // mfem::Array thetaCenterDofs, phiCenterDofs; + // std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement(); + // mfem::SparseMatrix& D_mat = forms->D->SpMat(); + // + // for (int i = 0; i < phiCenterDofs.Size(); ++i) + // { + // const int dof_idx = phiCenterDofs[i]; + // if (dof_idx >= 0 && dof_idx < D_mat.Height()) { + // D_mat(dof_idx, dof_idx) += penalty_param; + // } + // } + // --- Build the BlockOperator --- m_polytropOperator = std::make_unique( std::move(forms->M), @@ -173,8 +187,9 @@ void PolySolver::assembleAndFinalizeForm(auto &f) { void PolySolver::solve() const { // --- Set the initial guess for the solution --- setInitialGuess(); - setOperatorEssentialTrueDofs(); + + // --- Cast the GridFunctions to mfem::Vector --- const auto thetaVec = static_cast(*m_theta); // NOLINT(*-slicing) const auto phiVec = static_cast(*m_phi); // NOLINT(*-slicing) @@ -182,23 +197,31 @@ void PolySolver::solve() const { // Finalize with the initial state of theta for the initial jacobian calculation m_polytropOperator->finalize(thetaVec); - // It's safer to get the offsets directly from the operator after finalization - const mfem::Array& block_offsets = m_polytropOperator->get_reduced_block_offsets(); // Assuming a getter exists or accessing member if public/friend - mfem::BlockVector state_vector(block_offsets); - state_vector.GetBlock(0) = thetaVec; // NOLINT(*-slicing) - state_vector.GetBlock(1) = phiVec; // NOLINT(*-slicing) + // --- Broadcast initial condition to the full state vector --- + const mfem::Array& full_block_offsets = m_polytropOperator->get_block_offsets(); + mfem::Vector x_full(full_block_offsets.Last()); + mfem::BlockVector x_full_block(x_full, full_block_offsets); + x_full_block.GetBlock(0) = thetaVec; // NOLINT(*-slicing) + x_full_block.GetBlock(1) = phiVec; // NOLINT(*-slicing) + // --- Extract only the free DOFs from the full state vector --- + const mfem::Array& freeDofs = m_polytropOperator->get_free_dofs(); + mfem::Vector x_free(m_polytropOperator->get_reduced_system_size()); + x_full.GetSubVector(freeDofs, x_free); // Extract the free DOFs from the full vector + + // --- Initialize RHS --- mfem::Vector zero_rhs(m_polytropOperator->get_reduced_system_size()); zero_rhs = 0.0; + // --- Setup and run the Newton solver --- const solverBundle sb = setupNewtonSolver(); - sb.newton.Mult(zero_rhs, state_vector); + sb.newton.Mult(zero_rhs, x_free); - // TODO: Need some system here to reconstruct the full system vector from what MFEM's Newton Solver - // will return, since that will end up being the vector for the reduced system. + // --- Reconstruct the full state vector from the reduced solution --- + mfem::BlockVector solution = m_polytropOperator->reconstruct_full_block_state_vector(x_free); // --- Save and view an approximate 1D solution --- - saveAndViewSolution(state_vector); + saveAndViewSolution(solution); } SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const { @@ -206,10 +229,12 @@ SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const { mfem::Array phi_ess_tdof_list; mfem::Array thetaCenterDofs, phiCenterDofs; // phiCenterDofs are not used - mfem::Array thetaCenterVals; + mfem::Array thetaCenterVals, phiCenterVals; std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement(); thetaCenterVals.SetSize(thetaCenterDofs.Size()); + // phiCenterVals.SetSize(phiCenterDofs.Size()); + // phiCenterVals = 0.0; thetaCenterVals = 1.0; mfem::Array ess_brd(m_mesh.bdr_attributes.Max()); @@ -228,6 +253,9 @@ SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const { theta_ess_tdof_list.Append(thetaCenterDofs); thetaSurfaceVals.Append(thetaCenterVals); + // phi_ess_tdof_list.Append(phiCenterDofs); + // phiSurfaceVals.Append(phiCenterVals); + SSE::MFEMArrayPair thetaPair = std::make_pair(theta_ess_tdof_list, thetaSurfaceVals); SSE::MFEMArrayPair phiPair = std::make_pair(phi_ess_tdof_list, phiSurfaceVals); SSE::MFEMArrayPairSet pairSet = std::make_pair(thetaPair, phiPair); @@ -238,21 +266,59 @@ SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const { std::pair, mfem::Array> PolySolver::findCenterElement() const { mfem::Array thetaCenterDofs; mfem::Array phiCenterDofs; - mfem::DenseMatrix centerPoint(m_mesh.SpaceDimension(), 1); - centerPoint(0, 0) = 0.0; - centerPoint(1, 0) = 0.0; - centerPoint(2, 0) = 0.0; - mfem::Array elementIDs; - mfem::Array ips; - m_mesh.FindPoints(centerPoint, elementIDs, ips); - mfem::Array tempDofs; - for (int i = 0; i < elementIDs.Size(); i++) { - m_feTheta->GetElementDofs(elementIDs[i], tempDofs); - thetaCenterDofs.Append(tempDofs); - m_fePhi->GetElementDofs(elementIDs[i], tempDofs); - phiCenterDofs.Append(tempDofs); + // --- 1. Find the index of the single mesh vertex at the origin --- + int center_vertex_idx = -1; + constexpr double tol = 1e-9; // A small tolerance for floating point comparison + + for (int i = 0; i < m_mesh.GetNV(); ++i) { + const double* vertex_coords = m_mesh.GetVertex(i); + if (std::abs(vertex_coords[0]) < tol && + std::abs(vertex_coords[1]) < tol && + std::abs(vertex_coords[2]) < tol) { + + center_vertex_idx = i; + break; // Found it, assume there's only one. + } } + + if (center_vertex_idx == -1) { + MFEM_ABORT("Could not find the center vertex at [0,0,0]. Check mesh construction."); + } + + // --- 2. Get Theta (H1) DoFs associated ONLY with that vertex --- + // CORRECTED: Use GetVertexDofs, not GetVDofs. + m_feTheta->GetVertexDofs(center_vertex_idx, thetaCenterDofs); + + + mfem::Array central_element_ids; + + // PERF: could probably move this to a member variable and populate during construction + mfem::Table* vertex_to_elements_table = m_mesh.GetVertexToElementTable(); + vertex_to_elements_table->Finalize(); + mfem::Array element_ids; + vertex_to_elements_table->GetRow(center_vertex_idx, element_ids); + delete vertex_to_elements_table; + + for (int i = 0; i < element_ids.Size(); ++i) { + int element_id = element_ids[i]; + mfem::Array tempDofs; + m_fePhi->GetElementDofs(element_id, tempDofs); + + // decode negative dofs to their true, physical, dof indices + for (int j = 0; j < tempDofs.Size(); ++j) { + int dof = tempDofs[j]; + if (dof < 0) { + dof = -dof - 1; // Convert to positive index + } + phiCenterDofs.Append(dof); + } + } + + + phiCenterDofs.Sort(); + phiCenterDofs.Unique(); + return std::make_pair(thetaCenterDofs, phiCenterDofs); } @@ -302,10 +368,18 @@ void PolySolver::setInitialGuess() const { // solver towards a more consistent answer (x_φ - target) m_phi->ProjectBdrCoefficientNormal(phiSurfaceVectors, ess_brd); + auto [thetaCenterDofs, phiCenterDofs] = findCenterElement(); + + for (int i = 0; i < phiCenterDofs.Size(); ++i) + { + (*m_phi)(phiCenterDofs[i]) = 0.0; + } + if (m_config.get("Poly:Solver:ViewInitialGuess", false)) { Probe::glVisView(*m_theta, m_mesh, "θ init"); Probe::glVisView(*m_phi, m_mesh, "φ init"); } + std::cout << "HERE" << std::endl; } diff --git a/src/poly/utils/private/polytropeOperator.cpp b/src/poly/utils/private/polytropeOperator.cpp index 1be8402..898c03c 100644 --- a/src/poly/utils/private/polytropeOperator.cpp +++ b/src/poly/utils/private/polytropeOperator.cpp @@ -26,7 +26,109 @@ #include "mfem_smout.h" #include +// --- SchurCompliment Class Implementation --- +// SchurCompliment: Constructors +SchurCompliment::SchurCompliment( + const mfem::SparseMatrix &QOp, + const mfem::SparseMatrix &DOp, + const mfem::SparseMatrix &MOp, + const mfem::Solver &GradInvOp) : +mfem::Operator(DOp.Height(), DOp.Width()) +{ + SetOperator(QOp, DOp, MOp, GradInvOp); + m_nPhi = m_DOp->Height(); + m_nTheta = m_MOp->Height(); +} + +SchurCompliment::SchurCompliment( + const mfem::SparseMatrix &QOp, + const mfem::SparseMatrix &DOp, + const mfem::SparseMatrix &MOp) : +mfem::Operator(DOp.Height(), DOp.Width()) +{ + updateConstantTerms(QOp, DOp, MOp); + m_nPhi = m_DOp->Height(); + m_nTheta = m_MOp->Height(); +} + +// SchurCompliment: Public Interface Methods +void SchurCompliment::Mult(const mfem::Vector &x, mfem::Vector &y) const { + // Check that the input vector is the correct size + if (x.Size() != m_nPhi) { + MFEM_ABORT("Input vector x has size " + std::to_string(x.Size()) + ", expected " + std::to_string(m_nPhi)); + } + if (y.Size() != m_nPhi) { + MFEM_ABORT("Output vector y has size " + std::to_string(y.Size()) + ", expected " + std::to_string(m_nPhi)); + } + + // Check that the operators are set + if (m_QOp == nullptr) { + MFEM_ABORT("QOp is null in SchurCompliment::Mult"); + } + if (m_DOp == nullptr) { + MFEM_ABORT("DOp is null in SchurCompliment::Mult"); + } + if (m_MOp == nullptr) { + MFEM_ABORT("MOp is null in SchurCompliment::Mult"); + } + if (m_GradInvOp == nullptr) { + MFEM_ABORT("GradInvOp is null in SchurCompliment::Mult"); + } + + mfem::Vector v1(m_nTheta); // M * x + m_MOp -> Mult(x, v1); // M * x + + mfem::Vector v2(m_nTheta); // GradInv * M * x + m_GradInvOp -> Mult(v1, v2); // GradInv * M * x + + mfem::Vector v3(m_nPhi); // Q * GradInv * M * x + m_QOp -> Mult(v2, v3); // Q * GradInv * M * x + + mfem::Vector v4(m_nPhi); // D * x + m_DOp -> Mult(x, v4); // D * x + + subtract(v4, v3, y); // (D - Q * GradInv * M) * x +} + +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); +} + +void SchurCompliment::updateInverseNonlinearJacobian(const mfem::Solver &gradInv) { + m_GradInvOp = &gradInv; +} + +// SchurCompliment: Private Helper Methods +void SchurCompliment::updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp) { + m_QOp = &QOp; + m_DOp = &DOp; + m_MOp = &MOp; +} + + +// --- GMRESInverter Class Implementation --- + +// GMRESInverter: Constructor +GMRESInverter::GMRESInverter(const SchurCompliment &op) : +mfem::Operator(op.Height(), op.Width()), +m_op(op) { + m_solver.SetOperator(m_op); + m_solver.SetMaxIter(100); + m_solver.SetRelTol(1e-1); + m_solver.SetAbsTol(1e-1); +} + +// GMRESInverter: Public Interface Methods +void GMRESInverter::Mult(const mfem::Vector &x, mfem::Vector &y) const { + m_solver.Mult(x, y); // Approximates m_op^-1 * x +} + + +// --- PolytropeOperator Class Implementation --- + +// PolytropeOperator: Constructor PolytropeOperator::PolytropeOperator( std::unique_ptr M, @@ -51,128 +153,7 @@ PolytropeOperator::PolytropeOperator( m_invNonlinearJacobian = std::make_unique(0, 3); } -void PolytropeOperator::populate_free_dof_array() { - m_freeDofs.SetSize(0); - for (int i = 0; i < m_blockOffsets.Last(); i++) { - const int thetaSearchIndex = i; - const int phiSearchIndex = i - m_blockOffsets[1]; - if (phiSearchIndex < 0){ - if (m_theta_ess_tdofs.first.Find(thetaSearchIndex) == -1) { - m_freeDofs.Append(i); - } - } else { - if (m_phi_ess_tdofs.first.Find(phiSearchIndex) == -1) { - m_freeDofs.Append(i); - } - } - } -} - -void PolytropeOperator::scatter_boundary_conditions() { - mfem::Vector thetaStateValues(m_theta_ess_tdofs.first.Size()); - for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) { - thetaStateValues[i] = m_theta_ess_tdofs.second[i]; - } - mfem::Vector phiStateValues(m_phi_ess_tdofs.first.Size()); - for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) { - phiStateValues[i] = m_phi_ess_tdofs.second[i]; - } - - mfem::Array phiDofIndices(m_phi_ess_tdofs.first.Size()); - for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) { - phiDofIndices[i] = m_phi_ess_tdofs.first[i] + m_blockOffsets[1]; - } - - m_state.SetSize(m_blockOffsets.Last()); - m_state = 0.0; - m_state.SetSubVector(m_theta_ess_tdofs.first, thetaStateValues); - m_state.SetSubVector(phiDofIndices, phiStateValues); - -} - -void PolytropeOperator::construct_matrix_representations() { - 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_MReduced = std::make_unique(serif::utilities::buildReducedMatrix(*m_Mmat, m_phi_ess_tdofs.first, m_theta_ess_tdofs.first)); - m_QReduced = std::make_unique(serif::utilities::buildReducedMatrix(*m_Qmat, m_theta_ess_tdofs.first, m_phi_ess_tdofs.first)); - m_DReduced = std::make_unique(serif::utilities::buildReducedMatrix(*m_Dmat, m_phi_ess_tdofs.first, m_phi_ess_tdofs.first)); - - m_negQ_mat = std::make_unique(m_QReduced.get(), -1.0); -} - -void PolytropeOperator::construct_reduced_block_offsets() { - m_reducedBlockOffsets.SetSize(3); - m_reducedBlockOffsets[0] = 0; // R0 block (theta) - m_reducedBlockOffsets[1] = m_MReduced->Height(); // R1 block (theta) - m_reducedBlockOffsets[2] = m_QReduced->Height() + m_reducedBlockOffsets[1]; // R2 block (phi) -} - -void PolytropeOperator::construct_jacobian_constant_terms() { - m_jacobian = std::make_unique(m_reducedBlockOffsets); - m_jacobian->SetBlock(0, 1, m_MReduced.get()); //<- M (constant) - m_jacobian->SetBlock(1, 0, m_negQ_mat.get()); //<- -Q (constant) - m_jacobian->SetBlock(1, 1, m_DReduced.get()); //<- D (constant) -} - -void PolytropeOperator::finalize(const mfem::Vector &initTheta) { - using serif::utilities::buildReducedMatrix; - - if (m_isFinalized) { - return; - } - - // These functions must be called in this order since they depend on each others post state - // TODO: Refactor this so that either there are explicit checks to make sure the order is correct or make - // them pure functions - construct_matrix_representations(); - construct_reduced_block_offsets(); - construct_jacobian_constant_terms(); - scatter_boundary_conditions(); - populate_free_dof_array(); - - // Override the size based on the reduced system - height = m_reducedBlockOffsets.Last(); - width = m_reducedBlockOffsets.Last(); - - m_isFinalized = true; - -} - -const mfem::BlockOperator &PolytropeOperator::get_jacobian_operator() const { - 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::get_preconditioner() const { - if (m_schurPreconditioner == nullptr) { - MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call."); - } - if (!m_isFinalized) { - MFEM_ABORT("PolytropeOperator not finalized prior to call to GetPreconditioner()."); - } - return *m_schurPreconditioner; -} - -int PolytropeOperator::get_reduced_system_size() const { - if (!m_isFinalized) { - MFEM_ABORT("PolytropeOperator not finalized prior to call to GetReducedSystemSize()."); - } - return m_reducedBlockOffsets.Last(); -} - -const mfem::Vector &PolytropeOperator::reconstruct_full_state_vector(const mfem::Vector &reducedState) const { - m_state.SetSubVector(m_freeDofs, reducedState); // Scatter the reduced state vector into the full state vector - return m_state; -} - +// PolytropeOperator: Core Operator Overrides void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const { if (!m_isFinalized) { MFEM_ABORT("PolytropeOperator::Mult called before finalize"); @@ -225,7 +206,194 @@ void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) con } } +mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const { + //TODO: This now needs to be updated to deal with the reduced system size + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator::GetGradient called before finalize"); + } + m_state.SetSubVector(m_freeDofs, xFree); // Scatter the free dofs from the input vector xFree into the state vector + // --- Get the gradient of f --- + mfem::BlockVector x_block(const_cast(m_state), m_blockOffsets); + const mfem::Vector& x_theta = x_block.GetBlock(0); + // PERF: There are a lot of copies and loops here, probably performance could be gained by flattering some of these. + auto &grad = m_f->GetGradient(x_theta); + // updatePreconditioner(grad); + const auto gradMatrix = dynamic_cast(&grad); + + if (gradMatrix == nullptr) { + MFEM_ABORT("PolytropeOperator::GetGradient: Gradient is not a SparseMatrix."); + } + + m_gradReduced = std::make_unique ( + serif::utilities::build_reduced_matrix( + *gradMatrix, + m_theta_ess_tdofs.first, + m_theta_ess_tdofs.first + ) + ); + + m_jacobian->SetBlock(0, 0, m_gradReduced.get()); + + return *m_jacobian; +} + + +// PolytropeOperator: Setup and Finalization +void PolytropeOperator::finalize(const mfem::Vector &initTheta) { + using serif::utilities::build_reduced_matrix; + + if (m_isFinalized) { + return; + } + + // These functions must be called in this order since they depend on each others post state + // TODO: Refactor this so that either there are explicit checks to make sure the order is correct or make + // them pure functions + construct_matrix_representations(); + construct_reduced_block_offsets(); + construct_jacobian_constant_terms(); + scatter_boundary_conditions(); + populate_free_dof_array(); + + // Override the size based on the reduced system + height = m_reducedBlockOffsets.Last(); + width = m_reducedBlockOffsets.Last(); + + m_isFinalized = true; + +} + +// PolytropeOperator: Essential True DOF Management +void PolytropeOperator::set_essential_true_dofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) { + m_isFinalized = false; + m_theta_ess_tdofs = theta_ess_tdofs; + m_phi_ess_tdofs = phi_ess_tdofs; + + if (m_f) { + m_f->SetEssentialTrueDofs(theta_ess_tdofs.first); + } else { + MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs"); + } +} + +void PolytropeOperator::set_essential_true_dofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) { + set_essential_true_dofs(ess_tdof_pair_set.first, ess_tdof_pair_set.second); +} + +SSE::MFEMArrayPairSet PolytropeOperator::get_essential_true_dofs() const { + return std::make_pair(m_theta_ess_tdofs, m_phi_ess_tdofs); +} + +// PolytropeOperator: Getter Methods +const mfem::BlockOperator &PolytropeOperator::get_jacobian_operator() const { + 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::get_preconditioner() const { + if (m_schurPreconditioner == nullptr) { + MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call."); + } + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator not finalized prior to call to GetPreconditioner()."); + } + return *m_schurPreconditioner; +} + +int PolytropeOperator::get_reduced_system_size() const { + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator not finalized prior to call to GetReducedSystemSize()."); + } + return m_reducedBlockOffsets.Last(); +} + +// PolytropeOperator: State Reconstruction +const mfem::Vector &PolytropeOperator::reconstruct_full_state_vector(const mfem::Vector &reducedState) const { + m_state.SetSubVector(m_freeDofs, reducedState); // Scatter the reduced state vector into the full state vector + return m_state; +} + +const mfem::BlockVector PolytropeOperator::reconstruct_full_block_state_vector(const mfem::Vector &reducedState) const { + m_state.SetSubVector(m_freeDofs, reducedState); // Scatter the reduced state vector into the full state vector + mfem::BlockVector x_block(m_state, m_blockOffsets); + return x_block; +} + +// PolytropeOperator: DOF Population +void PolytropeOperator::populate_free_dof_array() { + m_freeDofs.SetSize(0); + for (int i = 0; i < m_blockOffsets.Last(); i++) { + const int thetaSearchIndex = i; + const int phiSearchIndex = i - m_blockOffsets[1]; + if (phiSearchIndex < 0){ + if (m_theta_ess_tdofs.first.Find(thetaSearchIndex) == -1) { + m_freeDofs.Append(i); + } + } else { + if (m_phi_ess_tdofs.first.Find(phiSearchIndex) == -1) { + m_freeDofs.Append(i); + } + } + } +} + +// PolytropeOperator: Private Helper Methods - Construction and Setup +void PolytropeOperator::construct_matrix_representations() { + 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_MReduced = std::make_unique(serif::utilities::build_reduced_matrix(*m_Mmat, m_phi_ess_tdofs.first, m_theta_ess_tdofs.first)); + m_QReduced = std::make_unique(serif::utilities::build_reduced_matrix(*m_Qmat, m_theta_ess_tdofs.first, m_phi_ess_tdofs.first)); + m_DReduced = std::make_unique(serif::utilities::build_reduced_matrix(*m_Dmat, m_phi_ess_tdofs.first, m_phi_ess_tdofs.first)); + + m_negQ_mat = std::make_unique(m_QReduced.get(), -1.0); +} + +void PolytropeOperator::construct_reduced_block_offsets() { + m_reducedBlockOffsets.SetSize(3); + m_reducedBlockOffsets[0] = 0; // R0 block (theta) + m_reducedBlockOffsets[1] = m_MReduced->Height(); // R1 block (theta) + m_reducedBlockOffsets[2] = m_QReduced->Height() + m_reducedBlockOffsets[1]; // R2 block (phi) +} + +void PolytropeOperator::construct_jacobian_constant_terms() { + m_jacobian = std::make_unique(m_reducedBlockOffsets); + m_jacobian->SetBlock(0, 1, m_MReduced.get()); //<- M (constant) + m_jacobian->SetBlock(1, 0, m_negQ_mat.get()); //<- -Q (constant) + m_jacobian->SetBlock(1, 1, m_DReduced.get()); //<- D (constant) +} + +void PolytropeOperator::scatter_boundary_conditions() { + mfem::Vector thetaStateValues(m_theta_ess_tdofs.first.Size()); + for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) { + thetaStateValues[i] = m_theta_ess_tdofs.second[i]; + } + mfem::Vector phiStateValues(m_phi_ess_tdofs.first.Size()); + for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) { + phiStateValues[i] = m_phi_ess_tdofs.second[i]; // TODO: figure out if this needs to be normalized + } + + mfem::Array phiDofIndices(m_phi_ess_tdofs.first.Size()); + for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) { + phiDofIndices[i] = m_phi_ess_tdofs.first[i] + m_blockOffsets[1]; + } + + m_state.SetSize(m_blockOffsets.Last()); + m_state = 0.0; + m_state.SetSubVector(m_theta_ess_tdofs.first, thetaStateValues); + m_state.SetSubVector(phiDofIndices, phiStateValues); + +} + +// PolytropeOperator: Private Helper Methods - Jacobian and Preconditioner Updates void PolytropeOperator::update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const { m_invNonlinearJacobian->SetOperator(grad); } @@ -257,137 +425,3 @@ void PolytropeOperator::update_preconditioner(const mfem::Operator &grad) const update_inverse_nonlinear_jacobian(grad); update_inverse_schur_compliment(); } - -mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const { - //TODO: This now needs to be updated to deal with the reduced system size - if (!m_isFinalized) { - MFEM_ABORT("PolytropeOperator::GetGradient called before finalize"); - } - m_state.SetSubVector(m_freeDofs, xFree); // Scatter the free dofs from the input vector xFree into the state vector - // --- Get the gradient of f --- - mfem::BlockVector x_block(const_cast(m_state), m_blockOffsets); - const mfem::Vector& x_theta = x_block.GetBlock(0); - - // PERF: There are a lot of copies and loops here, probably performance could be gained by flattering some of these. - auto &grad = m_f->GetGradient(x_theta); - // updatePreconditioner(grad); - const auto gradMatrix = dynamic_cast(&grad); - - if (gradMatrix == nullptr) { - MFEM_ABORT("PolytropeOperator::GetGradient: Gradient is not a SparseMatrix."); - } - - mfem::SparseMatrix reducedGrad = serif::utilities::buildReducedMatrix(*gradMatrix, m_theta_ess_tdofs.first, m_theta_ess_tdofs.first); - - m_jacobian->SetBlock(0, 0, &reducedGrad); - - return *m_jacobian; -} -void PolytropeOperator::set_essential_true_dofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) { - m_isFinalized = false; - m_theta_ess_tdofs = theta_ess_tdofs; - m_phi_ess_tdofs = phi_ess_tdofs; - - if (m_f) { - m_f->SetEssentialTrueDofs(theta_ess_tdofs.first); - } else { - MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs"); - } -} - -void PolytropeOperator::set_essential_true_dofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) { - set_essential_true_dofs(ess_tdof_pair_set.first, ess_tdof_pair_set.second); -} - -SSE::MFEMArrayPairSet PolytropeOperator::get_essential_true_dofs() const { - return std::make_pair(m_theta_ess_tdofs, m_phi_ess_tdofs); -} -GMRESInverter::GMRESInverter(const SchurCompliment &op) : -mfem::Operator(op.Height(), op.Width()), -m_op(op) { - m_solver.SetOperator(m_op); - m_solver.SetMaxIter(100); - m_solver.SetRelTol(1e-1); - m_solver.SetAbsTol(1e-1); -} - -void GMRESInverter::Mult(const mfem::Vector &x, mfem::Vector &y) const { - m_solver.Mult(x, y); // Approximates m_op^-1 * x -} - - -SchurCompliment::SchurCompliment( - const mfem::SparseMatrix &QOp, - const mfem::SparseMatrix &DOp, - const mfem::SparseMatrix &MOp, - const mfem::Solver &GradInvOp) : -mfem::Operator(DOp.Height(), DOp.Width()) -{ - SetOperator(QOp, DOp, MOp, GradInvOp); - m_nPhi = m_DOp->Height(); - m_nTheta = m_MOp->Height(); -} - -SchurCompliment::SchurCompliment( - const mfem::SparseMatrix &QOp, - const mfem::SparseMatrix &DOp, - const mfem::SparseMatrix &MOp) : -mfem::Operator(DOp.Height(), DOp.Width()) -{ - updateConstantTerms(QOp, DOp, MOp); - m_nPhi = m_DOp->Height(); - m_nTheta = m_MOp->Height(); -} - -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); -} - -void SchurCompliment::updateInverseNonlinearJacobian(const mfem::Solver &gradInv) { - m_GradInvOp = &gradInv; -} - -void SchurCompliment::updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp) { - m_QOp = &QOp; - m_DOp = &DOp; - m_MOp = &MOp; -} - -void SchurCompliment::Mult(const mfem::Vector &x, mfem::Vector &y) const { - // Check that the input vector is the correct size - if (x.Size() != m_nPhi) { - MFEM_ABORT("Input vector x has size " + std::to_string(x.Size()) + ", expected " + std::to_string(m_nPhi)); - } - if (y.Size() != m_nPhi) { - MFEM_ABORT("Output vector y has size " + std::to_string(y.Size()) + ", expected " + std::to_string(m_nPhi)); - } - - // Check that the operators are set - if (m_QOp == nullptr) { - MFEM_ABORT("QOp is null in SchurCompliment::Mult"); - } - if (m_DOp == nullptr) { - MFEM_ABORT("DOp is null in SchurCompliment::Mult"); - } - if (m_MOp == nullptr) { - MFEM_ABORT("MOp is null in SchurCompliment::Mult"); - } - if (m_GradInvOp == nullptr) { - MFEM_ABORT("GradInvOp is null in SchurCompliment::Mult"); - } - - mfem::Vector v1(m_nTheta); // M * x - m_MOp -> Mult(x, v1); // M * x - - mfem::Vector v2(m_nTheta); // GradInv * M * x - m_GradInvOp -> Mult(v1, v2); // GradInv * M * x - - mfem::Vector v3(m_nPhi); // Q * GradInv * M * x - m_QOp -> Mult(v2, v3); // Q * GradInv * M * x - - mfem::Vector v4(m_nPhi); // D * x - m_DOp -> Mult(x, v4); // D * x - - subtract(v4, v3, y); // (D - Q * GradInv * M) * x -} diff --git a/src/poly/utils/private/utilities.cpp b/src/poly/utils/private/utilities.cpp index 5eb02e7..f0ffe08 100644 --- a/src/poly/utils/private/utilities.cpp +++ b/src/poly/utils/private/utilities.cpp @@ -2,7 +2,7 @@ #include "mfem.hpp" namespace serif::utilities { - mfem::SparseMatrix buildReducedMatrix( + mfem::SparseMatrix build_reduced_matrix( const mfem::SparseMatrix& matrix, const mfem::Array& trialEssentialDofs, const mfem::Array& testEssentialDofs @@ -93,4 +93,11 @@ namespace serif::utilities { return A_new; } + + mfem::Vector build_dof_identification_vector(const mfem::Array& allDofs, const::mfem::Array& highlightDofs) { + mfem::Vector v(allDofs.Size()); + v = 0.0; // Initialize the vector to zero + v.SetSubVector(highlightDofs, 1.0); // Set the highlighted dofs to 1.0 + return v; + } } diff --git a/src/poly/utils/public/polytropeOperator.h b/src/poly/utils/public/polytropeOperator.h index 0714c64..db332df 100644 --- a/src/poly/utils/public/polytropeOperator.h +++ b/src/poly/utils/public/polytropeOperator.h @@ -208,6 +208,57 @@ public: */ void set_essential_true_dofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set); + + /** + * @brief Reconstructs the full state vector (including essential DOFs) from a reduced state vector (free DOFs). + * @param reducedState The vector containing only the free degrees of freedom. + * @return Constant reference to the internal full state vector, updated with the reducedState. + */ + [[nodiscard]] const mfem::Vector &reconstruct_full_state_vector(const mfem::Vector &reducedState) const; + + /** + * @breif Reconstruct the full state vector (including essential DOFs) from a reduced state vector (free DOFs) as well as the block offsets. + * @param reducedState The vector containing only the free degrees of freedom. + * @return Constant reference to the internal full state vector, updated with the reducedState as a block vector. + */ + [[nodiscard]] const mfem::BlockVector reconstruct_full_block_state_vector(const mfem::Vector &reducedState) const; + + /** + * @brief Populates the internal array of free (non-essential) degree of freedom indices. + * This is called during finalize(). + */ + void populate_free_dof_array(); + + /// --- Getters for key internal state and operators --- + /** + * @brief Gets the Jacobian operator. + * Asserts that the operator is finalized and the Jacobian has been computed. + * @return Constant reference to the block Jacobian operator. + */ + const mfem::BlockOperator &get_jacobian_operator() const; + + /** + * @brief Gets the block diagonal preconditioner for the Schur complement system. + * Asserts that the operator is finalized and the preconditioner has been computed. + * @return Reference to the block diagonal preconditioner. + */ + mfem::BlockDiagonalPreconditioner &get_preconditioner() const; + + + /// --- Getters for information on internal state --- + /** + * @brief Gets the full state vector, including essential DOFs. + * @return Constant reference to the internal state vector. + */ + const mfem::Array& get_free_dofs() const { return m_freeDofs; } ///< Getter for the free DOFs array. + + /** + * @brief Gets the size of the reduced system (number of free DOFs). + * Asserts that the operator is finalized. + * @return The total number of free degrees of freedom. + */ + int get_reduced_system_size() const; + /** * @brief Gets the currently set essential true degrees of freedom. * @return A pair containing the essential TDOF pairs for theta and phi. @@ -226,40 +277,6 @@ public: */ const mfem::Array& get_reduced_block_offsets() const {return m_reducedBlockOffsets; } - /** - * @brief Gets the Jacobian operator. - * Asserts that the operator is finalized and the Jacobian has been computed. - * @return Constant reference to the block Jacobian operator. - */ - const mfem::BlockOperator &get_jacobian_operator() const; - - /** - * @brief Gets the block diagonal preconditioner for the Schur complement system. - * Asserts that the operator is finalized and the preconditioner has been computed. - * @return Reference to the block diagonal preconditioner. - */ - mfem::BlockDiagonalPreconditioner &get_preconditioner() const; - - /** - * @brief Gets the size of the reduced system (number of free DOFs). - * Asserts that the operator is finalized. - * @return The total number of free degrees of freedom. - */ - int get_reduced_system_size() const; - - /** - * @brief Reconstructs the full state vector (including essential DOFs) from a reduced state vector (free DOFs). - * @param reducedState The vector containing only the free degrees of freedom. - * @return Constant reference to the internal full state vector, updated with the reducedState. - */ - [[nodiscard]] const mfem::Vector &reconstruct_full_state_vector(const mfem::Vector &reducedState) const; - - /** - * @brief Populates the internal array of free (non-essential) degree of freedom indices. - * This is called during finalize(). - */ - void populate_free_dof_array(); - private: // --- Logging --- Probe::LogManager& m_logManager = Probe::LogManager::getInstance(); ///< Reference to the global log manager. @@ -280,6 +297,7 @@ private: std::unique_ptr m_MReduced; ///< Reduced M matrix (free DOFs only). std::unique_ptr m_QReduced; ///< Reduced Q matrix (free DOFs only). std::unique_ptr m_DReduced; ///< Reduced D matrix (free DOFs only). + mutable std::unique_ptr m_gradReduced; ///< Reduced gradient operator (G) for the nonlinear part f(θ). // --- State Vectors and DOF Management --- mutable mfem::Vector m_state; ///< Full state vector [θ, φ]^T, including essential DOFs. diff --git a/src/poly/utils/public/utilities.h b/src/poly/utils/public/utilities.h index c7e529b..a1a1e42 100644 --- a/src/poly/utils/public/utilities.h +++ b/src/poly/utils/public/utilities.h @@ -3,9 +3,41 @@ #include "mfem.hpp" namespace serif::utilities { - mfem::SparseMatrix buildReducedMatrix( + [[nodiscard]] mfem::SparseMatrix build_reduced_matrix( const mfem::SparseMatrix& matrix, const mfem::Array& trialEssentialDofs, const mfem::Array& testEssentialDofs ); + + /** + * @brief Generate a vector of 1s and 0s where 1 elemetns cooresponds to queried dofs. Useful for degugging + * @param allDofs array, counding from 0, of all dofs in the system + * @param highlightDofs the dofs that you want to identify + * @return + * + * *Example Usage:* + * One could use this to identify, for example, which dofs are being identified as the central dofs + * @code + * ... + * mfem::Array phiDofs, thetaDofs; + * phiDofs.SetSize(m_fePhi->GetNDofs()); + * thetaDofs.SetSize(m_feTheta->GetNDofs()); + * const mfem::Vector phiHighlightVector = serif::utilities::build_dof_identification_vector(phiDofs, phiCenterDofs); + * const mfem::Vector thetaHighlightVector = serif::utilities::build_dof_identification_vector(thetaDofs, thetaCenterDofs); + * Probe::glVisView( + * const_cast(phiHighlightVector), + * *m_fePhi, + * "Phi Center Dofs" + * ); + * Probe::glVisView( + * const_cast(thetaHighlightVector), + * *m_feTheta, + * "Theta Center Dofs" + * ); + * @endcode + */ + mfem::Vector build_dof_identification_vector( + const mfem::Array& allDofs, + const::mfem::Array& highlightDofs + ); }