diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 74830fc..5ae30f1 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -31,7 +31,7 @@ #include "config.h" #include "integrators.h" #include "mfem.hpp" -#include "operator.h" +#include "polytropeOperator.h" #include "polyCoeff.h" #include "probe.h" #include "resourceManager.h" @@ -119,8 +119,7 @@ void PolySolver::assembleBlockSystem() { std::move(forms->Q), std::move(forms->D), std::move(forms->f), - blockOffsets, - m_polytropicIndex); + blockOffsets); } mfem::Array PolySolver::computeBlockOffsets() const { @@ -149,7 +148,7 @@ std::unique_ptr PolySolver::buildIndividualForms(const mfem::Array(negOneVec); // --- Add the integrators to the forms --- - forms->M->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator(*m_negationCoeff)); + forms->M->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator()); forms->Q->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); forms->D->AddDomainIntegrator(new mfem::VectorFEMassIntegrator()); @@ -184,12 +183,12 @@ void PolySolver::solve() const { m_polytropOperator->finalize(thetaVec); // 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 + 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) - mfem::Vector zero_rhs(block_offsets.Last()); + mfem::Vector zero_rhs(m_polytropOperator->get_reduced_system_size()); zero_rhs = 0.0; const solverBundle sb = setupNewtonSolver(); @@ -308,7 +307,7 @@ void PolySolver::setInitialGuess() const { } void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) const { - mfem::BlockVector x_block(const_cast(state_vector), m_polytropOperator->GetBlockOffsets()); + mfem::BlockVector x_block(const_cast(state_vector), m_polytropOperator->get_block_offsets()); mfem::Vector& x_theta = x_block.GetBlock(0); mfem::Vector& x_phi = x_block.GetBlock(1); @@ -335,7 +334,7 @@ void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) cons void PolySolver::setOperatorEssentialTrueDofs() const { const SSE::MFEMArrayPairSet ess_tdof_pair_set = getEssentialTrueDof(); - m_polytropOperator->SetEssentialTrueDofs(ess_tdof_pair_set); + m_polytropOperator->set_essential_true_dofs(ess_tdof_pair_set); } void PolySolver::LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel, double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const { diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 6448c9a..11c4071 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -26,7 +26,7 @@ #include "integrators.h" #include "4DSTARTypes.h" -#include "operator.h" +#include "polytropeOperator.h" #include "config.h" #include "meshIO.h" #include "probe.h" diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build index 8fc7d68..cd668cc 100644 --- a/src/poly/utils/meson.build +++ b/src/poly/utils/meson.build @@ -20,7 +20,8 @@ # *********************************************************************** # polyutils_sources = files( 'private/integrators.cpp', - 'private/operator.cpp' + 'private/polytropeOperator.cpp', + 'private/utilities.cpp', ) dependencies = [ diff --git a/src/poly/utils/private/polytropeOperator.cpp b/src/poly/utils/private/polytropeOperator.cpp index e81a032..6eab1e3 100644 --- a/src/poly/utils/private/polytropeOperator.cpp +++ b/src/poly/utils/private/polytropeOperator.cpp @@ -18,8 +18,10 @@ // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // // *********************************************************************** */ -#include "operator.h" +#include "polytropeOperator.h" #include "4DSTARTypes.h" +#include "utilities.h" + #include "mfem.hpp" #include "mfem_smout.h" #include @@ -31,9 +33,9 @@ PolytropeOperator::PolytropeOperator( std::unique_ptr Q, std::unique_ptr D, std::unique_ptr f, - const mfem::Array &blockOffsets, - const double index) : + const mfem::Array &blockOffsets) : + // TODO: Need to update this so that the size is that of the reduced system operator mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector m_blockOffsets(blockOffsets), m_jacobian(nullptr) { @@ -49,63 +51,96 @@ PolytropeOperator::PolytropeOperator( m_invNonlinearJacobian = std::make_unique(0, 3); } -void PolytropeOperator::finalize(const mfem::Vector &initTheta) { - if (m_isFinalized) { - return; +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()); - saveSparseMatrixBinary(*m_Mmat, "MmatRawO2.bin"); - saveSparseMatrixBinary(*m_Qmat, "QmatRawO2.bin"); - saveSparseMatrixBinary(*m_Dmat, "DmatRawO2.bin"); + 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)); - // Remove the essential dofs from the constant matrices - for (const auto& dof: m_theta_ess_tdofs.first) { - std::cout << "Eliminating dof: " << dof << std::endl; - m_Mmat->EliminateRow(dof); - m_Qmat->EliminateCol(dof); + 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; } - saveSparseMatrixBinary(*m_Mmat, "MmatBCThetaO2.bin"); - saveSparseMatrixBinary(*m_Qmat, "QmatBCThetaO2.bin"); - saveSparseMatrixBinary(*m_Dmat, "DmatBCThetaO2.bin"); + // 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(); - for (const auto& dof: m_phi_ess_tdofs.first) { - m_Mmat->EliminateCol(dof); - m_Qmat->EliminateRow(dof); - m_Dmat->EliminateRowCol(dof); - } - saveSparseMatrixBinary(*m_Mmat, "MmatBCO2.bin"); - saveSparseMatrixBinary(*m_Qmat, "QmatBCO2.bin"); - saveSparseMatrixBinary(*m_Dmat, "DmatBCO2.bin"); - - 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_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_Mmat.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); + // Override the size based on the reduced system + height = m_reducedBlockOffsets.Last(); + width = m_reducedBlockOffsets.Last(); m_isFinalized = true; - // Build the initial preconditioner based on some initial guess - const auto &grad = m_f->GetGradient(initTheta); - // updatePreconditioner(grad); - } -const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const { +const mfem::BlockOperator &PolytropeOperator::get_jacobian_operator() const { if (m_jacobian == nullptr) { MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call."); } @@ -116,7 +151,7 @@ const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const { return *m_jacobian; } -mfem::BlockDiagonalPreconditioner& PolytropeOperator::GetPreconditioner() const { +mfem::BlockDiagonalPreconditioner& PolytropeOperator::get_preconditioner() const { if (m_schurPreconditioner == nullptr) { MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call."); } @@ -126,12 +161,25 @@ mfem::BlockDiagonalPreconditioner& PolytropeOperator::GetPreconditioner() const return *m_schurPreconditioner; } -void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { +int PolytropeOperator::get_reduced_system_size() const { + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator not finalized prior to call to GetReducedSystemSize()."); + } + return m_reducedBlockOffsets.Last(); +} + +void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const { if (!m_isFinalized) { MFEM_ABORT("PolytropeOperator::Mult called before finalize"); } + + // TODO: confirm that the vectors xFree and m_freeDofs are always parallel + m_state.SetSubVector(m_freeDofs, xFree); // Scatter the free dofs from the input vector xFree into the state vector + mfem::Vector y; + y.SetSize(m_blockOffsets.Last()); + // -- Create BlockVector views for input x and output y - mfem::BlockVector x_block(const_cast(x), m_blockOffsets); + mfem::BlockVector x_block(const_cast(m_state), m_blockOffsets); mfem::BlockVector y_block(y, m_blockOffsets); // -- Get Vector views for individual blocks @@ -162,31 +210,22 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { add(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++) { - 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] = x_theta(idx) - targetValue; // inhomogenous essential bc. + yFree.SetSize(m_reducedBlockOffsets.Last()); + MFEM_ASSERT(m_freeDofs.Size() == m_reducedBlockOffsets.Last(), "PolytropeOperator::Mult: Size of free dofs does not match reduced block offsets size."); + for (int i = 0, j = 0; i < y.Size(); ++i) { + if (m_freeDofs.Find(i) != -1) { + yFree[j] = y[i]; + j++; } } - - 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()) { - const double &targetValue = m_phi_ess_tdofs.second[i]; - y_block.GetBlock(1)[idx] = x_phi(idx) - targetValue; // inhomogenous essential bc. - } - } - - std::cout << "||r_θ|| = " << y_block.GetBlock(0).Norml2(); - std::cout << ", ||r_φ|| = " << y_block.GetBlock(1).Norml2() << std::endl; } -void PolytropeOperator::updateInverseNonlinearJacobian(const mfem::Operator &grad) const { +void PolytropeOperator::update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const { m_invNonlinearJacobian->SetOperator(grad); } -void PolytropeOperator::updateInverseSchurCompliment() const { +void PolytropeOperator::update_inverse_schur_compliment() const { // TODO: This entire function could probably be refactored out if (!m_isFinalized) { MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before finalize"); @@ -209,27 +248,37 @@ void PolytropeOperator::updateInverseSchurCompliment() const { } -void PolytropeOperator::updatePreconditioner(const mfem::Operator &grad) const { - updateInverseNonlinearJacobian(grad); - updateInverseSchurCompliment(); +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 &x) const { +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(x), m_blockOffsets); + 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); - m_jacobian->SetBlock(0, 0, &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::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) { +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; @@ -241,11 +290,11 @@ void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess } } -void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) { - SetEssentialTrueDofs(ess_tdof_pair_set.first, ess_tdof_pair_set.second); +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::GetEssentialTrueDofs() const { +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) : diff --git a/src/poly/utils/private/utilities.cpp b/src/poly/utils/private/utilities.cpp index 227b549..5eb02e7 100644 --- a/src/poly/utils/private/utilities.cpp +++ b/src/poly/utils/private/utilities.cpp @@ -1,5 +1,96 @@ -// -// Created by Emily Boudreaux on 6/4/25. -// - #include "utilities.h" +#include "mfem.hpp" + +namespace serif::utilities { + mfem::SparseMatrix buildReducedMatrix( + const mfem::SparseMatrix& matrix, + const mfem::Array& trialEssentialDofs, + const mfem::Array& testEssentialDofs + ) { + int M_orig = matrix.Height(); + int N_orig = matrix.Width(); + + // 1. Create boolean lookup tables for eliminated rows/columns + // These tables help quickly check if an original row/column index is eliminated. + mfem::Array row_is_eliminated(M_orig); + row_is_eliminated = false; // Initialize all to false (no rows eliminated yet) + for (int i = 0; i < testEssentialDofs.Size(); ++i) { + int r_idx = testEssentialDofs[i]; + if (r_idx >= 0 && r_idx < M_orig) { // Check for valid index + row_is_eliminated[r_idx] = true; + } + } + + mfem::Array col_is_eliminated(N_orig); + col_is_eliminated = false; // Initialize all to false (no columns eliminated yet) + for (int i = 0; i < trialEssentialDofs.Size(); ++i) { + int c_idx = trialEssentialDofs[i]; + if (c_idx >= 0 && c_idx < N_orig) { // Check for valid index + col_is_eliminated[c_idx] = true; + } + } + + // 2. Create mappings from old (original) indices to new indices + // Also, count the number of rows and columns in the new, reduced matrix. + mfem::Array old_row_to_new_row(M_orig); + int num_new_rows = 0; + for (int i = 0; i < M_orig; ++i) { + if (row_is_eliminated[i]) { + old_row_to_new_row[i] = -1; // Mark as eliminated (no corresponding new row) + } else { + old_row_to_new_row[i] = num_new_rows++; // Assign the next available new row index + } + } + + mfem::Array old_col_to_new_col(N_orig); + int num_new_cols = 0; + for (int i = 0; i < N_orig; ++i) { + if (col_is_eliminated[i]) { + old_col_to_new_col[i] = -1; // Mark as eliminated (no corresponding new column) + } else { + old_col_to_new_col[i] = num_new_cols++; // Assign the next available new column index + } + } + + // 3. Create the new sparse matrix with the calculated reduced dimensions. + // It's initially empty (no non-zero entries). + mfem::SparseMatrix A_new(num_new_rows, num_new_cols); + + // 4. Iterate through the non-zero entries of the original matrix (A_orig). + // A_orig is typically stored in Compressed Sparse Row (CSR) format. + // GetI() returns row pointers, GetJ() returns column indices, GetData() returns values. + const int* I_orig = matrix.GetI(); + const int* J_orig = matrix.GetJ(); + const double* V_orig = matrix.GetData(); // Assuming double, can be templated if needed + + for (int r_orig = 0; r_orig < M_orig; ++r_orig) { + // If the original row is marked for elimination, skip all its entries. + if (row_is_eliminated[r_orig]) { + continue; + } + + // Get the new row index for the current original row. + int r_new = old_row_to_new_row[r_orig]; + + // Iterate through non-zero entries in the current original row r_orig. + // I_orig[r_orig] is the start index in J_orig and V_orig for row r_orig. + // I_orig[r_orig+1]-1 is the end index. + for (int k = I_orig[r_orig]; k < I_orig[r_orig + 1]; ++k) { + int c_orig = J_orig[k]; // Original column index of the non-zero entry + double val = V_orig[k]; // Value of the non-zero entry + + if (col_is_eliminated[c_orig]) { + continue; + } + + int c_new = old_col_to_new_col[c_orig]; + + A_new.Add(r_new, c_new, val); + } + } + + A_new.Finalize(); + + return A_new; + } +} diff --git a/src/poly/utils/public/polytropeOperator.h b/src/poly/utils/public/polytropeOperator.h index 051845f..6b905ee 100644 --- a/src/poly/utils/public/polytropeOperator.h +++ b/src/poly/utils/public/polytropeOperator.h @@ -70,28 +70,34 @@ public: std::unique_ptr Q, std::unique_ptr D, std::unique_ptr f, - const mfem::Array &blockOffsets, - const double index); + const mfem::Array &blockOffsets); + + void populate_free_dof_array(); + + ~PolytropeOperator() override = default; - void Mult(const mfem::Vector &x, mfem::Vector &y) const override; + void Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const override; - mfem::Operator& GetGradient(const mfem::Vector &x) const override; + mfem::Operator& GetGradient(const mfem::Vector &xFree) const override; - void SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs); - void SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set); + void set_essential_true_dofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs); + void set_essential_true_dofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set); - SSE::MFEMArrayPairSet GetEssentialTrueDofs() const; + SSE::MFEMArrayPairSet get_essential_true_dofs() const; bool isFinalized() const { return m_isFinalized; } void finalize(const mfem::Vector &initTheta); - const mfem::Array& GetBlockOffsets() const { return m_blockOffsets; } + const mfem::Array& get_block_offsets() const { return m_blockOffsets; } + const mfem::Array& get_reduced_block_offsets() const {return m_reducedBlockOffsets; } - const mfem::BlockOperator &GetJacobianOperator() const; + const mfem::BlockOperator &get_jacobian_operator() const; - mfem::BlockDiagonalPreconditioner &GetPreconditioner() const; + mfem::BlockDiagonalPreconditioner &get_preconditioner() const; + + int get_reduced_system_size() const; private: @@ -102,26 +108,24 @@ 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; + mutable mfem::Vector m_state; + mfem::Array m_freeDofs; + + // These are used to store the reduced system matrices after essential dofs have been eliminated + std::unique_ptr m_MReduced; + std::unique_ptr m_QReduced; + std::unique_ptr m_DReduced; + const mfem::Array m_blockOffsets; + mfem::Array m_reducedBlockOffsets; SSE::MFEMArrayPair m_theta_ess_tdofs; SSE::MFEMArrayPair m_phi_ess_tdofs; - std::unique_ptr m_negM_mat; std::unique_ptr m_negQ_mat; mutable std::unique_ptr m_jacobian; @@ -134,8 +138,12 @@ private: bool m_isFinalized = false; private: - void updateInverseNonlinearJacobian(const mfem::Operator &grad) const; - void updateInverseSchurCompliment() const; - void updatePreconditioner(const mfem::Operator &grad) const; + void update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const; + void update_inverse_schur_compliment() const; + void update_preconditioner(const mfem::Operator &grad) const; + void scatter_boundary_conditions(); + void construct_matrix_representations(); + void construct_reduced_block_offsets(); + void construct_jacobian_constant_terms(); }; diff --git a/src/poly/utils/public/utilities.h b/src/poly/utils/public/utilities.h new file mode 100644 index 0000000..c7e529b --- /dev/null +++ b/src/poly/utils/public/utilities.h @@ -0,0 +1,11 @@ +#pragma once + +#include "mfem.hpp" + +namespace serif::utilities { + mfem::SparseMatrix buildReducedMatrix( + const mfem::SparseMatrix& matrix, + const mfem::Array& trialEssentialDofs, + const mfem::Array& testEssentialDofs + ); +}