tboudreaux/SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

fix(poly): fixed numerous bugs related to inconsistent system sizing with the reduced operator

Browse Source 
this has restored the symmetry which we relied on before.
This commit is contained in:
Emily Boudreaux
2025-06-09 10:19:18 -04:00
parent 6e1453cf6e
commit 2a91d57ad7
5 changed files with 481 additions and 316 deletions
Download Patch File Download Diff File Expand all files Collapse all files

122
src/poly/solver/private/polySolver.cpp
View File

@@ -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<formBundle> forms = buildIndividualForms(blockOffsets);
// const double penalty_param = m_config.get<double>("Poly::Solver::ZeroDerivativePenalty", 1e10);
// mfem::Array<int> 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<PolytropeOperator>(
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<mfem::Vector>(*m_theta); // NOLINT(*-slicing)
const auto phiVec = static_cast<mfem::Vector>(*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<int>& 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<int>& 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<int>& 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<int> phi_ess_tdof_list;
mfem::Array<int> thetaCenterDofs, phiCenterDofs; // phiCenterDofs are not used
mfem::Array<double> thetaCenterVals;
mfem::Array<double> thetaCenterVals, phiCenterVals;
std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
thetaCenterVals.SetSize(thetaCenterDofs.Size());
// phiCenterVals.SetSize(phiCenterDofs.Size());
// phiCenterVals = 0.0;
thetaCenterVals = 1.0;
mfem::Array<int> 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<int>, mfem::Array<int>> PolySolver::findCenterElement() const {
mfem::Array<int> thetaCenterDofs;
mfem::Array<int> 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<int> elementIDs;
mfem::Array<mfem::IntegrationPoint> ips;
m_mesh.FindPoints(centerPoint, elementIDs, ips);
mfem::Array<int> 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<int> 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<int> 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<int> 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<bool>("Poly:Solver:ViewInitialGuess", false)) {
Probe::glVisView(*m_theta, m_mesh, "θ init");
Probe::glVisView(*m_phi, m_mesh, "φ init");
}
std::cout << "HERE" << std::endl;
}

546
src/poly/utils/private/polytropeOperator.cpp
View File

@@ -26,7 +26,109 @@
#include "mfem_smout.h"
#include <memory>
// --- 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<mfem::MixedBilinearForm> M,
@@ -51,128 +153,7 @@ PolytropeOperator::PolytropeOperator(
m_invNonlinearJacobian = std::make_unique<mfem::GSSmoother>(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<int> 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<mfem::SparseMatrix>(m_M->SpMat());
m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
m_MReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::buildReducedMatrix(*m_Mmat, m_phi_ess_tdofs.first, m_theta_ess_tdofs.first));
m_QReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::buildReducedMatrix(*m_Qmat, m_theta_ess_tdofs.first, m_phi_ess_tdofs.first));
m_DReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::buildReducedMatrix(*m_Dmat, m_phi_ess_tdofs.first, m_phi_ess_tdofs.first));
m_negQ_mat = std::make_unique<mfem::ScaledOperator>(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<mfem::BlockOperator>(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<mfem::Vector&>(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<mfem::SparseMatrix*>(&grad);
if (gradMatrix == nullptr) {
MFEM_ABORT("PolytropeOperator::GetGradient: Gradient is not a SparseMatrix.");
}
m_gradReduced = std::make_unique<mfem::SparseMatrix> (
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<mfem::SparseMatrix>(m_M->SpMat());
m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
m_MReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::build_reduced_matrix(*m_Mmat, m_phi_ess_tdofs.first, m_theta_ess_tdofs.first));
m_QReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::build_reduced_matrix(*m_Qmat, m_theta_ess_tdofs.first, m_phi_ess_tdofs.first));
m_DReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::build_reduced_matrix(*m_Dmat, m_phi_ess_tdofs.first, m_phi_ess_tdofs.first));
m_negQ_mat = std::make_unique<mfem::ScaledOperator>(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<mfem::BlockOperator>(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<int> 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<mfem::Vector&>(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<mfem::SparseMatrix*>(&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
}

9
src/poly/utils/private/utilities.cpp
View File

@@ -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<int>& trialEssentialDofs,
const mfem::Array<int>& testEssentialDofs
@@ -93,4 +93,11 @@ namespace serif::utilities {
return A_new;
}
mfem::Vector build_dof_identification_vector(const mfem::Array<int>& allDofs, const::mfem::Array<int>& 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;
}
}

86
src/poly/utils/public/polytropeOperator.h
View File

@@ -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<int>& 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<int>& 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<mfem::SparseMatrix> m_MReduced; ///< Reduced M matrix (free DOFs only).
std::unique_ptr<mfem::SparseMatrix> m_QReduced; ///< Reduced Q matrix (free DOFs only).
std::unique_ptr<mfem::SparseMatrix> m_DReduced; ///< Reduced D matrix (free DOFs only).
mutable std::unique_ptr<mfem::SparseMatrix> 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.

34
src/poly/utils/public/utilities.h
View File

@@ -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<int>& trialEssentialDofs,
const mfem::Array<int>& 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<int> 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<mfem::Vector&>(phiHighlightVector),
* *m_fePhi,
* "Phi Center Dofs"
* );
* Probe::glVisView(
* const_cast<mfem::Vector&>(thetaHighlightVector),
* *m_feTheta,
* "Theta Center Dofs"
* );
* @endcode
*/
mfem::Vector build_dof_identification_vector(
const mfem::Array<int>& allDofs,
const::mfem::Array<int>& highlightDofs
);
}