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

Add certified polynomial enthalpy bounds and feasible direction repair

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 09:23:19 -04:00
parent 7a07caeee3
commit 2c084b1c62
3 changed files with 713 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
experiments/POLYNOMIAL_ENTHALPY_ADMISSIBILITY_2026-09-17.md Normal file
View File

@@ -0,0 +1,17 @@
# Polynomial enthalpy admissibility and feasible direction repair
This branch adds a sufficient polynomial nonnegativity certificate and an interchangeable enthalpy direction repair. It does not alter the EOS or declare a constrained residual minimum an equilibrium.
`TensorBernsteinPolynomial<Dimension>` encodes spatial dimension in its type and accepts polynomial order through the mesh's runtime nodal points. A cached one-dimensional Lagrange-to-Bernstein transformation is applied by tensor contractions. Every nontrivial arithmetic operation in the transformation and de Casteljau subdivision encloses its result using `nextafter`; a polynomial is certified only when all coefficient lower bounds in every leaf are nonnegative. The implementation assumes ordinary IEEE binary64 arithmetic, gradual underflow, and no unsafe fast-math transformations. The underlying mathematical FE is the Lagrange polynomial at the supplied floating-point nodal coordinates.
A negative control coefficient is not a negative polynomial witness. Such cells are refined. A negative upper interval bound at a tensor corner proves a negative polynomial value. Cells that remain ambiguous at the work budget return `unresolved`, including some truly nonnegative polynomials with tangential interior zeros. No negative tolerance is accepted as nonnegative. Exact endpoint interpolation preserves identically zero tensor faces without manufacturing a small negative error interval.
The production adapter supports conforming scalar nodal H1 hexahedra with closed Gauss-Lobatto, closed uniform, or closed Gauss-Legendre points. Unsupported element types are rejected explicitly. Its reduced field map and MFEM prolongation preserve shared true DOFs; no elementwise coefficient clipping occurs. The prepared object assumes the finite-element space and topology remain unchanged. Mesh adaptation requires a new object.
Given a certified base `h`, a certified anchor `a`, and proposed direction `p`, the repaired full candidate is `(1-beta)*(h+p)+beta*a`. Thus the altered enthalpy direction is `(1-beta)*p+beta*(a-h)`. Zero-trace reduced rows are fixed to zero. The adapter constructs a sufficient convex cone from a subdivision on which the anchor has nonnegative Bernstein bounds, and computes the largest required fraction across all coefficients and MPI ranks. This is the minimum sufficient fraction for those conservative coefficient inequalities, not the least fraction for the complete nonnegative polynomial set. A caller can choose any larger fraction through one, then verify the actual rounded full candidate. The original direction is left unchanged when verification fails.
This changes the enthalpy direction while allowing the stellar strategy to retain the other field directions. It is a feasible direction repair, not a constrained Newton optimizer. Full residual acceptance, actual descent checks, and the original physical convergence requirements remain necessary. Invalid or unresolved bases are rejected explicitly; this mechanism is not a feasibility-restoration solver. A certified checkpoint or initial seed is required as anchor. A frozen repaired full candidate and certified base define a nonnegative exact-arithmetic segment for backtracking; every floating-point candidate must nevertheless pass the same check.
Collective adapter calls must occur in the same order on all ranks. Sign status, polynomial bounds, and repair fraction are globally combined. Witness element indices, coordinates, and work counters remain rank-local. Unequal repair fractions and work budgets are rejected collectively.
Tests cover nodal-positive/interior-negative cubics, a positive polynomial with negative coarse Bernstein coefficients, exact zero faces, very small negative pockets, unresolved tangential zeros, arbitrary runtime degree, nonfinite input, finite-element ordering, true-DOF sharing, invalid bases, immutable failure output, and prescribed trace handling. Integration tests are intended for both one and two MPI ranks.

508
libmeanfield/interface/solver/enthalpy_admissibility.cppm Normal file
View File

