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

Reject unsupported pressure traces before polynomial state certification

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 10:15:57 -04:00
parent 6e812aeca2
commit 071417f0c1
3 changed files with 61 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

42
libmeanfield/interface/solver/stellar_enthalpy_strategy.cppm
View File

@@ -18,6 +18,44 @@ export import :solver.enthalpy_admissibility;
export import :normalization.stellar_equilibrium;
export import :utils.domain;
export namespace mean_field::solver::detail {
template <class SurfaceConstraint>
concept PolynomialEnthalpySurfaceConstraint = requires {
typename SurfaceConstraint::CarrierQuantity;
typename SurfaceConstraint::Relation;
requires std::same_as<typename SurfaceConstraint::CarrierQuantity,dimensions::quantity::SpecificEnthalpy>;
requires std::same_as<typename SurfaceConstraint::Relation,eos::SpecificEnthalpyFromPressure>;
};
// The polynomial engine fixes homogeneous trace rows. Check the physical
// prescription before inspecting coefficients, so an unsupported surface
// condition is never reported as an invalid numerical starting state.
template <PolynomialEnthalpySurfaceConstraint SurfaceConstraint, class EquationOfState>
void RequirePolynomialEnthalpySurface(const SurfaceConstraint &surface,
const EquationOfState &equationOfState, MPI_Comm communicator) {
int localUnsupported = surface.targetPressure().value() == 0.0 ? 0 : 1;
int globalUnsupported = 0;
if (MPI_Allreduce(&localUnsupported,&globalUnsupported,1,MPI_INT,MPI_MAX,communicator) != MPI_SUCCESS)
throw std::runtime_error("Could not validate the polynomial enthalpy surface prescription.");
if (globalUnsupported)
throw std::invalid_argument("PolynomialEnthalpy supports only zero prescribed surface pressure (Psurf=0).");
bool homogeneous = false;
try {
const double trace = eos::evaluate<dimensions::quantity::SpecificEnthalpy>(
equationOfState,surface.targetPressure()).value();
homogeneous = std::isfinite(trace) && trace == 0.0;
} catch (...) {
// A rank-local EOS exception must not strand other ranks in the
// collective compatibility check that follows.
}
localUnsupported = homogeneous ? 0 : 1;
if (MPI_Allreduce(&localUnsupported,&globalUnsupported,1,MPI_INT,MPI_MAX,communicator) != MPI_SUCCESS)
throw std::runtime_error("Could not validate the zero-pressure enthalpy reference.");
if (globalUnsupported)
throw std::invalid_argument("PolynomialEnthalpy requires the EOS to prescribe h(Psurf=0)=0.");
}
}
export namespace mean_field::solver {
/*
* The solver owns this workspace only when a positivity policy is selected.
@@ -28,6 +66,8 @@ export namespace mean_field::solver {
class PreparedStellarEnthalpyStrategy final {
using Problem = std::remove_cvref_t<decltype(std::declval<State &>().Problem())>;
using Form = typename Problem::FormType;
static_assert(detail::PolynomialEnthalpySurfaceConstraint<typename Problem::CompiledSurfaceConstraintType>,
"PolynomialEnthalpy requires a state-independent specific-enthalpy pressure surface relation.");
using Schema = utils::domain::CoreEnvelopeVacuumDomainSchema;
using Domain = field::FieldDomainT<field::Enthalpy>;
static constexpr int enthalpyIndex = utils::blocks::type_index_v<
@@ -39,6 +79,8 @@ export namespace mean_field::solver {
: m_state(&state),
m_offset(state.Problem().GetManifest().layout().offset(enthalpyBlock)),
m_size(state.Problem().GetManifest().layout().size(enthalpyBlock)) {
detail::RequirePolynomialEnthalpySurface(state.Problem().GetCompiledSurfaceConstraint(),
state.Problem().GetStellarModel().equationOfState(),state.Problem().GetCommunicator());
const auto &fem = equilibrium::detail::StellarEquilibriumProblemFactory::FiniteElementModel(state.Problem());
if (!fem.mesh || !fem.enthalpyFes)
throw std::invalid_argument("Stellar enthalpy certification requires a complete enthalpy FE space.");