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

feat(solver/continuation): system to auto continue a solver along some number of parameters

Browse Source 
A representative example is angular momentum. It may be useful to not spin a model up all the way initially (due to convergence issues) rather the continuation system allows for a model to be spun up over a few solves so that it can more reliably reach a stable state
This commit is contained in:
Emily Boudreaux
2026-11-09 13:24:42 -05:00
parent 87a1026be6
commit 64cc306a2c
10 changed files with 1340 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/include/serif/solver/continuation/continuation.hpp Normal file
View File

@@ -0,0 +1,6 @@
#pragma once
#include "serif/solver/continuation/family.hpp"
#include "serif/solver/continuation/identified_equations.hpp"
#include "serif/solver/continuation/polytropic.hpp"
#include "serif/solver/continuation/driver.hpp"

255
src/include/serif/solver/continuation/controller.hpp Normal file
View File

@@ -0,0 +1,255 @@
#pragma once
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <limits>
#include <optional>
#include "serif/utils/error/errors.hpp"
namespace serif::solver::continuation {
enum class CoordinateOrder { automatic, preferred, strict };
enum class AttemptOutcome { success, recoverable_rejection, fatal_failure };
enum class ControllerStatus { ready, pending, complete, attempt_limit, no_progress, fatal_failure };
/** @brief Numerical increments are fractions of each coordinate's complete physical span. */
template <std::size_t Count>
struct CoordinateOptions final {
static_assert(Count > 0, "Continuation requires at least one physical coordinate.");
std::array<double, Count> initial_step_fraction = filled(0.25);
std::array<double, Count> minimum_step_fraction = filled(1.0e-6);
std::array<double, Count> maximum_step_fraction = filled(1.0);
std::array<std::size_t, Count> axis_order = []() {
std::array<std::size_t, Count> result{};
for (std::size_t axis = 0; axis < Count; ++axis) {
result[axis] = axis;
}
return result;
}();
CoordinateOrder order{CoordinateOrder::automatic};
double growth_factor{1.5};
double shrink_factor{0.5};
std::uint64_t maximum_attempts{128};
std::uint64_t maximum_rejections_per_axis{16};
/** @brief Initialize every coordinate with the same numerical control. */
[[nodiscard]] static constexpr std::array<double, Count> filled(const double value) {
std::array<double, Count> result{};
result.fill(value);
return result;
}
};
template <std::size_t Count>
struct CoordinateProposal final {
std::uint64_t id{};
std::size_t axis{};
std::array<double, Count> coordinates{};
double step_fraction{};
};
struct CoordinateCounters final {
std::uint64_t attempts{};
std::uint64_t successes{};
std::uint64_t rejections{};
std::uint64_t fatal_failures{};
std::uint64_t work{};
};
/** @brief Own only the continuation path; a caller promotes physical states after successful correction. */
template <std::size_t Count>
class CoordinateController final {
public:
using Coordinates = std::array<double, Count>;
using Proposal = CoordinateProposal<Count>;
/** @brief Validate finite endpoints and controls before any proposal can consume work. */
CoordinateController(Coordinates start, Coordinates target, CoordinateOptions<Count> options = {}) :
m_start(start), m_target(target), m_options(options), m_accepted(start), m_step(options.initial_step_fraction) {
validate_options();
for (std::size_t axis = 0; axis < Count; ++axis) {
m_span[axis] = target[axis] - start[axis];
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
std::isfinite(start[axis]) && std::isfinite(target[axis]) && std::isfinite(m_span[axis]),
"Continuation endpoints and their differences must be finite; invalid coordinate {}.", axis);
}
if (m_accepted == m_target) {
m_status = ControllerStatus::complete;
}
}
/** @brief Read the original physical coordinates. */
[[nodiscard]] const Coordinates& start() const noexcept { return m_start; }
/** @brief Read the exact requested endpoint. */
[[nodiscard]] const Coordinates& target() const noexcept { return m_target; }
/** @brief Read only successfully corrected coordinates. */
[[nodiscard]] const Coordinates& accepted() const noexcept { return m_accepted; }
/** @brief Inspect the selected step controls. */
[[nodiscard]] const CoordinateOptions<Count>& options() const noexcept { return m_options; }
/** @brief Include failed attempts in work accounting. */
[[nodiscard]] const CoordinateCounters& counters() const noexcept { return m_counters; }
/** @brief Explain whether another proposal may be requested. */
[[nodiscard]] ControllerStatus status() const noexcept { return m_status; }
/** @brief Inspect the outstanding proposal without accepting its coordinates. */
[[nodiscard]] const std::optional<Proposal>& pending() const noexcept { return m_pending; }
/** @brief Propose one bounded representable change without modifying accepted coordinates. */
[[nodiscard]] std::optional<Proposal> propose() {
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT>(!m_pending,
"Resolve the outstanding continuation proposal before requesting another.");
if (m_status != ControllerStatus::ready) {
return std::nullopt;
}
if (m_counters.attempts == m_options.maximum_attempts) {
m_status = ControllerStatus::attempt_limit;
return std::nullopt;
}
for (std::size_t attempt = 0; attempt < Count; ++attempt) {
const std::optional<std::size_t> selected = select_axis();
if (!selected) {
break;
}
const std::size_t axis = *selected;
Coordinates coordinates = m_accepted;
const double distance = m_step[axis] * std::abs(m_span[axis]);
if (distance >= std::abs(m_target[axis] - m_accepted[axis])) {
coordinates[axis] = m_target[axis];
} else {
coordinates[axis] += std::copysign(distance, m_span[axis]);
}
const double actual = coordinates[axis] - m_accepted[axis];
if (!std::isfinite(coordinates[axis]) || actual == 0.0 || std::signbit(actual) != std::signbit(m_span[axis])) {
m_blocked[axis] = true;
continue;
}
if ((m_span[axis] > 0.0 && coordinates[axis] > m_target[axis]) || (m_span[axis] < 0.0 && coordinates[axis] < m_target[axis])) {
coordinates[axis] = m_target[axis];
}
++m_counters.attempts;
m_pending = Proposal{.id = m_counters.attempts, .axis = axis, .coordinates = coordinates,
.step_fraction = std::abs((coordinates[axis] - m_accepted[axis]) / m_span[axis])};
m_status = ControllerStatus::pending;
return m_pending;
}
m_status = ControllerStatus::no_progress;
return std::nullopt;
}
/** @brief Validate a resolution before the caller commits any corresponding physical state. */
void validate_resolution(std::uint64_t id, AttemptOutcome outcome, std::uint64_t work = 0) const {
using namespace utils::errors;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(m_pending && m_status == ControllerStatus::pending && m_pending->id == id,
"Continuation resolution must identify its outstanding proposal.");
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
outcome == AttemptOutcome::success || outcome == AttemptOutcome::recoverable_rejection || outcome == AttemptOutcome::fatal_failure,
"Unknown continuation attempt outcome.");
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::overflow_error>(work <= std::numeric_limits<std::uint64_t>::max() - m_counters.work,
"Continuation work counter would overflow.");
}
/** @brief Resolve exactly one proposal; clipped rejected increments shrink from the actual distance attempted. */
void resolve(std::uint64_t id, AttemptOutcome outcome, std::uint64_t work = 0) {
validate_resolution(id, outcome, work);
const std::size_t axis = m_pending->axis;
const double actual_fraction = m_pending->step_fraction;
m_counters.work += work;
if (outcome == AttemptOutcome::success) {
m_accepted = m_pending->coordinates;
++m_counters.successes;
m_blocked.fill(false);
m_rejections_at_base.fill(0);
if (m_step[axis] >= m_options.maximum_step_fraction[axis] / m_options.growth_factor) {
m_step[axis] = m_options.maximum_step_fraction[axis];
} else {
m_step[axis] *= m_options.growth_factor;
}
m_status = m_accepted == m_target ? ControllerStatus::complete : ControllerStatus::ready;
} else if (outcome == AttemptOutcome::fatal_failure) {
++m_counters.fatal_failures;
m_status = ControllerStatus::fatal_failure;
} else {
++m_counters.rejections;
++m_rejections_at_base[axis];
const double previous = std::min(m_step[axis], actual_fraction);
const double next = std::max(m_options.minimum_step_fraction[axis], previous * m_options.shrink_factor);
if (m_rejections_at_base[axis] >= m_options.maximum_rejections_per_axis || actual_fraction <= m_options.minimum_step_fraction[axis] || next >= previous) {
m_blocked[axis] = true;
} else {
m_step[axis] = next;
}
m_status = ControllerStatus::ready;
}
m_pending.reset();
if (m_status == ControllerStatus::ready) {
if (m_counters.attempts == m_options.maximum_attempts) {
m_status = ControllerStatus::attempt_limit;
} else if (!select_axis()) {
m_status = ControllerStatus::no_progress;
}
}
}
private:
/** @brief Reject invalid adaptation, ordering, and work limits as configuration errors. */
void validate_options() const {
using namespace utils::errors;
const auto& options = m_options;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
options.order == CoordinateOrder::automatic || options.order == CoordinateOrder::preferred || options.order == CoordinateOrder::strict,
"Unknown continuation coordinate ordering.");
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
std::isfinite(options.growth_factor) && options.growth_factor >= 1.0 && std::isfinite(options.shrink_factor) &&
options.shrink_factor > 0.0 && options.shrink_factor < 1.0 && options.maximum_attempts > 0 && options.maximum_rejections_per_axis > 0,
"Continuation growth must be at least one, shrinkage between zero and one, and work limits positive.");
std::array<bool, Count> seen{};
for (std::size_t axis = 0; axis < Count; ++axis) {
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(options.axis_order[axis] < Count && !seen[options.axis_order[axis]],
"Continuation axis ordering must be a permutation of the selected coordinates.");
seen[options.axis_order[axis]] = true;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
std::isfinite(options.minimum_step_fraction[axis]) && std::isfinite(options.initial_step_fraction[axis]) &&
std::isfinite(options.maximum_step_fraction[axis]) && options.minimum_step_fraction[axis] > 0.0 &&
options.minimum_step_fraction[axis] <= options.initial_step_fraction[axis] &&
options.initial_step_fraction[axis] <= options.maximum_step_fraction[axis] && options.maximum_step_fraction[axis] <= 1.0,
"Continuation fractions must satisfy 0 < minimum <= initial <= maximum <= 1, including coordinate {}.", axis);
}
}
/** @brief Prefer incomplete coordinates by policy, allowing blocked preferred axes to yield to another axis. */
[[nodiscard]] std::optional<std::size_t> select_axis() const noexcept {
std::optional<std::size_t> selected;
double largest = -1.0;
for (const std::size_t axis : m_options.axis_order) {
if (m_accepted[axis] == m_target[axis]) {
continue;
}
if (m_options.order == CoordinateOrder::strict) {
return m_blocked[axis] ? std::nullopt : std::optional<std::size_t>{axis};
}
if (m_blocked[axis]) {
continue;
}
if (m_options.order == CoordinateOrder::preferred) {
return axis;
}
const double remaining = std::abs((m_target[axis] - m_accepted[axis]) / m_span[axis]);
if (remaining > largest) {
largest = remaining;
selected = axis;
}
}
return selected;
}
Coordinates m_start, m_target;
CoordinateOptions<Count> m_options;
Coordinates m_accepted{}, m_span{}, m_step{};
std::array<bool, Count> m_blocked{};
std::array<std::uint64_t, Count> m_rejections_at_base{};
CoordinateCounters m_counters{};
ControllerStatus m_status{ControllerStatus::ready};
std::optional<Proposal> m_pending;
};
}

