From 9f67fccd4599f3ae66ad46c0d1e1020aeb8e848e Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 7 Oct 2026 11:14:14 -0400 Subject: [PATCH] feat(equations/angular_momentum): added angular momentum constraint equation --- .../serif/equations/angular_momentum.hpp | 206 ++++++++++++++++++ src/include/serif/equations/equations.hpp | 0 2 files changed, 206 insertions(+) create mode 100644 src/include/serif/equations/angular_momentum.hpp create mode 100644 src/include/serif/equations/equations.hpp diff --git a/src/include/serif/equations/angular_momentum.hpp b/src/include/serif/equations/angular_momentum.hpp new file mode 100644 index 00000000..c25063f5 --- /dev/null +++ b/src/include/serif/equations/angular_momentum.hpp @@ -0,0 +1,206 @@ +#pragma once + +/** + * @file angular_momentum.hpp + * @author Emily M. Boudreaux + * + * @brief File containing the angular momentum constraint equation, which enforces a global constraint on the total angular momentum of a system. + * + * @note The angular momentum constraint is a global constraint that is enforced by integrating the point-wise contributions of the angular momentum density over the + * entire domain. The integrand is given by \f$ j_\parallel = \rho \Omega s^2\f$, where \f$s^2 = |\hat{n} × (r − r_{c})|^2\f$ is the cylindrical radius squared, + * and \f$\hat{n}\f$ is the unit rotation axis. The solver will integrate this over the domain and compare it to the target angular momentum to enforce the constraint. + * + * A few important notes about the implementation + * + * - The way the equation system plugs into the mfem backend is such that equations should only ever return point wise contributions to the residual, not + * a global integrated residual. This is why we return \f$ j_\parallel = \rho \Omega s^2\f$ as the residual, rather than \f$ j_\parallel - J_{\text{target}}\f$. However + * we still need to subtract the target angular momentum from the integrated residual to enforce the constraint. This is done by implementing the method + * ::residual_offset() which returns \f$ -J_{\text{target}}\f$. The backend will automatically detect and, if present, subtract this from the integrated residual to + * enforce the constraint. + * + * - The solver can choose to validate the angular momentum constraint by integrating the ::validation_integrand() over the domain and then passing + * this result to ::validate_integral(double result). Note that this is not specific to this equation but a feature of the overall system + */ + +#include +#include + +#include "serif/equations/contributions.hpp" +#include "serif/equations/point.hpp" +#include "serif/equations/rigid_rotation.hpp" + +namespace serif::equations +{ + /** + * @brief Angular momentum constraint equation, which enforces a global constraint on the total angular momentum of a system. + * + * The point-wise integrand is given by \f$ j_\parallel = \rho \Omega s^2\f$, where \f$s^2 = |\hat{n} × (r − r_{c})|^2\f$ is the cylindrical radius squared, + * and \f$\hat{n}\f$ is the unit rotation axis. The solver will integrate this over the domain and compare it to the target angular momentum to enforce the constraint. + * + * Use this equation in conjunction with a RigidRotation object to constrain the angular momentum of a system to some fixed value. Note that this constraint + * is weakly applied and therefore should be checked for convergence in the final solution. + */ + class AngularMomentumConstraint final : public Contribution> + { + public: + using Density = discretization::blocks::density::mass; + using ResidualQuantity = dimensions::AngularMomentum; + using AngularVelocity = discretization::blocks::rotation::angular_velocity; + using DensityValue = discretization::forms::Operand; + using AngularVelocityValue = discretization::forms::Operand; + using Operands = utils::types::TypeList; + + static constexpr std::string_view name = "angular momentum"; + + /// Prescribe Jₜ ≥ 0 and the axis through rₐ; Ω remains an unknown and may be signed in trial states. + explicit AngularMomentumConstraint( + const double target_axial_angular_momentum, + const RigidRotation::Vector &rotation_axis = {0.0, 0.0, 1.0}, + const RigidRotation::Vector &axis_origin = {} + ) : + m_target_axial_angular_momentum(target_axial_angular_momentum), + m_unit_rotation(RigidRotation::unit_rotation(rotation_axis, axis_origin)) + { + using namespace utils::errors; + serif_ensure(std::isfinite(target_axial_angular_momentum) && target_axial_angular_momentum >= 0.0, "Target angular momentum must be finite and nonnegative."); + } + + /** + * @brief Evaluate the angular-momentum constraints point-wise integrand and its partial derivatives. + * + * The point-wise integrand is calculated as \f$ j_\parallel = \rho \Omega s^2\f$ where \f$s^2 = |\hat{n} × (r − r_{c})|^2\f$ is the cylindrical radius squared. The partial derivatives are calculated as: + * - \f$ \frac{\partial j_\parallel}{\partial \rho} = \Omega s^2\f$ + * - \f$ \frac{\partial j_\parallel}{\partial \Omega} = \rho s^2\f$ + * - \f$ \nabla j_\parallel = 2 \rho\Omega\lambda_{\perp}\f$ where \f$\lambda_{\perp} = \hat{n} × (r − r_{c})\f$ is the displacement perpendicular to the rotation axis. + * + * @note There is a subtly here which may trick the careful programmer up. Note that we set the residual + * to \f$R = \rho\Omega s^{2}\f$; rather that the more obvious choice of \f$R = \rho\Omega s^{2} - J_{\text{target}}\f$. + * The important thing to remember is that these equations are evaluated point-wise. Put another way if you were + * to evaluate the physical residual as the local angular momentum contribution - the global angular momentum at some single point the meaning would be ill-defined. There is no local + * angular momentum constraint, we only have a global one. What we actually need to do is give the backend someone of knowing what the total angular + * momentum of the system is. Therefore what we return here is integrand of + * the residual which is then automatically integrated over the volume by the backend (so when integrating all these point wise contributions over the domain the backend will get the total angular momentum). Further, if an equation + * supplies the method ::residual_offset() (as this one does) then the backend will automatically subtract that + * from the integrated residual (e.g. the backend has computed the total angular momentum and then subtracts residual_offset to find the actual residual, which will approach 0 as the total angular momentum approaches the supplied value). A heuristic which may be useful to keep in mind is that global constraints should + * supply a residual offset whereas local equations often will not need to. + * + * @param point The point at which to evaluate the integrand and its partial derivatives. + */ + [[nodiscard]] PointLinearization evaluate(const EquationPoint& point) const + { + const double mass_density = point.fields.get(); + const double angular_velocity = point.fields.get(); + + // Unit rotation: s² = 2Ψ₁ = |n̂ × (r − rₐ)|²; ∇Ψ₁ = λ⊥. + const double cylindrical_radius_squared = 2.0 * m_unit_rotation.centrifugal_potential(point.position); + + // dI∥/dV = ρs²; λ⊥ is the displacement perpendicular to the axis. + const double moment_of_inertia_density = mass_density * cylindrical_radius_squared; + const auto perpendicular_displacement = m_unit_rotation.centrifugal_potential_gradient(point.position); + + PointLinearization result; + result.residual = angular_velocity * moment_of_inertia_density; // j∥ = Ωρs². + + result.field_derivatives.get() = angular_velocity * cylindrical_radius_squared; // ∂j∥/∂ρ = Ωs². + result.field_derivatives.get() = moment_of_inertia_density; // ∂j∥/∂Ω = ρs². + + // ∇j∥ = 2Ωρλ⊥, holding ρ and Ω fixed. + for (int component = 0; component < 3; ++component) + { + result.position_derivative[component] = 2.0 * angular_velocity * mass_density * perpendicular_displacement[component]; + } + return result; + } + + /** + * @brief Return the additional quadrature degree for the density coupling. + * + * At fixed position and angular velocity, \f$\delta j_{\parallel} = \Omega s^{2} \delta \rho\f$. + * The backend already counts the density trial degree \f$p_{\rho}\f$. Since \f$\Omega\f$ is + * spatially constant and \f$s^{2}\f$ is quadratic in position, the remaining + * coefficient \f$\Omega s^{2}\f$ has polynomial degree at most \f$2p_{r}\f$. + * + * @param position_order Polynomial degree bound for physical position + * expressed in reference-element coordinates. + * @return The coefficient degree bound \f$2p_{r}\f$, excluding the geometry weight. + */ + [[nodiscard]] static int additional_order(DensityValue, int position_order) + { + return 2 * position_order; + } + + /** + * @brief Return the additional quadrature degree for the angular-velocity coupling. + * + * At fixed density and position, \f$\delta j_{\parallel} = \rho s^{2}\delta \Omega\f$. + * The angular-velocity trial function is spatially constant. The remaining + * coefficient \f$\rho s^{2}\f$ has polynomial degree at most \f$p_{\rho} + 2p_{r}\f$. + * + * @param position_order Polynomial degree bound for physical position + * expressed in reference-element coordinates. + * @return The coefficient degree bound \f$p_{\rho} + 2p_{r}\f$, excluding the geometry weight. + */ + [[nodiscard]] static int additional_order(AngularVelocityValue, int position_order) + { + return Density::FunctionalSpace::order + 2 * position_order; + } + + /** + * @brief Return the residual offset for the angular-momentum constraint. + * + * The backend integrates the point-wise residual contributions over the domain. The global constraint + * is satisfied when the total angular momentum equals the target value. Therefore we subtract the target + * from the integrated residual to get the actual residual. + * + * @note The backend automatically looks for this method and calls it if it exists. + * + * @return The negative of the target axial angular momentum, $-J_{\text{target}}$. + */ + [[nodiscard]] double residual_offset() const + { + return -m_target_axial_angular_momentum; + } + + + /** + * @brief Return the already computed inertia density, \f$\rho s^{2}\f$, for validation. + * + * When an equation provides validate_integral(), the operator backend integrates + * this quantity over the equation's domain and passes the result to that method. + * Here the integral is the axial moment of inertia, + * \f$I_{\parallel} = \int_\star \rho s^{2}\,dV\f$. + * + * The equation defines both the quantity to integrate and its admissibility + * requirements; the backend performs the integration. Here we require finite, + * nonnegative inertia as a necessary, but not sufficient, condition for a + * physically admissible mass distribution. + * + * @param physics The point-wise linearization, whose angular-velocity derivative + * already contains \f$\partial j_{\parallel}/\partial\Omega = \rho s^{2}\f$. + */ + [[nodiscard]] static double validation_integrand(const EquationPoint&, const PointLinearization& physics) + { + return physics.field_derivatives.get(); + } + + /** + * @brief Require the integrated axial moment of inertia to be finite and nonnegative. + * + * The operator backend supplies the integral of validation_integrand(). + * Failure rejects preparation of the candidate state. Passing this check + * does not establish nonnegative density everywhere or satisfaction of + * \f$\Omega I_{\parallel} - J_{\mathrm{target}} = 0\f$. + * + * @param axial_moment_of_inertia The integrated axial moment of inertia, + * \f$I_{\parallel} = \int_\star \rho s^{2}\,dV\f$. + */ + static void validate_integral(const double axial_moment_of_inertia) { + using namespace utils::errors; + serif_ensure(std::isfinite(axial_moment_of_inertia) && axial_moment_of_inertia >= 0.0, "The moment of inertia must be finite and nonnegative."); + } + + private: + double m_target_axial_angular_momentum; ///< The prescribed total angular momentum along the rotation axis, \f$ J_{\text{target}} \geq 0 \f$. + RigidRotation m_unit_rotation; ///< The unit rotation, storing the normalized axis and its origin, used to compute \f$s^{2}\f$ and its gradient. + }; +} // namespace serif::equations diff --git a/src/include/serif/equations/equations.hpp b/src/include/serif/equations/equations.hpp new file mode 100644 index 00000000..e69de29b