From 64cc306a2ccb28aa29f3fa074c9ac372d0b8fd1e Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 9 Nov 2026 13:24:42 -0500 Subject: [PATCH] feat(solver/continuation): system to auto continue a solver along some number of parameters 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 --- .../solver/continuation/continuation.hpp | 6 + .../serif/solver/continuation/controller.hpp | 255 ++++++++++++ .../serif/solver/continuation/driver.hpp | 232 +++++++++++ .../serif/solver/continuation/equilibrium.hpp | 88 +++++ .../serif/solver/continuation/family.hpp | 372 ++++++++++++++++++ .../continuation/identified_equations.hpp | 39 ++ .../serif/solver/continuation/options.hpp | 106 +++++ .../serif/solver/continuation/parameters.hpp | 104 +++++ .../serif/solver/continuation/polytropic.hpp | 46 +++ .../continuation/stellar_parameters.hpp | 92 +++++ 10 files changed, 1340 insertions(+) create mode 100644 src/include/serif/solver/continuation/continuation.hpp create mode 100644 src/include/serif/solver/continuation/controller.hpp create mode 100644 src/include/serif/solver/continuation/driver.hpp create mode 100644 src/include/serif/solver/continuation/equilibrium.hpp create mode 100644 src/include/serif/solver/continuation/family.hpp create mode 100644 src/include/serif/solver/continuation/identified_equations.hpp create mode 100644 src/include/serif/solver/continuation/options.hpp create mode 100644 src/include/serif/solver/continuation/parameters.hpp create mode 100644 src/include/serif/solver/continuation/polytropic.hpp create mode 100644 src/include/serif/solver/continuation/stellar_parameters.hpp diff --git a/src/include/serif/solver/continuation/continuation.hpp b/src/include/serif/solver/continuation/continuation.hpp new file mode 100644 index 00000000..1c1e360a --- /dev/null +++ b/src/include/serif/solver/continuation/continuation.hpp @@ -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" diff --git a/src/include/serif/solver/continuation/controller.hpp b/src/include/serif/solver/continuation/controller.hpp new file mode 100644 index 00000000..19b4831e --- /dev/null +++ b/src/include/serif/solver/continuation/controller.hpp @@ -0,0 +1,255 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#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 + struct CoordinateOptions final { + static_assert(Count > 0, "Continuation requires at least one physical coordinate."); + std::array initial_step_fraction = filled(0.25); + std::array minimum_step_fraction = filled(1.0e-6); + std::array maximum_step_fraction = filled(1.0); + std::array axis_order = []() { + std::array 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 filled(const double value) { + std::array result{}; + result.fill(value); + return result; + } + }; + + template + struct CoordinateProposal final { + std::uint64_t id{}; + std::size_t axis{}; + std::array 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 + class CoordinateController final { + public: + using Coordinates = std::array; + using Proposal = CoordinateProposal; + + /** @brief Validate finite endpoints and controls before any proposal can consume work. */ + CoordinateController(Coordinates start, Coordinates target, CoordinateOptions 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( + 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& 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& pending() const noexcept { return m_pending; } + + /** @brief Propose one bounded representable change without modifying accepted coordinates. */ + [[nodiscard]] std::optional propose() { + utils::errors::serif_ensure(!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 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(m_pending && m_status == ControllerStatus::pending && m_pending->id == id, + "Continuation resolution must identify its outstanding proposal."); + serif_ensure( + outcome == AttemptOutcome::success || outcome == AttemptOutcome::recoverable_rejection || outcome == AttemptOutcome::fatal_failure, + "Unknown continuation attempt outcome."); + serif_ensure(work <= std::numeric_limits::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( + options.order == CoordinateOrder::automatic || options.order == CoordinateOrder::preferred || options.order == CoordinateOrder::strict, + "Unknown continuation coordinate ordering."); + serif_ensure( + 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 seen{}; + for (std::size_t axis = 0; axis < Count; ++axis) { + serif_ensure(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( + 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 select_axis() const noexcept { + std::optional 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{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 m_options; + Coordinates m_accepted{}, m_span{}, m_step{}; + std::array m_blocked{}; + std::array m_rejections_at_base{}; + CoordinateCounters m_counters{}; + ControllerStatus m_status{ControllerStatus::ready}; + std::optional m_pending; + }; +} diff --git a/src/include/serif/solver/continuation/driver.hpp b/src/include/serif/solver/continuation/driver.hpp new file mode 100644 index 00000000..f520b8cd --- /dev/null +++ b/src/include/serif/solver/continuation/driver.hpp @@ -0,0 +1,232 @@ +#pragma once + +#include +#include + +#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 + struct StageRecord final { + std::array coordinates{}; + std::optional changed_axis; + bool starting_star{false}; + StageDisposition disposition{StageDisposition::not_converged}; + std::uint64_t work{}; + std::optional residual_norm; + std::string message; + }; + + /** @brief Retain the last converged physical state separately from the requested endpoint and failed candidates. */ + template + struct Result final { + Status status{Status::starting_star_failed}; + std::string message; + std::array target{}; + std::optional> accepted_coordinates; + std::unique_ptr accepted; + std::unique_ptr failed; + std::vector> 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 + struct StartingStar final { + Model model; + State state; + }; + + namespace detail { + /** @brief Try bounded simpler starting models, unless the user explicitly fixed the starting coordinates. */ + template + [[nodiscard]] std::vector starting_candidates(const Family& family, const Options& options) { + std::vector 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 + void record(Result& result, StageRecord 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 + [[nodiscard]] auto run(const Family& family, const State& initial_state, Corrector&& correct, Options options, double absolute_tolerance) { + using namespace utils::errors; + using Report = std::remove_cvref_t>; + using Record = StageRecord; + constexpr std::size_t count = Family::coordinate_count; + serif_ensure(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 validation(options.start, family.target_coordinates(), options.steps); + Result result{.target = family.target_coordinates(), .absolute_tolerance = absolute_tolerance}; + result.stages.reserve(std::min(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 axis) { + Record stage{.coordinates = coordinates, .changed_axis = axis, .starting_star = starting}; + auto model = family.try_model_at(coordinates); + std::unique_ptr 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(std::invoke(correct, std::move(*model), State(state), CorrectionControls{absolute_tolerance, remaining})); + stage.work = static_cast(report->diagnostics().attempted_iterations); + serif_ensure(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 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 + [[nodiscard]] auto run(const Family& family, const StartingStar& initial, Corrector&& correct, + Options options, double absolute_tolerance) { + utils::errors::serif_ensure(!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(correct), std::move(options), absolute_tolerance); + } +} diff --git a/src/include/serif/solver/continuation/equilibrium.hpp b/src/include/serif/solver/continuation/equilibrium.hpp new file mode 100644 index 00000000..7aeb0863 --- /dev/null +++ b/src/include/serif/solver/continuation/equilibrium.hpp @@ -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 + [[nodiscard]] auto corrector(const Family& family, Factory& configure, SolverOptions options, Solve solve) { + options.method.validate(); + utils::errors::serif_ensure( + 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; + using FiniteElements = std::remove_cvref_t().finite_elements)>; + using TransferGuard = decltype(backend::mfem::physical_transfer_guard(std::declval())); + return [&family, &configure, options = std::move(options), solve = std::move(solve), transfer = std::optional{}] + (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( + 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( + (*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(std::min( + 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 > + requires std::copy_constructible + [[nodiscard]] auto continue_equilibrium(const Family& family, Factory configure, const State& initial_state, + continuation::Options 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 + [[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 > + requires std::copy_constructible + [[nodiscard]] auto continue_equilibrium(const Family& family, Factory configure, const State& initial_state, const Vector& displacement, + continuation::Options 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 + [[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)); + } +} diff --git a/src/include/serif/solver/continuation/family.hpp b/src/include/serif/solver/continuation/family.hpp new file mode 100644 index 00000000..9a3a6ee3 --- /dev/null +++ b/src/include/serif/solver/continuation/family.hpp @@ -0,0 +1,372 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#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 + [[nodiscard]] constexpr auto depends_on(Parameters...) noexcept { return ParameterList{}; } + + template + 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 + [[nodiscard]] bool callable_valid(const Function& value) { + if constexpr (std::is_pointer_v || std::is_member_pointer_v) { + 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::epsilon() * std::max(std::abs(first), std::abs(second))); + } + + template + struct Contained; + template + struct Contained, Superset> : std::bool_constant<(utils::dependencies::ContainsKey && ...)> {}; + + template + concept LawDefinition = requires { + typename Candidate::OutputParameter; + typename Candidate::DependencyParameters; + requires PhysicalParameter; + requires utils::dependencies::IsKeyList; + }; + + template + struct FamilyTraits; + template + struct FamilyTraits, std::tuple> { + using Derived = ParameterList; + static constexpr bool valid = sizeof...(Selected) > 0 && utils::dependencies::IsKeyList> && + utils::dependencies::IsKeyList && + (!utils::dependencies::ContainsKey> && ...) && + (Contained>::value && ...); + static constexpr bool complete = valid && (Contained::Type, Derived>::value && ...); + template + static constexpr bool owners_present = (Model::template contains_specification && ...) && + (Model::template contains_specification && ...); + template + static constexpr bool affects = (std::same_as || ...) || + (std::same_as || ...); + }; + + template + struct AllIdentified; + template + struct AllIdentified> : std::bool_constant<(IdentifiedSpecification && ...)> {}; + + template + concept BindableFamily = model::StellarModelType && FamilyTraits::complete && + FamilyTraits::template owners_present && IdentifiedSpecification && + AllIdentified::value && + (std::is_empty_v || IdentifiedSpecification); + + struct CopyEquations final { + /** @brief A copyable equation definition needs no separate construction policy. */ + template + requires std::copy_constructible + [[nodiscard]] Equations operator()(const Equations& value) const { return value; } + }; + + /** @brief Bind target-dependent laws once, owning all resulting captures. */ + template + [[nodiscard]] auto bind_law(Law&& value, const Model& target) { + if constexpr (requires { std::forward(value).bind(target); }) { + return std::forward(value).bind(target); + } else { + return std::remove_cvref_t(std::forward(value)); + } + } + } + + /** @brief Own a typed law, including move-only captures; its inputs are selected coordinates rather than other laws. */ + template + requires utils::dependencies::IsKeyList> && + std::invocable&, typename Inputs::Value...> && + std::same_as&, typename Inputs::Value...>>, typename Output::Value> + [[nodiscard]] auto law(Output, ParameterList, Function&& function) { + utils::errors::serif_ensure(detail::callable_valid(function), + "A continuation law cannot use a null callable pointer."); + return ParameterLaw, std::remove_cvref_t>{.function = std::forward(function)}; + } + + template + class BoundFamily; + + /** @brief Own target physics and construct fresh family members without a mesh, evaluator, or mutable solve context. */ + template + class BoundFamily, std::tuple, EquationFactory> final { + using Selection = ParameterList; + using Traits = detail::FamilyTraits>; + using DerivedValues = std::tuple; + public: + using ModelType = Model; + static constexpr std::size_t coordinate_count = sizeof...(Selected); + using Coordinates = std::array; + template + static constexpr bool has_coordinate = utils::dependencies::ContainsKey; + + /** @brief Locate one selected physical parameter without runtime type dispatch. */ + template + requires has_coordinate + [[nodiscard]] static consteval std::size_t coordinate_index(Parameter = {}) { + return utils::dependencies::key_index(); + } + + /** @brief Expose stable parameter names for diagnostics and persisted coordinates. */ + [[nodiscard]] static constexpr std::array 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, 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>); + const std::vector identity = model_signature(m_target); + const bool finite = std::ranges::all_of(identity, [](double value) { return std::isfinite(value); }); + utils::errors::serif_ensure(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(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 = [&]() { + auto value = Parameter::read(m_target.template specification()); + if constexpr (requires { { Parameter::starting_value(value) } -> std::same_as; }) { + value = Parameter::starting_value(value); + } + utils::errors::serif_ensure(Parameter::valid(value), + "The default starting {} lies outside its physical domain.", Parameter::name); + return coordinate_value(value); + }; + return {start.template operator()()...}; + } + + /** @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()))...}; + } + + /** @brief Include unchanged EOS, constraints, and stateful equation definitions in family identity. */ + [[nodiscard]] std::vector model_signature(const Model& candidate) const { + std::vector signature; + const auto append = [&](const auto& value) { + const auto identity = SpecificationTraits>::identity(value); + signature.push_back(static_cast(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) { + append(candidate.equations()); + } + return signature; + } + + /** @brief Reject unrelated physical anchors even if their selected coordinates match. */ + [[nodiscard]] std::expected 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(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 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{}); + 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; + if constexpr (Traits::template affects) { + if (!endpoint) { + return SpecificationTraits::rebuild(specification, values); + } + } + return SpecificationTraits::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{})) { + 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(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 + [[nodiscard]] typename Parameter::Value get(Parameter) const { + if constexpr (has_coordinate) { + return typename Parameter::Value{coordinates[coordinate_index(Parameter{})]}; + } else if constexpr (utils::dependencies::ContainsKey) { + return std::get()>(derived); + } else { + return Parameter::read(target.template specification()); + } + } + }; + + /** @brief Read only independent coordinates into a typed dependent law. */ + template + [[nodiscard]] auto evaluate_law(const Law& value, const Coordinates& coordinates, ParameterList) const { + using Output = typename Law::OutputParameter; + using Result = std::expected; + 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 + [[nodiscard]] std::expected evaluate_laws(const Coordinates& coordinates, std::index_sequence) const { + auto results = std::tuple{evaluate_law(std::get(m_laws), coordinates, typename Laws::DependencyParameters{})...}; + std::optional failure; + const auto check = [&](const auto& result) { + if (!result && !failure) { + failure = result.error(); + } + }; + (check(std::get(results)), ...); + if (failure) { + return std::unexpected(std::move(*failure)); + } + return DerivedValues{std::move(*std::get(results))...}; + } + + /** @brief Read back every installed parameter; exact selected coordinates cannot be silently ignored. */ + template + [[nodiscard]] bool realized(const Model& result, const Coordinates& coordinates, const DerivedValues& derived, std::index_sequence) const { + return ((coordinate_value(Selected::read(result.template specification())) == coordinates[coordinate_index(Selected{})]) && ...) && + (detail::values_agree(coordinate_value(Laws::OutputParameter::read(result.template specification())), + coordinate_value(std::get(derived))) && ...); + } + + Model m_target; + std::tuple m_laws; + EquationFactory m_equations; + Coordinates m_target_coordinates; + }; + + /** @brief Own a reusable physical family definition before selecting a target model. */ + template + 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 + requires detail::FamilyTraits(), std::declval>()))>::valid + [[nodiscard]] auto with(Added... laws) const& requires std::copy_constructible { + return Family{std::move(laws)...}))>{ + std::tuple_cat(m_laws, std::tuple{std::move(laws)...})}; + } + + /** @brief Transfer move-only law captures when extending a temporary definition. */ + template + requires detail::FamilyTraits(), std::declval>()))>::valid + [[nodiscard]] auto with(Added... laws) && { + auto joined = std::tuple_cat(std::move(m_laws), std::tuple{std::move(laws)...}); + return Family{std::move(joined)}; + } + + /** @brief Bind a reproducible model, using an explicit factory for move-only equation definitions when necessary. */ + template + requires detail::BindableFamily && std::copy_constructible && + std::invocable + [[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 + requires detail::BindableFamily && std::invocable + [[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 + [[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(value), target)...}; + }, std::forward(values)); + return BoundFamily{std::move(target), std::move(bound), std::move(factory)}; + } + Laws m_laws; + }; + + /** @brief Select independent physical coordinates with duplicate detection at compile time. */ + template + requires (sizeof...(Parameters) > 0) && utils::dependencies::IsKeyList> + [[nodiscard]] auto in(Parameters...) { return Family, std::tuple<>>{}; } +} diff --git a/src/include/serif/solver/continuation/identified_equations.hpp b/src/include/serif/solver/continuation/identified_equations.hpp new file mode 100644 index 00000000..6d1af027 --- /dev/null +++ b/src/include/serif/solver/continuation/identified_equations.hpp @@ -0,0 +1,39 @@ +#pragma once + +#include +#include +#include +#include + +#include "serif/solver/continuation/parameters.hpp" + +namespace serif::solver::continuation { + /** @brief Give a deferred equation definition an explicit physical identity without runtime polymorphism. */ + template + struct IdentifiedEquations final { + Definition definition; + std::tuple 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{coordinate_value(values)...}; }, identity); + } + + /** @brief Rebuild equations from the stage's own model while retaining the supplied callable. */ + template + requires std::invocable + [[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 + [[nodiscard]] auto identified_equations(Definition definition, Values... values) { + return IdentifiedEquations{std::move(definition), std::tuple{std::move(values)...}}; + } +} diff --git a/src/include/serif/solver/continuation/options.hpp b/src/include/serif/solver/continuation/options.hpp new file mode 100644 index 00000000..994057ca --- /dev/null +++ b/src/include/serif/solver/continuation/options.hpp @@ -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 + struct Options final { + CoordinateOptions steps{}; + typename Family::Coordinates start{}; + std::array 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 + struct StartingAt final { + typename Parameter::Value value; + /** @brief Change one typed starting coordinate with strong rejection of duplicate or invalid input. */ + template + requires Family::template has_coordinate + void apply(Options& options) const { + constexpr std::size_t axis = Family::coordinate_index(Parameter{}); + utils::errors::serif_ensure(!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 + [[nodiscard]] auto starting_at(Parameter, typename Parameter::Value value) { return StartingAt{value}; } + + template + requires (sizeof...(Parameters) > 0) && utils::dependencies::IsKeyList> + struct Ordering final { + /** @brief A preference may yield to other axes; strict stages require every coordinate exactly once. */ + template + requires (Family::template has_coordinate && ...) && (!Strict || sizeof...(Parameters) == Family::coordinate_count) + void apply(Options& options) const { + utils::errors::serif_ensure(!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 + requires (sizeof...(Parameters) > 0) && utils::dependencies::IsKeyList> + [[nodiscard]] auto prefer(Parameters...) { return Ordering{}; } + + /** @brief Complete each physical stage before attempting the next. */ + template + requires (sizeof...(Parameters) > 0) && utils::dependencies::IsKeyList> + [[nodiscard]] auto stages(Parameters...) { return Ordering{}; } + + template + struct NumericalControls final { + CoordinateOptions options; + /** @brief Preserve explicit physical ordering independently of option argument order. */ + template + requires (Family::coordinate_count == Count) + void apply(Options& target) const { + utils::errors::serif_ensure(!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 + [[nodiscard]] auto numerics(CoordinateOptions options) { return NumericalControls{options}; } + + /** @brief Resolve physical settings from family defaults before any corrector runs. */ + template + requires (requires(const Settings& setting, Options& options) { setting.apply(options); } && ...) + [[nodiscard]] auto make_options(const Family& family, const Settings&... settings) { + Options result{.start = family.starting_coordinates()}; + (settings.apply(result), ...); + return result; + } +} diff --git a/src/include/serif/solver/continuation/parameters.hpp b/src/include/serif/solver/continuation/parameters.hpp new file mode 100644 index 00000000..04843fba --- /dev/null +++ b/src/include/serif/solver/continuation/parameters.hpp @@ -0,0 +1,104 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "serif/dimensions/quantities.hpp" +#include "serif/utils/dependencies/key_list_operations.hpp" + +namespace serif::solver::continuation { + template + using ParameterList = utils::dependencies::KeyList; + + 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 + 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 + [[nodiscard]] static Specification rebuild(const Specification& value, const Values& parameters) + requires requires { { value.with_continuation_parameters(parameters) } -> std::same_as; } { + 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; }) { + return value.clone_for_continuation(); + } else { + return value; + } + } + }; + + template + concept CoordinateValue = std::same_as || dimensions::QuantityValueType; + + /** @brief Convert physical units only at the coordinate-controller boundary. */ + template + [[nodiscard]] constexpr double coordinate_value(Value value) noexcept { + if constexpr (std::same_as) { + return value; + } else { + return value.value(); + } + } + + namespace detail { + template + struct ParameterCatalog { using Type = ParameterList<>; }; + template + struct ParameterCatalog> { using Type = typename Owner::ContinuationParameters; }; + + template + struct Catalog : ParameterCatalog {}; + template + struct Catalog::Parameters>> { using Type = typename SpecificationTraits::Parameters; }; + + template + struct RequiredDependents { using Type = ParameterList<>; }; + template + struct RequiredDependents> { + using Type = typename Parameter::RequiredDependentParameters; + }; + } + + /** @brief A typed physical parameter owns its domain predicate, units, owner, and stable identity. */ + template + concept PhysicalParameter = requires(const typename Parameter::Specification& specification, typename Parameter::Value value) { + requires CoordinateValue; + { Parameter::name } -> std::convertible_to; + { Parameter::stable_id } -> std::convertible_to; + { Parameter::read(specification) } -> std::same_as; + { Parameter::valid(value) } -> std::same_as; + requires utils::dependencies::ContainsKey::Type>; + }; + + /** @brief Require explicit identity data, never an object address or incidental memory representation. */ + template + concept IdentifiedSpecification = requires(const Specification& value) { + SpecificationTraits::identity(value); + requires std::ranges::forward_range::identity(value))>; + requires std::ranges::sized_range::identity(value))>; + requires std::same_as::identity(value))>, double>; + }; +} diff --git a/src/include/serif/solver/continuation/polytropic.hpp b/src/include/serif/solver/continuation/polytropic.hpp new file mode 100644 index 00000000..3d6d2897 --- /dev/null +++ b/src/include/serif/solver/continuation/polytropic.hpp @@ -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; + + /** @brief Validate a positive finite thermodynamic pressure anchor. */ + explicit PreservePressureAtReferenceDensity(dimensions::DensityValue density) : m_density(density) { + utils::errors::serif_ensure( + 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 + requires Model::template contains_specification + [[nodiscard]] auto bind(const Model& target) const { + const auto& eos = target.template specification(); + 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( + 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}; + } +} diff --git a/src/include/serif/solver/continuation/stellar_parameters.hpp b/src/include/serif/solver/continuation/stellar_parameters.hpp new file mode 100644 index 00000000..a25e5ed8 --- /dev/null +++ b/src/include/serif/solver/continuation/stellar_parameters.hpp @@ -0,0 +1,92 @@ +#pragma once + +#include + +#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; + 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 { + using Parameters = ParameterList; + /** @brief Include both EOS parameters in family membership. */ + [[nodiscard]] static std::array 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 + [[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 { + using Parameters = ParameterList; + /** @brief Rotation axes and centers are immutable family inputs as well as angular momentum. */ + [[nodiscard]] static std::array 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 + [[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 + requires std::same_as + struct SpecificationTraits { + /** @brief Read the fixed physical target. */ + [[nodiscard]] static std::array identity(const Specification& value) noexcept { return {value.target().value()}; } + /** @brief Retain unchanged prescribed physics. */ + [[nodiscard]] static Specification clone(const Specification& value) { return value; } + }; +}