232
src/include/serif/solver/continuation/driver.hpp Normal file
View File

@@ -0,0 +1,232 @@
#pragma once
#include <memory>
#include <vector>
#include "serif/solver/continuation/options.hpp"
namespace serif::solver::continuation {
enum class Status { target_reached, starting_star_failed, attempt_limit, work_limit, no_progress, invalid_model };
enum class StageDisposition { converged, not_converged, invalid_model };
struct CorrectionControls final {
double absolute_tolerance{};
std::uint64_t maximum_work{};
};
template <std::size_t Count>
struct StageRecord final {
std::array<double, Count> coordinates{};
std::optional<std::size_t> changed_axis;
bool starting_star{false};
StageDisposition disposition{StageDisposition::not_converged};
std::uint64_t work{};
std::optional<double> residual_norm;
std::string message;
};
/** @brief Retain the last converged physical state separately from the requested endpoint and failed candidates. */
template <typename Report, std::size_t Count>
struct Result final {
Status status{Status::starting_star_failed};
std::string message;
std::array<double, Count> target{};
std::optional<std::array<double, Count>> accepted_coordinates;
std::unique_ptr<Report> accepted;
std::unique_ptr<Report> failed;
std::vector<StageRecord<Count>> stages;
CoordinateCounters counters{};
std::uint64_t starting_attempts{};
std::uint64_t starting_work{};
std::uint64_t omitted_records{};
double absolute_tolerance{};
/** @brief A converged intermediate star does not certify the requested target. */
[[nodiscard]] bool converged() const noexcept { return status == Status::target_reached; }
/** @brief Inspect the owning last successful solve even when the target was not reached. */
[[nodiscard]] const Report* last_converged() const noexcept { return accepted.get(); }
/** @brief Count starting qualification and every continuation correction, including failed attempts. */
[[nodiscard]] std::uint64_t total_attempts() const noexcept { return starting_attempts + counters.attempts; }
/** @brief Include all recorded nonlinear work rather than only successful stages. */
[[nodiscard]] std::uint64_t total_work() const noexcept { return starting_work + counters.work; }
};
/** @brief A supplied physical anchor still requires fresh qualification at the continuation tolerance. */
template <typename Model, typename State>
struct StartingStar final {
Model model;
State state;
};
namespace detail {
/** @brief Try bounded simpler starting models, unless the user explicitly fixed the starting coordinates. */
template <typename Family>
[[nodiscard]] std::vector<typename Family::Coordinates> starting_candidates(const Family& family, const Options<Family>& options) {
std::vector<typename Family::Coordinates> result{options.start};
if (!options.explicit_start) {
const auto target = family.target_coordinates();
for (std::size_t axis = 0; axis < Family::coordinate_count; ++axis) {
auto candidate = options.start;
candidate[axis] = target[axis];
if (std::ranges::find(result, candidate) == result.end()) {
result.push_back(candidate);
}
}
if (std::ranges::find(result, target) == result.end()) {
result.push_back(target);
}
}
return result;
}
/** @brief Bound diagnostic storage while retaining total work and the newest stage. */
template <typename Report, std::size_t Count>
void record(Result<Report, Count>& result, StageRecord<Count> stage, std::size_t maximum) {
if (maximum == 0) {
++result.omitted_records;
} else if (result.stages.size() < maximum) {
result.stages.push_back(std::move(stage));
} else {
result.stages.back() = std::move(stage);
++result.omitted_records;
}
}
}
/**
* @brief Follow a physical family using owning correction results and physical state transfer.
* @param correct A callable accepting a fresh model, a physical state value, and CorrectionControls. Its result must expose
* converged(), has_state(), state(), diagnostics().attempted_iterations/final_residual_norm, and message().
* @note The corrector must honor the supplied work cap and one fixed positive absolute tolerance. Failed states are never
* transferred into retries. Numerical nonconvergence is recoverable; malformed family models stop the path.
*/
template <typename Family, typename State, typename Corrector>
[[nodiscard]] auto run(const Family& family, const State& initial_state, Corrector&& correct, Options<Family> options, double absolute_tolerance) {
using namespace utils::errors;
using Report = std::remove_cvref_t<std::invoke_result_t<Corrector&, typename Family::ModelType, State, CorrectionControls>>;
using Record = StageRecord<Family::coordinate_count>;
constexpr std::size_t count = Family::coordinate_count;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(std::isfinite(absolute_tolerance) && absolute_tolerance > 0.0 && options.maximum_work > 0,
"Continuation requires a fixed positive absolute tolerance and a positive total work budget.");
const CoordinateController<count> validation(options.start, family.target_coordinates(), options.steps);
Result<Report, count> result{.target = family.target_coordinates(), .absolute_tolerance = absolute_tolerance};
result.stages.reserve(std::min<std::uint64_t>(options.maximum_records, options.steps.maximum_attempts));
// Both initial qualification and later trials use this same certification rule and physical-state boundary.
const auto attempt = [&](const typename Family::Coordinates& coordinates, const State& state, bool starting, std::optional<std::size_t> axis) {
Record stage{.coordinates = coordinates, .changed_axis = axis, .starting_star = starting};
auto model = family.try_model_at(coordinates);
std::unique_ptr<Report> report;
if (!model) {
stage.disposition = StageDisposition::invalid_model;
stage.message = model.error().message;
} else {
const std::uint64_t remaining = options.maximum_work - result.total_work();
report = std::make_unique<Report>(std::invoke(correct, std::move(*model), State(state), CorrectionControls{absolute_tolerance, remaining}));
stage.work = static_cast<std::uint64_t>(report->diagnostics().attempted_iterations);
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(stage.work <= remaining, "A continuation corrector exceeded its supplied work budget.");
stage.residual_norm = report->diagnostics().final_residual_norm;
const bool certified = report->converged() && report->has_state() && stage.residual_norm &&
std::isfinite(*stage.residual_norm) && *stage.residual_norm >= 0.0 && *stage.residual_norm <= absolute_tolerance;
stage.disposition = certified ? StageDisposition::converged : StageDisposition::not_converged;
stage.message = report->message();
}
return std::pair{std::move(report), std::move(stage)};
};
for (const auto& coordinates : detail::starting_candidates(family, options)) {
if (result.starting_attempts == options.steps.maximum_attempts || result.total_work() == options.maximum_work) {
break;
}
++result.starting_attempts;
auto [report, stage] = attempt(coordinates, initial_state, true, std::nullopt);
result.starting_work += stage.work;
const StageDisposition disposition = stage.disposition;
const std::string message = stage.message;
detail::record(result, std::move(stage), options.maximum_records);
if (disposition == StageDisposition::invalid_model) {
result.status = Status::invalid_model;
result.message = message;
return result;
}
if (disposition == StageDisposition::converged) {
result.accepted.swap(report);
result.accepted_coordinates = coordinates;
break;
}
result.failed.swap(report);
}
if (!result.accepted) {
result.status = Status::starting_star_failed;
result.message = "No permitted starting model produced a star certified at the fixed continuation tolerance.";
return result;
}
if (*result.accepted_coordinates == result.target) {
result.status = Status::target_reached;
return result;
}
if (result.total_attempts() == options.steps.maximum_attempts) {
result.status = Status::attempt_limit;
return result;
}
auto steps = options.steps;
steps.maximum_attempts -= result.starting_attempts;
CoordinateController<count> controller(*result.accepted_coordinates, result.target, steps);
while (result.total_work() < options.maximum_work) {
const auto proposal = controller.propose();
if (!proposal) {
break;
}
auto [report, stage] = attempt(proposal->coordinates, result.accepted->state(), false, proposal->axis);
const StageDisposition disposition = stage.disposition;
AttemptOutcome outcome = AttemptOutcome::recoverable_rejection;
if (disposition == StageDisposition::converged) {
outcome = AttemptOutcome::success;
} else if (disposition == StageDisposition::invalid_model) {
outcome = AttemptOutcome::fatal_failure;
}
const std::uint64_t work = stage.work;
const std::string message = stage.message;
controller.validate_resolution(proposal->id, outcome, work);
detail::record(result, std::move(stage), options.maximum_records);
controller.resolve(proposal->id, outcome, work);
if (outcome == AttemptOutcome::success) {
result.accepted.swap(report);
result.accepted_coordinates = controller.accepted();
} else {
result.failed.swap(report);
}
result.counters = controller.counters();
if (outcome == AttemptOutcome::fatal_failure) {
result.status = Status::invalid_model;
result.message = message;
return result;
}
if (controller.status() == ControllerStatus::complete) {
result.status = Status::target_reached;
return result;
}
}
result.counters = controller.counters();
if (result.total_work() == options.maximum_work) {
result.status = Status::work_limit;
} else if (controller.status() == ControllerStatus::attempt_limit) {
result.status = Status::attempt_limit;
} else {
result.status = Status::no_progress;
}
return result;
}
/** @brief Qualify a supplied star afresh after checking every physical family input. */
template <typename Family, typename Model, typename State, typename Corrector>
[[nodiscard]] auto run(const Family& family, const StartingStar<Model, State>& initial, Corrector&& correct,
Options<Family> options, double absolute_tolerance) {
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(!options.explicit_start,
"A supplied starting star determines its own physical coordinates. Remove starting_at settings when continuing from that star.");
family.validate_member(initial.model);
options.start = family.coordinates_of(initial.model);
options.explicit_start = true;
return run(family, initial.state, std::forward<Corrector>(correct), std::move(options), absolute_tolerance);
}
}

