From 9bd77a661fd265c9695e07bd5ca4de7f20871d07 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 8 Oct 2026 13:11:09 -0400 Subject: [PATCH] feat(equations): added hydrostatic_equilibrium equation --- .../equations/hydrostatic_equilibrium.hpp | 190 ++++++++++++++++++ 1 file changed, 190 insertions(+) create mode 100644 src/include/serif/equations/hydrostatic_equilibrium.hpp diff --git a/src/include/serif/equations/hydrostatic_equilibrium.hpp b/src/include/serif/equations/hydrostatic_equilibrium.hpp new file mode 100644 index 00000000..db4bcbda --- /dev/null +++ b/src/include/serif/equations/hydrostatic_equilibrium.hpp @@ -0,0 +1,190 @@ +#pragma once + +#include +#include +#include + +#include "serif/discretization/blocks/fields.hpp" +#include "serif/equations/point.hpp" +#include "serif/equations/rigid_rotation.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::equations { + + // So we had a previous version of this code which had roughly 20k lines of bespoke operators for each equation + // down that route lie madness eh? + // We have tried to make this production version far more generic. HydrostaticEquilibrium here is a useful + // test case since its on the simpler side (being a scalar valued thing) + + /* + * + * @class BarotropicHydrostaticEquilibrium + * @brief The constraint balancing the pressure gradient, gravitational potential, and centrifugal potential. + * + * @todo This may make sense to break into smaller contributions in the same way as the displacement force equation + * does already. + * + * If you are a physicist, astronomer, or general science goober looking to modify physics in SERiF + * this is likely a good place to start. Here we lay out a simple case of how to implement an + * equation. + * + * Consider the strong form of this equation, generally + * + * \f[ + * \theta = -\left(\frac{1}{\rho}\right)\nabla P - \nabla \phi + \omega^2 s + * \f] + * + * where + * \f$s = \sqrt{x^2 + y^2}\f$ and is the cylindrical radius + * + * In this particular equation we specialize this for a barotrope, for a barotrope + * + * \f[ + * dh = \frac{dP}{\rho} + * \f] + * therefore + * \f[ + * \nabla h = \frac{1}{\rho} \nabla P + * \f] + * + * Then + * + * \f[ + * \theta = -\nabla h - \nabla \phi + \omega^2 s + * \f] + * or equivalently + * \f[ + * \theta = \nabla h + \nabla \phi - \omega^2 s + * \f] + * + * Then, in the corotating frame, we can define an effective potential \f$\psi = \frac{1}{2}\omega^{2}\left(x^{2} + y^{2}\right)\f$ + * + * therefore + * + * \f[ + * \theta = \nabla h + \nabla \phi - \nabla \psi = \nabla (h + \phi - \psi) + * \f] + * + * hence + * + * + * \f[ + * h + \phi - \psi = C + * \f] + * or + * \f[ + * h + \phi - \psi - C = 0 + * \f] + * + * This is a pointwise residual equation. Now broadly for finite element methods we need the weak form of this equation; this system is + * broadly and roughly intended to reduce friction. Look at the above equation then investigate HydrostaticEquilibrium::evaluate + * where we set the residual in the return value. The backend system integrates the produce of these residuals and a test function over + * a domain. The goal here is that the above step is where you as a research get to stop. + * + * Now there are a few other things to specify here. For one which term does this residual contribute too. You must decide this + * however there are natural choices. In the case of HydrostaticEquilibrium the residual is a scalar and so it contributes to the enthalpy equation. + * You must alias the residual term to a Test member. + * + * In addition to enthalpy computing the residual also depends on the gravitational potential, though the residual does not + * contribute to the gravitational potential equation. You must specify this dependency by listing the gravitational potential as an operand in the Operands member. + * + * In addition to the residual there are a few additional derivatives that should be implemented or explicitly set to 0, broadly these + * should be a field derivative for each operand and, if the term has a spatial component, a position derivative. + * + * The implementation below shows these trivially + * + * Other Notes + * =========== + * - This is a particularly simple case as we start with a scalar valued operation. + * - This equation does not use the Contributions system (mostly this is because I wrote this equation first before I wrote the contribution system). + * There is a benefit to this in so far as I think it makes it a bit more approachable initially as it is more self contained; however, + * it does break the pattern that other equations use. + * - This is very much a place where the polytropic assumption has slipped in. This equation will need to be reworked or replaced to move beyond polytropes. + */ + class BarotropicHydrostaticEquilibrium final { + public: + using Test = discretization::blocks::enthalpy::specific; ///< The test type is the model specific enthalpy, which is the quantity being constrained by the hydrostatic equilibrium equation. + using ResidualQuantity = dimensions::SpecificEnthalpy; ///< The residual quantity is the specific enthalpy, which is the difference between the model specific enthalpy and the specific enthalpy given by the gravitational and centrifugal potentials. + using Potential = discretization::blocks::gravity::potential; ///< We alias the gravitational potential to this type for convenience. + using Constant = discretization::blocks::barotropic_constant::mass_normalization; ///< We alias the barotropic constant to this type for convenience. + using AngularVelocity = discretization::blocks::rotation::angular_velocity; ///< We alias the angular velocity to this type for convenience. + + /** + * @brief The operand here is a pack of the model specific enthalpy, the gravitational potential, the barotropic constant, and the angular velocity. + * + * Think of this as saying that the residual for this equation depends on all of the things in this type list. + */ + using Operands = utils::types::TypeList< + discretization::forms::Operand, + discretization::forms::Operand, + discretization::forms::Operand, + discretization::forms::Operand + >; + + static constexpr std::string_view name = "hydrostatic equilibrium"; ///< Compile time name for the hydrostatic equilibrium equation, used for debugging and logging. + + /** + * @brief Construct a barotropic hydrostatic equilibrium equation. + * + * @param axis The axis of rotation for the centrifugal potential. This is a vector in physical space. + * @param center The center of rotation for the centrifugal potential. This is a vector in physical space. + * + * @todo Currently there is nothing constraining this axis or position to agree with the axis and center or rotation set + * in the force equation. The defaults are the same; however, we should introduce some system to ensure that these agree + * and that if they dont agree an error is thrown (since I cannot imagine any world in which a reasonable model + * would have these disagree). + */ + explicit BarotropicHydrostaticEquilibrium( + const RigidRotation::Vector &axis = {0.0, 0.0, 1.0}, + const RigidRotation::Vector ¢er = {} + ) : + m_rotation(RigidRotation::unit_rotation(axis, center)) {} + + /// Evaluate h + Phi - Psi - C and its partial derivatives at a physical point. + [[nodiscard]] PointLinearization evaluate(const EquationPoint& point) const { + PointLinearization result; + const double omega = point.fields.get(); + const double omega_squared = omega * omega; + + // Below there is some level of maybe wierd (though not awful in the context of the language) C++ syntax. We need to get access to the correct + // field values in a type safe manner, they are tuples so we use the get method with the + // types we aliased at the top of this class. At other points in the code you may see something like + // point.fields.template get() which is a bit more verbose but is the same thing. The template keyword is required in that case because + // there are cases where there is an ambiguity between a template and a comparator operator (e.g. <) (I think this is the reason, frankly I've never + // really sorted it out, I generally write the .template version and them my IDE tells me the template is not needed sometimes...) + + // Here is where we implement the residual, this is directly from the derivation above + result.residual = point.fields.template get() + point.fields.get() - omega_squared * m_rotation.centrifugal_potential(point.position) - point.fields.get(); + + result.field_derivatives.get() = 1.0; + result.field_derivatives.get() = 1.0; + result.field_derivatives.get() = -1.0; + result.field_derivatives.get() = -2.0 * omega * m_rotation.centrifugal_potential(point.position); + + result.position_derivative = m_rotation.centrifugal_potential_gradient(point.position); + for (double& component : result.position_derivative) { component *= -omega_squared; } + + return result; + } + + /** + * @brief Return the additional quadrature order required for the operands of this equation. + * If the operand has a functional space then the additional order is 2 * position_order - functional_space_order, otherwise it is 2 * position_order. + * + * @param operand The operand type for which to determine the additional quadrature order. + * @param position_order The order of the position basis functions. + * @return The additional quadrature order required for the given operand type. + */ + template + [[nodiscard]] int additional_order(OperandT, int position_order) const { + if constexpr (discretization::blocks::HasFunctionalSpace) { + return std::max(0, 2 * position_order - OperandT::term_type::FunctionalSpace::order); + } else { + return 2 * position_order; + } + } + + private: + RigidRotation m_rotation; ///< The rigid rotation object used to compute the centrifugal potential and its gradient. This is constructed from the axis and center of rotation provided to the constructor. + }; +} // namespace serif::equations