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

feat(equations): added hydrostatic_equilibrium equation

Browse Source
This commit is contained in:
Emily Boudreaux
2026-10-08 13:11:09 -04:00
parent 6a273fdd9d
commit 9bd77a661f
1 changed files with 190 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

190
src/include/serif/equations/hydrostatic_equilibrium.hpp Normal file
View File

@@ -0,0 +1,190 @@
#pragma once
#include <algorithm>
#include <cmath>
#include <string_view>
#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<Test>,
discretization::forms::Operand<Potential>,
discretization::forms::Operand<Constant>,
discretization::forms::Operand<AngularVelocity>
>;
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 &center = {}
) :
m_rotation(RigidRotation::unit_rotation(axis, center)) {}
/// Evaluate h + Phi - Psi - C and its partial derivatives at a physical point.
[[nodiscard]] PointLinearization<Operands> evaluate(const EquationPoint<Operands>& point) const {
PointLinearization<Operands> result;
const double omega = point.fields.get<AngularVelocity>();
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<Test>() 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<Test>() + point.fields.get<Potential>() - omega_squared * m_rotation.centrifugal_potential(point.position) - point.fields.get<Constant>();
result.field_derivatives.get<Test>() = 1.0;
result.field_derivatives.get<Potential>() = 1.0;
result.field_derivatives.get<Constant>() = -1.0;
result.field_derivatives.get<AngularVelocity>() = -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 <typename OperandT>
[[nodiscard]] int additional_order(OperandT, int position_order) const {
if constexpr (discretization::blocks::HasFunctionalSpace<typename OperandT::term_type>) {
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