88
src/include/serif/solver/continuation/equilibrium.hpp Normal file
View File

@@ -0,0 +1,88 @@
#pragma once
#include "serif/solver/continuation/driver.hpp"
#include "serif/backend/mfem/solver/continuation/detail/physical_transfer.hpp"
#include "serif/solver/equilibrium.hpp"
namespace serif::solver::continuation {
/** @brief Continuation uses one absolute qualification threshold across all family members. */
[[nodiscard]] inline SolverOptions<> default_solver_options() {
return {.method = {.relative_tolerance = 0.0, .absolute_tolerance = 1.0e-10}};
}
namespace detail {
/** @brief Adapt one fresh configured discretization per stage to the same public equilibrium solve. */
template <typename Family, typename Factory, typename SolverOptions, typename Solve>
[[nodiscard]] auto corrector(const Family& family, Factory& configure, SolverOptions options, Solve solve) {
options.method.validate();
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
options.method.relative_tolerance == 0.0 && options.method.absolute_tolerance > 0.0,
"Continuation requires zero relative tolerance and a fixed positive absolute Newton tolerance.");
using Configuration = std::invoke_result_t<Factory&, typename Family::ModelType>;
using FiniteElements = std::remove_cvref_t<decltype(std::declval<Configuration>().finite_elements)>;
using TransferGuard = decltype(backend::mfem::physical_transfer_guard(std::declval<const FiniteElements&>()));
return [&family, &configure, options = std::move(options), solve = std::move(solve), transfer = std::optional<TransferGuard>{}]
(auto model, auto state, CorrectionControls controls) mutable {
const auto expected = family.model_signature(model);
auto configuration = std::invoke(configure, std::move(model));
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
family.model_signature(configuration.model) == expected,
"The continuation configuration factory changed the requested physical model. Preserve its prescriptions and equation identity.");
if (!transfer) {
transfer.emplace(backend::mfem::physical_transfer_guard(configuration.finite_elements));
}
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
(*transfer)(configuration.finite_elements),
"Continuation cannot transfer physical coefficients to changed finite-element bindings, reference geometry, or deformation coordinates.");
auto stage_options = options;
stage_options.method.absolute_tolerance = controls.absolute_tolerance;
stage_options.method.maximum_iterations = static_cast<int>(std::min<std::uint64_t>(
stage_options.method.maximum_iterations, controls.maximum_work));
return std::invoke(solve, std::move(configuration), std::move(state), std::move(stage_options));
};
}
}
}
namespace serif::solver {
/**
* @brief Continue moving-star equilibria through an explicitly declared physical family.
* @param configure Own or borrow a callable that returns a configured discretization for its supplied model.
* It must preserve the physical coordinate layout and finite-element reference of transferred states; fresh evaluators are implicit.
* @note The explicit factory supports move-only resource policies without inventing a mutable clone/context hierarchy.
*/
template <typename Family, typename Factory, typename State, typename Options = SolverOptions<>>
requires std::copy_constructible<Options>
[[nodiscard]] auto continue_equilibrium(const Family& family, Factory configure, const State& initial_state,
continuation::Options<Family> continuation_options, Options solver_options = continuation::default_solver_options()) {
const double tolerance = solver_options.method.absolute_tolerance;
auto correct = continuation::detail::corrector(family, configure, std::move(solver_options), [](auto configuration, auto state, auto options) {
return solve_equilibrium(std::move(configuration), std::move(state), std::move(options));
});
return continuation::run(family, initial_state, std::move(correct), std::move(continuation_options), tolerance);
}
/** @brief Start with the family's physics-owned default coordinates and bounded fallback route. */
template <typename Family, typename Factory, typename State>
[[nodiscard]] auto continue_equilibrium(const Family& family, Factory configure, const State& initial_state) {
return continue_equilibrium(family, std::move(configure), initial_state, serif::solver::continuation::make_options(family));
}
/** @brief Continue a fixed-geometry model while preserving its prescribed displacement. */
template <typename Family, typename Factory, typename State, typename Options = SolverOptions<>>
requires std::copy_constructible<Options>
[[nodiscard]] auto continue_equilibrium(const Family& family, Factory configure, const State& initial_state, const Vector& displacement,
continuation::Options<Family> continuation_options, Options solver_options = continuation::default_solver_options()) {
const double tolerance = solver_options.method.absolute_tolerance;
auto correct = continuation::detail::corrector(family, configure, std::move(solver_options), [&displacement](auto configuration, auto state, auto options) {
return solve_equilibrium(std::move(configuration), std::move(state), displacement, std::move(options));
});
return continuation::run(family, initial_state, std::move(correct), std::move(continuation_options), tolerance);
}
/** @brief Use the same default route and absolute qualification tolerance for fixed geometry. */
template <typename Family, typename Factory, typename State>
[[nodiscard]] auto continue_equilibrium(const Family& family, Factory configure, const State& initial_state, const Vector& displacement) {
return continue_equilibrium(family, std::move(configure), initial_state, displacement, serif::solver::continuation::make_options(family));
}
}

