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

feat(quadrature): quadrature system brought over

Browse Source 
Ported and dramatically cleaned up the quadrature system. This includes centralizing all field definitions
This commit is contained in:
Emily Boudreaux
2026-09-19 07:35:57 -04:00
parent 7b4f30d945
commit c44584c1e7
75 changed files with 8936 additions and 172 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
src/lib/serif/dimensions/runtime/runtime.cpp
View File

@@ -1,13 +0,0 @@
#include "serif/dimensions/runtime/runtime.hpp"
namespace serif::dimensions::runtime {
constexpr ThermodynamicQuantityID::ThermodynamicQuantityID(const std::string_view name) noexcept : m_name(name) {}
[[nodiscard]] constexpr std::string_view ThermodynamicQuantityID::name() const noexcept {
return m_name;
}
[[nodiscard]] constexpr bool RuntimeRelationDescriptor::hasPartialDerivative(const std::size_t inputIndex) const noexcept {
return inputIndex < inputQuantities.size() && (partialDerivativeMask & (std::uint64_t{1} << inputIndex)) != 0;
}
}

103
src/lib/serif/discretization/quadrature/policy.cpp Normal file
View File

@@ -0,0 +1,103 @@
#include <stdexcept>
#include <format>
#include <print>
#include "../../../../include/serif/discretization/quadrature/policy.hpp"
#include "serif/utils/warnings/warning.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::discretization::quadrature {
using namespace serif::utils::warnings;
using namespace serif::utils::errors;
Policy::Policy(RuleSet ruleSet) : m_ruleSet(std::move(ruleSet)) {
}
int Policy::compute_base_order(const Query &query) {
if (query.base_order.has_value()) {
if (*query.base_order < 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_QUADRATURE_ORDER>("Quadrature base_order must be non-negative. Received: {}. This is likely a configuration error", *query.base_order);
}
return *query.base_order;
}
if (query.trial_order < 0 || query.test_order < 0 || query.coefficient_order < 0 || query.geometry_weight_order < 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_QUADRATURE_ORDER>("Quadrature orders must be non-negative. Received: trial_order={}, test_order={}, coefficient_order={}, geometry_weight_order={}. This is likely a configuration error", query.trial_order, query.test_order, query.coefficient_order, query.geometry_weight_order);
}
if (query.trial_order_reduction < 0 || query.test_order_reduction < 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_QUADRATURE_ORDER>("Quadrature order reductions must be non-negative. Received: trial_order_reduction={}, test_order_reduction={}. Note these are reductions, they are subtracted from the corresponding order. You do not need to include a negative sign with them. This is likely a configuration error", query.trial_order_reduction, query.test_order_reduction);
}
const int reduced_trial_order = query.trial_order - query.trial_order_reduction;
const int reduced_test_order = query.test_order - query.test_order_reduction;
if (reduced_trial_order < 0) {
serif_warning("Quadrature trial order reduction ({}) is greater than the trial order ({}). Clamping this value to 0", query.trial_order_reduction, query.trial_order);
}
if (reduced_test_order < 0) {
serif_warning("Quadrature test order reduction ({}) is greater than the test order ({}). Clamping this value to 0", query.test_order_reduction, query.test_order);
}
const int trial_order = std::max(0, reduced_trial_order);
const int test_order = std::max(0, reduced_test_order);
return trial_order + test_order + query.coefficient_order + query.geometry_weight_order;
}
const RuleControl &Policy::get_role_control(QuadratureRole role) const {
return m_ruleSet.role_control(role);
}
Resolution Policy::resolve(const Query &query) const {
const int base_order = compute_base_order(query);
const RuleControl *form_control = m_ruleSet.find_form_control(query.form);
const RuleControl &role_control = get_role_control(query.role);
const RuleControl &fallback_control = m_ruleSet.fallback_control();
std::optional<int> fixed_order;
if (form_control != nullptr && form_control->fixed_order.has_value()) {
fixed_order = form_control->fixed_order;
} else if (role_control.fixed_order.has_value()) {
fixed_order = role_control.fixed_order;
} else if (fallback_control.fixed_order.has_value()) {
fixed_order = fallback_control.fixed_order;
}
if (fixed_order.has_value()) {
if (*fixed_order < 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_QUADRATURE_ORDER>("Quadrature fixed_order must be non-negative. Received: {}. This is likely a configuration error", *fixed_order);
}
return Resolution{
.base_order = base_order,
.boost = 0,
.order = *fixed_order,
.used_fixed_order = true
};
}
int boost = fallback_control.boost + role_control.boost;
if (form_control != nullptr) {
boost += form_control->boost;
}
const int order = base_order + boost;
if (order < 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_QUADRATURE_ORDER>("Quadrature order must be non-negative. Received: {}. Realistically this should not be able to happen and represents a likely bug in the code. Please contact the SERiF developers.", order);
};
return Resolution{
.base_order = base_order,
.boost = boost,
.order = order,
.used_fixed_order = false
};
}
}

