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

feat(poly): refactoring PolytropeOperator to work on the reduced system so as to avoid rank deficiencies

Browse Source
This commit is contained in:
Emily Boudreaux
2025-06-05 12:37:00 -04:00
parent 4eb8b71271
commit a31c966146
7 changed files with 272 additions and 113 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -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 <memory>
@@ -31,9 +33,9 @@ PolytropeOperator::PolytropeOperator(
std::unique_ptr<mfem::MixedBilinearForm> Q,
std::unique_ptr<mfem::BilinearForm> D,
std::unique_ptr<mfem::NonlinearForm> f,
const mfem::Array<int> &blockOffsets,
const double index) :
const mfem::Array<int> &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<mfem::GSSmoother>(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<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());
saveSparseMatrixBinary(*m_Mmat, "MmatRawO2.bin");
saveSparseMatrixBinary(*m_Qmat, "QmatRawO2.bin");
saveSparseMatrixBinary(*m_Dmat, "DmatRawO2.bin");
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));
// 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<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;
}
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<mfem::ScaledOperator>(m_Mmat.get(), -1.0);
m_negQ_mat = std::make_unique<mfem::ScaledOperator>(m_Qmat.get(), -1.0);
// m_schurCompliment = std::make_unique<SchurCompliment>(*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<mfem::BlockOperator>(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<GMRESInverter>(*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<mfem::Vector&>(x), m_blockOffsets);
mfem::BlockVector x_block(const_cast<mfem::Vector&>(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<mfem::Vector&>(x), m_blockOffsets);
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);
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) :