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

feat(PolytropeOperator): added polytropic index as a member variable

Browse Source
This commit is contained in:
Emily Boudreaux
2025-05-12 14:27:01 -04:00
parent 14eb21bd31
commit 1ee919a4a9
2 changed files with 24 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -55,8 +55,8 @@ void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr<mfem::Spa
} }
} }
// static int s_newtonStepGrad = 0; static int s_newtonStepGrad = 0;
// static int s_newtonStepMult = 0; static int s_newtonStepMult = 0;
// //
PolytropeOperator::PolytropeOperator( PolytropeOperator::PolytropeOperator(
@@ -64,11 +64,13 @@ PolytropeOperator::PolytropeOperator(
std::unique_ptr<mfem::MixedBilinearForm> Q, std::unique_ptr<mfem::MixedBilinearForm> Q,
std::unique_ptr<mfem::BilinearForm> D, std::unique_ptr<mfem::BilinearForm> D,
std::unique_ptr<mfem::NonlinearForm> f, std::unique_ptr<mfem::NonlinearForm> f,
const mfem::Array<int> &blockOffsets) : const mfem::Array<int> &blockOffsets,
const double index) :
mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector
m_blockOffsets(blockOffsets), m_blockOffsets(blockOffsets),
m_jacobian(nullptr) { m_jacobian(nullptr),
m_polytropicIndex(index){
m_M = std::move(M); m_M = std::move(M);
m_Q = std::move(Q); m_Q = std::move(Q);
@@ -181,7 +183,6 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
y_block.GetBlock(0)[idx] = x_theta(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilities arising 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;
for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) { 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()) { if (int idx = m_phi_ess_tdofs.first[i]; idx >= 0 && idx < y_R1.Size()) {
@@ -190,10 +191,11 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
} }
} }
// std::cout << "φ Norm: " << y_block.GetBlock(1).Norml2() << std::endl; std::cout << "||r_θ|| = " << y_block.GetBlock(0).Norml2();
std::cout << ", ||r_φ|| = " << y_block.GetBlock(1).Norml2() << std::endl;
// //
// std::cout << "Writing out residuals to file..." << std::endl; // std::cout << "Writing out residuals to file..." << std::endl;
// std::ofstream outfile("Residuals_" + std::to_string(s_newtonStepMult) + ".csv", std::ios::trunc); // std::ofstream outfile("Residuals_n" + std::to_string(m_polytropicIndex) + "_" + std::to_string(s_newtonStepMult) + ".csv", std::ios::trunc);
// if (!outfile.is_open()) { // if (!outfile.is_open()) {
// MFEM_ABORT("Could not open file for writing residuals"); // MFEM_ABORT("Could not open file for writing residuals");
// } // }
@@ -267,24 +269,23 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const {
// auto *gradPtr = &grad; // auto *gradPtr = &grad;
// updatePreconditioner(grad); // updatePreconditioner(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); m_jacobian->SetBlock(0, 0, &grad);
// // TODO: Save the jacobian here
// std::cout << "Writing out jacobian to file..." << std::endl; std::cout << "Writing out jacobian to file..." << std::endl;
// std::vector<mfem::Operator*> ops; std::vector<mfem::Operator*> ops;
// ops.push_back(&m_jacobian->GetBlock(0, 0)); ops.push_back(&m_jacobian->GetBlock(0, 0));
// ops.push_back(&m_jacobian->GetBlock(0, 1)); ops.push_back(&m_jacobian->GetBlock(0, 1));
// ops.push_back(&m_jacobian->GetBlock(1, 0)); ops.push_back(&m_jacobian->GetBlock(1, 0));
// ops.push_back(&m_jacobian->GetBlock(1, 1)); 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"); saveBlockFormToBinary(ops, {{0, 0}, {0, 1}, {1, 0}, {1, 1}}, {false, false, false, false}, "jacobian_" + std::to_string(m_polytropicIndex) + "_" + std::to_string(s_newtonStepGrad) + ".bin");
// s_newtonStepGrad++; s_newtonStepGrad++;
// std::cout << "Done writing out jacobian to file." << std::endl; std::cout << "Done writing out jacobian to file." << std::endl;
// // Exit the code here
// // throw std::runtime_error("Done writing out jacobian to file");
return *m_jacobian; return *m_jacobian;
} }
void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) { void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) {
m_isFinalized = false; m_isFinalized = false;
m_theta_ess_tdofs = theta_ess_tdofs; m_theta_ess_tdofs = theta_ess_tdofs;

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

@@ -33,7 +33,8 @@ public:
std::unique_ptr<mfem::MixedBilinearForm> Q, std::unique_ptr<mfem::MixedBilinearForm> Q,
std::unique_ptr<mfem::BilinearForm> D, std::unique_ptr<mfem::BilinearForm> D,
std::unique_ptr<mfem::NonlinearForm> f, std::unique_ptr<mfem::NonlinearForm> f,
const mfem::Array<int> &blockOffsets); const mfem::Array<int> &blockOffsets,
const double index);
~PolytropeOperator() override = default; ~PolytropeOperator() override = default;
void Mult(const mfem::Vector &x, mfem::Vector &y) const override; void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
@@ -90,6 +91,7 @@ private:
mutable std::unique_ptr<mfem::BlockDiagonalPreconditioner> m_schurPreconditioner; mutable std::unique_ptr<mfem::BlockDiagonalPreconditioner> m_schurPreconditioner;
bool m_isFinalized = false; bool m_isFinalized = false;
double m_polytropicIndex;
private: private:
void updateInverseNonlinearJacobian(const mfem::Operator &grad) const; void updateInverseNonlinearJacobian(const mfem::Operator &grad) const;