372
src/include/serif/solver/continuation/family.hpp Normal file
View File

@@ -0,0 +1,372 @@
#pragma once
#include <algorithm>
#include <cmath>
#include <functional>
#include <limits>
#include <optional>
#include <tuple>
#include <utility>
#include <vector>
#include "serif/solver/continuation/parameters.hpp"
#include "serif/model/stellar_model.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::solver::continuation {
/** @brief Name only independent physical inputs to a dependent law. */
template <PhysicalParameter... Parameters>
[[nodiscard]] constexpr auto depends_on(Parameters...) noexcept { return ParameterList<Parameters...>{}; }
template <PhysicalParameter Output, typename Dependencies, typename Function>
struct ParameterLaw final {
using OutputParameter = Output;
using DependencyParameters = Dependencies;
Function function;
};
namespace detail {
/** @brief Reject null callable pointers even when a public aggregate bypasses its factory. */
template <typename Function>
[[nodiscard]] bool callable_valid(const Function& value) {
if constexpr (std::is_pointer_v<Function> || std::is_member_pointer_v<Function>) {
return value != nullptr;
} else {
return true;
}
}
/** @brief Compare independently evaluated physical laws without an absolute unit-dependent tolerance. */
[[nodiscard]] inline bool values_agree(double first, double second) noexcept {
return std::isfinite(first) && std::isfinite(second) && (first == second ||
std::abs(first - second) <= 64.0 * std::numeric_limits<double>::epsilon() * std::max(std::abs(first), std::abs(second)));
}
template <typename List, typename Superset>
struct Contained;
template <typename... Parameters, typename Superset>
struct Contained<ParameterList<Parameters...>, Superset> : std::bool_constant<(utils::dependencies::ContainsKey<Parameters, Superset> && ...)> {};
template <typename Candidate>
concept LawDefinition = requires {
typename Candidate::OutputParameter;
typename Candidate::DependencyParameters;
requires PhysicalParameter<typename Candidate::OutputParameter>;
requires utils::dependencies::IsKeyList<typename Candidate::DependencyParameters>;
};
template <typename Selected, typename Laws>
struct FamilyTraits;
template <typename... Selected, typename... Laws>
struct FamilyTraits<ParameterList<Selected...>, std::tuple<Laws...>> {
using Derived = ParameterList<typename Laws::OutputParameter...>;
static constexpr bool valid = sizeof...(Selected) > 0 && utils::dependencies::IsKeyList<ParameterList<Selected...>> &&
utils::dependencies::IsKeyList<Derived> &&
(!utils::dependencies::ContainsKey<typename Laws::OutputParameter, ParameterList<Selected...>> && ...) &&
(Contained<typename Laws::DependencyParameters, ParameterList<Selected...>>::value && ...);
static constexpr bool complete = valid && (Contained<typename RequiredDependents<Selected>::Type, Derived>::value && ...);
template <typename Model>
static constexpr bool owners_present = (Model::template contains_specification<typename Selected::Specification> && ...) &&
(Model::template contains_specification<typename Laws::OutputParameter::Specification> && ...);
template <typename Specification>
static constexpr bool affects = (std::same_as<Specification, typename Selected::Specification> || ...) ||
(std::same_as<Specification, typename Laws::OutputParameter::Specification> || ...);
};
template <typename Specifications>
struct AllIdentified;
template <typename... Specifications>
struct AllIdentified<utils::types::TypeList<Specifications...>> : std::bool_constant<(IdentifiedSpecification<Specifications> && ...)> {};
template <typename Model, typename Selected, typename Laws>
concept BindableFamily = model::StellarModelType<Model> && FamilyTraits<Selected, Laws>::complete &&
FamilyTraits<Selected, Laws>::template owners_present<Model> && IdentifiedSpecification<typename Model::EquationOfState> &&
AllIdentified<typename Model::SpecificationTypes>::value &&
(std::is_empty_v<typename Model::Equations> || IdentifiedSpecification<typename Model::Equations>);
struct CopyEquations final {
/** @brief A copyable equation definition needs no separate construction policy. */
template <typename Equations>
requires std::copy_constructible<Equations>
[[nodiscard]] Equations operator()(const Equations& value) const { return value; }
};
/** @brief Bind target-dependent laws once, owning all resulting captures. */
template <typename Law, typename Model>
[[nodiscard]] auto bind_law(Law&& value, const Model& target) {
if constexpr (requires { std::forward<Law>(value).bind(target); }) {
return std::forward<Law>(value).bind(target);
} else {
return std::remove_cvref_t<Law>(std::forward<Law>(value));
}
}
}
/** @brief Own a typed law, including move-only captures; its inputs are selected coordinates rather than other laws. */
template <PhysicalParameter Output, PhysicalParameter... Inputs, typename Function>
requires utils::dependencies::IsKeyList<ParameterList<Inputs...>> &&
std::invocable<const std::remove_cvref_t<Function>&, typename Inputs::Value...> &&
std::same_as<std::remove_cvref_t<std::invoke_result_t<const std::remove_cvref_t<Function>&, typename Inputs::Value...>>, typename Output::Value>
[[nodiscard]] auto law(Output, ParameterList<Inputs...>, Function&& function) {
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(detail::callable_valid(function),
"A continuation law cannot use a null callable pointer.");
return ParameterLaw<Output, ParameterList<Inputs...>, std::remove_cvref_t<Function>>{.function = std::forward<Function>(function)};
}
template <typename Model, typename Selected, typename Laws, typename EquationFactory>
class BoundFamily;
/** @brief Own target physics and construct fresh family members without a mesh, evaluator, or mutable solve context. */
template <model::StellarModelType Model, PhysicalParameter... Selected, typename... Laws, typename EquationFactory>
class BoundFamily<Model, ParameterList<Selected...>, std::tuple<Laws...>, EquationFactory> final {
using Selection = ParameterList<Selected...>;
using Traits = detail::FamilyTraits<Selection, std::tuple<Laws...>>;
using DerivedValues = std::tuple<typename Laws::OutputParameter::Value...>;
public:
using ModelType = Model;
static constexpr std::size_t coordinate_count = sizeof...(Selected);
using Coordinates = std::array<double, coordinate_count>;
template <typename Parameter>
static constexpr bool has_coordinate = utils::dependencies::ContainsKey<Parameter, Selection>;
/** @brief Locate one selected physical parameter without runtime type dispatch. */
template <PhysicalParameter Parameter>
requires has_coordinate<Parameter>
[[nodiscard]] static consteval std::size_t coordinate_index(Parameter = {}) {
return utils::dependencies::key_index<Parameter, Selection>();
}
/** @brief Expose stable parameter names for diagnostics and persisted coordinates. */
[[nodiscard]] static constexpr std::array<ParameterDescriptor, coordinate_count> parameter_descriptors() noexcept {
return {{{Selected::name, Selected::stable_id}...}};
}
/** @brief Reject inconsistent laws and unverifiable target identities before starting continuation. */
BoundFamily(Model target, std::tuple<Laws...> laws, EquationFactory equations) :
m_target(std::move(target)), m_laws(std::move(laws)), m_equations(std::move(equations)), m_target_coordinates(coordinates_of(m_target)) {
static_assert(detail::BindableFamily<Model, Selection, std::tuple<Laws...>>);
const std::vector<double> identity = model_signature(m_target);
const bool finite = std::ranges::all_of(identity, [](double value) { return std::isfinite(value); });
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(finite,
"Every physical input in a continuation target identity must be finite.");
const auto endpoint = try_model_at(m_target_coordinates);
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(endpoint.has_value(),
"The continuation target is inconsistent: {}", endpoint ? std::string{} : endpoint.error().message);
}
/** @brief Inspect the immutable requested physical model. */
[[nodiscard]] const Model& target_model() const noexcept { return m_target; }
/** @brief Preserve the exact endpoint independently of adaptive step sizes. */
[[nodiscard]] Coordinates target_coordinates() const noexcept { return m_target_coordinates; }
/** @brief Use physics-owned simpler coordinates where a parameter supplies them. */
[[nodiscard]] Coordinates starting_coordinates() const {
const auto start = [&]<PhysicalParameter Parameter>() {
auto value = Parameter::read(m_target.template specification<typename Parameter::Specification>());
if constexpr (requires { { Parameter::starting_value(value) } -> std::same_as<typename Parameter::Value>; }) {
value = Parameter::starting_value(value);
}
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(Parameter::valid(value),
"The default starting {} lies outside its physical domain.", Parameter::name);
return coordinate_value(value);
};
return {start.template operator()<Selected>()...};
}
/** @brief Extract coordinates; extraction alone does not certify family membership. */
[[nodiscard]] Coordinates coordinates_of(const Model& candidate) const {
return {coordinate_value(Selected::read(candidate.template specification<typename Selected::Specification>()))...};
}
/** @brief Include unchanged EOS, constraints, and stateful equation definitions in family identity. */
[[nodiscard]] std::vector<double> model_signature(const Model& candidate) const {
std::vector<double> signature;
const auto append = [&](const auto& value) {
const auto identity = SpecificationTraits<std::remove_cvref_t<decltype(value)>>::identity(value);
signature.push_back(static_cast<double>(std::ranges::size(identity)));
signature.insert(signature.end(), std::ranges::begin(identity), std::ranges::end(identity));
};
candidate.visit_specifications(append);
if constexpr (!std::is_empty_v<typename Model::Equations>) {
append(candidate.equations());
}
return signature;
}
/** @brief Reject unrelated physical anchors even if their selected coordinates match. */
[[nodiscard]] std::expected<void, Failure> try_validate_member(const Model& candidate) const {
const auto expected = try_model_at(coordinates_of(candidate));
if (!expected) {
return std::unexpected(expected.error());
}
const auto actual_identity = model_signature(candidate);
const auto expected_identity = model_signature(*expected);
if (actual_identity.size() != expected_identity.size() || !std::ranges::equal(actual_identity, expected_identity, detail::values_agree)) {
return std::unexpected(Failure{.message = "The starting model does not belong to the declared physical family."});
}
return {};
}
/** @brief Validate a user-supplied anchor at the configuration boundary. */
void validate_member(const Model& candidate) const {
const auto result = try_validate_member(candidate);
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(result.has_value(),
"{}", result ? std::string{} : result.error().message);
}
/** @brief Construct a candidate and report out-of-domain laws as ordinary values for the continuation driver. */
[[nodiscard]] std::expected<Model, Failure> try_model_at(const Coordinates& coordinates) const {
if (!(Selected::valid(typename Selected::Value{coordinates[coordinate_index(Selected{})]}) && ...)) {
return std::unexpected(Failure{.message = "A continuation coordinate lies outside its physical parameter domain."});
}
auto derived = evaluate_laws(coordinates, std::index_sequence_for<Laws...>{});
if (!derived) {
return std::unexpected(derived.error());
}
const bool endpoint = coordinates == m_target_coordinates;
const Values values{m_target, coordinates, *derived};
auto candidate = m_target.transform_specifications([&](const auto& specification) {
using Specification = std::remove_cvref_t<decltype(specification)>;
if constexpr (Traits::template affects<Specification>) {
if (!endpoint) {
return SpecificationTraits<Specification>::rebuild(specification, values);
}
}
return SpecificationTraits<Specification>::clone(specification);
}, std::invoke(m_equations, m_target.equations()));
if (endpoint && model_signature(candidate) != model_signature(m_target)) {
return std::unexpected(Failure{.message = "A continuation clone changed the exact target physical identity."});
}
if (!realized(candidate, coordinates, *derived, std::index_sequence_for<Laws...>{})) {
return std::unexpected(Failure{.message = "A continuation law or prescription did not reproduce its selected and derived physical parameters."});
}
return candidate;
}
/** @brief Construct a family member when invalid coordinates are a caller configuration error. */
[[nodiscard]] Model model_at(const Coordinates& coordinates) const {
auto result = try_model_at(coordinates);
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(result.has_value(),
"{}", result ? std::string{} : result.error().message);
return std::move(*result);
}
private:
struct Values final {
const Model& target;
const Coordinates& coordinates;
const DerivedValues& derived;
/** @brief Read a selected, derived, or unchanged physical value during synchronous reconstruction. */
template <PhysicalParameter Parameter>
[[nodiscard]] typename Parameter::Value get(Parameter) const {
if constexpr (has_coordinate<Parameter>) {
return typename Parameter::Value{coordinates[coordinate_index(Parameter{})]};
} else if constexpr (utils::dependencies::ContainsKey<Parameter, typename Traits::Derived>) {
return std::get<utils::dependencies::key_index<Parameter, typename Traits::Derived>()>(derived);
} else {
return Parameter::read(target.template specification<typename Parameter::Specification>());
}
}
};
/** @brief Read only independent coordinates into a typed dependent law. */
template <typename Law, typename... Inputs>
[[nodiscard]] auto evaluate_law(const Law& value, const Coordinates& coordinates, ParameterList<Inputs...>) const {
using Output = typename Law::OutputParameter;
using Result = std::expected<typename Output::Value, Failure>;
if (!detail::callable_valid(value.function)) {
return Result{std::unexpected(Failure{.message = "A continuation law has a null callable pointer."})};
}
auto result = std::invoke(value.function, typename Inputs::Value{coordinates[coordinate_index(Inputs{})]}...);
if (!Output::valid(result)) {
return Result{std::unexpected(Failure{.message = "A dependent continuation law returned an invalid physical value."})};
}
return Result{std::move(result)};
}
/** @brief Evaluate every independent law before rebuilding any physical prescription. */
template <std::size_t... Indices>
[[nodiscard]] std::expected<DerivedValues, Failure> evaluate_laws(const Coordinates& coordinates, std::index_sequence<Indices...>) const {
auto results = std::tuple{evaluate_law(std::get<Indices>(m_laws), coordinates, typename Laws::DependencyParameters{})...};
std::optional<Failure> failure;
const auto check = [&](const auto& result) {
if (!result && !failure) {
failure = result.error();
}
};
(check(std::get<Indices>(results)), ...);
if (failure) {
return std::unexpected(std::move(*failure));
}
return DerivedValues{std::move(*std::get<Indices>(results))...};
}
/** @brief Read back every installed parameter; exact selected coordinates cannot be silently ignored. */
template <std::size_t... Indices>
[[nodiscard]] bool realized(const Model& result, const Coordinates& coordinates, const DerivedValues& derived, std::index_sequence<Indices...>) const {
return ((coordinate_value(Selected::read(result.template specification<typename Selected::Specification>())) == coordinates[coordinate_index(Selected{})]) && ...) &&
(detail::values_agree(coordinate_value(Laws::OutputParameter::read(result.template specification<typename Laws::OutputParameter::Specification>())),
coordinate_value(std::get<Indices>(derived))) && ...);
}
Model m_target;
std::tuple<Laws...> m_laws;
EquationFactory m_equations;
Coordinates m_target_coordinates;
};
/** @brief Own a reusable physical family definition before selecting a target model. */
template <typename Selected, typename Laws>
class Family final {
public:
/** @brief Capture dependent laws by value. */
explicit Family(Laws laws = {}) : m_laws(std::move(laws)) {}
/** @brief Extend a copyable family without changing its independent coordinates. */
template <detail::LawDefinition... Added>
requires detail::FamilyTraits<Selected, decltype(std::tuple_cat(std::declval<Laws>(), std::declval<std::tuple<Added...>>()))>::valid
[[nodiscard]] auto with(Added... laws) const& requires std::copy_constructible<Laws> {
return Family<Selected, decltype(std::tuple_cat(m_laws, std::tuple<Added...>{std::move(laws)...}))>{
std::tuple_cat(m_laws, std::tuple<Added...>{std::move(laws)...})};
}
/** @brief Transfer move-only law captures when extending a temporary definition. */
template <detail::LawDefinition... Added>
requires detail::FamilyTraits<Selected, decltype(std::tuple_cat(std::declval<Laws>(), std::declval<std::tuple<Added...>>()))>::valid
[[nodiscard]] auto with(Added... laws) && {
auto joined = std::tuple_cat(std::move(m_laws), std::tuple<Added...>{std::move(laws)...});
return Family<Selected, decltype(joined)>{std::move(joined)};
}
/** @brief Bind a reproducible model, using an explicit factory for move-only equation definitions when necessary. */
template <typename Model, typename EquationFactory = detail::CopyEquations>
requires detail::BindableFamily<Model, Selected, Laws> && std::copy_constructible<Laws> &&
std::invocable<const EquationFactory&, const typename Model::Equations&>
[[nodiscard]] auto bind(Model target, EquationFactory factory = {}) const& {
return bind_impl(std::move(target), m_laws, std::move(factory));
}
/** @brief Bind a target while transferring ownership of move-only laws. */
template <typename Model, typename EquationFactory = detail::CopyEquations>
requires detail::BindableFamily<Model, Selected, Laws> && std::invocable<const EquationFactory&, const typename Model::Equations&>
[[nodiscard]] auto bind(Model target, EquationFactory factory = {}) && {
return bind_impl(std::move(target), std::move(m_laws), std::move(factory));
}
private:
/** @brief Resolve target-dependent laws before moving their physical anchor into the family. */
template <typename Model, typename LawValues, typename EquationFactory>
[[nodiscard]] static auto bind_impl(Model target, LawValues&& values, EquationFactory factory) {
auto bound = std::apply([&](auto&&... value) {
return std::tuple{detail::bind_law(std::forward<decltype(value)>(value), target)...};
}, std::forward<LawValues>(values));
return BoundFamily<Model, Selected, decltype(bound), EquationFactory>{std::move(target), std::move(bound), std::move(factory)};
}
Laws m_laws;
};
/** @brief Select independent physical coordinates with duplicate detection at compile time. */
template <PhysicalParameter... Parameters>
requires (sizeof...(Parameters) > 0) && utils::dependencies::IsKeyList<ParameterList<Parameters...>>
[[nodiscard]] auto in(Parameters...) { return Family<ParameterList<Parameters...>, std::tuple<>>{}; }
}

