From cb2ac127aad41d149390ee9925c3481174382f33 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 12 Oct 2026 09:15:32 -0400 Subject: [PATCH] feat(equations): added pressure_force equation this equation describes the contribution of pressure to the mechanical balance of the system --- .../serif/equations/pressure_force.hpp | 154 ++++++++++++++++++ 1 file changed, 154 insertions(+) create mode 100644 src/include/serif/equations/pressure_force.hpp diff --git a/src/include/serif/equations/pressure_force.hpp b/src/include/serif/equations/pressure_force.hpp new file mode 100644 index 00000000..63f21ff2 --- /dev/null +++ b/src/include/serif/equations/pressure_force.hpp @@ -0,0 +1,154 @@ +#pragma once + +/** + * @file pressure_force.hpp + * @author Emily M. Boudreaux + * + * @brief File containing the PressureForce contribution to the displacement balance 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/contributions.hpp" +#include "serif/equations/point.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::equations { + /** + * @brief The pressure force contribution to the displacement balance equation. This equation + * is really the only one that makes heavy use of the eos. Broadly, when we switch to a non polytropic + * case this will likely need some major refactoring. + * + * @tparam EOST The EOS model type to use for evaluating the pressure from the specific enthalpy. + * + * @requires The EOS model must support the relation PressureFromSpecificEnthalpy and its partial derivative with respect to specific enthalpy. + * + * @contribution Residuals to the displacement vector residual block where only on the parts of that block which are part of the stellar domain. This contribution is through the divergence of the test function. + */ + template + class PressureForce final : public Contribution< + discretization::forms::Operand, + discretization::domain::StellarDomains + > { + public: + using Enthalpy = discretization::blocks::enthalpy::specific; ///< Alias for the specific enthalpy + using EnthalpyValue = discretization::forms::Operand; ///< Alias for the specific enthalpy operand + + using Operands = utils::types::TypeList; ///< This equation only depends on the specific enthalpy operand + + static constexpr std::string_view name = "pressure force"; ///< Compile time name for the pressure force equation, used for debugging and logging. + + /** + * @brief Construct a pressure force contribution with a given EOS model and enthalpy power. + * + * @param equation_of_state The EOS model to use for evaluating the pressure from the specific enthalpy. Note that the PressureForce will take ownership of the EOS model. + * @param enthalpy_power The effective power of the pressure as a function of specific enthalpy, which must be positive and finite. This is used for quadrature selection. + * + * @throws std::invalid_argument If the enthalpy power is not finite and strictly positive. + */ + explicit PressureForce( + 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 pressure quadrature power must be positive and finite."); + } + + /** + * @brief Construct a pressure force contribution with a given EOS model. This constructor is only available for polytropic EOS models, where the enthalpy power is inferred from the polytropic index. + * In this case we call the other constructor with the polytropic index + 1 as the enthalpy power. + * + * @param equation_of_state The polytropic EOS model to use for evaluating the pressure from the specific enthalpy. Note that the PressureForce will take ownership of the EOS model. + * + * @requires The EOS model must be a polytropic EOS model, otherwise this constructor is not available. + */ + explicit PressureForce( + EOST equation_of_state + ) + requires std::same_as + : PressureForce(equation_of_state, equation_of_state.polytropic_index() + 1.0) {} + + /** + * @brief Evaluate the pressure force contribution at a given point. + * + * @param point The point at which to evaluate the pressure force contribution. + * @return A PointLinearization object containing the residual and field derivatives for the pressure force contribution at the given point. + * + * + * Here we let the point-wise residual return as \f$-P\f$, recall that since we set the contribution operand to be the divergence and supported on + * the stellar domain this will be integrated as (for some trial function \f$\phi\f$) + * + * \f[ + * R = \int_{\Omega_{\text{stellar}}} -\left(\nabla \mathbf{\phi} \right)P dV + * \f] + * + * Further, this implements the partial derivative of this residual with respect to the + * specific enthalpy as + * + * - \f$ \frac{\partial R}{\partial h} = -\frac{\partial P}{\partial h} \f$ + * + * @throws eos::EOSEvaluationError If the EOS model returns a non-finite pressure or pressure derivative for the given specific enthalpy. + * @throws eos::EOSEvaluationError If the specific enthalpy is not finite or is negative, which is outside the domain of the EOS model. + */ + [[nodiscard]] PointLinearization evaluate( + const EquationPoint& point + ) const { + using namespace utils::errors; + + // TODO: Here is a place where the assumptions about a polytrope have slipped into my code. + // As we move to a more general case we cannot rely on the enthalpy to be a general carrier field + const double enthalpy = point.fields.get(); + + serif_ensure(std::isfinite(enthalpy), "Pressure force requires finite enthalpy."); + serif_ensure(enthalpy >= 0.0, "Stellar pressure force requires non-negative enthalpy."); + + const dimensions::SpecificEnthalpyValue input{enthalpy}; + + const double pressure = eos::evaluate(m_eos, input).value(); + const double pressure_derivative = eos::partial_derivative(m_eos, input).value(); + + serif_ensure(std::isfinite(pressure) && std::isfinite(pressure_derivative), "Pressure force received a non-finite EOS result."); + + PointLinearization result; + result.residual = -pressure; + result.field_derivatives.get() = -pressure_derivative; + + return result; + } + + /** + * @brief Return the additional quadrature order requires for integrating the pressure force contribution. + * This is calculated as \f$\lceil \max(0, n-1) p \rceil\f$ where \f$n\f$ is the enthalpy power and \f$p\f$ is the order of the specific enthalpy functional space. + * + * @param operand-type The operand type for which to determine the additional quadrature order. This is not used in this implementation since the pressure force contribution only depends on the specific enthalpy operand. + * @param position_order The order of the position operand, which is not used in this implementation. + * + * @throws std::invalid_argument If the calculated additional quadrature order exceeds the maximum value of an int. + * + * @return The additional quadrature order required for integrating the pressure force contribution. + */ + [[nodiscard]] int additional_order(EnthalpyValue /*operand-type*/, int /*position_order*/) const { + using namespace utils::errors; + + 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 pressure quadrature order exceeds the supported range."); + + return static_cast(extra); + } + + private: + EOST m_eos; ///< The EOS model to use for evaluating the pressure from the specific enthalpy. + double m_enthalpy_power; ///< The effective power of the pressure as a function of specific enthalpy, which must be positive and finite. This is used for quadrature selection. + }; +} // namespace serif::equations