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

feat(equations): added central_density

Browse Source 
This equation provides a constraint on central density through specific enthalpy (I think this is general so long as the eos can convert, though I wrote this for the barotropic case). This may need to be changed latter to be more general (though also generally we should not constrain central density). this is really there to deal with the specific non-rotating n=3 case where the central density is ill-defined due to a singularity.
This commit is contained in:
Emily Boudreaux
2026-10-08 10:44:52 -04:00
parent ba6264e22c
commit 0571c63fad
1 changed files with 135 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

135
src/include/serif/equations/central_density.hpp Normal file
View File

@@ -0,0 +1,135 @@
#pragma once
/**
* @file central_density.hpp
* @author Emily M. Boudreaux
*
* @brief Prescribe central density through the continuous enthalpy: h(0) - h_EOS(rho_c) = 0. Generally this is not needed; however, in the particular
* case of an n=3 polytrope with no rotation the central density becomes ill-defined and must be prescribed to avoid a singularity in the solution.
*/
#include <array>
#include <cmath>
#include <string_view>
#include <utility>
#include "serif/dimensions/type_alias.hpp"
#include "serif/eos/evaluation.hpp"
#include "serif/equations/contributions.hpp"
#include "serif/equations/point.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::equations {
/// Prescribe central density through the continuous enthalpy: h(0) - h_EOS(rho_c) = 0.
/**
* @brief Prescribe central density through the continuous enthalpy: \f$h(0) - h_{EOS}(\rho_c) = 0\f$. Generally this is not needed; however, in the particular
* case of an n=3 polytrope with no rotation the central density becomes ill-defined and must be prescribed to avoid a singularity in the solution. Note that
* this equation is not responsible for sorting out what point is the center of the star; rather, this is provided to the
* evaluate method.
*
* @tparam EOST The EOS model type to use for evaluating the enthalpy from the prescribed density.
*
* @requires The EOS model must support the relation SpecificEnthalpyFromDensity
*
* @contribution Residuals to the specific enthalpy row of the system
*
* @throws std::invalid_argument If the prescribed density is not finite and strictly positive.
* @throws std::domain_error If the EOS returns a non-finite enthalpy for the prescribed density.
*
*/
template <eos::models::EOSModel EOST>
requires eos::models::EOSSupportingRelation<EOST, eos::relations::SpecificEnthalpyFromDensity>
class CentralDensityConstraint final : public Contribution<discretization::forms::Operand<discretization::blocks::enthalpy::specific>> {
public:
using Operands = utils::types::TypeList<TestOperand>; ///< The operand types of the equation, that is the types of the fields which are used to evaluate the equation.
using Border = discretization::blocks::central_density::border; ///< The border type of the equation, that is what row of the global system this equation contributes to.
using ResidualQuantity = dimensions::SpecificEnthalpy; ///< 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 = "central density constraint"; ///< Compile time name for the equation, used for debugging and logging.
static constexpr bool position_independent = true; ///< Whether the equation is independent of the position of the point being evaluated, that is whether the equation is invariant under translations and rotations of the point. This is only used for compile-time checks
static constexpr bool affine = true; ///< Whether the equation is affine, that is whether the equation is linear in the fields being evaluated. This is only used for compile-time checks
static constexpr std::array<double, 3> reference_position{}; ///< The reference position of the equation, this is only used for compile-time checks and is not used in the evaluation of the equation.
/// Own the EOS and resolve the strictly positive prescribed density once.
/**
* @brief Construct a central density constraint.
*
* @param equation_of_state The EOS model to use for evaluating the enthalpy from the prescribed density.
* @param target_density The prescribed central density, which must be finite and strictly positive.
*
* @throws std::invalid_argument If the prescribed density is not finite and strictly positive.
* @throws std::domain_error If the EOS returns a non-finite enthalpy for the prescribed density.
*
* @todo We should consider whether we actually need to own the EOS, this choice made lifetime management easier but introduces copy overhead. Moreover,
* it means that the EOS must be copyable and, more importantly in my view, that the caller must know that the EOS is being
* copied.
*/
CentralDensityConstraint(EOST equation_of_state, dimensions::DensityValue target_density) :
m_eos(std::move(equation_of_state)),
m_target_density(target_density) {
using namespace utils::errors;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(std::isfinite(target_density.value()) && target_density.value() > 0.0, "Central density must be finite and strictly positive.");
m_target_enthalpy = eos::evaluate<dimensions::SpecificEnthalpy>(m_eos, target_density).value();
serif_ensure<SERiFErrorCode::NON_FINITE_VALUE_ERROR, std::domain_error>(std::isfinite(m_target_enthalpy), "The EOS returned a non-finite central enthalpy.");
}
/// The border backend places this scalar condition in its global phase row.
/**
* @brief Evaluate the central density constraint at a given point.
*
* @param point The point at which to evaluate the equation, which must be the center of the star.
*
* @note It is the responsibility of the caller (and generally of the border equation backend) to ensure that this residual
* is applied to the correct row of the global system, the row defined by the ::Border type alias.
*
* @return The linearization of the equation at the given point, which contains the residual and the derivatives of the residual with respect to the fields being evaluated.
*/
[[nodiscard]] PointLinearization<Operands> evaluate(const EquationPoint<Operands>& point) const {
PointLinearization<Operands> result;
result.residual = point.fields.get<Test>() - m_target_enthalpy;
result.field_derivatives.get<Test>() = 1.0;
return result;
}
/**
* @brief Return the prescribed central density.
*
* @return The prescribed central density.
*/
[[nodiscard]] dimensions::DensityValue target_density() const noexcept {
return m_target_density;
}
/**
* @brief Return the target enthalpy corresponding to the prescribed central density.
*
* @return The target enthalpy corresponding to the prescribed central density.
*/
[[nodiscard]] dimensions::SpecificEnthalpyValue target_enthalpy() const noexcept {
return dimensions::SpecificEnthalpyValue{m_target_enthalpy};
}
/**
* @brief Convert some enthalpy to a density for physical diagnostics, when the inverse relation exists.
* @param enthalpy The specific enthalpy to convert.
*
* @return The corresponding density.
*
* @requires The EOS model must support the relation DensityFromSpecificEnthalpy, otherwise this method is not available.
*
* @note Generally this method is only for debugging, if you want to use the density you should really call the EOS directly.
*/
[[nodiscard]] dimensions::DensityValue density_from_enthalpy(dimensions::SpecificEnthalpyValue enthalpy) const
requires eos::models::EOSSupportingRelation<EOST, eos::relations::DensityFromSpecificEnthalpy> {
return eos::evaluate<dimensions::MassDensity>(m_eos, enthalpy);
}
private:
EOST m_eos; ///< The EOS model to use for evaluating the enthalpy from the prescribed density.
dimensions::DensityValue m_target_density; ///< The prescribed central density, which must be finite and strictly positive.
double m_target_enthalpy; ///< The target enthalpy corresponding to the prescribed central density, which is computed from the EOS and must be finite.
};
} // namespace serif::equations