39
src/include/serif/solver/continuation/identified_equations.hpp Normal file
View File

@@ -0,0 +1,39 @@
#pragma once
#include <array>
#include <functional>
#include <tuple>
#include <utility>
#include "serif/solver/continuation/parameters.hpp"
namespace serif::solver::continuation {
/** @brief Give a deferred equation definition an explicit physical identity without runtime polymorphism. */
template <typename Definition, CoordinateValue... Values>
struct IdentifiedEquations final {
Definition definition;
std::tuple<Values...> identity;
/** @brief Read fixed captured quantities in their declared physical units for family membership. */
[[nodiscard]] auto continuation_identity() const {
return std::apply([](const auto&... values) { return std::array<double, sizeof...(Values)>{coordinate_value(values)...}; }, identity);
}
/** @brief Rebuild equations from the stage's own model while retaining the supplied callable. */
template <typename Model>
requires std::invocable<const Definition&, const Model&>
[[nodiscard]] decltype(auto) operator()(const Model& model) const { return std::invoke(definition, model); }
};
/**
* @brief Declare every fixed physical constant captured by a deferred equation definition.
* @note Include all captured physical inputs, such as G and surface pressure. Values read from
* the supplied model belong to its specification identity already. The declaration is a user
* contract: C++ cannot inspect arbitrary captures or prove that a callback uses only these inputs.
* Move-only definitions remain supported with an explicit equation factory when binding a family.
*/
template <typename Definition, CoordinateValue... Values>
[[nodiscard]] auto identified_equations(Definition definition, Values... values) {
return IdentifiedEquations<Definition, Values...>{std::move(definition), std::tuple<Values...>{std::move(values)...}};
}
}

