From 1a222ef4b1c357ed168fa6e3ceb4526029c92836 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 9 Oct 2026 16:15:26 -0400 Subject: [PATCH] feat(equations): Added mixed formulation gravity equations we solve gravity with a mixed formulation to avoid the loss of one order of accuracy while differentiating the potential. --- src/include/serif/equations/mixed_gravity.hpp | 313 ++++++++++++++++++ 1 file changed, 313 insertions(+) create mode 100644 src/include/serif/equations/mixed_gravity.hpp diff --git a/src/include/serif/equations/mixed_gravity.hpp b/src/include/serif/equations/mixed_gravity.hpp new file mode 100644 index 00000000..5526c645 --- /dev/null +++ b/src/include/serif/equations/mixed_gravity.hpp @@ -0,0 +1,313 @@ +#pragma once + +/** + * @file mixed_gravity.hpp + * @author Emily M. Boudreaux + * + * @brief File containing the mixed gravity equation and its contributions. + * + * In SERiF we used a mixed formulation when solving for the gravitational potential and gravitational field. + * This is a pragmatic choice to try to avoid loss of accuracy when differentiating the potential to get the field. + * + * More specifically, if we were to just solve for the potential and then differentiate to get the field + * we would lose one order of accuracy in the field when compared to the potential. Due to + * the required accuracy in virialization (which is fundamentally related to the gravitational field) + * we cannot afford to loose that order and maintain the low-ish resolution mesh we use. I was already + * solving for potential when I went this route so I wanted to keep it. I suppose there is some world where you + * could solve for the field and then integrate it to get the potential; however, I have + * not thought deeply about the effects that would have on either the accuracy or run-time + * performance. + * + * As it stands we solve for the potential and the field separately but simultaneously and we weakly + * constrain them to agree. + */ + +#include +#include +#include +#include +#include "serif/equations/contributions.hpp" +#include "serif/equations/point.hpp" +#include "serif/utils/error/errors.hpp" + + +namespace serif::equations { + + // TODO: I know that these deeply nested namespaces are really a C++ no-no (or as we in the biz call them an anit-pattern). I think + // my background as a python dev makes me like them, but perhaps at some point we should go through and flatten this out so we + // avoid having these long alias definitions at the starts of classes and namespaces (which is also bad in a header to begin with)? + // Something to think about + + namespace gravity_contributions { + using Gradient = discretization::blocks::gravity::gradient; ///< Alias for the gravitational field type, which is the gradient of the gravitational potential. + using Potential = discretization::blocks::gravity::potential; ///< Alias for the gravitational potential type, which is the scalar potential field that generates the gravitational field. + using Density = discretization::blocks::density::mass; ///< Alias for the mass density type, which is the source term in the Poisson equation for gravity. + + using Value = discretization::forms::Operation::Value; ///< Alias for the value operation, which is used to evaluate the point-wise value of a field. + using Divergence = discretization::forms::Operation::Divergence; ///< Alias for the divergence operation, which is used to evaluate the point-wise divergence of a vector field. + + /** + * @tparam Term The type of the term to create an operand for, which must be a valid discretization block. + * @tparam Operation The type of the operation to create an operand for, which must + * + * @brief Alias for an @ref Operand of type Operand. This is a convenience alias to avoid having to write out the full type name. See + * my above (likely unparsed in doxygen comment about my bad habit of using deeply nested namespaces and types. I should probably fix that at some point). + */ + template + using Operand = discretization::forms::Operand; + + /** + * @class GradientValue + * + * @brief The point wise contribution to the gravity-gradient relation. Note that these contributions + * are just the point-wise field values. You can get to this by going through the + * weak form derivation. The derivative is just the identity matrix since the residual is just the field value. + * + * @contribution Residuals to the gravity-gradient relation block of the system. + */ + struct GradientValue : Contribution> { + using Operands = utils::types::TypeList>; ///< This equation only depends on the gradient field + + /** + * @brief Evaluate the gravity-gradient relation at a given point. Fundamentally all this this is doing is + * returning the gradient of the gravitational potential since to satisfy the gravity-gradient relation that should be zero. + * The derivative of the residual with respect to the gradient along each axis is then just the identity matrix. + * + * @param point The point at which to evaluate the gravity-gradient relation. + * @return A PointLinearization object containing the residual and field derivatives for the gravity-gradient relation at the given point. + */ + static PointLinearization evaluate( + const EquationPoint& point + ) { + PointLinearization result; + + result.residual = point.fields.get(); + + for (int axis = 0; axis < 3; ++axis) { + result.field_derivatives.get()[axis][axis] = 1.0; + } + + return result; + } + + /** + * @brief Return the additional order of the gravity-gradient relation with respect to the gradient operand. + * + * @param operand-type The type of the operand being evaluated, which must be the gradient operand. + * @param position-order The order of the position being evaluated, which is not used in this implementation. + * @return The additional order of the gravity-gradient relation with respect to the gradient operand, which is always 0. + */ + static int additional_order(Operand /*operand-type*/, int /*position-order*/) { + return 0; + } + }; + + /** + * @class PotentialValue + * + * @brief The point wise contribution to the Poisson balance equation. Note that these contributions + * are just the point-wise potential values. You can get to this by going through the + * weak form derivation. The derivative is just one since the residual is just the variable. + * + * @contribution Residuals to the gravity-gradient relation block of the system. Note how this is declares the divergence operator. Therefore whatever is returned from here + * will be integrated against the divergence of the test function. + */ + struct PotentialValue : Contribution> { + using Operands = utils::types::TypeList>; ///< This equation only depends on the potential field + + /** + * @brief Evaluate the Poisson balance equation at a given point. Fundamentally all this this is doing is + * returning the gravitational potential since to satisfy the Poisson balance equation that should be zero. + * The derivative of the residual with respect to the potential is then just one. + * + * When you go through the weak form derivation and the integration by parts you find that + * the point wise contribution to the gravity gradient row of the system is + * the divergence of the test function times the potential. \f$ \left(\nabla \cdot v\right) \phi \f$. + * The derivative is then trivially one. + * + * @param point The point at which to evaluate the Poisson balance equation. + * @return A PointLinearization object containing the residual and field derivatives for the Poisson balance equation at the given point. + */ + static PointLinearization evaluate(const EquationPoint& point) { + PointLinearization result; + + result.residual = point.fields.get(); + result.field_derivatives.get() = 1.0; + + return result; + } + + /** + * @brief Return the additional order of the Poisson balance equation with respect to the potential operand. The contribution + * is linear in potential so the additional order is always 0. + * + * @param operand-type The type of the operand being evaluated, which must be the potential operand. + * @param position-order The order of the position being evaluated, which is not used in this implementation. + * @return The additional order of the Poisson balance equation with respect to the potential operand, which is always 0. + */ + static int additional_order(Operand /*operand-type*/, int /*position-order*/) { + return 0; + } + }; + + /** + * @class GradientDivergence + * + * @brief The point wise contribution to the Poisson balance equation. Note that these contributions + * are just the point-wise divergence of the gradient field. You can get to this by going through the + * weak form derivation. The derivative is just the identity matrix since the residual is just the divergence of the gradient field. + * + * @contribution Residuals to the gravitational potential block of the system. + */ + struct GradientDivergence : Contribution> { + using Operands = utils::types::TypeList>; ///< This equation only depends on the divergence of the gradient field + + /** + * @brief Evaluate the point-wise residual of the potential. Specifically, this uses the + * gets the divergence of the gradient. Note how directly related this is to the name of this + * class + * + * @param point The point at which to evaluate the Poisson balance equation. + * @return A PointLinearization object containing the residual and field derivatives for the Poisson balance equation at the given point. + */ + static PointLinearization evaluate( + const EquationPoint& point + ) { + PointLinearization result; + + result.residual = point.fields.get>(); + result.field_derivatives.get>() = 1.0; + + return result; + } + + /** + * @brief Return the additional order of the Poisson balance equation with respect to the divergence of the gradient operand. The contribution + * is linear in the divergence of the gradient so the additional order is always 0. + * + * @param operand-type The type of the operand being evaluated, which must be the divergence of the gradient operand. + * @param position-order The order of the position being evaluated, which is not used in this implementation. + * @return The additional order of the Poisson balance equation with respect to the divergence of the gradient operand, which is always 0. + */ + static int additional_order(Operand /*operand-type*/, int /*position-order*/) { + return 0; + } + }; + + /** + * @class DensitySource + * + * @brief The point wise contribution to the Poisson balance equation. Note that these contributions + * are just the point-wise density values multiplied by \f$-4\pi G\f#. You can get to this by going through the + * weak form derivation. The derivative is just \f$-4\pi G\f$. + * + * @contribution Residuals to the gravitational potential block of the system where that block is supported by the stellar domain (e.g. mass is only defined on the stellar domain). + */ + struct DensitySource : Contribution, discretization::domain::StellarDomains> { + using Operands = utils::types::TypeList>; ///< This equation only depends on the density field + double four_pi_G; ///< The coefficient to multiply the density by, which is typically \f$-4\pi G\f$. (we let this vary so that the user can select the value of G they want to use in their simulation, primarily since I am writing this with natural units where G=1 but we may eventually switch internally to cgs) + + /** + * @brief Construct a DensitySource contribution with a given coefficient. + * + * @param coefficient The coefficient to multiply the density by, which is typically \f$-4\pi G\f$. + */ + explicit DensitySource(const double coefficient) : four_pi_G(coefficient) {} + + /** + * @brief Evaluate the point-wise residual of the potential. Specifically, this uses the + * gets the density multiplied by \f$-4\pi G\f$. Note how directly related this is to the name of this + * class + * + * @param point The point at which to evaluate the Poisson balance equation. + * @return A PointLinearization object containing the residual and field derivatives for the Poisson balance equation at the given point. + */ + PointLinearization evaluate( + const EquationPoint& point + ) const { + PointLinearization result; + + result.residual = -four_pi_G * point.fields.get(); + result.field_derivatives.get() = -four_pi_G; + + return result; + } + + /** + * @brief Return the additional order of the Poisson balance equation with respect to the density operand. The contribution + * is linear in density so the additional order is always 0. + * + * @param operand-type The type of the operand being evaluated, which must be the density operand. + * @param position-order The order of the position being evaluated, which is not used in this implementation. + * @return The additional order of the Poisson balance equation with respect to the density operand, which is always 0. + */ + static int additional_order(Operand, int) { + return 0; + } + }; + } // namespace gravity_contributions + + /** + * @class GravityGradientRelation + * + * @brief The gravity-gradient relation equation, which weakly constrains the gravitational potential and gravitational field to agree. + * + * @contribution Residuals to the gravity-gradient relation block of the system. + */ + struct GravityGradientRelation { + using Test = discretization::blocks::gravity::gradient; ///< The test function for the gravity-gradient relation equation, which is the gravitational field. + using ResidualQuantity = dimensions::Acceleration; ///< The physical quantity which is being constrained by this equation, that is the quantity which is being set to a target value. + + static constexpr std::string_view name = "gravity-gradient relation"; ///< Compile time name for the equation, used for debugging and logging. + + /** + * @brief Return the contributions to the gravity-gradient relation equation. + * + * @return A tuple containing the contributions to the gravity-gradient relation equation, which are the point-wise residuals of the gradient field and the divergence of the test function times the potential field. + */ + static auto contributions() { + return std::tuple{gravity_contributions::GradientValue{}, gravity_contributions::PotentialValue{}}; + } + + }; + + /** + * @class PoissonBalance + * + * @brief The Poisson balance equation, which weakly constrains the gravitational potential to satisfy the Poisson equation with the mass density as a source term. + * + * @contribution Residuals to the gravitational potential block of the system where that block is supported by the stellar domain (e.g. mass is only defined on the stellar domain). + */ + class PoissonBalance { + public: + using Test = discretization::blocks::gravity::potential; ///< The test function for the Poisson balance equation, which is the gravitational potential. + using ResidualQuantity = dimensions::InverseTimeSquared; ///< The physical quantity which is being constrained by this equation, that is the quantity which is being set to a target value. + + static constexpr std::string_view name = "Poisson balance"; ///< Compile time name for the equation, used for debugging and logging. + + /** + * @brief Construct a Poisson balance equation with a given coefficient. + * + * @param four_pi_G The coefficient to multiply the density by, which is typically \f$-4\pi G\f$. + * + * @throws std::invalid_argument If the coefficient is not finite or not positive. + */ + explicit PoissonBalance(const double four_pi_G) + : m_four_pi_G(four_pi_G) { + using namespace utils::errors; + serif_ensure(std::isfinite(four_pi_G) && four_pi_G > 0.0, "Poisson balance requires positive finite 4 pi G."); + } + + /** + * @brief Return the contributions to the Poisson balance equation. + * + * @return A tuple containing the contributions to the Poisson balance equation, which are the point-wise residuals of the divergence of the gradient field and the density field multiplied by \f$-4\pi G\f$. + */ + auto contributions() const { + return std::tuple{gravity_contributions::GradientDivergence{}, gravity_contributions::DensitySource{m_four_pi_G}}; + } + + private: + double m_four_pi_G; ///< The coefficient to multiply the density by, which is typically \f$-4\pi G\f$. + }; +} // namespace serif::equations