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

refactor(jacobiInvert): moved all jacobi inverting logic into dedicated function

Browse Source
This commit is contained in:
Emily Boudreaux
2025-04-21 07:33:39 -04:00
parent ae317fda5e
commit e7ad73c0f9
1 changed files with 26 additions and 55 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -125,50 +125,6 @@ const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const {
return *m_jacobian;
}
mfem::Vector PolytropeOperator::GetJ00Diag() const {
}
mfem::Vector PolytropeOperator::GetJ01Diag() const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::Get01Diag -> GetJ01Diag called before finalization of PolytropeOperator");
}
if (m_Mmat == nullptr) {
MFEM_ABORT("PolytropeOperator::Get01Diag -> M sparse matrix has not been initialized before GetJ01Diag() call.");
}
mfem::Vector J01diag;
m_Mmat->GetDiag(J01diag);
J01diag *= -1;
return J01diag;
}
mfem::Vector PolytropeOperator::GetJ10diag() const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::Get10Diag -> GetJ10Diag called before finalization of PolytropeOperator");
}
if (m_Qmat == nullptr) {
MFEM_ABORT("PolytropeOperator::Get10Diag -> Q sparse matrix has not been initialized before GetJ10Diag() call.");
}
mfem::Vector J10diag;
m_Qmat->GetDiag(J10diag);
J10diag *= -1;
return J10diag;
}
mfem::Vector PolytropeOperator::GetJ11diag() const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::Get11Diag -> GetJ11Diag called before finalization of PolytropeOperator");
}
if (m_Dmat == nullptr) {
MFEM_ABORT("PolytropeOperator::Get11Diag -> D sparse matrix has not been initialized before GetJ11Diag() call.");
}
mfem::Vector J11diag;
m_Dmat->GetDiag(J11diag);
J11diag *= -1;
return J11diag;
}
void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::Mult called before finalize");
@@ -209,7 +165,6 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
subtract(Dphi_term, Qtheta_term, y_R1);
// -- Apply essential boundary conditions --
for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) {
int idx = m_theta_ess_tdofs.first[i];
@@ -232,16 +187,35 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
}
// TODO: I was *very* stupid and I accidently deleted a lot of the code
// which I had written to find and use the preconditioner. This needs
// to be reimplemented. Once that is working you can get back to
// Trying to understand the multimodal hump in the residuals vector.
// There is a jupyter notebook about this. I was thinking that this was
// perhapse related to the non consistent application of boundary conditions.
void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr<mfem::SparseMatrix>& invMat, std::string name="matrix") {
// TODO: This likely can be made much more efficient and will probably be called in tight loops, a good
// place for some easy optimization might be here.
// Get the diagonal of the matrix
mfem::Vector diag;
mat.GetDiag(diag);
// Invert the diagonal
for (int i = 0; i < diag.Size(); i++) {
if (diag(i) == 0.0) {
MFEM_ABORT("Diagonal element (" + std::to_string(i) +") in " + name + " is zero, cannot invert.");
}
diag(i) = 1.0 / diag(i);
}
invMat = std::make_unique<mfem::SparseMatrix>(diag);
}
void PolytropeOperator::updateInverseNonlinearJacobian(const mfem::Operator &grad) const {
if (const auto *sparse_mat = dynamic_cast<const mfem::SparseMatrix*>(&grad); sparse_mat != nullptr) {
mfem::Vector gradDiag;
sparse_mat->GetDiag(gradDiag);
for (int i = 0; i < gradDiag.Size(); i++) {
gradDiag(i) = 1.0/gradDiag(i); // Invert the diagonals of the jacobian
}
m_invNonlinearJacobian = std::make_unique<mfem::SparseMatrix>(gradDiag);
approxJacobiInvert(*sparse_mat, m_invNonlinearJacobian, "Nonlinear Jacobian");
} else {
MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear Jacobian");
MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear jacobian where nonlinear jacobian is not trivially castable to a sparse matrix");
}
}
@@ -258,9 +232,6 @@ void PolytropeOperator::updateInverseSchurCompliment() const {
schurCompliment->Add(-1, *temp_QxdfInvxM); // D - Q * df^{-1} * M
writeOperatorToCSV(*schurCompliment, "/Users/tboudreaux/Desktop/SchursCompliment.csv");
std::cout << "Wrote" << std::endl;
// Note: EMB (April 9, 2025) Where I left for the day, so I don't pull my hair out tomorrow:
/*
* Need to calculate the inverse of Schur's compliment to precondition the jacobian