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.
This commit is contained in:
313
src/include/serif/equations/mixed_gravity.hpp
Normal file
313
src/include/serif/equations/mixed_gravity.hpp
Normal file
@@ -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 <cmath>
|
||||
#include <tuple>
|
||||
#include <string_view>
|
||||
#include <stdexcept>
|
||||
#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<Term, Operation>. 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 <typename Term, typename Operation = Value>
|
||||
using Operand = discretization::forms::Operand<Term, Operation>;
|
||||
|
||||
/**
|
||||
* @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<Operand<Gradient>> {
|
||||
using Operands = utils::types::TypeList<Operand<Gradient>>; ///< 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<Operands, discretization::forms::VectorValue> evaluate(
|
||||
const EquationPoint<Operands>& point
|
||||
) {
|
||||
PointLinearization<Operands, discretization::forms::VectorValue> result;
|
||||
|
||||
result.residual = point.fields.get<Gradient>();
|
||||
|
||||
for (int axis = 0; axis < 3; ++axis) {
|
||||
result.field_derivatives.get<Gradient>()[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<Gradient> /*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<Operand<Gradient, Divergence>> {
|
||||
using Operands = utils::types::TypeList<Operand<Potential>>; ///< 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<Operands> evaluate(const EquationPoint<Operands>& point) {
|
||||
PointLinearization<Operands> result;
|
||||
|
||||
result.residual = point.fields.get<Potential>();
|
||||
result.field_derivatives.get<Potential>() = 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<Potential> /*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<Operand<Potential>> {
|
||||
using Operands = utils::types::TypeList<Operand<Gradient, Divergence>>; ///< 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<Operands> evaluate(
|
||||
const EquationPoint<Operands>& point
|
||||
) {
|
||||
PointLinearization<Operands> result;
|
||||
|
||||
result.residual = point.fields.get<Operand<Gradient, Divergence>>();
|
||||
result.field_derivatives.get<Operand<Gradient, Divergence>>() = 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<Gradient, Divergence> /*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<Operand<Potential>, discretization::domain::StellarDomains> {
|
||||
using Operands = utils::types::TypeList<Operand<Density>>; ///< 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<Operands> evaluate(
|
||||
const EquationPoint<Operands>& point
|
||||
) const {
|
||||
PointLinearization<Operands> result;
|
||||
|
||||
result.residual = -four_pi_G * point.fields.get<Density>();
|
||||
result.field_derivatives.get<Density>() = -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<Density>, 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<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(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
|
||||
Reference in New Issue
Block a user