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

feat(solver/seed): Initial guess system

Browse Source 
the seed system has two jobs, get the lane-emden solution for a reasonable n and project that onto the mesh. Note that the seed system does not solve lane-emden, that is in utils/laneemden (so that it can be reused by other systems). Currently SERiF only solves polytrope so the lane-emden initial guess is generally very good; however, as we move to non polytropic EOS we will still use a projected lane-emden as our initial guess for the starting model
This commit is contained in:
Emily Boudreaux
2026-11-09 13:35:24 -05:00
parent 9f4207bc6a
commit ac2cccda29
5 changed files with 260 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

87
src/lib/serif/solver/seed/fitted.cpp Normal file
View File

@@ -0,0 +1,87 @@
#include "serif/solver/seed/fitted.hpp"
#include <cmath>
namespace serif::solver::seed {
/** @brief Validate a regular source and own its explicit, thermodynamically unchanged radial stretch. */
std::expected<FittedRadialSeed, std::string> fit_radial_seed(const RadialSeed& seed, dimensions::LengthValue radius) {
const auto& samples = seed.solution.samples;
if (samples.size() < 2 || samples.front().xi != 0.0 || samples.front().theta != 1.0 || samples.front().derivative != 0.0 ||
!seed.solution.surface || samples.back().xi != seed.solution.surface->xi || samples.back().theta != 0.0 ||
samples.back().derivative != seed.solution.surface->derivative) {
return std::unexpected("Fitting a radial guess requires a complete regular Lane–Emden profile ending at its finite surface.");
}
double previous_radius = -1.0;
for (const auto& sample : samples) {
if (!std::isfinite(sample.xi) || sample.xi <= previous_radius || !std::isfinite(sample.theta) || sample.theta < 0.0 || !std::isfinite(sample.derivative)) {
return std::unexpected("Fitting a radial guess requires finite nonnegative samples at strictly increasing dimensionless radii.");
}
previous_radius = sample.xi;
}
const auto center = evaluate_radial_seed(seed, dimensions::LengthValue{0.0});
if (!center) {
return std::unexpected(center.error());
}
if (!std::isfinite(radius.value()) || radius.value() <= 0.0 || !std::isfinite(seed.moment_of_inertia) ||
seed.moment_of_inertia <= 0.0 || !std::isfinite(seed.angular_velocity)) {
return std::unexpected("Fitting a radial guess requires a finite positive radius and inertia and a finite angular velocity.");
}
for (const double coordinate : seed.center) {
if (!std::isfinite(coordinate)) {
return std::unexpected("Fitting a radial guess requires a finite physical center.");
}
}
const long double stretch = static_cast<long double>(radius.value()) / seed.scales.radius.value();
// Apply each length factor to its physical quantity: s^5 alone can overflow while I*s^5 remains representable.
const double mass = static_cast<double>(seed.scales.mass.value() * stretch * stretch * stretch);
const double inertia = static_cast<double>(seed.moment_of_inertia * stretch * stretch * stretch * stretch * stretch);
const double rotation = static_cast<double>(seed.angular_velocity / stretch / stretch / stretch / stretch / stretch);
const double potential = static_cast<double>(center->potential.value() * stretch * stretch);
if (!std::isfinite(mass) || mass <= 0.0 || !std::isfinite(inertia) || inertia <= 0.0 ||
!std::isfinite(rotation) || (seed.angular_velocity != 0.0 && rotation == 0.0) || !std::isfinite(potential)) {
return std::unexpected("Fitting the radial guess exceeds representable mass, gravity or rotation scales. Use consistent physical units and a less extreme radius ratio.");
}
return FittedRadialSeed{.source = seed, .radius = radius, .mass = dimensions::MassValue{mass},
.moment_of_inertia = inertia, .angular_velocity = rotation};
}
/** @brief Evaluate unchanged material fields with the gravity generated by their stretched density. */
std::expected<utils::laneemden::PhysicalSample, std::string> evaluate_radial_seed(
const FittedRadialSeed& seed, dimensions::LengthValue radius) {
if (!std::isfinite(radius.value()) || radius.value() < 0.0 || !std::isfinite(seed.radius.value()) || seed.radius.value() <= 0.0 ||
!std::isfinite(seed.mass.value()) || seed.mass.value() <= 0.0 || !std::isfinite(seed.source.scales.radius.value()) || seed.source.scales.radius.value() <= 0.0 ||
!std::isfinite(seed.source.scales.gravitational_constant.value()) || seed.source.scales.gravitational_constant.value() <= 0.0) {
return std::unexpected("A fitted radial guess requires a finite nonnegative evaluation radius and finite positive source radius, fitted radius, mass and G.");
}
if (radius.value() >= seed.radius.value()) {
const long double gravitational_mass = static_cast<long double>(seed.source.scales.gravitational_constant.value()) * seed.mass.value();
const double potential = static_cast<double>(-gravitational_mass / radius.value());
const double gravity = static_cast<double>(gravitational_mass / radius.value() / radius.value());
if (!std::isfinite(potential) || !std::isfinite(gravity)) {
return std::unexpected("The fitted exterior gravity fields overflow in the supplied physical units.");
}
return utils::laneemden::PhysicalSample{
.density = dimensions::DensityValue{0.0}, .pressure = dimensions::PressureValue{0.0},
.enthalpy = dimensions::QuantityValue<dimensions::SpecificEnthalpy>{0.0},
.potential = dimensions::QuantityValue<dimensions::GravitationalPotential>{potential},
.enclosed_mass = seed.mass, .gravity = dimensions::QuantityValue<dimensions::Acceleration>{gravity}
};
}
const long double stretch = static_cast<long double>(seed.radius.value()) / seed.source.scales.radius.value();
const double source_radius = static_cast<double>((static_cast<long double>(radius.value()) / seed.radius.value()) * seed.source.scales.radius.value());
auto sample = evaluate_radial_seed(seed.source, dimensions::LengthValue{source_radius});
if (!sample) {
return sample;
}
const double mass = static_cast<double>(sample->enclosed_mass.value() * stretch * stretch * stretch);
const double potential = static_cast<double>(sample->potential.value() * stretch * stretch);
const double gravity = static_cast<double>(sample->gravity.value() * stretch);
if (!std::isfinite(mass) || !std::isfinite(potential) || !std::isfinite(gravity)) {
return std::unexpected("The fitted radial gravity fields overflow in the supplied physical units.");
}
sample->enclosed_mass = dimensions::MassValue{mass};
sample->potential = dimensions::QuantityValue<dimensions::GravitationalPotential>{potential};
sample->gravity = dimensions::QuantityValue<dimensions::Acceleration>{gravity};
return sample;
}
}