106
src/include/serif/solver/continuation/options.hpp Normal file
View File

@@ -0,0 +1,106 @@
#pragma once
#include "serif/solver/continuation/family.hpp"
#include "serif/solver/continuation/controller.hpp"
namespace serif::solver::continuation {
/** @brief Bound continuation work and storage independently of Newton's per-correction controls. */
template <typename Family>
struct Options final {
CoordinateOptions<Family::coordinate_count> steps{};
typename Family::Coordinates start{};
std::array<bool, Family::coordinate_count> specified{};
bool explicit_start{false};
bool has_ordering{false};
bool has_numerics{false};
std::uint64_t maximum_work{100000};
std::size_t maximum_records{128};
};
template <PhysicalParameter Parameter>
struct StartingAt final {
typename Parameter::Value value;
/** @brief Change one typed starting coordinate with strong rejection of duplicate or invalid input. */
template <typename Family>
requires Family::template has_coordinate<Parameter>
void apply(Options<Family>& options) const {
constexpr std::size_t axis = Family::coordinate_index(Parameter{});
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(!options.specified[axis] && Parameter::valid(value),
"The starting {} must be specified once and lie within its physical domain.", Parameter::name);
options.start[axis] = coordinate_value(value);
options.specified[axis] = true;
options.explicit_start = true;
}
};
/** @brief Select a physical starting value with the parameter's required units. */
template <PhysicalParameter Parameter>
[[nodiscard]] auto starting_at(Parameter, typename Parameter::Value value) { return StartingAt<Parameter>{value}; }
template <bool Strict, PhysicalParameter... Parameters>
requires (sizeof...(Parameters) > 0) && utils::dependencies::IsKeyList<ParameterList<Parameters...>>
struct Ordering final {
/** @brief A preference may yield to other axes; strict stages require every coordinate exactly once. */
template <typename Family>
requires (Family::template has_coordinate<Parameters> && ...) && (!Strict || sizeof...(Parameters) == Family::coordinate_count)
void apply(Options<Family>& options) const {
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(!options.has_ordering,
"Choose only one continuation ordering policy.");
constexpr std::array selected{Family::coordinate_index(Parameters{})...};
std::size_t index = 0;
for (const std::size_t axis : selected) {
options.steps.axis_order[index++] = axis;
}
for (std::size_t axis = 0; axis < Family::coordinate_count; ++axis) {
if (std::ranges::find(selected, axis) == selected.end()) {
options.steps.axis_order[index++] = axis;
}
}
options.steps.order = Strict ? CoordinateOrder::strict : CoordinateOrder::preferred;
options.has_ordering = true;
}
};
/** @brief Try these physical coordinates first, allowing blocked coordinates to yield. */
template <PhysicalParameter... Parameters>
requires (sizeof...(Parameters) > 0) && utils::dependencies::IsKeyList<ParameterList<Parameters...>>
[[nodiscard]] auto prefer(Parameters...) { return Ordering<false, Parameters...>{}; }
/** @brief Complete each physical stage before attempting the next. */
template <PhysicalParameter... Parameters>
requires (sizeof...(Parameters) > 0) && utils::dependencies::IsKeyList<ParameterList<Parameters...>>
[[nodiscard]] auto stages(Parameters...) { return Ordering<true, Parameters...>{}; }
template <std::size_t Count>
struct NumericalControls final {
CoordinateOptions<Count> options;
/** @brief Preserve explicit physical ordering independently of option argument order. */
template <typename Family>
requires (Family::coordinate_count == Count)
void apply(Options<Family>& target) const {
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(!target.has_numerics,
"Continuation numerical controls may be supplied only once.");
const auto order = target.steps.order;
const auto axes = target.steps.axis_order;
target.steps = options;
if (target.has_ordering) {
target.steps.order = order;
target.steps.axis_order = axes;
}
target.has_numerics = true;
}
};
/** @brief Expose advanced numerical coordinate stepping without mixing it into physical prescriptions. */
template <std::size_t Count>
[[nodiscard]] auto numerics(CoordinateOptions<Count> options) { return NumericalControls<Count>{options}; }
/** @brief Resolve physical settings from family defaults before any corrector runs. */
template <typename Family, typename... Settings>
requires (requires(const Settings& setting, Options<Family>& options) { setting.apply(options); } && ...)
[[nodiscard]] auto make_options(const Family& family, const Settings&... settings) {
Options<Family> result{.start = family.starting_coordinates()};
(settings.apply(result), ...);
return result;
}
}