@@ -0,0 +1,508 @@
module;
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <limits>
#include <map>
#include <memory>
#include <span>
#include <stdexcept>
#include <utility>
#include <vector>
#include <mfem.hpp>
#include <mpi.h>
export module mean_field:solver.enthalpy_admissibility;
export import :field.mfem;
export namespace mean_field::solver {
enum class PolynomialSignStatus { certified_nonnegative, unresolved, negative, non_finite };
struct EnthalpyAdmissibilityOptions final {
int maximumSubdivisionDepth{24};
int maximumSubcells{8192};
};
struct PolynomialSignReport final {
PolynomialSignStatus status{PolynomialSignStatus::certified_nonnegative};
double lowerBound{std::numeric_limits<double>::infinity()};
double upperBound{-std::numeric_limits<double>::infinity()};
int visitedSubcells{0};
int maximumDepth{0};
int witnessElement{-1};
std::array<double, 3> witnessReferencePoint{};
};
enum class EnthalpyDirectionRepairStatus {
feasible, invalid_base, invalid_anchor, invalid_direction, unresolved
};
struct EnthalpyDirectionRepairReport final {
EnthalpyDirectionRepairStatus status{EnthalpyDirectionRepairStatus::unresolved};
double minimumSafeFraction{1.0};
PolynomialSignReport base;
PolynomialSignReport anchor;
PolynomialSignReport candidate;
};
}
namespace mean_field::solver::enthalpy_detail {
// Bounds concern the represented polynomial, not a tolerance-relaxed EOS.
// Exact zero is retained algebraically so a prescribed zero trace does not
// become an unresolved negative interval merely through outward rounding.
struct Interval final {
double lo{0.0}, hi{0.0};
explicit Interval(double value = 0.0) : lo(value), hi(value) {}
Interval(double lower, double upper) : lo(lower), hi(upper) {}
bool zero() const { return lo == 0.0 && hi == 0.0; }
bool finite() const { return std::isfinite(lo) && std::isfinite(hi); }
};
inline double Down(double value) { return std::nextafter(value, -std::numeric_limits<double>::infinity()); }
inline double Up(double value) { return std::nextafter(value, std::numeric_limits<double>::infinity()); }
inline Interval Add(Interval a, Interval b) {
if (a.zero()) return b;
if (b.zero()) return a;
return {Down(a.lo + b.lo), Up(a.hi + b.hi)};
}
inline Interval Negate(Interval a) { return {-a.hi, -a.lo}; }
inline Interval Multiply(Interval a, Interval b) {
if (a.zero() || b.zero()) return Interval{};
if (a.lo == 1.0 && a.hi == 1.0) return b;
if (b.lo == 1.0 && b.hi == 1.0) return a;
const std::array products{a.lo*b.lo, a.lo*b.hi, a.hi*b.lo, a.hi*b.hi};
return {Down(*std::min_element(products.begin(), products.end())),
Up(*std::max_element(products.begin(), products.end()))};
}
inline Interval Divide(Interval a, Interval b) {
if (b.lo <= 0.0 && b.hi >= 0.0) throw std::invalid_argument("A tensor nodal basis has repeated nodes.");
if (a.zero()) return Interval{};
return Multiply(a, {Down(1.0/b.hi), Up(1.0/b.lo)});
}
inline void ValidateOptions(const EnthalpyAdmissibilityOptions &options) {
if (options.maximumSubdivisionDepth < 0 || options.maximumSubdivisionDepth > 60 ||
options.maximumSubcells < 1)
throw std::invalid_argument("Polynomial certification requires a finite positive work budget.");
}
inline int TensorSize(int points, int dimension) {
std::size_t count = 1;
for (int d = 0; d < dimension; ++d) {
count *= static_cast<std::size_t>(points);
if (count > static_cast<std::size_t>(std::numeric_limits<int>::max()))
throw std::invalid_argument("Tensor polynomial is too large.");
}
return static_cast<int>(count);
}
inline int CollectiveMaximum(int value, MPI_Comm comm) {
int result = 0;
if (MPI_Allreduce(&value, &result, 1, MPI_INT, MPI_MAX, comm) != MPI_SUCCESS)
throw std::runtime_error("Enthalpy admissibility could not combine rank-local status.");
return result;
}
}
export namespace mean_field::solver {
/*
* Dimension belongs to the type; polynomial order belongs to the mesh.
* Coefficients use lexicographic ordering with the first coordinate fastest.
* Only tensor Lagrange polynomials on [0,1]^Dimension are represented.
*/
template <int Dimension> requires (Dimension >= 1 && Dimension <= 3)
class TensorBernsteinPolynomial final {
using Interval = enthalpy_detail::Interval;
public:
explicit TensorBernsteinPolynomial(std::span<const double> nodes)
: m_nodes(nodes.begin(), nodes.end()), m_points(static_cast<int>(nodes.size())),
m_size(enthalpy_detail::TensorSize(m_points, Dimension)) {
if (m_points < 2 || m_nodes.front() != 0.0 || m_nodes.back() != 1.0)
throw std::invalid_argument("Tensor certification requires closed nodal points including 0 and 1.");
for (int i = 0; i < m_points; ++i)
if (!std::isfinite(m_nodes[i]) || (i && m_nodes[i] <= m_nodes[i-1]))
throw std::invalid_argument("Tensor nodes must be finite and strictly increasing.");
BuildTransform();
}
[[nodiscard]] int Size() const noexcept { return m_size; }
[[nodiscard]] int Order() const noexcept { return m_points-1; }
[[nodiscard]] PolynomialSignReport Check(std::span<const double> nodalValues,
EnthalpyAdmissibilityOptions options = {}) const {
enthalpy_detail::ValidateOptions(options);
auto coefficients = Transform(nodalValues);
struct Cell { std::vector<Interval> coefficients; int depth; std::array<double,3> lo, hi; };
std::vector<Cell> pending;
pending.push_back({std::move(coefficients), 0, {0,0,0}, {1,1,1}});
PolynomialSignReport report;
bool unresolved = false;
while (!pending.empty()) {
Cell cell = std::move(pending.back());
pending.pop_back();
++report.visitedSubcells;
report.maximumDepth = std::max(report.maximumDepth, cell.depth);
double lower = std::numeric_limits<double>::infinity(), upper = -lower;
bool finite = true;
for (const auto &coefficient : cell.coefficients) {
finite = finite && coefficient.finite();
lower = std::min(lower, coefficient.lo);
upper = std::max(upper, coefficient.hi);
}
if (!finite) { report.status = PolynomialSignStatus::non_finite; return report; }
// Bernstein coefficients at tensor corners are actual values;
// an arbitrary negative control coefficient is not a witness.
for (int corner = 0; corner < (1<<Dimension); ++corner) {
int index = 0, stride = 1;
for (int axis = 0; axis < Dimension; ++axis) {
const bool high = (corner & (1<<axis)) != 0;
index += high ? (m_points-1)*stride : 0;
stride *= m_points;
}
if (cell.coefficients[index].hi < 0.0) {
report.status = PolynomialSignStatus::negative;
report.lowerBound = lower;
report.upperBound = upper;
for (int axis = 0; axis < Dimension; ++axis)
report.witnessReferencePoint[axis] = (corner & (1<<axis)) ? cell.hi[axis] : cell.lo[axis];
return report;
}
}
if (lower >= 0.0) {
report.lowerBound = std::min(report.lowerBound, lower);
report.upperBound = std::max(report.upperBound, upper);
continue;
}
if (cell.depth >= options.maximumSubdivisionDepth ||
report.visitedSubcells + static_cast<int>(pending.size()) + 2 > options.maximumSubcells) {
unresolved = true;
report.lowerBound = std::min(report.lowerBound, lower);
report.upperBound = std::max(report.upperBound, upper);
continue;
}
const int axis = cell.depth % Dimension;
auto children = Split(cell.coefficients, axis);
const double midpoint = 0.5*(cell.lo[axis]+cell.hi[axis]);
auto leftHi = cell.hi, rightLo = cell.lo;
leftHi[axis] = midpoint;
rightLo[axis] = midpoint;
pending.push_back({std::move(children.second), cell.depth+1, rightLo, cell.hi});
pending.push_back({std::move(children.first), cell.depth+1, cell.lo, leftHi});
}
report.status = unresolved ? PolynomialSignStatus::unresolved : PolynomialSignStatus::certified_nonnegative;
return report;
}
// This computes a sufficient fraction in a convex inner approximation,
// not the minimum over all nonnegative polynomials or a Newton optimum.
[[nodiscard]] double MinimumRepairFraction(std::span<const double> anchor,
std::span<const double> candidate,
EnthalpyAdmissibilityOptions options = {}) const {
enthalpy_detail::ValidateOptions(options);
struct Cell { std::vector<Interval> anchor, candidate; int depth; };
std::vector<Cell> pending;
pending.push_back({Transform(anchor), Transform(candidate), 0});
double fraction = 0.0;
int visited = 0;
while (!pending.empty()) {
Cell cell = std::move(pending.back()); pending.pop_back(); ++visited;
bool anchorCertified = true;
for (const auto &a : cell.anchor) {
if (!a.finite()) return std::numeric_limits<double>::quiet_NaN();
anchorCertified = anchorCertified && a.lo >= 0.0;
}
if (!anchorCertified) {
if (cell.depth >= options.maximumSubdivisionDepth ||
visited + static_cast<int>(pending.size()) + 2 > options.maximumSubcells)
return std::numeric_limits<double>::quiet_NaN();
const int axis = cell.depth % Dimension;
auto a = Split(cell.anchor, axis), c = Split(cell.candidate, axis);
pending.push_back({std::move(a.second), std::move(c.second), cell.depth+1});
pending.push_back({std::move(a.first), std::move(c.first), cell.depth+1});
continue;
}
for (int i = 0; i < m_size; ++i) {
const double c = cell.candidate[i].lo, a = cell.anchor[i].lo;
if (!cell.candidate[i].finite()) return std::numeric_limits<double>::quiet_NaN();
if (c < 0.0) {
// Outward rounding covers both subtraction and division.
const double denominator = enthalpy_detail::Down(a-c);
const double needed = a == 0.0 ? 1.0 : enthalpy_detail::Up((-c)/denominator);
fraction = std::max(fraction, std::min(1.0, needed));
}
}
}
return fraction;
}
private:
void BuildTransform() {
using namespace enthalpy_detail;
m_transform.resize(static_cast<std::size_t>(m_points*m_points));
const int degree = m_points-1;
for (int node = 0; node < m_points; ++node) {
std::vector<Interval> power(1, Interval(1.0));
for (int other = 0; other < m_points; ++other) {
if (other == node) continue;
const Interval denominator = Add(Interval(m_nodes[node]), Negate(Interval(m_nodes[other])));
std::vector<Interval> next(power.size()+1);
for (std::size_t k = 0; k < power.size(); ++k) {
next[k] = Add(next[k], Divide(Multiply(power[k], Interval(-m_nodes[other])), denominator));
next[k+1] = Add(next[k+1], Divide(power[k], denominator));
}
power = std::move(next);
}
for (int coefficient = 0; coefficient < m_points; ++coefficient) {
Interval value;
Interval ratio(1.0);
for (int k = 0; k <= coefficient; ++k) {
value = Add(value, Multiply(power[k], ratio));
if (k < coefficient)
ratio = Multiply(ratio, Divide(Interval(coefficient-k), Interval(degree-k)));
}
m_transform[coefficient*m_points+node] = value;
}
}
// Endpoint interpolation is exact and is also the algebraic reason
// a face whose nodal values vanish has an exactly zero trace.
for (int node = 0; node < m_points; ++node) {
m_transform[node] = Interval(node == 0 ? 1.0 : 0.0);
m_transform[(m_points-1)*m_points+node] = Interval(node == m_points-1 ? 1.0 : 0.0);
}
}
[[nodiscard]] std::vector<Interval> Transform(std::span<const double> values) const {
using namespace enthalpy_detail;
if (values.size() != static_cast<std::size_t>(m_size))
throw std::invalid_argument("Tensor polynomial has incompatible nodal coefficients.");
std::vector<Interval> current; current.reserve(values.size());
for (double value : values) current.emplace_back(value);
int stride = 1;
for (int axis = 0; axis < Dimension; ++axis) {
std::vector<Interval> next(m_size);
for (int index = 0; index < m_size; ++index) {
const int row = (index/stride)%m_points, line = index-row*stride;
Interval sum;
for (int column = 0; column < m_points; ++column)
sum = Add(sum, Multiply(m_transform[row*m_points+column], current[line+column*stride]));
next[index] = sum;
}
current = std::move(next); stride *= m_points;
}
return current;
}
[[nodiscard]] std::pair<std::vector<Interval>,std::vector<Interval>>
Split(const std::vector<Interval> &values, int axis) const {
using namespace enthalpy_detail;
std::vector<Interval> left(m_size), right(m_size), row(m_points);
int stride = 1;
for (int d = 0; d < axis; ++d) stride *= m_points;
for (int base = 0; base < m_size; ++base) {
if ((base/stride)%m_points != 0) continue;
for (int j = 0; j < m_points; ++j) row[j] = values[base+j*stride];
left[base] = row[0]; right[base+(m_points-1)*stride] = row[m_points-1];
for (int level = 1; level < m_points; ++level) {
for (int j = 0; j < m_points-level; ++j)
row[j] = Multiply(Add(row[j], row[j+1]), Interval(0.5));
left[base+level*stride] = row[0];
right[base+(m_points-1-level)*stride] = row[m_points-1-level];
}
}
return {std::move(left),std::move(right)};
}
std::vector<double> m_nodes;
int m_points, m_size;
std::vector<Interval> m_transform;
};
/* Prepared against an immutable conforming scalar H1 space. All mutations
* remain in reduced true-DOF coordinates; elementwise coefficient clipping
* would break continuity and is deliberately absent. Collective methods
* must be called in the same order on all ranks of the supplied space. */
class PreparedEnthalpyAdmissibility final {
struct ElementData final {
int element;
const TensorBernsteinPolynomial<3> *basis;
std::vector<int> lexicographicToNative;
};
public:
PreparedEnthalpyAdmissibility(mfem::ParFiniteElementSpace &space, field::FieldDofMap dofMap,
std::span<const int> supportedElementIds,
std::span<const int> zeroTraceReducedDofs = {},
EnthalpyAdmissibilityOptions options = {})
: m_space(&space), m_adapter(std::move(dofMap), space),
m_zeroTrace(zeroTraceReducedDofs.begin(),zeroTraceReducedDofs.end()), m_options(options) {
const bool validOptions = options.maximumSubdivisionDepth >= 0 &&
options.maximumSubdivisionDepth <= 60 && options.maximumSubcells >= 1;
if (enthalpy_detail::CollectiveMaximum(validOptions ? 0 : 1,space.GetComm()))
throw std::invalid_argument("Polynomial certification requires a finite positive work budget.");
bool eligible = space.GetVDim() == 1 && !space.Nonconforming() &&
dynamic_cast<const mfem::H1_FECollection *>(space.FEColl()) != nullptr;
for (int e : supportedElementIds) {
if (e < 0 || e >= space.GetNE()) { eligible = false; continue; }
const auto *fe = dynamic_cast<const mfem::NodalTensorFiniteElement *>(space.GetFE(e));
eligible = eligible && fe && fe->GetDim() == 3 && fe->GetGeomType() == mfem::Geometry::CUBE &&
fe->GetMapType() == mfem::FiniteElement::VALUE &&
(fe->GetBasisType() == mfem::BasisType::GaussLobatto ||
fe->GetBasisType() == mfem::BasisType::ClosedUniform ||
fe->GetBasisType() == mfem::BasisType::ClosedGL) && fe->GetDof() ==
enthalpy_detail::TensorSize(fe->GetOrder()+1,3);
}
for (int row : m_zeroTrace) eligible = eligible && row >= 0 && row < m_adapter.dof_map().reduced_size();
if (enthalpy_detail::CollectiveMaximum(eligible ? 0 : 1, space.GetComm()))
throw std::invalid_argument("Enthalpy certification requires conforming scalar tensor H1 hexahedra.");
const std::array localOptions{options.maximumSubdivisionDepth,options.maximumSubcells};
std::array<int,2> minimum{}, maximum{};
if (MPI_Allreduce(localOptions.data(),minimum.data(),2,MPI_INT,MPI_MIN,space.GetComm()) != MPI_SUCCESS ||
MPI_Allreduce(localOptions.data(),maximum.data(),2,MPI_INT,MPI_MAX,space.GetComm()) != MPI_SUCCESS)
throw std::runtime_error("Could not validate collective polynomial-certification options.");
if (minimum != maximum) throw std::invalid_argument("Polynomial-certification budgets must agree on all ranks.");
for (int e : supportedElementIds) {
const auto &fe = dynamic_cast<const mfem::NodalTensorFiniteElement &>(*space.GetFE(e));
const auto key = std::pair{fe.GetOrder(),fe.GetBasisType()};
const auto &map = fe.GetDofMap();
std::vector<int> ordering(fe.GetDof());
for (int i = 0; i < fe.GetDof(); ++i) ordering[i] = map.Size() ? map[i] : i;
if (!m_bases.contains(key)) {
std::vector<double> nodes(fe.GetOrder()+1);
for (int i = 0; i <= fe.GetOrder(); ++i) nodes[i] = fe.GetNodes().IntPoint(ordering[i]).x;
m_bases.emplace(key,std::make_unique<TensorBernsteinPolynomial<3>>(nodes));
}
m_elements.push_back({e,m_bases.at(key).get(),std::move(ordering)});
}
}
[[nodiscard]] int Size() const noexcept { return m_adapter.dof_map().reduced_size(); }
[[nodiscard]] PolynomialSignReport Check(const mfem::Vector &reducedEnthalpy) const {
RequireSize(reducedEnthalpy);
PolynomialSignReport local;
bool finite = true;
for (int i = 0; i < reducedEnthalpy.Size(); ++i) finite = finite && std::isfinite(reducedEnthalpy(i));
if (enthalpy_detail::CollectiveMaximum(finite ? 0 : 1,m_space->GetComm())) {
local.status = PolynomialSignStatus::non_finite; return local;
}
mfem::ParGridFunction grid(m_space); m_adapter.scatter(reducedEnthalpy,grid);
for (const auto &element : m_elements) {
const auto values = ElementValues(grid,element);
auto report = element.basis->Check(values,m_options);
local.visitedSubcells += report.visitedSubcells;
local.maximumDepth = std::max(local.maximumDepth,report.maximumDepth);
local.lowerBound = std::min(local.lowerBound,report.lowerBound);
local.upperBound = std::max(local.upperBound,report.upperBound);
if (static_cast<int>(report.status) > static_cast<int>(local.status)) {
local.status = report.status;
local.witnessElement = element.element;
local.witnessReferencePoint = report.witnessReferencePoint;
}
}
local.status = static_cast<PolynomialSignStatus>(enthalpy_detail::CollectiveMaximum(
static_cast<int>(local.status),m_space->GetComm()));
double lower = local.lowerBound, upper = local.upperBound;
if (MPI_Allreduce(&lower,&local.lowerBound,1,MPI_DOUBLE,MPI_MIN,m_space->GetComm()) != MPI_SUCCESS ||
MPI_Allreduce(&upper,&local.upperBound,1,MPI_DOUBLE,MPI_MAX,m_space->GetComm()) != MPI_SUCCESS)
throw std::runtime_error("Could not reduce distributed polynomial bounds.");
// Witness indices and work counts intentionally remain rank-local.
return local;
}
[[nodiscard]] EnthalpyDirectionRepairReport MinimumRepairFraction(
const mfem::Vector &base, const mfem::Vector &anchor, const mfem::Vector &direction) const {
RequireSize(base); RequireSize(anchor); RequireSize(direction);
EnthalpyDirectionRepairReport report;
report.base = Check(base);
if (report.base.status != PolynomialSignStatus::certified_nonnegative || !ZeroTrace(base)) {
report.status = EnthalpyDirectionRepairStatus::invalid_base; return report;
}
report.anchor = Check(anchor);
if (report.anchor.status != PolynomialSignStatus::certified_nonnegative || !ZeroTrace(anchor)) {
report.status = EnthalpyDirectionRepairStatus::invalid_anchor; return report;
}
mfem::Vector candidate(base); candidate += direction;
for (int row : m_zeroTrace) candidate(row) = 0.0;
report.candidate = Check(candidate);
if (report.candidate.status == PolynomialSignStatus::non_finite) {
report.status = EnthalpyDirectionRepairStatus::invalid_direction; return report;
}
if (report.candidate.status == PolynomialSignStatus::certified_nonnegative) {
report.minimumSafeFraction = 0.0;
report.status = EnthalpyDirectionRepairStatus::feasible; return report;
}
mfem::ParGridFunction anchorGrid(m_space), candidateGrid(m_space);
m_adapter.scatter(anchor,anchorGrid); m_adapter.scatter(candidate,candidateGrid);
double localFraction = 0.0;
bool resolved = true;
for (const auto &element : m_elements) {
const double fraction = element.basis->MinimumRepairFraction(
ElementValues(anchorGrid,element),ElementValues(candidateGrid,element),m_options);
resolved = resolved && std::isfinite(fraction);
if (std::isfinite(fraction)) localFraction = std::max(localFraction,fraction);
}
if (enthalpy_detail::CollectiveMaximum(resolved ? 0 : 1,m_space->GetComm())) {
report.status = EnthalpyDirectionRepairStatus::unresolved; return report;
}
if (MPI_Allreduce(&localFraction,&report.minimumSafeFraction,1,MPI_DOUBLE,MPI_MAX,m_space->GetComm()) != MPI_SUCCESS)
throw std::runtime_error("Could not combine feasible enthalpy repair fractions.");
// A small inward move covers rounding in the later state addition;
// ApplyRepairFraction still checks the actual represented candidate.
if (report.minimumSafeFraction > 0.0 && report.minimumSafeFraction < 1.0)
report.minimumSafeFraction = std::min(1.0, report.minimumSafeFraction +
64.0*std::numeric_limits<double>::epsilon()*(1.0-report.minimumSafeFraction));
report.status = EnthalpyDirectionRepairStatus::feasible;
return report;
}
[[nodiscard]] PolynomialSignReport ApplyRepairFraction(
const mfem::Vector &base, const mfem::Vector &anchor, const mfem::Vector &originalDirection,
double fraction, mfem::Vector &repairedDirection) const {
RequireSize(base); RequireSize(anchor); RequireSize(originalDirection);
const bool validFraction = std::isfinite(fraction) && fraction >= 0.0 && fraction <= 1.0;
if (enthalpy_detail::CollectiveMaximum(validFraction ? 0 : 1,m_space->GetComm()))
throw std::invalid_argument("An enthalpy repair fraction must lie in [0,1].");
double minimum = 0.0, maximum = 0.0;
if (MPI_Allreduce(&fraction,&minimum,1,MPI_DOUBLE,MPI_MIN,m_space->GetComm()) != MPI_SUCCESS ||
MPI_Allreduce(&fraction,&maximum,1,MPI_DOUBLE,MPI_MAX,m_space->GetComm()) != MPI_SUCCESS)
throw std::runtime_error("Could not validate the collective enthalpy repair fraction.");
if (minimum != maximum) throw std::invalid_argument("The enthalpy repair fraction must agree across ranks.");
mfem::Vector repaired(Size()), candidate(Size());
for (int i = 0; i < Size(); ++i) {
repaired(i) = fraction == 0.0 ? originalDirection(i) :
fraction == 1.0 ? anchor(i)-base(i) :
(1.0-fraction)*originalDirection(i)+fraction*(anchor(i)-base(i));
}
for (int row : m_zeroTrace) repaired(row) = 0.0;
candidate = base; candidate += repaired;
auto report = Check(candidate);
if (report.status == PolynomialSignStatus::certified_nonnegative) repairedDirection = repaired;
return report;
}
private:
void RequireSize(const mfem::Vector &vector) const {
if (enthalpy_detail::CollectiveMaximum(vector.Size() == Size() ? 0 : 1,m_space->GetComm()))
throw std::invalid_argument("Enthalpy feasibility received incompatible reduced true-DOF vectors.");
}
[[nodiscard]] bool ZeroTrace(const mfem::Vector &values) const {
bool zero = true;
for (int row : m_zeroTrace) zero = zero && values(row) == 0.0;
return enthalpy_detail::CollectiveMaximum(zero ? 0 : 1,m_space->GetComm()) == 0;
}
[[nodiscard]] std::vector<double> ElementValues(const mfem::ParGridFunction &grid,
const ElementData &element) const {
mfem::Array<int> dofs; mfem::Vector native;
auto *transformation = m_space->GetElementDofs(element.element,dofs);
grid.GetSubVector(dofs,native);
if (transformation) transformation->InvTransformPrimal(native);
std::vector<double> values(element.lexicographicToNative.size());
for (std::size_t i = 0; i < values.size(); ++i) values[i] = native(element.lexicographicToNative[i]);
return values;
}
mfem::ParFiniteElementSpace *m_space;
field::FieldDofGridFunctionAdapter m_adapter;
std::vector<int> m_zeroTrace;
EnthalpyAdmissibilityOptions m_options;
std::map<std::pair<int,int>,std::unique_ptr<TensorBernsteinPolynomial<3>>> m_bases;
std::vector<ElementData> m_elements;
};
}