60
src/lib/serif/solver/seed/polytrope.cpp Normal file
View File

@@ -0,0 +1,60 @@
#include "serif/solver/seed/polytrope.hpp"
#include <cmath>
#include <numbers>
#include <utility>
namespace serif::solver::seed {
std::expected<RadialSeed, std::string> make_polytropic_seed(
const eos::models::Polytrope& equation_of_state, dimensions::MassValue mass, PolytropeOptions options) {
const double n = equation_of_state.polytropic_index();
if (!std::isfinite(mass.value()) || mass.value() <= 0.0 || !std::isfinite(options.gravitational_constant.value()) || options.gravitational_constant.value() <= 0.0) {
return std::unexpected("A polytropic seed requires a positive finite target mass and gravitational constant.");
}
auto solution = utils::laneemden::solve(n, options.integration);
if (!solution) {
return std::unexpected(solution.error());
}
if (!solution->surface) {
return std::unexpected("A finite-radius stellar seed requires a Lane-Emden first zero; choose a finite-radius index and sufficient integration extent.");
}
if (!options.central_density) {
if (std::abs(n - 3.0) < 1.0e-12) {
return std::unexpected("An n=3 polytrope's mass is independent of central density. Specify central_density to choose its initial radius.");
}
const auto unit = utils::laneemden::scale_from_polytropic_constant(*solution, equation_of_state.polytropic_constant(),
dimensions::DensityValue{1.0}, options.gravitational_constant);
if (!unit) {
return std::unexpected(unit.error());
}
const double density = std::exp((2.0 * n / (3.0 - n)) * std::log(mass.value() / unit->mass.value()));
options.central_density = dimensions::DensityValue{density};
}
auto scales = utils::laneemden::scale_from_polytropic_constant(*solution, equation_of_state.polytropic_constant(),
*options.central_density, options.gravitational_constant);
if (!scales) {
return std::unexpected(scales.error());
}
// An explicit density defines an initial guess, including at n=3. Its physical mass
// need not equal the model's target; that prescription belongs to the equilibrium equations.
// Simpson integration of the regular dimensionless profile estimates I for the initial rigid rotation.
constexpr int intervals = 2048;
const double step = solution->surface->xi / intervals;
double integral = 0.0;
for (int index = 1; index < intervals; ++index) {
const double xi = index * step;
const auto sample = utils::laneemden::evaluate(*solution, xi);
if (!sample) {
return std::unexpected(sample.error());
}
const double weight = index % 2 == 0 ? 2.0 : 4.0;
integral += weight * std::pow(sample->theta, n) * std::pow(xi, 4);
}
const double inertia = (8.0 * std::numbers::pi / 3.0) * scales->central_density.value() * std::pow(scales->radial_scale, 5) * integral * step / 3.0;
if (!std::isfinite(inertia) || inertia <= 0.0) {
return std::unexpected("The spherical seed's moment of inertia is not finite and positive. Check the physical scales.");
}
return RadialSeed{.solution = std::move(*solution), .scales = *scales, .moment_of_inertia = inertia};
}
}