fix(poly): phi boundary values now coorespond to theta flux through polytrope surface
This commit is contained in:
@@ -38,6 +38,8 @@ namespace polycoeff{
|
||||
// v = -std::pow(r, 2);
|
||||
}
|
||||
|
||||
//PERF: One area of future optimization might be turning these into matrix operations since they are
|
||||
//PERF: fundamentally just the dot product of [a b c ...] * [n^0 n^1 n^2 ...]
|
||||
double x1(const double n)
|
||||
{
|
||||
double r = 0;
|
||||
@@ -47,4 +49,12 @@ namespace polycoeff{
|
||||
return r;
|
||||
}
|
||||
|
||||
double thetaSurfaceFlux(const double n) {
|
||||
double surfaceFlux = 0;
|
||||
for (int i = 0; i < dThetaInterpCoeff::numTerms; i++) {
|
||||
surfaceFlux += dThetaInterpCoeff::coeff[i] * std::pow(n, dThetaInterpCoeff::power[i]);
|
||||
}
|
||||
return surfaceFlux;
|
||||
}
|
||||
|
||||
}
|
||||
@@ -42,15 +42,25 @@ namespace polycoeff
|
||||
void diffusionCoeff(const mfem::Vector &x, mfem::Vector &v);
|
||||
|
||||
double x1(const double n);
|
||||
/**
|
||||
* @brief Coefficients for the interpolations of the surface location of a polytrope
|
||||
* @param numTerms Number of terms in the polynomial interpolator
|
||||
* @param power Array of the powers of the polynomial interpolator
|
||||
* @param coeff Array of the coefficients of the polynomial interpolator
|
||||
*/
|
||||
struct x1InterpCoeff {
|
||||
constexpr static int numTerms = 51;
|
||||
constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50};
|
||||
constexpr static double coeff[51] = {2.449260049071003e+00, 4.340987403815444e-01, 2.592462347475860e+00, -2.283794512491552e+01, 1.135535208066129e+02, -3.460093673587010e+02, 6.948139716983651e+02, -9.564368645158952e+02, 9.184546798891624e+02, -6.130308987803503e+02, 2.747228735318193e+02, -7.416795903196430e+01, 7.318638538580859e+00, 1.749441306797260e+00, -4.240148582456829e-01, -5.818809544982156e-02, 1.514877217199105e-02, 3.228634707578998e-03, -2.862524323980516e-04, -1.622486968261819e-04, -1.253644717104076e-05, 4.334945141292894e-06, 1.296452565229763e-06, 8.634802209209870e-08, -3.337511676486084e-08, -1.094796628367775e-08, -1.228178760540410e-09, 1.416744125622751e-10, 8.513777351265677e-11, 1.624582561364811e-11, 8.377207519041114e-13, -4.363812865112836e-13, -1.535862757816461e-13, -2.485085045037669e-14, -6.566281276491033e-16, 8.405047965853478e-16, 2.673804441025638e-16, 4.176337890285142e-17, 8.150073570140493e-19, -1.531673805257016e-18, -4.746996933716653e-19, -6.976825828390195e-20, 8.807513368604331e-22, 3.307739310180278e-21, 8.593260940093030e-22, 7.385969061093440e-23, -2.211130577977291e-23, -8.557291048388455e-24, -4.169359901215994e-25, 4.609379358875657e-25, -3.590870335035984e-26};
|
||||
};
|
||||
|
||||
double thetaSurfaceFlux(const double n);
|
||||
|
||||
/**
|
||||
* @brief Coefficients for the interpolations of the surface location of a polytrope
|
||||
* @param numTerms Number of terms in the polynomial interpolator
|
||||
* @param power Array of the powers of the polynomial interpolator
|
||||
* @param coeff Array of the coefficients of the polynomial interpolator
|
||||
*/
|
||||
const struct x1InterpCoeff {
|
||||
constexpr static int numTerms = 51;
|
||||
constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50};
|
||||
constexpr static double coeff[51] = {2.449260049071003e+00, 4.340987403815444e-01, 2.592462347475860e+00, -2.283794512491552e+01, 1.135535208066129e+02, -3.460093673587010e+02, 6.948139716983651e+02, -9.564368645158952e+02, 9.184546798891624e+02, -6.130308987803503e+02, 2.747228735318193e+02, -7.416795903196430e+01, 7.318638538580859e+00, 1.749441306797260e+00, -4.240148582456829e-01, -5.818809544982156e-02, 1.514877217199105e-02, 3.228634707578998e-03, -2.862524323980516e-04, -1.622486968261819e-04, -1.253644717104076e-05, 4.334945141292894e-06, 1.296452565229763e-06, 8.634802209209870e-08, -3.337511676486084e-08, -1.094796628367775e-08, -1.228178760540410e-09, 1.416744125622751e-10, 8.513777351265677e-11, 1.624582561364811e-11, 8.377207519041114e-13, -4.363812865112836e-13, -1.535862757816461e-13, -2.485085045037669e-14, -6.566281276491033e-16, 8.405047965853478e-16, 2.673804441025638e-16, 4.176337890285142e-17, 8.150073570140493e-19, -1.531673805257016e-18, -4.746996933716653e-19, -6.976825828390195e-20, 8.807513368604331e-22, 3.307739310180278e-21, 8.593260940093030e-22, 7.385969061093440e-23, -2.211130577977291e-23, -8.557291048388455e-24, -4.169359901215994e-25, 4.609379358875657e-25, -3.590870335035984e-26};
|
||||
};
|
||||
|
||||
const struct dThetaInterpCoeff {
|
||||
constexpr static int numTerms = 51;
|
||||
constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50};
|
||||
constexpr static double coeff[51] = {-8.164634353536671e-01, 8.651804355981230e-01, -6.391939109867470e-01, 4.882118355278944e-01, -3.898475553686603e-01, 3.102560043891329e-01, -2.365565113144737e-01, 1.639832755778444e-01, -9.645271947133248e-02, 4.450140636696046e-02, -1.460885194751961e-02, 2.816293786806504e-03, -8.178583546493117e-05, -8.584494556484958e-05, 8.476252127593868e-06, 2.923593421315422e-06, -2.206768214995963e-07, -1.203227957690371e-07, -3.381181730985542e-09, 4.022824706790907e-09, 7.041049107708875e-10, -3.562885681170365e-11, -3.281525407784209e-11, -5.031807464141896e-12, 2.034136401832885e-13, 2.361284283178230e-13, 4.602774507763180e-14, 1.809170850970874e-15, -1.333813332262995e-15, -4.045891156434286e-16, -5.197949512114809e-17, 2.222220713310119e-18, 2.625223897130583e-18, 6.226001466529447e-19, 6.571419077089260e-20, -6.672159423054950e-21, -4.476242224056620e-21, -9.790792477821165e-22, -9.222211318122281e-23, 1.427942034536028e-23, 7.759197090219954e-24, 1.546886518887300e-24, 9.585471984274525e-26, -4.005276449706623e-26, -1.459299762834743e-26, -1.870491354620814e-27, 2.271573838802745e-28, 1.360979028415734e-28, 1.102172718357361e-29, -6.919347646474293e-30, 4.875282352118995e-31};
|
||||
};
|
||||
|
||||
} // namespace polyCoeff
|
||||
@@ -139,7 +139,7 @@ std::unique_ptr<formBundle> PolySolver::buildIndividualForms(const mfem::Array<i
|
||||
);
|
||||
|
||||
// --- -M negation -> M ---
|
||||
mfem::Vector negOneVec(m_fePhi->GetVDim());
|
||||
mfem::Vector negOneVec(m_mesh->SpaceDimension());
|
||||
negOneVec = -1.0;
|
||||
|
||||
m_negationCoeff = std::make_unique<mfem::VectorConstantCoefficient>(negOneVec);
|
||||
@@ -184,9 +184,6 @@ void PolySolver::solve() const {
|
||||
throw std::runtime_error("PolytropeOperator is not finalized. Cannot solve.");
|
||||
}
|
||||
|
||||
GetDofCoordinates(*m_feTheta, "dof2posTheta.csv");
|
||||
GetDofCoordinates(*m_fePhi, "dof2posPhi.csv");
|
||||
|
||||
// It's safer to get the offsets directly from the operator after finalization
|
||||
const mfem::Array<int>& block_offsets = m_polytropOperator->GetBlockOffsets(); // Assuming a getter exists or accessing member if public/friend
|
||||
mfem::BlockVector state_vector(block_offsets);
|
||||
@@ -194,23 +191,14 @@ void PolySolver::solve() const {
|
||||
state_vector.GetBlock(1) = phiVec; // NOLINT(*-slicing)
|
||||
|
||||
|
||||
mfem::Vector zero_rhs(block_offsets.Last());
|
||||
mfem::Vector zero_rhs(block_offsets.Last()); // TODO: Is this right?
|
||||
zero_rhs = 0.0;
|
||||
|
||||
const solverBundle sb = setupNewtonSolver();
|
||||
|
||||
// EMB 2025: Calling Mult gets the gradient of the operator for the NewtonSolver
|
||||
// This then is set as the operator for the solver for NewtonSolver
|
||||
// The solver (assuming it is an iterative solver) then sets the
|
||||
// operator for its preconditioner to this.
|
||||
// What this means is that there is no need to manually deal
|
||||
// with updating the preconditioner at every newton step as the
|
||||
// changes to the jacobian are automatically propagated through the
|
||||
// solving chain. This is at least true with MFEM 4.8-rc0
|
||||
std::string windowTitle = "testWindow";
|
||||
std::string keyset = "";
|
||||
sb.newton.Mult(zero_rhs, state_vector);
|
||||
|
||||
// Ax = b for x
|
||||
|
||||
// --- Save and view the solution ---
|
||||
saveAndViewSolution(state_vector);
|
||||
}
|
||||
@@ -219,28 +207,31 @@ SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const {
|
||||
mfem::Array<int> theta_ess_tdof_list;
|
||||
mfem::Array<int> phi_ess_tdof_list;
|
||||
|
||||
mfem::Array<int> thetaCenterDofs, phiCenterDofs;
|
||||
mfem::Array<double> thetaCenterVals, phiCenterVals;
|
||||
std::tie(thetaCenterDofs, phi_ess_tdof_list) = findCenterElement();
|
||||
mfem::Array<int> thetaCenterDofs, phiCenterDofs; // phiCenterDofs are not used
|
||||
mfem::Array<double> thetaCenterVals;
|
||||
std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
|
||||
thetaCenterVals.SetSize(thetaCenterDofs.Size());
|
||||
phiCenterVals.SetSize(phi_ess_tdof_list.Size());
|
||||
|
||||
thetaCenterVals = 1.0;
|
||||
phiCenterVals = 0.0;
|
||||
|
||||
mfem::Array<int> ess_brd(m_mesh->bdr_attributes.Max());
|
||||
ess_brd = 1;
|
||||
|
||||
mfem::Array<double> thetaSurfaceVals;
|
||||
mfem::Array<double> thetaSurfaceVals, phiSurfaceVals;
|
||||
m_feTheta->GetEssentialTrueDofs(ess_brd, theta_ess_tdof_list);
|
||||
m_fePhi->GetEssentialTrueDofs(ess_brd, phi_ess_tdof_list);
|
||||
|
||||
thetaSurfaceVals.SetSize(theta_ess_tdof_list.Size());
|
||||
thetaSurfaceVals = 0.0;
|
||||
phiSurfaceVals.SetSize(phi_ess_tdof_list.Size());
|
||||
phiSurfaceVals = polycoeff::thetaSurfaceFlux(m_polytropicIndex);
|
||||
|
||||
// combine the essential dofs with the center dofs
|
||||
theta_ess_tdof_list.Append(thetaCenterDofs);
|
||||
thetaSurfaceVals.Append(thetaCenterVals);
|
||||
|
||||
SSE::MFEMArrayPair thetaPair = std::make_pair(theta_ess_tdof_list, thetaSurfaceVals);
|
||||
SSE::MFEMArrayPair phiPair = std::make_pair(phi_ess_tdof_list, phiCenterVals);
|
||||
SSE::MFEMArrayPair phiPair = std::make_pair(phi_ess_tdof_list, phiSurfaceVals);
|
||||
SSE::MFEMArrayPairSet pairSet = std::make_pair(thetaPair, phiPair);
|
||||
|
||||
return pairSet;
|
||||
@@ -275,16 +266,39 @@ void PolySolver::setInitialGuess() const {
|
||||
const double radius = Probe::getMeshRadius(*m_mesh);
|
||||
const double u = 1/radius;
|
||||
|
||||
return (-1.0/radius) * r + 1;
|
||||
// return (-1.0/radius) * r + 1;
|
||||
// return -std::pow((u*r), 2)+1.0; // The series expansion is a better guess; however, this is cheaper and ensures that the value at the surface is very close to zero in a way that the series expansion does not
|
||||
return laneEmden::thetaSeriesExpansion(r, m_polytropicIndex, 5);
|
||||
}
|
||||
);
|
||||
|
||||
mfem::VectorFunctionCoefficient phiSurfaceVectors (m_mesh->SpaceDimension(),
|
||||
[this](const mfem::Vector &x, mfem::Vector &y) {
|
||||
const double r = x.Norml2();
|
||||
mfem::Vector xh(x);
|
||||
xh /= r; // Normalize the vector
|
||||
y.SetSize(m_mesh->SpaceDimension());
|
||||
y = xh;
|
||||
y *= polycoeff::thetaSurfaceFlux(m_polytropicIndex);
|
||||
}
|
||||
);
|
||||
// We want to apply specific boundary conditions to the surface
|
||||
mfem::Array<int> ess_brd(m_mesh->bdr_attributes.Max());
|
||||
ess_brd = 1;
|
||||
|
||||
// θ = 0 at surface
|
||||
mfem::ConstantCoefficient surfacePotential(0);
|
||||
|
||||
m_theta->ProjectCoefficient(thetaInitGuess);
|
||||
m_theta->ProjectBdrCoefficient(surfacePotential, ess_brd);
|
||||
|
||||
mfem::GradientGridFunctionCoefficient phiInitGuess (m_theta.get());
|
||||
m_phi->ProjectCoefficient(phiInitGuess);
|
||||
m_phi->ProjectBdrCoefficientNormal(phiSurfaceVectors, ess_brd);
|
||||
|
||||
if (m_config.get<bool>("Poly:Solver:ViewInitialGuess", false)) {
|
||||
Probe::glVisView(*m_theta, *m_mesh, "θ init");
|
||||
Probe::glVisView(*m_phi, *m_mesh, "ɸ init");
|
||||
Probe::glVisView(*m_phi, *m_mesh, "φ init");
|
||||
}
|
||||
|
||||
}
|
||||
@@ -335,7 +349,7 @@ void PolySolver::LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol
|
||||
LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel);
|
||||
}
|
||||
|
||||
void PolySolver::GetDofCoordinates(mfem::FiniteElementSpace &fes, const std::string& filename) const {
|
||||
void PolySolver::GetDofCoordinates(const mfem::FiniteElementSpace &fes, const std::string& filename) {
|
||||
mfem::Mesh *mesh = fes.GetMesh();
|
||||
double r = Probe::getMeshRadius(*mesh);
|
||||
std::ofstream outputFile(filename, std::ios::out | std::ios::trunc);
|
||||
|
||||
@@ -111,6 +111,8 @@ private: // Private methods
|
||||
* trying to account for the true size of the system. VSize on the other hand refers to the total number of dofs.
|
||||
*
|
||||
* @return blockOffsets The offsets for the blocks in the operator
|
||||
*
|
||||
* @pre m_feTheta and m_fePhi must be valid, populated FiniteElementSpaces.
|
||||
*/
|
||||
mfem::Array<int> computeBlockOffsets() const;
|
||||
|
||||
|
||||
@@ -97,6 +97,8 @@ namespace polyMFEMUtils {
|
||||
}
|
||||
|
||||
// Calculate the Jacobian
|
||||
|
||||
// TODO: Check for when theta ~ 0?
|
||||
const double u_safe = std::max(u_val, 0.0);
|
||||
const double d_u_nl = m_polytropicIndex * std::pow(u_safe, m_polytropicIndex - 1);
|
||||
const double x2_d_u_nl = d_u_nl;
|
||||
|
||||
@@ -28,6 +28,8 @@
|
||||
|
||||
#include "linalg/tensor.hpp"
|
||||
|
||||
static int s_PolyOperatorMultCalledCount = 0;
|
||||
|
||||
void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfem::DenseMatrix *mat) {
|
||||
if (!mat) {
|
||||
throw std::runtime_error("The operator is not a SparseMatrix.");
|
||||
@@ -78,7 +80,7 @@ void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfe
|
||||
*/
|
||||
void writeOperatorToCSV(const mfem::Operator &op,
|
||||
const std::string &filename,
|
||||
int precision = 6) // Add precision argument
|
||||
const int precision = 6) // Add precision argument
|
||||
{
|
||||
// Attempt to cast the Operator to a SparseMatrix
|
||||
const auto *sparse_mat = dynamic_cast<const mfem::SparseMatrix*>(&op);
|
||||
@@ -138,34 +140,19 @@ PolytropeOperator::PolytropeOperator(
|
||||
|
||||
void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
|
||||
if (m_isFinalized) {
|
||||
return;
|
||||
return; // do nothing if already finalized
|
||||
}
|
||||
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());
|
||||
|
||||
for (const int thetaDof : m_theta_ess_tdofs.first) {
|
||||
m_Mmat->EliminateRow(thetaDof);
|
||||
m_Qmat->EliminateCol(thetaDof);
|
||||
}
|
||||
// These are commented out because they theoretically are wrong (need to think more about how to apply essential dofs to a vector div field)
|
||||
// for (const int phiDof : m_phi_ess_tdofs.first) {
|
||||
// if (phiDof >=0 && phiDof < m_Dmat->Height()) {
|
||||
// m_Dmat->EliminateRowCol(phiDof);
|
||||
// m_Qmat->EliminateRow(phiDof);
|
||||
// m_Mmat->EliminateCol(phiDof);
|
||||
// }
|
||||
// }
|
||||
|
||||
m_negM_op = std::make_unique<mfem::ScaledOperator>(m_Mmat.get(), -1.0);
|
||||
m_negQ_op = std::make_unique<mfem::ScaledOperator>(m_Qmat.get(), -1.0);
|
||||
m_negM_op = std::make_unique<mfem::ScaledOperator>(m_M.get(), -1.0);
|
||||
m_negQ_op = std::make_unique<mfem::ScaledOperator>(m_Q.get(), -1.0);
|
||||
|
||||
// Set up the constant parts of the jacobian now
|
||||
m_jacobian = std::make_unique<mfem::BlockOperator>(m_blockOffsets);
|
||||
m_jacobian->SetBlock(0, 1, m_negM_op.get()); // -M (constant)
|
||||
m_jacobian->SetBlock(1, 0, m_negQ_op.get()); // -Q (constant)
|
||||
m_jacobian->SetBlock(1, 1, m_Dmat.get()); // D (constant)
|
||||
m_jacobian->SetBlock(0, 1, m_negM_op.get()); //<- -M (constant)
|
||||
m_jacobian->SetBlock(1, 0, m_negQ_op.get()); //<- -Q (constant)
|
||||
m_jacobian->SetBlock(1, 1, m_D.get()); //<- D (constant)
|
||||
|
||||
// Build the initial preconditioner based on some initial guess
|
||||
const auto &grad = m_f->GetGradient(initTheta);
|
||||
updatePreconditioner(grad);
|
||||
|
||||
@@ -173,15 +160,22 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
|
||||
}
|
||||
|
||||
const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const {
|
||||
if (m_jacobian == nullptr) {
|
||||
MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call.");
|
||||
}
|
||||
return *m_jacobian;
|
||||
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::GetPreconditioner() const {
|
||||
if (m_schurPreconditioner == nullptr) {
|
||||
MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call.");
|
||||
}
|
||||
if (!m_isFinalized) {
|
||||
MFEM_ABORT("PolytropeOperator::GetPreconditioner called before finalize");
|
||||
MFEM_ABORT("PolytropeOperator not finalized prior to call to GetPreconditioner().");
|
||||
}
|
||||
return *m_schurPreconditioner;
|
||||
}
|
||||
@@ -213,38 +207,39 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
|
||||
// R1 = Dɸ - Qθ
|
||||
|
||||
MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult");
|
||||
MFEM_ASSERT(m_Mmat.get() != nullptr, "SparseMatrix m_Mmat is null in PolytropeOperator::Mult");
|
||||
MFEM_ASSERT(m_Dmat.get() != nullptr, "SparseMatrix m_Dmat is null in PolytropeOperator::Mult");
|
||||
MFEM_ASSERT(m_Qmat.get() != nullptr, "SparseMatrix m_Qmat is null in PolytropeOperator::Mult");
|
||||
|
||||
m_f->Mult(x_theta, f_term);
|
||||
m_Mmat->Mult(x_phi, Mphi_term);
|
||||
m_Dmat->Mult(x_phi, Dphi_term);
|
||||
m_Qmat->Mult(x_theta, Qtheta_term);
|
||||
m_f->Mult(x_theta, f_term); // fixme: this may be the wrong way to assemble m_f?
|
||||
m_M->Mult(x_phi, Mphi_term);
|
||||
m_D->Mult(x_phi, Dphi_term);
|
||||
m_Q->Mult(x_theta, Qtheta_term);
|
||||
|
||||
subtract(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++) {
|
||||
int idx = m_theta_ess_tdofs.first[i];
|
||||
if (idx >= 0 && idx < y_R0.Size()) {
|
||||
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] = targetValue - x_theta(idx); // Zero out the essential theta dofs in the bi-linear form
|
||||
// y_block.GetBlock(0)[idx] = 0; // Zero out the essential theta dofs in the bi-linear form
|
||||
// y_block.GetBlock(0)[idx] = targetValue - x_theta(idx); // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising
|
||||
y_block.GetBlock(0)[idx] = 0; // Zero out the essential theta dofs in the bi-linear form
|
||||
}
|
||||
}
|
||||
|
||||
// TODO look into how the true dof -> vector component works
|
||||
// for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) {
|
||||
// int idx = m_phi_ess_tdofs.first[i];
|
||||
// if (idx >= 0 && idx < y_R1.Size()) {
|
||||
// // const double &targetValue = m_phi_ess_tdofs.second[i];
|
||||
// // y_block.GetBlock(1)[idx] = targetValue - x_phi(idx); // Zero out the essential phi dofs in the bi-linear form
|
||||
// y_block.GetBlock(1)[idx] = 0; // Zero out the essential phi dofs in the bi-linear form
|
||||
// }
|
||||
// }
|
||||
// TODO: Are the residuals for φ being calculated correctly?
|
||||
|
||||
std::ofstream outfileθ("PolyOperatorMultCallTheta_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc);
|
||||
outfileθ << "dof,R_θ\n";
|
||||
for (int i = 0; i < y_R0.Size(); i++) {
|
||||
outfileθ << i << "," << y_R0(i) << "\n";
|
||||
}
|
||||
outfileθ.close();
|
||||
std::ofstream outfileφ("PolyOperatorMultCallPhi_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc);
|
||||
outfileφ << "dof,R_ɸ\n";
|
||||
for (int i = 0; i < y_R1.Size(); i++) {
|
||||
outfileφ << i << "," << y_R1(i) << "\n";
|
||||
}
|
||||
outfileφ.close();
|
||||
++s_PolyOperatorMultCalledCount;
|
||||
|
||||
}
|
||||
|
||||
@@ -253,7 +248,7 @@ void PolytropeOperator::updateInverseNonlinearJacobian(const mfem::Operator &gra
|
||||
if (const auto *sparse_mat = dynamic_cast<const mfem::SparseMatrix*>(&grad); sparse_mat != nullptr) {
|
||||
approxJacobiInvert(*sparse_mat, m_invNonlinearJacobian, "Nonlinear Jacobian");
|
||||
} else {
|
||||
MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear jacobian where nonlinear jacobian is not trivially castable to a sparse matrix");
|
||||
MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear jacobian where nonlinear jacobian is not dynamically castable to a sparse matrix");
|
||||
}
|
||||
}
|
||||
|
||||
@@ -262,11 +257,11 @@ void PolytropeOperator::updateInverseSchurCompliment() const {
|
||||
if (m_invNonlinearJacobian == nullptr) {
|
||||
MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian");
|
||||
}
|
||||
mfem::SparseMatrix* schurCompliment(m_Dmat.get()); // This is now a copy of D
|
||||
mfem::SparseMatrix* schurCompliment(&m_D->SpMat()); // This is now a copy of D
|
||||
|
||||
// Calculate S = D - Q df^{-1} M
|
||||
mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(*m_Qmat, *m_invNonlinearJacobian); // Q * df^{-1}
|
||||
mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, *m_Mmat); // Q * df^{-1} * M
|
||||
mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(m_Q->SpMat(), *m_invNonlinearJacobian); // Q * df^{-1}
|
||||
const mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, m_M->SpMat()); // Q * df^{-1} * M
|
||||
|
||||
// PERF: Could potentially add some caching in here to only update the preconditioner when some condition has been met
|
||||
schurCompliment->Add(-1, *temp_QxdfInvxM); // D - Q * df^{-1} * M
|
||||
@@ -275,6 +270,9 @@ void PolytropeOperator::updateInverseSchurCompliment() const {
|
||||
if (m_schurPreconditioner == nullptr) {
|
||||
m_schurPreconditioner = std::make_unique<mfem::BlockDiagonalPreconditioner>(m_blockOffsets);
|
||||
}
|
||||
|
||||
// ⎡ḟ(θ)^-1 0⎤
|
||||
// ⎣0 S^-1⎦
|
||||
m_schurPreconditioner->SetDiagonalBlock(0, m_invNonlinearJacobian.get());
|
||||
m_schurPreconditioner->SetDiagonalBlock(1, m_invSchurCompliment.get());
|
||||
}
|
||||
@@ -312,6 +310,26 @@ void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess
|
||||
} else {
|
||||
MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs");
|
||||
}
|
||||
|
||||
if (m_M) {
|
||||
m_M->EliminateTestEssentialBC(theta_ess_tdofs.first);
|
||||
m_M->EliminateTrialEssentialBC(phi_ess_tdofs.first);
|
||||
} else {
|
||||
MFEM_ABORT("m_M is null in PolytropeOperator::SetEssentialTrueDofs");
|
||||
}
|
||||
|
||||
if (m_Q) {
|
||||
m_Q->EliminateTrialEssentialBC(theta_ess_tdofs.first);
|
||||
m_Q->EliminateTestEssentialBC(phi_ess_tdofs.first);
|
||||
} else {
|
||||
MFEM_ABORT("m_Q is null in PolytropeOperator::SetEssentialTrueDofs");
|
||||
}
|
||||
|
||||
if (m_D) {
|
||||
m_D->EliminateEssentialBC(phi_ess_tdofs.first);
|
||||
} else {
|
||||
MFEM_ABORT("m_D is null in PolytropeOperator::SetEssentialTrueDofs");
|
||||
}
|
||||
}
|
||||
|
||||
void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) {
|
||||
|
||||
@@ -70,10 +70,6 @@ private:
|
||||
SSE::MFEMArrayPair m_theta_ess_tdofs;
|
||||
SSE::MFEMArrayPair m_phi_ess_tdofs;
|
||||
|
||||
std::unique_ptr<mfem::SparseMatrix> m_Mmat;
|
||||
std::unique_ptr<mfem::SparseMatrix> m_Qmat;
|
||||
std::unique_ptr<mfem::SparseMatrix> m_Dmat;
|
||||
|
||||
std::unique_ptr<mfem::ScaledOperator> m_negM_op;
|
||||
std::unique_ptr<mfem::ScaledOperator> m_negQ_op;
|
||||
mutable std::unique_ptr<mfem::BlockOperator> m_jacobian;
|
||||
|
||||
Reference in New Issue
Block a user