104
src/include/serif/solver/continuation/parameters.hpp Normal file
View File

@@ -0,0 +1,104 @@
#pragma once
#include <array>
#include <concepts>
#include <expected>
#include <ranges>
#include <string>
#include <string_view>
#include <type_traits>
#include "serif/dimensions/quantities.hpp"
#include "serif/utils/dependencies/key_list_operations.hpp"
namespace serif::solver::continuation {
template <typename... Parameters>
using ParameterList = utils::dependencies::KeyList<Parameters...>;
struct Failure final {
std::string message;
};
struct ParameterDescriptor final {
std::string_view name;
std::string_view stable_id;
};
/** @brief Add continuation access without imposing numerical dependencies on a physical prescription. */
template <typename Specification>
struct SpecificationTraits {
/** @brief Custom physics explicitly identifies every immutable physical input. */
[[nodiscard]] static auto identity(const Specification& value)
requires requires { value.continuation_identity(); } {
return value.continuation_identity();
}
/** @brief Rebuild custom physics from selected and derived parameters. */
template <typename Values>
[[nodiscard]] static Specification rebuild(const Specification& value, const Values& parameters)
requires requires { { value.with_continuation_parameters(parameters) } -> std::same_as<Specification>; } {
return value.with_continuation_parameters(parameters);
}
/** @brief Use an explicit clone hook when physics cannot be copied as an ordinary value. */
[[nodiscard]] static Specification clone(const Specification& value) {
if constexpr (requires { { value.clone_for_continuation() } -> std::same_as<Specification>; }) {
return value.clone_for_continuation();
} else {
return value;
}
}
};
template <typename Value>
concept CoordinateValue = std::same_as<Value, double> || dimensions::QuantityValueType<Value>;
/** @brief Convert physical units only at the coordinate-controller boundary. */
template <CoordinateValue Value>
[[nodiscard]] constexpr double coordinate_value(Value value) noexcept {
if constexpr (std::same_as<Value, double>) {
return value;
} else {
return value.value();
}
}
namespace detail {
template <typename Owner, typename = void>
struct ParameterCatalog { using Type = ParameterList<>; };
template <typename Owner>
struct ParameterCatalog<Owner, std::void_t<typename Owner::ContinuationParameters>> { using Type = typename Owner::ContinuationParameters; };
template <typename Owner, typename = void>
struct Catalog : ParameterCatalog<Owner> {};
template <typename Owner>
struct Catalog<Owner, std::void_t<typename SpecificationTraits<Owner>::Parameters>> { using Type = typename SpecificationTraits<Owner>::Parameters; };
template <typename Parameter, typename = void>
struct RequiredDependents { using Type = ParameterList<>; };
template <typename Parameter>
struct RequiredDependents<Parameter, std::void_t<typename Parameter::RequiredDependentParameters>> {
using Type = typename Parameter::RequiredDependentParameters;
};
}
/** @brief A typed physical parameter owns its domain predicate, units, owner, and stable identity. */
template <typename Parameter>
concept PhysicalParameter = requires(const typename Parameter::Specification& specification, typename Parameter::Value value) {
requires CoordinateValue<typename Parameter::Value>;
{ Parameter::name } -> std::convertible_to<std::string_view>;
{ Parameter::stable_id } -> std::convertible_to<std::string_view>;
{ Parameter::read(specification) } -> std::same_as<typename Parameter::Value>;
{ Parameter::valid(value) } -> std::same_as<bool>;
requires utils::dependencies::ContainsKey<Parameter, typename detail::Catalog<typename Parameter::Specification>::Type>;
};
/** @brief Require explicit identity data, never an object address or incidental memory representation. */
template <typename Specification>
concept IdentifiedSpecification = requires(const Specification& value) {
SpecificationTraits<Specification>::identity(value);
requires std::ranges::forward_range<decltype(SpecificationTraits<Specification>::identity(value))>;
requires std::ranges::sized_range<decltype(SpecificationTraits<Specification>::identity(value))>;
requires std::same_as<std::ranges::range_value_t<decltype(SpecificationTraits<Specification>::identity(value))>, double>;
};
}

