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

feat(NonlinearPowerIntergrator): increased robustness to theta ~ 0 and theta < 0

Browse Source
This commit is contained in:
Emily Boudreaux
2025-05-11 14:38:22 -04:00
parent 454d49c3d3
commit 3d33839028
2 changed files with 61 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

68
src/poly/utils/private/integrators.cpp
View File

@@ -22,14 +22,23 @@
#include <cmath>
#include "integrators.h"
#include "config.h"
#include <string>
// static std::ofstream debugOut("gradient.csv", std::ios::trunc);
namespace polyMFEMUtils {
NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) :
m_polytropicIndex(n) {}
m_polytropicIndex(n),
m_epsilon(Config::getInstance().get<double>("Poly:Solver:Epsilon", 1.0e-8)) {
if (m_polytropicIndex < 0.0) {
throw std::invalid_argument("Polytropic index must be non-negative.");
}
if (m_polytropicIndex > 5.0) {
throw std::invalid_argument("Polytropic index must be less than 5.0.");
}
if (m_epsilon <= 0.0) {
throw std::invalid_argument("Epsilon (Poly:Solver:Epsilon) must be positive non-zero.");
}
}
void NonlinearPowerIntegrator::AssembleElementVector(
const mfem::FiniteElement &el,
@@ -54,13 +63,20 @@ namespace polyMFEMUtils {
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
const double u_safe = std::max(u_val, 0.0);
const double u_nl = std::pow(u_safe, m_polytropicIndex);
const double x2_u_nl = u_nl;
double u_nl;
if (u_val < m_epsilon) {
u_nl = fmod(u_val, m_polytropicIndex, m_epsilon);
} else {
u_nl = std::pow(u_val, m_polytropicIndex);
}
// const double u_safe = std::max(u_val, 0.0);
// const double u_nl = std::pow(u_safe, m_polytropicIndex);
//
// const double x2_u_nl = u_nl;
for (int i = 0; i < dof; i++){
elvect(i) += shape(i) * x2_u_nl * weight;
elvect(i) += shape(i) * u_nl * weight;
}
}
}
@@ -99,37 +115,23 @@ 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;
double d_u_nl;
if (u_val < m_epsilon) {
d_u_nl = dfmod(m_epsilon, m_polytropicIndex);
} else {
d_u_nl = m_polytropicIndex * std::pow(u_val, m_polytropicIndex - 1.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;
for (int i = 0; i < dof; i++) {
for (int j = 0; j < dof; j++) {
elmat(i, j) += shape(i) * x2_d_u_nl * shape(j) * weight;
elmat(i, j) += shape(i) * d_u_nl * shape(j) * weight;
}
}
// // --- Debug Code to write out gradient ---
// Trans.Transform(ip,physCoord);
// el.CalcDShape(ip, dshape);
//
// mfem::CalcInverse(Trans.Jacobian(), invJ);
//
// mfem::DenseMatrix dshapePhys;
// dshapePhys.SetSize(dof, physCoord.Size());
// mfem::Mult(dshape, invJ, dshapePhys);
//
// gradPhys = 0.0;
// for (int j = 0; j < dof; j++) {
// for (int d = 0; d < gradPhys.Size(); d++) {
// gradPhys(d) += elfun(j)*dshapePhys(j, d);
// }
// }
//
// debugOut
// << physCoord(0) << ", " << physCoord(1) << ", " << physCoord(2)
// << ", " << gradPhys(0) << ", " << gradPhys(1) << ", " << gradPhys(2) << '\n';
}
// debugOut.flush();
}
} // namespace polyMFEMUtils

25
src/poly/utils/public/integrators.h
View File

@@ -73,5 +73,30 @@ namespace polyMFEMUtils {
Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
double m_polytropicIndex;
double m_epsilon;
};
inline double dfmod(const double epsilon, const double n) {
if (n == 0.0) {
return 0.0;
}
if (n == 1.0) {
return 1.0;
}
return n * std::pow(epsilon, n - 1.0);
}
inline double fmod(const double theta, const double n, const double epsilon) {
if (n == 0.0) {
return 1.0;
}
// For n != 0
const double y_prime_at_epsilon = dfmod(epsilon, n); // Uses the robust dfmod
const double y_at_epsilon = std::pow(epsilon, n); // epsilon^n
// f_mod(theta) = y_at_epsilon + y_prime_at_epsilon * (theta - epsilon)
return y_at_epsilon + y_prime_at_epsilon * (theta - epsilon);
}
} // namespace polyMFEMUtils