diff --git a/experiments/negative_density_projection.hpp b/experiments/negative_density_projection.hpp index b6bc646..8741808 100644 --- a/experiments/negative_density_projection.hpp +++ b/experiments/negative_density_projection.hpp @@ -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::make_query( quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array{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::make_query( + 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 { diff --git a/experiments/plateau_enthalpy_step_diagnostics.hpp b/experiments/plateau_enthalpy_step_diagnostics.hpp index 87a752a..ba971dc 100644 --- a/experiments/plateau_enthalpy_step_diagnostics.hpp +++ b/experiments/plateau_enthalpy_step_diagnostics.hpp @@ -127,7 +127,16 @@ std::map MeasureEnthalpyStep(State &state, const mean_field mean_field::quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array{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::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; diff --git a/experiments/plateau_pressure_obstruction_diagnostics.hpp b/experiments/plateau_pressure_obstruction_diagnostics.hpp index 153dada..54ce1c4 100644 --- a/experiments/plateau_pressure_obstruction_diagnostics.hpp +++ b/experiments/plateau_pressure_obstruction_diagnostics.hpp @@ -409,7 +409,16 @@ void ProbeSurfaceTrace(State &s, const mean_field::fem::FEM &f, double n, : EnthalpyField::template make_query( quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array{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::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); diff --git a/libmeanfield/impl/operators/kernels/barotropic_closure_kernels.cpp b/libmeanfield/impl/operators/kernels/barotropic_closure_kernels.cpp index c66fd7d..9c83e51 100644 --- a/libmeanfield/impl/operators/kernels/barotropic_closure_kernels.cpp +++ b/libmeanfield/impl/operators/kernels/barotropic_closure_kernels.cpp @@ -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::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 " diff --git a/libmeanfield/impl/operators/prepared_barotropic_closure.cpp b/libmeanfield/impl/operators/prepared_barotropic_closure.cpp index a57e60e..3b1739d 100644 --- a/libmeanfield/impl/operators/prepared_barotropic_closure.cpp +++ b/libmeanfield/impl/operators/prepared_barotropic_closure.cpp @@ -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::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, diff --git a/tests/field/field_mfem.cpp b/tests/field/field_mfem.cpp index 44929c4..7b556c5 100644 --- a/tests/field/field_mfem.cpp +++ b/tests/field/field_mfem.cpp @@ -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); } -} \ No newline at end of file +} diff --git a/tests/operators/prepared_barotropic_closure.cpp b/tests/operators/prepared_barotropic_closure.cpp index 6c64bf6..9475090 100644 --- a/tests/operators/prepared_barotropic_closure.cpp +++ b/tests/operators/prepared_barotropic_closure.cpp @@ -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(std::ceil( + (condition.polytropicIndex - 1.0) * field::Enthalpy::Scalar::familyOrder + )); + const auto sourceQuery = field::Field::make_query< + field::Enthalpy::Form::EosClosureSource>( + mean_field::quadrature::QuadratureRole::discretization, transformation.OrderW(), + std::array{extraOrder}, mean_field::utils::DOMAINS::STELLAR, + mean_field::quadrature::MappingKind::general + ); + const auto massQuery = field::Field::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); diff --git a/tests/quadrature/policy.cpp b/tests/quadrature/policy.cpp index 3edac14..d744b7d 100644 --- a/tests/quadrature/policy.cpp +++ b/tests/quadrature/policy.cpp @@ -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);