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

feat(equations): added rigid rotation helper

Browse Source 
This is a class which adds some helper utilities for computing the effective potential and the hessian of a body under rigid rotation. Note this is not an equation in and of itself
This commit is contained in:
Emily Boudreaux
2026-10-12 10:33:24 -04:00
parent fee05a621d
commit f8bc1a37a5
1 changed files with 236 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

236
src/include/serif/equations/rigid_rotation.hpp Normal file
View File

@@ -0,0 +1,236 @@
#pragma once
/**
* @file rigid_rotation.hpp
* @author Emily M. Boudreaux
*
* @brief An equation implementing rigid rotation
*/
#include <cmath>
#include <stdexcept>
#include "serif/discretization/forms/value.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::equations
{
/// Rigid rotation: ω = Ωn̂, λ = r − rₐ, v = ω × λ, and Ψ = ½|v|² = ½Ω²s².
/// Here n̂ is the unit axis, rₐ a point on it, and s = |n̂ × λ| the cylindrical radius.
/// Ψ is the positive centrifugal potential: ∇Ψ points outward and Φ_eff = Φ_grav − Ψ.
/**
* @class RigidRotation
*
* @brief A class that implements rigid rotation, providing methods to compute the centrifugal potential, its gradient, and its Hessian.
*
* @note This is not an equation in the same way that other parts of this module are; rather, this is a utility class used in other
* equations.
*
* @todo This should probably be refactored into its own submodule to clean up equations a bit
*
* Broadly this class handles a number of the details of rigid rotation. Note that this is really setup to be static physics; that is to say,
* there are a series of things you setup with the constructor (the angular velocity vector and the axis origin), then those
* are used to evaluate a series of relations. There is not however any way to update and re-evaluate the same
* object.
*
* The relations implemented are
* - Centrifugal potential: \f$\Psi = \frac{1}{2}|\omega \times \lambda|^2\f$ where \f$\lambda = r - r_a\f$ is the position relative to the axis origin. This can be further simplified to \f$\Psi = \frac{1}{2}\Omega^{2}s^{2}\f$ where \f$s\f$ is the
* cylindrical radius from the axis of rotation.
* - Centrifugal potential gradient: \f$\nabla \Psi = H\ lambda\f$ where \f$H = |\omega|^2 I - \omega \otimes \omega\f$ is the Hessian of the centrifugal potential.
*/
class RigidRotation final
{
public:
using Vector = discretization::forms::VectorValue; ///< Alias for some vector (rank 1 tensor) type
using Matrix = discretization::forms::MatrixValue; ///< Alias for some matrix (rank 2 tensor) type
/**
* @brief Construct a RigidRotation object with the given angular velocity vector and axis origin.
*
* @param angular_velocity_vector The angular velocity vector, which must be finite.
* @param axis_origin The origin of the rotation axis, which must be finite. (default to the origin, (0, 0, 0))
*
* @throws std::invalid_argument If the angular velocity vector or axis origin is not finite.
*
* The constructor also evaluates the Hessian of the centrifugal potential at construction time. This is defined as
*
* \f[
* H = \nabla\nabla \Psi = |\omega|^2 I - \omega \otimes \omega
* \f]
*
* @todo We make an assumption here about the dimensionality of the space (e.g that the hessian is a 3x3 matrix). This
* is likely fine; however, if we ever want to be able to model 1D or 2D systems we will need to refactor this to be more general.
*/
explicit RigidRotation(
const Vector &angular_velocity_vector,
const Vector &axis_origin = {}
) :
m_angular_velocity_vector(angular_velocity_vector),
m_axis_origin(axis_origin)
{
using namespace utils::errors;
double angular_speed_squared = 0.0;
bool finite_axis_origin = true;
// |ω|² = ωₓ² + ωᵧ² + ω_z².
for (int component = 0; component < 3; ++component)
{
finite_axis_origin = finite_axis_origin && std::isfinite(axis_origin[component]);
angular_speed_squared += angular_velocity_vector[component] * angular_velocity_vector[component];
}
serif_ensure<SERiFErrorCode::NON_FINITE_VALUE_ERROR, std::invalid_argument>(finite_axis_origin && std::isfinite(angular_speed_squared), "Rigid rotation requires finite center coordinates and finite squared angular speed; reduce overflowing parameters.");
// H = ∇∇Ψ = |ω|²I − ω⊗ω.
for (int row = 0; row < 3; ++row)
{
for (int column = 0; column < 3; ++column)
{
m_potential_hessian[row][column] = (row == column ? angular_speed_squared : 0.0) - angular_velocity_vector[row] * angular_velocity_vector[column];
}
}
}
/**
* @brief Construct a RigidRotation object with a unit rotation axis and an axis origin.
*
* @param rotation_axis The rotation axis, which must be finite and nonzero. (default to the z-axis, (0, 0, 1))
* @param axis_origin The origin of the rotation axis, which must be finite. (default to the origin, (0, 0, 0))
*
* @throws std::invalid_argument If the norm of the rotation axis is not finite and positive.
*
* This method normalizes the rotation axis to create a unit vector. It then constructs a RigidRotation object using this unit vector and the provided axis origin.
*
* @note This constructs a new RigidRotation object does not modify the original object in place.
*
* @todo This may make sense to turn into a free-function rather than a class method.
*/
[[nodiscard]] static RigidRotation unit_rotation(Vector rotation_axis = {0.0, 0.0, 1.0}, const Vector &axis_origin = {})
{
using namespace utils::errors;
const double axis_norm = std::hypot(rotation_axis[0], rotation_axis[1], rotation_axis[2]);
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(std::isfinite(axis_norm) && axis_norm > 0.0, "A rotation axis must be finite and nonzero.");
for (double& component : rotation_axis)
{
component /= axis_norm;
}
return RigidRotation(rotation_axis, axis_origin);
}
/// Ψ = ½|ω × λ|²; using the cross product avoids subtracting nearly equal squared lengths near the axis.
/**
* @brief Compute the centrifugal potential at a given position. This is computed as
* \f$ \Psi = \frac{1}{2}|\omega \times \lambda|^2\f$ where \f$\lambda = r - r_a\f$ is the position relative to the axis origin.
*
* @param position The position at which to compute the centrifugal potential, which must be finite.
* @return The centrifugal potential at the given position.
*
* @todo We should add some validation for position into this
*/
[[nodiscard]] double centrifugal_potential(const Vector& position) const
{
double linear_speed_squared = 0.0;
for (int component = 0; component < 3; ++component)
{
// These are allowed to cycle through the components of the vector.
// For some component c this will look like. Note this is
// just a trick to get the appropriate indices for
// the cross product.
/*
* | component | next | last |
* |===========|======|======|
* | 0 | 1 | 2 |
* | 1 | 2 | 0 |
* | 2 | 0 | 1 |
*/
const int next = (component + 1) % 3;
const int last = (component + 2) % 3;
// vᵢ = ωⱼλₖ − ωₖλⱼ, with (i, j, k) cyclic.
// Get to this by remembering
// v = ω x λ
// so then, trivially from the cross product, v has the form
// v_x = ω_y λ_z - ω_z λ_y
// v_y = ω_z λ_x - ω_x λ_z
// v_z = ω_x λ_y - ω_y λ_x
// The equation for the potential is
// φ = 1/2 |v|^2
// Since we can calculate each component of v using the cyclic permutation of the
// indices it is trivially to compute the squared magnitude of v by summing the squares
// of each component. Then we simply multiply by 1/2 (outside of the loop)
// to get the potential.
const double velocity_component = m_angular_velocity_vector[next] * (position[last] - m_axis_origin[last]) - m_angular_velocity_vector[last] * (position[next] - m_axis_origin[next]);
linear_speed_squared += velocity_component * velocity_component;
}
return 0.5 * linear_speed_squared;
}
/**
* @brief Compute the gradient of the centrifugal potential at a given position. This is computed as
* \f$ \nabla \Psi = H \lambda\f$ where \f$H = |\omega|^2 I - \omega \otimes \omega\f$ is the Hessian of the centrifugal potential.
*
* @param position The position at which to compute the gradient of the centrifugal potential, which must be finite.
* @return The gradient of the centrifugal potential at the given position.
*
* @todo We should add some validation for position into this
*/
[[nodiscard]] Vector centrifugal_potential_gradient(const Vector& position) const
{
Vector result{};
for (int row = 0; row < 3; ++row)
{
for (int column = 0; column < 3; ++column)
{
result[row] += m_potential_hessian[row][column] * (position[column] - m_axis_origin[column]);
}
}
return result;
}
/**
* @brief Get the Hessian of the centrifugal potential, which is independent of position for rigid rotation. This is computed as
* \f$ H = |\omega|^2 I - \omega \otimes \omega\f$.
*
* @return The Hessian of the centrifugal potential.
*/
[[nodiscard]] const Matrix& centrifugal_potential_hessian() const
{
return m_potential_hessian;
}
/**
* @brief Get the angular velocity vector of the rigid rotation.
*
* @return The angular velocity vector.
*/
[[nodiscard]] const Vector& angular_velocity_vector() const
{
return m_angular_velocity_vector;
}
/**
* @brief Get the origin of the rotation axis.
*
* @return The origin of the rotation axis.
*/
[[nodiscard]] const Vector& axis_origin() const
{
return m_axis_origin;
}
private:
Vector m_angular_velocity_vector; ///< The angular velocity vector of the rigid rotation, which must be finite.
Vector m_axis_origin; ///< The origin of the rotation axis, which must be finite.
Matrix m_potential_hessian{}; ///< The Hessian of the centrifugal potential, which is independent of position for rigid rotation and is computed as \f$ H = |\omega|^2 I - \omega \otimes \omega\f$.
};
} // namespace serif::equations