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

refactor(PolytropeOperator): commented out debug code

Browse Source
This commit is contained in:
Emily Boudreaux
2025-05-11 14:40:19 -04:00
parent 8b9d46e996
commit 412a3be2ec
2 changed files with 41 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
src/poly/utils/private/operator.cpp
View File

@@ -55,6 +55,9 @@ void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr<mfem::Spa
}
}
// static int s_newtonStepGrad = 0;
// static int s_newtonStepMult = 0;
//
PolytropeOperator::PolytropeOperator(
std::unique_ptr<mfem::MixedBilinearForm> M,
@@ -82,6 +85,7 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
// Remove the essential dofs from the constant matrices
for (const auto& dof: m_theta_ess_tdofs.first) {
m_Mmat->EliminateRow(dof);
m_Qmat->EliminateCol(dof);
@@ -104,7 +108,7 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
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_Qmat.get()); //<- -Q (constant)
m_jacobian->SetBlock(1, 1, m_Dmat.get()); //<- D (constant)
m_jacobian->SetBlock(1, 1, m_Dmat.get()); //<- D (constant)
// Build the initial preconditioner based on some initial guess
const auto &grad = m_f->GetGradient(initTheta);
@@ -174,7 +178,7 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
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. This is commented out since seems it results in dramatic instabilies arrising
y_block.GetBlock(0)[idx] = x_theta(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilities arising
}
}
std::cout << "θ Norm: " << y_block.GetBlock(0).Norml2() << std::endl;
@@ -182,11 +186,31 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
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. This is commented out since seems it results in dramatic instabilies arrising
y_block.GetBlock(1)[idx] = x_phi(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilities arising
}
}
std::cout << "φ Norm: " << y_block.GetBlock(1).Norml2() << std::endl;
// std::cout << "φ Norm: " << y_block.GetBlock(1).Norml2() << std::endl;
//
// std::cout << "Writing out residuals to file..." << std::endl;
// std::ofstream outfile("Residuals_" + std::to_string(s_newtonStepMult) + ".csv", std::ios::trunc);
// if (!outfile.is_open()) {
// MFEM_ABORT("Could not open file for writing residuals");
// }
// outfile << "# Residuals for Newton Step: " << s_newtonStepMult << '\n';
// outfile << "# theta size: " << theta_size << '\n';
// outfile << "# phi size: " << phi_size << '\n';
// outfile << "dof,R\n";
// for (int i = 0; i < y_block.GetBlock(0).Size(); i++) {
// outfile << i << "," << y_block.GetBlock(0)[i] << '\n';
// }
// for (int i = 0; i < y_block.GetBlock(1).Size(); i++) {
// outfile << y_block.GetBlock(0).Size() + i << "," << y_block.GetBlock(1)[i] << '\n';
// }
// outfile.close();
// s_newtonStepMult++;
// std::cout << "Done writing out residuals to file." << std::endl;
//
}
@@ -243,7 +267,19 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const {
// auto *gradPtr = &grad;
// updatePreconditioner(grad);
m_jacobian->SetBlock(0, 0, &grad);
// // TODO: Save the state vector somewhere so that I can inspect how the actual solution
// // to theta and phi are evolving
// m_jacobian->SetBlock(0, 0, &grad);
// // TODO: Save the jacobian here
// std::cout << "Writing out jacobian to file..." << std::endl;
// std::vector<mfem::Operator*> ops;
// ops.push_back(&m_jacobian->GetBlock(0, 0));
// ops.push_back(&m_jacobian->GetBlock(0, 1));
// ops.push_back(&m_jacobian->GetBlock(1, 0));
// ops.push_back(&m_jacobian->GetBlock(1, 1));
// saveBlockFormToBinary(ops, {{0, 0}, {0, 1}, {1, 0}, {1, 1}}, {false, false, false, false}, "jacobian_" + std::to_string(s_newtonStepGrad) + ".bin");
// s_newtonStepGrad++;
// std::cout << "Done writing out jacobian to file." << std::endl;
return *m_jacobian;
}

1
src/poly/utils/public/operator.h
View File

@@ -38,7 +38,6 @@ public:
void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
mfem::Operator& GetGradient(const mfem::Vector &x) const override;
void SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs);