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

Cover both density mass and EOS source in closure quadrature

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-19 10:47:58 -04:00
parent d6554cae6b
commit dd955abca5
8 changed files with 98 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
experiments/negative_density_projection.hpp
View File

@@ -238,12 +238,19 @@ inline NegativeDensityProjectionReport MeasureNegativeDensityProjection(
const mfem::IntegrationRule* rule = nullptr;
int resolvedOrder = control.quadratureOrder;
if (control.productionQuadrature) {
// PreparedBarotropicClosureOperator uses this source query for
// both rho*q and rho(h)*q, not the separate legacy mass query.
// Reproduce the shared production rule, including the density
// mass degree when it exceeds the nonlinear source degree.
const auto query = field::Field<field::Enthalpy>::make_query<field::Enthalpy::Form::EosClosureSource>(
quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array<int, 1>{extraOrder},
utils::DOMAINS::STELLAR, quadrature::MappingKind::general);
const auto resolution = fem.quadratureFactory->get(query, transformation->GetGeometryType());
const auto sourceResolution = fem.quadratureFactory->get(query, transformation->GetGeometryType());
const auto massQuery = field::Field<field::Density>::make_query<field::Density::Form::EosClosureMass>(
quadrature::QuadratureRole::discretization, transformation->OrderW(), {},
utils::DOMAINS::STELLAR, quadrature::MappingKind::general);
const auto massResolution = fem.quadratureFactory->get(massQuery, transformation->GetGeometryType());
const auto &resolution = sourceResolution.resolution.order >= massResolution.resolution.order
? sourceResolution
: massResolution;
rule = resolution.integration_rule;
resolvedOrder = resolution.resolution.order;
} else {

11
experiments/plateau_enthalpy_step_diagnostics.hpp
View File

@@ -127,7 +127,16 @@ std::map<std::string, double> MeasureEnthalpyStep(State &state, const mean_field
mean_field::quadrature::QuadratureRole::discretization, transformation->OrderW(),
std::array<int, 1>{extraOrder}, mean_field::utils::DOMAINS::STELLAR,
mean_field::quadrature::MappingKind::general);
const auto resolution = fem.quadratureFactory->get(query, transformation->GetGeometryType());
auto resolution = fem.quadratureFactory->get(query, transformation->GetGeometryType());
if (!pressure) {
const auto massQuery = mean_field::field::Field<mean_field::field::Density>::template make_query<
mean_field::field::Density::Form::EosClosureMass>(
mean_field::quadrature::QuadratureRole::discretization, transformation->OrderW(), {},
mean_field::utils::DOMAINS::STELLAR, mean_field::quadrature::MappingKind::general);
const auto massResolution = fem.quadratureFactory->get(massQuery, transformation->GetGeometryType());
if (massResolution.resolution.order > resolution.resolution.order)
resolution = massResolution;
}
if (!resolution.integration_rule)
throw std::runtime_error("Missing production EOS quadrature rule.");
const auto &rule = *resolution.integration_rule;

11
experiments/plateau_pressure_obstruction_diagnostics.hpp
View File

@@ -409,7 +409,16 @@ void ProbeSurfaceTrace(State &s, const mean_field::fem::FEM &f, double n,
: EnthalpyField::template make_query<field::Enthalpy::Form::EosClosureSource>(
quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array<int, 1>{extraOrder},
utils::DOMAINS::STELLAR, quadrature::MappingKind::general);
const auto rule = f.quadratureFactory->get(query, transformation->GetGeometryType());
auto rule = f.quadratureFactory->get(query, transformation->GetGeometryType());
if (!pressure) {
const auto massQuery = field::Field<field::Density>::template make_query<
field::Density::Form::EosClosureMass>(
quadrature::QuadratureRole::discretization, transformation->OrderW(), {},
utils::DOMAINS::STELLAR, quadrature::MappingKind::general);
const auto massResolution = f.quadratureFactory->get(massQuery, transformation->GetGeometryType());
if (massResolution.resolution.order > rule.resolution.order)
rule = massResolution;
}
if (!rule.integration_rule || point < 0 || point >= rule.integration_rule->GetNPoints())
throw std::runtime_error("Surface trace inspection could not resolve its production limiting point.");
const auto &ip = rule.integration_rule->IntPoint(point);

13
libmeanfield/impl/operators/kernels/barotropic_closure_kernels.cpp
View File

@@ -102,7 +102,18 @@ namespace {
mean_field::quadrature::MappingKind::general
);
const auto resolution = f.quadratureFactory->get(query, transformation.GetGeometryType());
const auto sourceResolution = f.quadratureFactory->get(query, transformation.GetGeometryType());
// The reference kernel must also cover the density mass degree when
// density is enriched independently of the EOS source field.
const auto massQuery = mean_field::field::Field<mean_field::field::Density>::make_query<
mean_field::field::Density::Form::EosClosureMass>(
mean_field::quadrature::QuadratureRole::discretization, transformation.OrderW(), {},
mean_field::utils::DOMAINS::STELLAR, mean_field::quadrature::MappingKind::general
);
const auto massResolution = f.quadratureFactory->get(massQuery, transformation.GetGeometryType());
const auto &resolution = sourceResolution.resolution.order >= massResolution.resolution.order
? sourceResolution
: massResolution;
MFEM_VERIFY(
resolution.integration_rule != nullptr, "The quadrature policy did not return an "

14
libmeanfield/impl/operators/prepared_barotropic_closure.cpp
View File

@@ -146,7 +146,19 @@ namespace {
mean_field::quadrature::MappingKind::general
);
const auto resolution = f.quadratureFactory->get(query, transformation.GetGeometryType());
const auto sourceResolution = f.quadratureFactory->get(query, transformation.GetGeometryType());
// Enriching density independently can make rho*q higher order than
// rho(h)*q. Resolve both through the policy so fixed-order overrides
// retain their meaning while the shared rule covers both forms.
const auto massQuery = mean_field::field::Field<mean_field::field::Density>::make_query<
mean_field::field::Density::Form::EosClosureMass>(
mean_field::quadrature::QuadratureRole::discretization, transformation.OrderW(), {},
mean_field::utils::DOMAINS::STELLAR, mean_field::quadrature::MappingKind::general
);
const auto massResolution = f.quadratureFactory->get(massQuery, transformation.GetGeometryType());
const auto &resolution = sourceResolution.resolution.order >= massResolution.resolution.order
? sourceResolution
: massResolution;
MFEM_VERIFY(
resolution.integration_rule != nullptr,

6
tests/field/field_mfem.cpp
View File

@@ -393,11 +393,11 @@ TEST_CASE(
REQUIRE(densitySource.base_order.has_value());
/*
* L2_2 value order 2
* Registered density value order
* + geometry order 3
* + dynamic coefficient order 4.
*/
CHECK(*densitySource.base_order == 9);
CHECK(*densitySource.base_order == field::Density::Scalar::familyOrder + 3 + 4);
CHECK(densitySource.term == quadrature::Term::density_projection);
@@ -775,4 +775,4 @@ TEST_CASE(
for (int index = 0; index < support.activeTrueDofMarker.Size(); ++index) {
CHECK(support.activeTrueDofMarker[index] == 1);
}
}
}

32
tests/operators/prepared_barotropic_closure.cpp
View File

@@ -839,6 +839,38 @@ namespace prepared_barotropic_closure_test_utils {
make_state_view(baseDensity, baseEnthalpy, baseDisplacement), make_dependencies()
);
// Agreement with the stateless kernel alone cannot expose a
// shared quadrature mistake. Check the actual prepared rule
// against both form orders, allowing MFEM's odd-order aliases.
int checkedRules = 0;
preparedOperator.VisitMappedGeometryRules([&](const int element, const mfem::IntegrationRule &rule) {
const auto &transformation = *f.mesh->GetElementTransformation(element);
const int extraOrder = static_cast<int>(std::ceil(
(condition.polytropicIndex - 1.0) * field::Enthalpy::Scalar::familyOrder
));
const auto sourceQuery = field::Field<field::Enthalpy>::make_query<
field::Enthalpy::Form::EosClosureSource>(
mean_field::quadrature::QuadratureRole::discretization, transformation.OrderW(),
std::array<int, 1>{extraOrder}, mean_field::utils::DOMAINS::STELLAR,
mean_field::quadrature::MappingKind::general
);
const auto massQuery = field::Field<field::Density>::make_query<
field::Density::Form::EosClosureMass>(
mean_field::quadrature::QuadratureRole::discretization, transformation.OrderW(), {},
mean_field::utils::DOMAINS::STELLAR, mean_field::quadrature::MappingKind::general
);
const auto sourceResolution =
f.quadratureFactory->get(sourceQuery, transformation.GetGeometryType());
const auto massResolution =
f.quadratureFactory->get(massQuery, transformation.GetGeometryType());
INFO("Prepared closure element = " << element);
INFO("Actual quadrature exactness = " << rule.GetOrder());
CHECK(rule.GetOrder() >= sourceResolution.resolution.order);
CHECK(rule.GetOrder() >= massResolution.resolution.order);
++checkedRules;
});
CHECK(checkedRules > 0);
mfem::Vector preparedResidual;
mfem::Vector preparedAction;
preparedOperator.BuildResidual(preparedResidual);

12
tests/quadrature/policy.cpp
View File

@@ -545,14 +545,18 @@ TEST_CASE(
STATIC_CHECK(hdiv_query.base_order.has_value());
STATIC_CHECK(*hdiv_query.base_order == 8);
STATIC_CHECK(*divergence_query.base_order == 6);
STATIC_CHECK(*source_query.base_order == 6);
STATIC_CHECK(*center_of_mass_query.base_order == 5);
STATIC_CHECK(*eos_closure_query.base_order == 13);
STATIC_CHECK(
*source_query.base_order == field::Density::Scalar::familyOrder + field::Gravity::Potential::familyOrder + 2
);
STATIC_CHECK(*center_of_mass_query.base_order == field::Density::Scalar::familyOrder + 1 + 2);
STATIC_CHECK(
*eos_closure_query.base_order == field::Enthalpy::Scalar::familyOrder + field::Density::Scalar::familyOrder + 6 + 2
);
STATIC_CHECK(*equilibrium_gravity_query.base_order == 7);
STATIC_CHECK(*rotation_query.base_order == 7);
STATIC_CHECK(*isobaric_surface_query.base_order == 8);
STATIC_CHECK(*mesh_extension_query.base_order == 6);
STATIC_CHECK(*mass_normalization_query.base_order == 4);
STATIC_CHECK(*mass_normalization_query.base_order == field::Density::Scalar::familyOrder + 2);
STATIC_CHECK(*pressure_integral_query.base_order == 14);
STATIC_CHECK(*equilibrium_constant_query.base_order == 5);