23
src/lib/serif/eos/models/polytropic.cpp
View File

@@ -2,11 +2,14 @@
#include "serif/utils/misc/finite.hpp"
#include "serif/eos/exceptions.hpp"
#include "serif/utils/error/codes.hpp"
namespace serif::eos::models {
using utils::misc::validate_finite;
using utils::misc::validate_greater_than;
using utils::misc::validate_less_than;
using utils::misc::validate_nonnegativity;
using utils::errors::SERiFErrorCode;
using namespace dimensions;
using namespace relations;
@@ -17,10 +20,10 @@ namespace serif::eos::models {
m_polytropic_index(n),
m_polytropic_constant(K),
m_enthalpy_scale((n + 1.0) * K) {
validate_finite<EOSEvaluationError>(m_polytropic_index, EOSEvaluationErrorCode::outside_domain);
validate_greater_than<EOSEvaluationError>(m_polytropic_index, 1.0, EOSEvaluationErrorCode::outside_domain);
validate_finite<EOSEvaluationError>(m_polytropic_constant, EOSEvaluationErrorCode::outside_domain);
validate_greater_than<EOSEvaluationError>(m_polytropic_constant, 0.0, EOSEvaluationErrorCode::outside_domain);
validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_OUTSIDE_DOMAIN>(m_polytropic_index);
validate_greater_than<EOSEvaluationError, SERiFErrorCode::EOS_OUTSIDE_DOMAIN>(m_polytropic_index, 1.0);
validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_OUTSIDE_DOMAIN>(m_polytropic_constant);
validate_greater_than<EOSEvaluationError, SERiFErrorCode::EOS_OUTSIDE_DOMAIN>(m_polytropic_constant, 0.0);
}
[[nodiscard]] double Polytrope::polytropic_index() const noexcept {
@@ -36,7 +39,7 @@ namespace serif::eos::models {
}
[[nodiscard]] PressureValue Polytrope::evaluate(PressureFromDensity, const DensityValue density) const {
validate_nonnegativity<EOSEvaluationError>(density.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(density.value());
if (density.value() == 0.0) {
// P = 0
@@ -50,7 +53,7 @@ namespace serif::eos::models {
}
[[nodiscard]] SpecificEnthalpyValue Polytrope::evaluate(SpecificEnthalpyFromDensity, const DensityValue density) const {
validate_nonnegativity<EOSEvaluationError>(density.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(density.value());
if (density.value() == 0.0) {
// h = 0
return SpecificEnthalpyValue{0.0};
@@ -63,7 +66,7 @@ namespace serif::eos::models {
}
[[nodiscard]] DensityValue Polytrope::evaluate(DensityFromSpecificEnthalpy, SpecificEnthalpyValue specific_enthalpy) const {
validate_finite<EOSEvaluationError>(specific_enthalpy.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(specific_enthalpy.value());
if (specific_enthalpy.value() <= 0.0) {
// ρ = 0
@@ -86,7 +89,7 @@ namespace serif::eos::models {
}
[[nodiscard]] SpecificEnthalpyValue Polytrope::evaluate(SpecificEnthalpyFromPressure, PressureValue pressure) const {
validate_nonnegativity<EOSEvaluationError>(pressure.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(pressure.value());
if (pressure.value() == 0.0) {
// h = 0
return SpecificEnthalpyValue{0.0};
@@ -102,7 +105,7 @@ namespace serif::eos::models {
}
[[nodiscard]] PartialDerivative<Density, SpecificEnthalpy> Polytrope::partial_derivative(DensityFromSpecificEnthalpy, WithRespectTo<SpecificEnthalpy>, const SpecificEnthalpyValue specific_enthalpy) const {
validate_finite<EOSEvaluationError>(specific_enthalpy.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(specific_enthalpy.value());
if (specific_enthalpy.value() < 0.0) {
// dρ/dh = 0
@@ -131,7 +134,7 @@ namespace serif::eos::models {
}
[[nodiscard]] PartialDerivative<Pressure, Density> Polytrope::partial_derivative(PressureFromDensity, WithRespectTo<Density>, const DensityValue density) const {
validate_nonnegativity<EOSEvaluationError>(density.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(density.value());
if (density.value() == 0.0) {
// dP/dρ = 0