From 08758315251fdfcde627d4ad481ae2d36f6c99b0 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 12 Oct 2026 10:48:08 -0400 Subject: [PATCH] feat(equations): added rotation_force. This constrains the rotation force based on ridig rotation. When we add differential rotation we will need to dramatically modify this equation --- .../serif/equations/rotational_force.hpp | 114 ++++++++++++++++++ 1 file changed, 114 insertions(+) create mode 100644 src/include/serif/equations/rotational_force.hpp diff --git a/src/include/serif/equations/rotational_force.hpp b/src/include/serif/equations/rotational_force.hpp new file mode 100644 index 00000000..4ba6c642 --- /dev/null +++ b/src/include/serif/equations/rotational_force.hpp @@ -0,0 +1,114 @@ +#pragma once +/** + * @file rotational_force.hpp + * @author Emily M. Boudreaux + * + * @brief File containing the RotationalForce contribution class, which computes the rotational force contribution to the displacement balance equation. + */ + +#include + +#include "serif/discretization/blocks/fields.hpp" +#include "serif/equations/contributions.hpp" +#include "serif/equations/point.hpp" +#include "serif/equations/rigid_rotation.hpp" + +namespace serif::equations { + /** + * @class RotationalForce + * + * @brief The rotational force contribution to the displacement balance equation. This equation computes the rotational force contribution + * to the displacement balance equation, which is given by \f$-\int_{\Omega_{\text{stellar}}}\rho\nabla\Psi \cdot w \, dV\f$ where \f$\Psi\f$ is the centrifugal potential + * and \f$w\f$ is the test function. + * + * @contribution Residuals to the displacement vector residual block where only on the parts of that block which are part of the stellar domain. + * + */ + class RotationalForce final : public Contribution, discretization::domain::StellarDomains> { + public: + using Density = discretization::blocks::density::mass; ///< Alias for mass density + using DensityValue = discretization::forms::Operand; ///< Operand type for mass density + + using AngularVelocity = discretization::blocks::rotation::angular_velocity; ///< Alias for angular velocity + using AngularVelocityValue = discretization::forms::Operand; ///< Operand type for angular velocity + + using Operands = utils::types::TypeList; ///< The rotational force depends on the mass density and the angular velocity operands. + using Result = discretization::forms::VectorValue; ///< The result type will be a vector value. + + static constexpr std::string_view name = "rotational force"; ///< Compile time name for the rotational force equation, used for debugging and logging. + + /** + * @brief Construct a rotational force contribution with a given axis and center of rotation. + * + * @param axis The axis of rotation, which must be finite. Defaults to the z-axis. + * @param center The center of rotation, which must be finite. Defaults to the origin. + */ + explicit RotationalForce( + const Result &axis = {0.0, 0.0, 1.0}, + const Result ¢er = {} + ) : + m_rotation(RigidRotation::unit_rotation(axis, center)) {} + + + /** + * @brief Evaluate the rotational force contribution at a given point. + * + * @param point The point at which to evaluate the rotational force contribution, which must contain the mass density and angular velocity operands. + * @return A PointLinearization containing the residual and derivatives of the rotational force contribution at the given point. + * + * Here the point-wise integrand is computed as a three-vector, where each component \f$r_{i}\f$ is given by \f$ r_{i} = \rho\left(-\omega^{2}\nabla\Psi_{i}\right) + * + * There are then derivatives of this given with respect to mass density, angular velocity, and position + * + * - \frac{\partial r_{i}}{\partial \rho} = -\omega^{2}\nabla\Psi_{i} + * - \frac{\partial r_{i}}{\partial \omega} = -2\omega\rho\nabla\Psi_{i} + * - \frac{\partial r_{i}}{\partial x_{j}} = -\rho \omega^{2} \frac{\partial^{2} \Psi_{i}}{\partial x_{i} \partial x_{j}} + */ + [[nodiscard]] PointLinearization evaluate( + const EquationPoint& point + ) const { + const double density = point.fields.get(); + const double omega = point.fields.get(); + const double omega_squared = omega * omega; + const auto gradient = m_rotation.centrifugal_potential_gradient(point.position); + + PointLinearization result; + + for (int row = 0; row < 3; ++row) { + result.residual[row] = density * (-omega_squared * gradient[row]); + result.field_derivatives.get()[row] = -omega_squared * gradient[row]; + result.field_derivatives.get()[row] = -2.0 * omega * density * gradient[row]; + for (int column = 0; column < 3; ++column) { + result.position_derivative[row][column] = -density * (omega_squared * m_rotation.centrifugal_potential_hessian()[row][column]); + } + } + + return result; + } + + /** + * @brief Get the additional quadrature order required for the operands of this equation when integrating density. + * + * @param operand-type The operand type for which to determine the additional quadrature order. Note this is only used to select the correct overload. + * @param position_order The quadrature order associated with the position. + * @return The additional quadrature order required for the given operand type. For density this just returns the position order; + */ + [[nodiscard]] static int additional_order(DensityValue /*operand-type*/, const int position_order) { + return position_order; + } + + /** + * @brief Get the additional quadrature order required for the operands of this equation when integrating angular velocity. + * + * @param operand-type The operand type for which to determine the additional quadrature order. Note this is only used to select the correct overload. + * @param position_order The quadrature order associated with the position. + * @return The additional quadrature order required for the given operand type. For angular velocity this returns the sum of the density functional space order and the position order. + */ + [[nodiscard]] static int additional_order(AngularVelocityValue /*operand-type*/, const int position_order) { + return Density::FunctionalSpace::order + position_order; + } + + private: + RigidRotation m_rotation; ///< The rigid rotation object used to compute the centrifugal potential and its derivatives, which is constructed from the axis and center of rotation provided to the constructor. + }; +} // namespace serif::equations