46
src/include/serif/solver/continuation/polytropic.hpp Normal file
View File

@@ -0,0 +1,46 @@
#pragma once
#include "serif/solver/continuation/family.hpp"
#include "serif/solver/continuation/stellar_parameters.hpp"
namespace serif::solver::continuation {
/** @brief Vary K with n to preserve pressure at one explicitly selected physical reference density. */
class PreservePressureAtReferenceDensity final {
public:
using OutputParameter = PolytropicConstant;
using DependencyParameters = ParameterList<PolytropicIndex>;
/** @brief Validate a positive finite thermodynamic pressure anchor. */
explicit PreservePressureAtReferenceDensity(dimensions::DensityValue density) : m_density(density) {
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
std::isfinite(density.value()) && density.value() > 0.0, "Pressure-reference density must be finite and positive.");
}
/** @brief Freeze target values into an owning law whose exact endpoint retains the original K. */
template <model::StellarModelType Model>
requires Model::template contains_specification<eos::models::Polytrope>
[[nodiscard]] auto bind(const Model& target) const {
const auto& eos = target.template specification<eos::models::Polytrope>();
const double index = eos.polytropic_index();
const double constant = eos.polytropic_constant();
const double density = m_density.value();
const double pressure = eos.evaluate(eos::relations::PressureFromDensity{}, m_density).value();
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
std::isfinite(pressure) && pressure > 0.0, "The target pressure at the reference density must be finite and positive.");
return law(PolytropicConstant{}, depends_on(PolytropicIndex{}), [index, constant, density](double candidate) {
if (candidate == index) {
return constant;
}
return std::exp(std::log(constant) + (1.0 / index - 1.0 / candidate) * std::log(density));
});
}
private:
dimensions::DensityValue m_density;
};
/** @brief Choose a pressure-preserving physical EOS family, independently of numerical normalization. */
[[nodiscard]] inline auto preserve_pressure_at_reference_density(dimensions::DensityValue density) {
return PreservePressureAtReferenceDensity{density};
}
}

92
src/include/serif/solver/continuation/stellar_parameters.hpp Normal file
View File

@@ -0,0 +1,92 @@
#pragma once
#include <cmath>
#include "serif/solver/continuation/parameters.hpp"
#include "serif/eos/models/polytropic.hpp"
#include "serif/model/specifications/fixed_angular_momentum.hpp"
#include "serif/model/specifications/fixed_total_mass.hpp"
namespace serif::solver::continuation {
struct PolytropicConstant;
struct PolytropicIndex final {
using Specification = eos::models::Polytrope;
using Value = double;
using RequiredDependentParameters = ParameterList<PolytropicConstant>;
static constexpr std::string_view name = "polytropic index";
static constexpr std::string_view stable_id = "polytrope.index";
/** @brief Read the physical EOS index. */
[[nodiscard]] static double read(const Specification& value) noexcept { return value.polytropic_index(); }
/** @brief Match the domain supported by SERiF's barotropic EOS. */
[[nodiscard]] static bool valid(double value) noexcept { return std::isfinite(value) && value >= 1.0; }
/** @brief Start with the analytically simple n=1 family member. */
[[nodiscard]] static double starting_value(double) noexcept { return 1.0; }
};
struct PolytropicConstant final {
using Specification = eos::models::Polytrope;
using Value = double;
static constexpr std::string_view name = "polytropic constant";
static constexpr std::string_view stable_id = "polytrope.constant";
/** @brief Read the dimensional numerical coefficient in the current physical units. */
[[nodiscard]] static double read(const Specification& value) noexcept { return value.polytropic_constant(); }
/** @brief A finite positive coefficient gives positive pressure at positive density. */
[[nodiscard]] static bool valid(double value) noexcept { return std::isfinite(value) && value > 0.0; }
};
struct AngularMomentum final {
using Specification = model::FixedAngularMomentum;
using Value = dimensions::AngularMomentumValue;
static constexpr std::string_view name = "angular momentum";
static constexpr std::string_view stable_id = "fixed_angular_momentum.target";
/** @brief Read the prescribed angular momentum, independently of solved angular velocity. */
[[nodiscard]] static Value read(const Specification& value) noexcept { return value.target(); }
/** @brief Match the nonnegative rotation prescription. */
[[nodiscard]] static bool valid(Value value) noexcept { return std::isfinite(value.value()) && value.value() >= 0.0; }
/** @brief Start from the nonrotating family member. */
[[nodiscard]] static Value starting_value(Value) noexcept { return Value{0.0}; }
};
template <>
struct SpecificationTraits<eos::models::Polytrope> {
using Parameters = ParameterList<PolytropicIndex, PolytropicConstant>;
/** @brief Include both EOS parameters in family membership. */
[[nodiscard]] static std::array<double, 2> identity(const eos::models::Polytrope& value) noexcept {
return {value.polytropic_index(), value.polytropic_constant()};
}
/** @brief Install n and K together so a dependent law cannot be lost. */
template <typename Values>
[[nodiscard]] static eos::models::Polytrope rebuild(const eos::models::Polytrope&, const Values& values) {
return eos::models::Polytrope{values.get(PolytropicIndex{}), values.get(PolytropicConstant{})};
}
/** @brief Preserve the exact target EOS rather than reevaluating a rounded law. */
[[nodiscard]] static eos::models::Polytrope clone(const eos::models::Polytrope& value) { return value; }
};
template <>
struct SpecificationTraits<model::FixedAngularMomentum> {
using Parameters = ParameterList<AngularMomentum>;
/** @brief Rotation axes and centers are immutable family inputs as well as angular momentum. */
[[nodiscard]] static std::array<double, 7> identity(const model::FixedAngularMomentum& value) noexcept {
return {value.target().value(), value.axis()[0], value.axis()[1], value.axis()[2], value.center()[0], value.center()[1], value.center()[2]};
}
/** @brief Rebuild the target while retaining its rotation axis and center. */
template <typename Values>
[[nodiscard]] static model::FixedAngularMomentum rebuild(const model::FixedAngularMomentum& value, const Values& values) {
return model::FixedAngularMomentum{{.angular_momentum = values.get(AngularMomentum{}), .axis = value.axis(), .center = value.center()}};
}
/** @brief Copy the exact endpoint prescription. */
[[nodiscard]] static model::FixedAngularMomentum clone(const model::FixedAngularMomentum& value) { return value; }
};
/** @brief Preserve unselected scalar prescriptions in the complete physical identity. */
template <typename Specification>
requires std::same_as<Specification, model::FixedTotalMass>
struct SpecificationTraits<Specification> {
/** @brief Read the fixed physical target. */
[[nodiscard]] static std::array<double, 1> identity(const Specification& value) noexcept { return {value.target().value()}; }
/** @brief Retain unchanged prescribed physics. */
[[nodiscard]] static Specification clone(const Specification& value) { return value; }
};
}