From c7eab38e1da1bb3c54905ed0ca373bb4eb4ac815 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 8 Oct 2026 09:32:19 -0400 Subject: [PATCH] feat(equations): Added barotropic closure This enforces that the model density and ethalpy are self consistent with the eos density from the model enthalpy --- docs/filters/requires.py | 0 docs/style/doxygen-extra.css | 0 .../serif/equations/barotropic_closure.hpp | 169 ++++++++++++++++++ 3 files changed, 169 insertions(+) create mode 100644 docs/filters/requires.py create mode 100644 docs/style/doxygen-extra.css create mode 100644 src/include/serif/equations/barotropic_closure.hpp diff --git a/docs/filters/requires.py b/docs/filters/requires.py new file mode 100644 index 00000000..e69de29b diff --git a/docs/style/doxygen-extra.css b/docs/style/doxygen-extra.css new file mode 100644 index 00000000..e69de29b diff --git a/src/include/serif/equations/barotropic_closure.hpp b/src/include/serif/equations/barotropic_closure.hpp new file mode 100644 index 00000000..7990a142 --- /dev/null +++ b/src/include/serif/equations/barotropic_closure.hpp @@ -0,0 +1,169 @@ +#pragma once + +/** + * @file barotropic_closure.hpp + * @author Emily M. Boudreaux + * + * @brief File containing the barotropic closure equation + */ + +#include +#include +#include +#include +#include + +#include "serif/discretization/blocks/fields.hpp" +#include "serif/eos/evaluation.hpp" +#include "serif/eos/exceptions.hpp" +#include "serif/eos/models/polytropic.hpp" +#include "serif/equations/point.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::equations { + // On stellar domains: integral [rho - rho_EOS(h)] q_rho dV = 0. + + /** + * @class BarotropicClosure + * @tparam EOST The equation of state type, this must satisfy the BarotropicClosureEOS concept. The requires section below lays out + * what the BarotriopicClosureEOS concept requires. + * + * @requires The equation of state type must publish a valid Relation catalog + * @requires The equation of state type must implement all relations required by the catalog it publishes + * @requires The equation of state type must support the relation density from specific enthalpy + * @requires The equation of state type must support the partial deriviative density from specific enthalpy with respect to specific enthalpy + * + * In a barotropic model surfaces of constant pressure coincide with surfaces of constant density and constant specific enthalpy. + * The barotropic closure equation enforces that the density of the model argrees with the density of the the equation of state + * given the model specific enthalpy. Put another way, for a barotropic model, this equation is what enforces + * that the ethalpy and density must be consistent with each other and with the equation of state. The residual + * is given by \f$ \rho - \rho_{\text{EOS}}(h) \f$, where \f$ \rho \f$ is the model density, \f$ h \f$ is the model specific enthalpy. + * + * The partial derivatives implemented are + * - \f$ \frac{\partial r}{\partial \rho} = 1 \f$ + * - \f$ \frac{\partial \rho}{\partial h} = -\frac{\partial \rho_{\text{EOS}}(h)}{\partial h} \f$ + * + * @note The eos derivative is implemented by the eos::partial_derivative function and is up to the eos model to implement correctly. + */ + template + class BarotropicClosure final { + public: + using Test = discretization::blocks::density::mass; ///< The test type is the model density, which is the quantity being constrained by the barotropic closure. + using ResidualQuantity = dimensions::MassDensity; ///< The residual quantity is the mass density, which is the difference between the model density and the density given by the equation of state. + using Enthalpy = discretization::blocks::enthalpy::specific; ///< The enthalpy type is the model specific enthalpy, which is the quantity being used to compute the density from the equation of state. Really this is just an alias. + + /** + * @brief The operand here is a pack of the model mass density and the model enthalpy. + */ + using Operands = utils::types::TypeList, discretization::forms::Operand>; + + static constexpr std::string_view name = "barotropic closure"; ///< Compile time name for the barotropic closure equation, used for debugging and logging. + + /** + * @brief Construct a barotropic closure equation. + * + * @param equation_of_state The equation of state to use for the closure. + * @param enthalpy_power The power of the enthalpy to use for the quadrature. Must be positive and finite. This is only used if the ::additional_order method is called. + * + * @throws std::invalid_argument If the enthalpy power is not positive and finite. + * + * @todo should we add a requires that the equation of state is the same as eos::models::Polytrope here? + */ + explicit BarotropicClosure( + EOST equation_of_state, + const double enthalpy_power + ) : + m_eos(std::move(equation_of_state)), m_enthalpy_power(enthalpy_power) { + using namespace utils::errors; + + serif_ensure(std::isfinite(enthalpy_power) && enthalpy_power > 0.0,"The density quadrature power must be positive and finite."); + } + + /** + * @brief Construct a barotropic closure equation for a polytropic equation of state. + * + * @param equation_of_state The polytropic equation of state to use for the closure. + * + * @requires The equation of state type must be a polytropic equation of state. + * + * @note This constructor is only available for polytropic equations of state. It uses the polytropic index as the enthalpy power. + * @note Generally this is the constructor that you should likely use. + */ + explicit BarotropicClosure( + EOST equation_of_state + ) + requires std::same_as + : BarotropicClosure(equation_of_state, equation_of_state.polytropic_index()) {} + + /** + * @brief Evaluate the barotropic closure equation point-wise integrand and its partial derivatives. + * + * The point-wise integrand is calculated as \f$ \rho - \rho_{\text{EOS}}(h) \f$, where \f$ \rho \f$ is the model density, and \f$ h \f$ is the model specific enthalpy. + * The partial derivatives are calculated as: + * - \f$ \frac{\partial r}{\partial \rho} = 1 \f$ + * - \f$ \frac{\partial \rho}{\partial h} = -\frac{\partial \rho_{\text{EOS}}(h)}{\partial h} \f$ + * + * @param point The point at which to evaluate the equation. + * @returns The point-wise linearization of the barotropic closure equation, which contains the residual and the partial derivatives with respect to the model density and enthalpy. + * + * @throws eos::EOSEvaluationError If the density or enthalpy are not finite, or if the enthalpy is negative, or if the equation of state returns a non-finite density or derivative. + */ + [[nodiscard]] PointLinearization evaluate( + const EquationPoint& point + ) const { + using namespace utils::errors; + + const double density = point.fields.template get(); + const double enthalpy = point.fields.template get(); + + serif_ensure(std::isfinite(density) && std::isfinite(enthalpy), "Closure requires finite density and enthalpy."); + serif_ensure(enthalpy >= 0.0, "Stellar closure requires non-negative enthalpy."); + + const dimensions::SpecificEnthalpyValue input{enthalpy}; + + const double eos_density = eos::evaluate(m_eos, input).value(); + const double density_derivative = eos::partial_derivative(m_eos, input).value(); + + serif_ensure(std::isfinite(eos_density) && std::isfinite(density_derivative), "Barotropic closure received a non-finite EOS result."); + + PointLinearization result; + result.residual = density - eos_density; + result.field_derivatives.template get() = 1.0; + result.field_derivatives.template get() = -density_derivative; + + return result; + } + + /** + * @brief Return the additional quadrature order required for the barotropic closure equation. + * + * The additional order is calculated as \f$ \text{order} = \text{density order} + 2 \times \text{position order} \f$, where the density order is the order of the functional space of the model density, and the position order is the order of the position functional space. + * + * @tparam OperandT The operand type for which to calculate the additional order. This is only used to determine if the operand is the enthalpy or not. If not enthalpy then the additional order is returned as 0. + * @param operand The operand type for which to calculate the additional order. This is only used to determine if the operand is the enthalpy or not. If not enthalpy then the additional order is returned as 0. + * @param position_order The order of the position functional space. + * + * @returns The additional quadrature order required for the barotropic closure equation. + * + * @throws std::invalid_argument If the extra order exceeds the supported range of int. + */ + template + [[nodiscard]] int additional_order(OperandT /*operand*/, int /*position_order*/) const { + using namespace utils::errors; + + if constexpr (std::same_as) { + const double extra = std::ceil(std::max(0.0, m_enthalpy_power - 1.0) * Enthalpy::FunctionalSpace::order); + + serif_ensure(extra <= std::numeric_limits::max(), "The EOS quadrature order exceeds the supported range."); + + return static_cast(extra); + } + + return 0; + } + + private: + EOST m_eos; ///< The equation of state to use for the barotropic closure. (@todo should this be a reference?) + double m_enthalpy_power; ///< The power of the enthalpy to use for the quadrature. Must be positive and finite. This is only used if the ::additional_order method is called. + }; +} // namespace serif::equations