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

Began porting the mapping system

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-20 14:37:11 -04:00
parent f69fac8f1e
commit 08a8f231d2
46 changed files with 1336 additions and 115 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
src/lib/serif/eos/models/polytropic.cpp
View File

@@ -104,24 +104,24 @@ namespace serif::eos::models {
};
}
[[nodiscard]] PartialDerivative<Density, SpecificEnthalpy> Polytrope::partial_derivative(DensityFromSpecificEnthalpy, WithRespectTo<SpecificEnthalpy>, const SpecificEnthalpyValue specific_enthalpy) const {
[[nodiscard]] PartialDerivative<MassDensity, SpecificEnthalpy> Polytrope::partial_derivative(DensityFromSpecificEnthalpy, WithRespectTo<SpecificEnthalpy>, const SpecificEnthalpyValue specific_enthalpy) const {
validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(specific_enthalpy.value());
if (specific_enthalpy.value() < 0.0) {
// dρ/dh = 0
return PartialDerivative<Density, SpecificEnthalpy>{0.0};
return PartialDerivative<MassDensity, SpecificEnthalpy>{0.0};
}
if (specific_enthalpy.value() == 0.0) {
// dρ/dh = n == 1 -> 1/hs
// dρ/dh = else -> 0
return PartialDerivative<Density, SpecificEnthalpy>{
return PartialDerivative<MassDensity, SpecificEnthalpy>{
m_polytropic_index == 1.0 ? 1.0 / m_enthalpy_scale : 0.0
};
}
// dρ/dh = (n/hs) * (h/hs)^(n-1)
return PartialDerivative<Density, SpecificEnthalpy>{
return PartialDerivative<MassDensity, SpecificEnthalpy>{
m_polytropic_index / m_enthalpy_scale * std::pow(specific_enthalpy.value() / m_enthalpy_scale, m_polytropic_index - 1.0)
};
}
@@ -133,16 +133,16 @@ namespace serif::eos::models {
return PartialDerivative<Pressure, SpecificEnthalpy>{density.value()};
}
[[nodiscard]] PartialDerivative<Pressure, Density> Polytrope::partial_derivative(PressureFromDensity, WithRespectTo<Density>, const DensityValue density) const {
[[nodiscard]] PartialDerivative<Pressure, MassDensity> Polytrope::partial_derivative(PressureFromDensity, WithRespectTo<MassDensity>, const DensityValue density) const {
validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(density.value());
if (density.value() == 0.0) {
// dP/dρ = 0
return PartialDerivative<Pressure, Density>{0.0};
return PartialDerivative<Pressure, MassDensity>{0.0};
}
// dP/dρ = K * (1 + 1/n) * (ρ)^(1/n)
return PartialDerivative<Pressure, Density>{
return PartialDerivative<Pressure, MassDensity>{
m_polytropic_constant * (1.0 + 1.0 / m_polytropic_index) * std::pow(density.value(), 1.0 / m_polytropic_index)
};
}

269
src/lib/serif/mapping/compactification/kelvin.cpp Normal file
View File

@@ -0,0 +1,269 @@
#include <cmath>
#include "serif/utils/checks/finite.hpp"
#include "serif/mapping/compactification/kelvin.hpp"
#include "serif/utils/user_options/mapping/compactification/kelvin.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::mapping::compactification {
using namespace serif::utils::user_options;
using namespace serif::utils::errors;
using namespace serif::utils::checks;
KelvinCompactification::KelvinCompactification(const KelvinCompactificationOptions& options) : m_options(options) {
if (!std::isfinite(m_options.r_star_ref)) {
serif_error<std::invalid_argument, SERiFErrorCode::NON_FINITE_MAPPING_INPUT>("Kelvin compactification requires a finite reference stellar radius. Received {} ", m_options.r_star_ref);
}
if (!std::isfinite(m_options.r_inf_ref)) {
serif_error<std::invalid_argument, SERiFErrorCode::NON_FINITE_MAPPING_INPUT>("Kelvin compactification requires a finite reference infinity radius. Received {} ", m_options.r_star_ref);
}
if (m_options.r_star_ref <= 0 || m_options.r_inf_ref < m_options.r_star_ref)
{
serif_error<std::invalid_argument, SERiFErrorCode::IMPROPERLY_ORDERED_REFERENCE_RADII>("Kelvin compactification requires 0 < r_star_ref < r_inf_ref. Got r_star_ref={} and r_inf_ref={}", m_options.r_star_ref, m_options.r_inf_ref);
}
if (!std::isfinite(m_options.coordinate_tolerance) || m_options.coordinate_tolerance < 0.0 || m_options.coordinate_tolerance > 1.0)
{
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COORDINATE_TOLERANCE>("Kelvin compactification coordinate tolerance must be finite and on [0, 1). Received: {}", m_options.coordinate_tolerance);
}
}
bool KelvinCompactification::is_affine_in_displacement() const noexcept {
return true; // Always true for kelvin compactification
}
MappingStatus KelvinCompactification::evaluate(const ExteriorMapInput& input, ExteriorMapResult& result) const {
const int dimension = input.reference_position.Size();
if (dimension <= 0 || input.displaced_position.Size() != dimension || input.compactification_coordinate_gradient.Size() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (input.displacement_jacobian.Height() != dimension || input.displacement_jacobian.Width() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
const bool reference_pos_is_finite = is_finite(input.reference_position);
const bool displaced_pos_is_finite = is_finite(input.displaced_position);
const bool coord_grad_is_finite = is_finite(input.compactification_coordinate_gradient);
const bool displacement_jacobian_is_finite = is_finite(input.displacement_jacobian);
if (!(reference_pos_is_finite && displaced_pos_is_finite && coord_grad_is_finite && displacement_jacobian_is_finite)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
RadialFactors factors{};
// The key radial factors we use are the scale (which stretches the finite computational domain into the infinite physical domain) and the scale derivative (which is used to compute the mapping jacobian)
const MappingStatus factor_status = compute_radial_factors(input.compactification_coordinate, factors);
if (factor_status != MappingStatus::VALID) {
return factor_status;
}
result.physical_position.SetSize(dimension);
result.mapping_jacobian.SetSize(dimension, dimension);
for (int i = 0; i < dimension; ++i) {
// Note how the physical position is just the product of the displaced position and the scale factor.
result.physical_position(i) = factors.scale * input.displaced_position(i);
// The mapping jacobian comes from trivial application of the product rule
// recall: r_i = scale * x_i where r is the physical position and x is the displaced position.
// then we can differentiate wrt. X_j holding nothing fixed. Note the capital X here, this is the mesh coordinate not the displaced position.
// Lets call this jacobian F
// F = \frac{\partial r_i}{\partial X_{j}}
// F then tells us how the physical position changes as we move along mesh coordinates
// Lets then apply this to the function we have for the kelvin compactification
// F = scale * \frac{\partial x_i}{\partial X_j} + x_i * \frac{\partial scale}{\partial X_j}
// Below you can see the displacement jacobian (\frac{\partial x_i}{\partial X_j}) and the scale derivative (\frac{\partial scale}{\partial X_j}) being applied to compute the mapping jacobian.
for (int j = 0; j < dimension; ++j) {
const double scale_gradient = factors.scale_derivative * input.compactification_coordinate_gradient(j);
result.mapping_jacobian(i, j) = factors.scale * input.displacement_jacobian(i, j) + input.displaced_position(i) * scale_gradient;
}
}
if (!is_finite(result.physical_position) || !is_finite(result.mapping_jacobian)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
const double mapping_determinant = result.mapping_jacobian.Det();
if (!std::isfinite(mapping_determinant)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
if (mapping_determinant <= 0.0) {
return MappingStatus::NON_POSITIVE_DETERMINATE;
}
return MappingStatus::VALID;
}
MappingStatus KelvinCompactification::evaluate_variation(const ExteriorMapInput& input, const ExteriorMapResult& result, const ExteriorMapDirection& direction, ExteriorMapVariation& variation) const {
const int dimension = input.reference_position.Size();
if (dimension <= 0 || input.displaced_position.Size() != dimension || input.compactification_coordinate_gradient.Size() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (input.displacement_jacobian.Height() != dimension || input.displacement_jacobian.Width() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (result.physical_position.Size() != dimension || result.mapping_jacobian.Height() != dimension || result.mapping_jacobian.Width() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (direction.displaced_position_variation.Size() != dimension || direction.displacement_jacobian_variation.Height() != dimension || direction.displacement_jacobian_variation.Width() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
const bool reference_position_is_finite = is_finite(input.reference_position);
const bool displaced_position_is_finite = is_finite(input.displaced_position);
const bool compactification_coord_grad_is_finite = is_finite(input.compactification_coordinate_gradient);
const bool displacement_jacobian_is_finite = is_finite(input.displacement_jacobian);
if (!(reference_position_is_finite && displaced_position_is_finite && compactification_coord_grad_is_finite && displacement_jacobian_is_finite)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
const bool result_physical_position_is_finite = is_finite(result.physical_position);
const bool result_mapping_jacobian_is_finite = is_finite(result.mapping_jacobian);
const bool direction_displaced_position_is_finite = is_finite(direction.displaced_position_variation);
const bool direction_displacement_jacobian_is_finite = is_finite(direction.displacement_jacobian_variation);
if (!(result_physical_position_is_finite && result_mapping_jacobian_is_finite && direction_displaced_position_is_finite && direction_displacement_jacobian_is_finite)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT; // Note we are checking a parameter called result here; however, this is the result from an evaluate call, it is an input to this function.
}
RadialFactors factors{};
const MappingStatus factor_status = compute_radial_factors(input.compactification_coordinate, factors);
if (factor_status != MappingStatus::VALID) {
return factor_status;
}
variation.physical_position_variation.SetSize(dimension);
variation.mapping_jacobian_variation.SetSize(dimension, dimension);
for (int i = 0; i < dimension; ++i) {
variation.physical_position_variation(i) = factors.scale * direction.displaced_position_variation(i);
for (int j = 0; j < dimension; ++j) {
const double scale_gradient = factors.scale_derivative * input.compactification_coordinate_gradient(j);
variation.mapping_jacobian_variation(i, j) = factors.scale * direction.displacement_jacobian_variation(i, j) + direction.displaced_position_variation(i) * scale_gradient;
}
}
if (!is_finite(variation.physical_position_variation) || !is_finite(variation.mapping_jacobian_variation)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
[[nodiscard]] std::string_view KelvinCompactification::get_name() const noexcept {
return m_name;
}
[[nodiscard]] double KelvinCompactification::get_reference_stellar_radius() const noexcept {
return m_options.r_star_ref;
}
[[nodiscard]] double KelvinCompactification::get_reference_infinity_radius() const noexcept {
return m_options.r_inf_ref;
}
[[nodiscard]] double KelvinCompactification::get_coordinate_tolerance() const noexcept {
return m_options.coordinate_tolerance;
}
[[nodiscard]] MappingStatus KelvinCompactification::evaluate_jacobian_variation_factors(double compactification_coordinate, const mfem::Vector& coordinateGradient, double& scale, mfem::Vector& scaleGradient) const {
if (coordinateGradient.Size() <= 0) {
return MappingStatus::INVALID_DIMENSION;
}
if (!is_finite(coordinateGradient)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
RadialFactors factors{};
const MappingStatus status = compute_radial_factors(compactification_coordinate, factors);
if (status != MappingStatus::VALID) {
return status;
}
scale = factors.scale;
scaleGradient = coordinateGradient;
scaleGradient *= factors.scale_derivative;
if (!is_finite(scaleGradient)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
MappingStatus KelvinCompactification::compute_radial_factors(double compactification_coordinate, RadialFactors& radial_factors) const {
if (!std::isfinite(compactification_coordinate)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
const double tolerance = m_options.coordinate_tolerance;
if (compactification_coordinate < -tolerance || compactification_coordinate > 1.0 + tolerance) {
return MappingStatus::OUTSIDE_REFERENCE_DOMAIN;
}
double coordinate = compactification_coordinate;
if (coordinate < 0.0) {
coordinate = 0.0; // Not sure this is the best behavior here...
}
if (coordinate >= 1.0 - tolerance) {
return MappingStatus::AT_COMPACTIFIED_INFINITY;
}
// Here we need to do some transformations from the options defined on the mesh to useful computational coordinates
// r_inf_ref is the computational / reference radius of the infinity surface (the edge of the entire domain) and r_star_ref is the radius of the spherical
// stellar model inscribed within. Therefore radial extent is the computational radial distance between the stellar surface and the infinity surface.
// Note that this is separate from the parameterize compactification coordinate.
const double radial_extent = m_options.r_inf_ref - m_options.r_star_ref;
// This places us at the correct spot in computational space given the current compactification coordinate. Say you have compactification = 0.5,
// an r_star_ref of 2 and a radial extent of 3, this this will place you at 2 + 0.5 * 3 = 3.5 in computational space, which is half way between the stellar surface and the infinity surface.
const double computational_radius = m_options.r_star_ref + coordinate * radial_extent;
if (!std::isfinite(computational_radius) || computational_radius <= 0.0) {
return MappingStatus::INVALID_REFERENCE_RADIUS;
}
// Invert the exterior coordinate so it runs from 0 at the star to 1 at compactified infinity
const double one_minus_coordinate = 1.0 - coordinate;
const double denominator = computational_radius * one_minus_coordinate;
if (!std::isfinite(denominator) || denominator <= 0.0) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
// The scale here is the factor which stretches the finite computational domain into the infinite physical domain. Properties we need this to have
// include that it should go to 1 at the stellar surface and go to infinity at the compactified infinity.
// Mathematically this is scale = |r|/|x| where r is the physical radius and x is the computational radius.
// Put another way, scale is the ratio of the target physical radius for the current compactification coordinate
// to the current mesh radius.
const double scale = m_options.r_star_ref / denominator;
const double scale_derivative = scale * (1.0 / one_minus_coordinate - radial_extent / computational_radius);
if (!std::isfinite(scale) || !std::isfinite(scale_derivative)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
radial_factors.compactification_coordinate = coordinate;
radial_factors.computational_radius = computational_radius;
radial_factors.scale = scale;
radial_factors.scale_derivative = scale_derivative;
return MappingStatus::VALID;
}
}