188
tests/solver/enthalpy_admissibility.cpp Normal file
View File

@@ -0,0 +1,188 @@
#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include <numeric>
#include <span>
#include <vector>
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <mfem.hpp>
#include <mpi.h>
import mean_field;
namespace {
using mean_field::solver::PolynomialSignStatus;
constexpr std::array cubicNodes{0.0, 1.0/3.0, 2.0/3.0, 1.0};
constexpr std::array quadraticNodes{0.0, 0.5, 1.0};
template <class Function>
std::vector<double> Values(std::span<const double> nodes, int dimension, Function function) {
int count = 1;
for (int d = 0; d < dimension; ++d) count *= static_cast<int>(nodes.size());
std::vector<double> values(count);
for (int i = 0; i < count; ++i) {
int index = i;
std::array<double,3> point{};
for (int d = 0; d < dimension; ++d) {
point[d] = nodes[index%nodes.size()]; index /= static_cast<int>(nodes.size());
}
values[i] = function(point);
}
return values;
}
struct ParallelFixture {
mfem::Mesh serial{mfem::Mesh::MakeCartesian3D(2,2,1,mfem::Element::HEXAHEDRON,1.0,1.0,1.0)};
mfem::ParMesh mesh{MPI_COMM_WORLD,serial};
mfem::H1_FECollection collection{3,3};
mfem::ParFiniteElementSpace space{&mesh,&collection};
std::vector<int> elements;
mean_field::field::FieldDofMap dofMap;
ParallelFixture() {
elements.resize(space.GetNE()); std::iota(elements.begin(),elements.end(),0);
mfem::Array<int> trueDofs(space.GetTrueVSize());
for (int i = 0; i < trueDofs.Size(); ++i) trueDofs[i] = i;
dofMap = mean_field::field::FieldDofMap(space.GetTrueVSize(),trueDofs);
}
template <class Function> mfem::Vector Project(Function function) {
mfem::FunctionCoefficient coefficient(function);
mfem::ParGridFunction grid(&space); grid.ProjectCoefficient(coefficient);
mfem::Vector values; grid.GetTrueDofs(values); return values;
}
};
}
TEST_CASE("Bernstein bounds detect a nodal-positive cubic with an interior negative pocket",
"[solver][enthalpy_admissibility][unit]") {
mean_field::solver::TensorBernsteinPolynomial<1> polynomial(cubicNodes);
const auto values = Values(cubicNodes,1,[](auto p) {
return ((p[0]-0.17)*(p[0]-0.17)-0.002)*(1.0+0.2*p[0]);
});
REQUIRE(std::all_of(values.begin(),values.end(),[](double value) { return value > 0.0; }));
const auto report = polynomial.Check(values);
REQUIRE(report.status == PolynomialSignStatus::negative);
const double x = report.witnessReferencePoint[0];
CHECK(((x-0.17)*(x-0.17)-0.002)*(1.0+0.2*x) < 0.0);
}
TEST_CASE("Negative coarse Bernstein coefficients cause refinement rather than a false rejection",
"[solver][enthalpy_admissibility][unit]") {
mean_field::solver::TensorBernsteinPolynomial<1> polynomial(quadraticNodes);
const auto values = Values(quadraticNodes,1,[](auto p) { return (p[0]-0.5)*(p[0]-0.5)+0.01; });
CHECK(polynomial.Check(values,{.maximumSubdivisionDepth=0}).status == PolynomialSignStatus::unresolved);
const auto refined = polynomial.Check(values);
CHECK(refined.status == PolynomialSignStatus::certified_nonnegative);
CHECK(refined.visitedSubcells > 1);
CHECK(refined.lowerBound >= 0.0);
}
TEST_CASE("Tensor certification retains structural zero faces without an admissibility epsilon",
"[solver][enthalpy_admissibility][unit]") {
mean_field::solver::TensorBernsteinPolynomial<3> polynomial(cubicNodes);
const auto values = Values(cubicNodes,3,[](auto p) {
return p[0]*(1.0-p[0])*p[1]*(1.0-p[1])*p[2]*(1.0-p[2]);
});
const auto report = polynomial.Check(values);
CHECK(report.status == PolynomialSignStatus::certified_nonnegative);
CHECK(report.lowerBound == 0.0);
const std::vector<double> zero(polynomial.Size(),0.0);
CHECK(polynomial.Check(zero).status == PolynomialSignStatus::certified_nonnegative);
}
TEST_CASE("A small negative pocket is not hidden by relative or absolute positivity tolerances",
"[solver][enthalpy_admissibility][unit]") {
mean_field::solver::TensorBernsteinPolynomial<1> polynomial(quadraticNodes);
const auto values = Values(quadraticNodes,1,[](auto p) { return (p[0]-0.37)*(p[0]-0.37)-1.0e-10; });
CHECK(polynomial.Check(values).status == PolynomialSignStatus::negative);
auto invalid = values; invalid[1] = std::numeric_limits<double>::quiet_NaN();
CHECK(polynomial.Check(invalid).status == PolynomialSignStatus::non_finite);
invalid[1] = std::numeric_limits<double>::infinity();
CHECK(polynomial.Check(invalid).status == PolynomialSignStatus::non_finite);
CHECK_THROWS_AS(polynomial.Check(std::span<const double>{}),std::invalid_argument);
}
TEST_CASE("Certification supports runtime polynomial degree and honest unresolved results at tangential zeros",
"[solver][enthalpy_admissibility][unit]") {
const std::array nodes{0.0,0.2,0.4,0.6,0.8,1.0};
mean_field::solver::TensorBernsteinPolynomial<2> polynomial(nodes);
const auto values = Values(nodes,2,[](auto p) { return 0.2+std::pow(p[0],5)+p[1]; });
CHECK(polynomial.Order() == 5);
CHECK(polynomial.Check(values).status == PolynomialSignStatus::certified_nonnegative);
mean_field::solver::TensorBernsteinPolynomial<1> quadratic(quadraticNodes);
const auto tangent = Values(quadraticNodes,1,[](auto p) { return (p[0]-0.37)*(p[0]-0.37); });
CHECK(quadratic.Check(tangent,{.maximumSubdivisionDepth=6}).status == PolynomialSignStatus::unresolved);
}
TEST_CASE("Feasible-anchor fractions describe a sufficient convex polynomial cone",
"[solver][enthalpy_admissibility][unit]") {
mean_field::solver::TensorBernsteinPolynomial<1> polynomial(cubicNodes);
const auto anchor = Values(cubicNodes,1,[](auto p) { return p[0]*(1.0-p[0]); });
const auto candidate = Values(cubicNodes,1,[](auto p) { return p[0]*(1.0-p[0])*(1.0-6.0*p[0]); });
const double fraction = polynomial.MinimumRepairFraction(anchor,candidate);
REQUIRE(std::isfinite(fraction)); REQUIRE(fraction > 0.0); REQUIRE(fraction <= 1.0);
std::vector<double> repaired(anchor.size());
const double safe = fraction+1.0e-10*(1.0-fraction);
for (std::size_t i = 0; i < repaired.size(); ++i) repaired[i] = (1.0-safe)*candidate[i]+safe*anchor[i];
CHECK(polynomial.Check(repaired).status == PolynomialSignStatus::certified_nonnegative);
CHECK(repaired.front() == 0.0); CHECK(repaired.back() == 0.0);
}
TEST_CASE("Prepared admissibility preserves shared true DOFs and rejects an invalid starting state",
"[solver][enthalpy_admissibility][integration][parallel]") {
using namespace mean_field::solver;
ParallelFixture fixture;
PreparedEnthalpyAdmissibility prepared(fixture.space,fixture.dofMap,fixture.elements);
const auto base = fixture.Project([](const mfem::Vector &) { return 1.0; });
const auto anchor = fixture.Project([](const mfem::Vector &p) { return 1.0+0.2*p(1); });
const auto direction = fixture.Project([](const mfem::Vector &p) { return -2.0-p(0)+0.5*p(2); });
const auto repair = prepared.MinimumRepairFraction(base,anchor,direction);
REQUIRE(repair.status == EnthalpyDirectionRepairStatus::feasible);
REQUIRE(repair.minimumSafeFraction > 0.0);
mfem::Vector repaired;
const double beta = repair.minimumSafeFraction+1.0e-8*(1.0-repair.minimumSafeFraction);
REQUIRE(prepared.ApplyRepairFraction(base,anchor,direction,beta,repaired).status ==
PolynomialSignStatus::certified_nonnegative);
mfem::Vector candidate(base); candidate += repaired;
CHECK(prepared.Check(candidate).status == PolynomialSignStatus::certified_nonnegative);
for (int i = 0; i < repaired.Size(); ++i)
CHECK(repaired(i) == Catch::Approx((1.0-beta)*direction(i)+beta*(anchor(i)-base(i))));
mfem::Vector invalid(base); invalid *= -1.0;
CHECK(prepared.MinimumRepairFraction(invalid,anchor,direction).status ==
EnthalpyDirectionRepairStatus::invalid_base);
int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD,&rank);
mfem::Vector nonfinite(base);
if (rank == 0 && nonfinite.Size()) nonfinite(0) = std::numeric_limits<double>::quiet_NaN();
CHECK(prepared.Check(nonfinite).status == PolynomialSignStatus::non_finite);
const auto saved = repaired;
CHECK(prepared.ApplyRepairFraction(base,anchor,direction,0.0,repaired).status == PolynomialSignStatus::negative);
for (int i = 0; i < saved.Size(); ++i) CHECK(repaired(i) == saved(i));
}
TEST_CASE("Feasible direction repair enforces exact homogeneous trace and rejects incompatible traces",
"[solver][enthalpy_admissibility][integration][parallel]") {
using namespace mean_field::solver;
ParallelFixture fixture;
mfem::Array<int> marker(fixture.mesh.bdr_attributes.Max()); marker = 1;
mfem::Array<int> boundary; fixture.space.GetEssentialTrueDofs(marker,boundary);
PreparedEnthalpyAdmissibility prepared(fixture.space,fixture.dofMap,fixture.elements,
std::span<const int>(boundary.GetData(),boundary.Size()));
const auto anchor = fixture.Project([](const mfem::Vector &p) {
return p(0)*(1.0-p(0))*p(1)*(1.0-p(1))*p(2)*(1.0-p(2));
});
const auto base = anchor;
mfem::Vector direction(anchor); direction *= -2.0;
for (int row : boundary) direction(row) = -0.25;
const auto report = prepared.MinimumRepairFraction(base,anchor,direction);
REQUIRE(report.status == EnthalpyDirectionRepairStatus::feasible);
mfem::Vector repaired;
REQUIRE(prepared.ApplyRepairFraction(base,anchor,direction,1.0,repaired).status ==
PolynomialSignStatus::certified_nonnegative);
for (int row : boundary) CHECK(repaired(row) == 0.0);
auto invalidTrace = base;
for (int row : boundary) invalidTrace(row) = 1.0e-20;
CHECK(prepared.MinimumRepairFraction(invalidTrace,anchor,direction).status ==
EnthalpyDirectionRepairStatus::invalid_base);
}