From 9da779417bb54503fc74cd21115bbcd1f2b2b865 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 9 Nov 2026 13:27:14 -0500 Subject: [PATCH] feat(solver/nonlinear): Implements a damped newton line-search with backtracking This uses a linear solver from solver/linear which then calls some linear solver from the backend. Note that MFEM's build in newton solver is very simplistic and not sufficient for our needs; therefore we implement a more robust one here (and this is why newton is in the frontend and not the backend) --- src/include/serif/solver/nonlinear/newton.hpp | 119 +++++++++++ .../serif/solver/nonlinear/observers.hpp | 107 ++++++++++ .../serif/solver/nonlinear/options.hpp | 106 ++++++++++ .../serif/solver/nonlinear/qualification.hpp | 55 +++++ src/include/serif/solver/nonlinear/result.hpp | 181 +++++++++++++++++ src/include/serif/solver/nonlinear/steps.hpp | 188 ++++++++++++++++++ src/lib/serif/solver/nonlinear/options.cpp | 67 +++++++ 7 files changed, 823 insertions(+) create mode 100644 src/include/serif/solver/nonlinear/newton.hpp create mode 100644 src/include/serif/solver/nonlinear/observers.hpp create mode 100644 src/include/serif/solver/nonlinear/options.hpp create mode 100644 src/include/serif/solver/nonlinear/qualification.hpp create mode 100644 src/include/serif/solver/nonlinear/result.hpp create mode 100644 src/include/serif/solver/nonlinear/steps.hpp create mode 100644 src/lib/serif/solver/nonlinear/options.cpp diff --git a/src/include/serif/solver/nonlinear/newton.hpp b/src/include/serif/solver/nonlinear/newton.hpp new file mode 100644 index 00000000..94a71d46 --- /dev/null +++ b/src/include/serif/solver/nonlinear/newton.hpp @@ -0,0 +1,119 @@ +#pragma once + +#include "serif/solver/nonlinear/observers.hpp" +#include "serif/solver/nonlinear/steps.hpp" + +namespace serif::solver { + namespace detail { + /** @brief A custom method must supply the exact correction contract; malformed overrides never select a fallback. */ + template + concept NewtonStrategyFor = requires(Problem& problem, const State& state, Options& options, + const typename decltype(problem.prepare(state))::value_type& prepared, typename Problem::Vector& vector) { + { options.strategies.solve(problem, state, prepared, vector, vector, options) } -> std::same_as; + }; + + /** + * @brief Find equilibrium while retaining the last fully qualified physical state. + * @note Preparation belongs to the problem; this driver alone publishes accepted steps. + * Scratch vectors use the problem's native storage. Timing belongs entirely to observers. + */ + template + [[nodiscard]] auto run_newton(Problem& problem, State initial_state, Options options, Observer observer, + std::optional restart) { + using Vector = typename Problem::Vector; + using Prepared = typename decltype(problem.prepare(initial_state))::value_type; + using Result = SolveResult; + + validate_start(options, restart, initial_state.values().size()); + Result result(problem.model_owner()); + SolveDiagnostics diagnostics = restart ? restart->diagnostics : SolveDiagnostics{}; + SolveObservation observation(observer, options.strategies, diagnostics); + auto initial = evaluate_start(problem, initial_state, options, restart, diagnostics); + if (!initial) { + return observation.finish(result, initial.error().status, initial.error().message); + } + MetricEvaluation accepted_metric = initial->metric; + result.accept(std::move(initial_state), std::move(initial->system)); + const int size = static_cast(result.state().values().size()); + Vector correction(size), physical_direction(size); + Vector rhs(static_cast(result.system()->residual().size())); + std::ranges::fill(correction, 0.0); + double next_step = options.line_search.initial_step; + if (restart) { + std::ranges::copy(restart->correction, correction.begin()); + next_step = restart->next_step; + } + save_progress(result, options, diagnostics, accepted_metric, std::span(correction), next_step); + observation.notify(restart ? EventKind::resumed : EventKind::initial, result.state(), result.system(), + diagnostics.attempted_iterations, 0, 0.0, accepted_metric); + + // Check terminal states before preparing an inverse: a root may have a singular Jacobian. + std::optional prepared; + while (true) { + if (accepted_metric.residual_norm <= *diagnostics.convergence_threshold) { + return observation.finish(result, SolveStatus::converged, "The star satisfies the equilibrium tolerance.", accepted_metric); + } + if (observation.stop_requested()) { + return observation.finish(result, SolveStatus::stopped, "Newton stopped at an accepted boundary at the observer's request.", accepted_metric); + } + if (diagnostics.attempted_iterations >= options.method.maximum_iterations) { + return observation.finish(result, SolveStatus::maximum_iterations, "Newton reached its iteration budget; the result retains the last accepted star.", accepted_metric); + } + if (!prepared) { + auto preparation = problem.prepare(result.state()); + if (const auto failure = preparation_rejection(problem, preparation)) { + return observation.finish(result, SolveStatus::preparation_failure, *failure, accepted_metric); + } + prepared = std::move(*preparation); + } + const int iteration = diagnostics.attempted_iterations++; + observation.notify(EventKind::before_iteration, result.state(), result.system(), iteration, 0, 0.0, accepted_metric); + const auto direction = compute_direction(problem, result, *prepared, options, diagnostics, + rhs, correction, std::span(physical_direction)); + if (!direction) { + return observation.finish(result, direction.error().status, direction.error().message, accepted_metric); + } + auto trial = find_acceptable_step(problem, result, accepted_metric, correction, std::span(physical_direction), + *direction, next_step, options, diagnostics, observation); + if (!trial) { + return observation.finish(result, trial.error().status, trial.error().message, accepted_metric); + } + + // Nothing becomes accepted until the candidate and its nonterminal inverse are both usable. + observation.notify(EventKind::before_accept, trial->state, trial->system, iteration, trial->trial, trial->step, trial->metric); + result.accept(std::move(trial->state), std::move(trial->system)); + prepared = std::move(trial->preparation); + accepted_metric = trial->metric; + ++diagnostics.accepted_iterations; + diagnostics.final_residual_norm = accepted_metric.residual_norm; + diagnostics.last_accepted_step = trial->step; + for (double& value : correction) { + value *= 1.0 - trial->step; + } + if (!finite_correction(std::span(correction), size)) { + std::ranges::fill(correction, 0.0); + } + next_step = std::min(options.line_search.initial_step, trial->step / options.line_search.contraction); + save_progress(result, options, diagnostics, accepted_metric, std::span(correction), next_step); + observation.notify(EventKind::accepted, result.state(), result.system(), iteration, trial->trial, trial->step, accepted_metric); + } + } + } + + /** + * @brief Solve a custom numerical problem using its native Vector, evaluate, prepare and solve_linear operations. + * @note Strategies here are executable numerical methods; solve_equilibrium uses the same strategies. + */ + template , typename Observer = NoObserver> + requires detail::NewtonStrategyFor + [[nodiscard]] auto solve_newton(Problem& problem, State initial_state, Options options = {}, Observer observer = {}) { + return detail::run_newton(problem, std::move(initial_state), std::move(options), std::move(observer), std::nullopt); + } + + /** @brief Resume the original convergence reference, correction reuse, and cumulative work budget. */ + template , typename Observer = NoObserver> + requires detail::NewtonStrategyFor + [[nodiscard]] auto resume_newton(Problem& problem, NewtonCheckpoint checkpoint, Options options = {}, Observer observer = {}) { + return detail::run_newton(problem, std::move(checkpoint.state), std::move(options), std::move(observer), std::move(checkpoint.progress)); + } +} diff --git a/src/include/serif/solver/nonlinear/observers.hpp b/src/include/serif/solver/nonlinear/observers.hpp new file mode 100644 index 00000000..502d2d6a --- /dev/null +++ b/src/include/serif/solver/nonlinear/observers.hpp @@ -0,0 +1,107 @@ +#pragma once + +#include +#include +#include + +#include "serif/solver/nonlinear/result.hpp" + +namespace serif::solver { + /** @brief Optional elapsed wall time, kept separately from numerical restart history. */ + struct SolveTiming { + double seconds{0.0}; + }; + + /** + * @brief Measure solves explicitly while retaining an ordinary progress or stop callback. + * @note Reusing the same timing record accumulates resumed solves. Timing starts before initial + * evaluation and ends on any returned numerical result, including a rejected initial state. + * Exceptions propagate normally; an interrupted call contributes no completed-solve duration. + */ + template + class TimedObserver { + public: + /** @brief Borrow the timing destination and own the user's callback. */ + explicit TimedObserver(SolveTiming& timing, Observer observer = {}) : + m_timing(timing), m_observer(std::move(observer)) {} + + /** @brief Begin measurement when the solver starts evaluating its initial state. */ + void started() { + m_start = std::chrono::steady_clock::now(); + if constexpr (requires { m_observer.started(); }) { + m_observer.started(); + } + } + + /** @brief Accumulate a normally completed call independently of its convergence status. */ + void finished(SolveStatus status, const SolveDiagnostics& diagnostics) { + if constexpr (requires { m_observer.finished(status, diagnostics); }) { + m_observer.finished(status, diagnostics); + } + m_timing.seconds += std::chrono::duration(std::chrono::steady_clock::now() - m_start).count(); + } + + /** @brief Preserve the callback's void or stop-request return type. */ + decltype(auto) operator()(SolveEvent event, const auto& state, const auto& system) { + return std::invoke(m_observer, std::move(event), state, system); + } + + private: + SolveTiming& m_timing; + Observer m_observer; + std::chrono::steady_clock::time_point m_start; + }; + + namespace detail { + /** @brief Adapt callbacks and latch requests; only the Newton driver decides when stopping is safe. */ + template + class SolveObservation { + public: + /** @brief Borrow the single diagnostics record and the selected numerical execution. */ + SolveObservation(Observer& observer, Strategies& strategies, SolveDiagnostics& diagnostics) : + m_observer(observer), m_strategies(strategies), m_diagnostics(diagnostics) { + if constexpr (requires { observer.started(); }) { + observer.started(); + } + } + + /** @brief Publish an immutable snapshot and remember a stop request until an accepted boundary. */ + void notify(EventKind kind, const auto& state, const auto& system, int iteration = 0, + int trial = 0, double step = 0.0, std::optional metric = {}, std::string message = {}) { + m_strategies.contribute_diagnostics(m_diagnostics); + SolveEvent event{.kind = kind, .iteration = iteration, .trial = trial, .step = step, + .metric = metric, .message = std::move(message), .diagnostics = m_diagnostics}; + if constexpr (std::same_as, SolveControl>) { + m_stop_requested |= std::invoke(m_observer, std::move(event), state, system) == SolveControl::stop; + } else { + std::invoke(m_observer, std::move(event), state, system); + } + } + + /** @brief Check the latched request without letting callbacks accept or reject a trial. */ + [[nodiscard]] bool stop_requested() const noexcept { + return m_stop_requested; + } + + /** @brief Finalize reporting even when the initial state was rejected and no snapshot exists. */ + [[nodiscard]] auto finish(auto& result, SolveStatus status, std::string message, MetricEvaluation metric = {}) { + m_strategies.contribute_diagnostics(m_diagnostics); + result.finish(status, std::move(message), m_diagnostics); + if (result.has_state()) { + notify(EventKind::finished, result.state(), result.system(), m_diagnostics.attempted_iterations, + 0, 0.0, metric, result.message()); + } + if constexpr (requires { m_observer.finished(status, m_diagnostics); }) { + m_observer.finished(status, m_diagnostics); + } + return std::move(result); + } + + private: + Observer& m_observer; + Strategies& m_strategies; + SolveDiagnostics& m_diagnostics; + bool m_stop_requested{false}; + }; + } +} diff --git a/src/include/serif/solver/nonlinear/options.hpp b/src/include/serif/solver/nonlinear/options.hpp new file mode 100644 index 00000000..342b841d --- /dev/null +++ b/src/include/serif/solver/nonlinear/options.hpp @@ -0,0 +1,106 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "serif/solver/linear/linear.hpp" + +namespace serif::solver { + /** @brief Accuracy and work limits for finding one equilibrium star. */ + struct Newton { + double relative_tolerance{1.0e-8}; + double absolute_tolerance{0.0}; + int maximum_iterations{50}; + // Inexact Newton needs a reliable direction, not an over-solved early correction. + linear::Control linear_solve{.relative_tolerance = 0.03, .maximum_iterations = 200}; + + /** @brief Check tolerances and budgets before evaluating any candidate. */ + void validate() const; + /** @brief Freeze the convergence threshold from this solve's initial residual. */ + [[nodiscard]] double convergence_threshold(double initial_norm) const; + }; + + /** @brief Shorten a Newton step until the complete residual merit decreases sufficiently. */ + struct Backtracking { + double initial_step{1.0}; + double contraction{0.5}; + double fraction_to_boundary{0.9}; + double sufficient_decrease{1.0e-4}; + double minimum_step{1.0e-8}; + int maximum_trials{20}; + bool geometry_preflight{true}; + + /** @brief Require a decreasing, finite sequence of positive trial lengths. */ + void validate() const; + }; + + /** @brief Separate the convergence norm from the line-search merit for custom metrics. */ + struct MetricEvaluation { + double residual_norm{0.0}; + double merit{0.0}; + + /** @brief A metric must give finite, nonnegative values to drive solver decisions. */ + [[nodiscard]] bool valid() const noexcept { + return std::isfinite(residual_norm) && residual_norm >= 0.0 && std::isfinite(merit) && merit >= 0.0; + } + }; + + /** @brief Use the Euclidean norm of SERiF's normalized residual and half its square. */ + struct NormalizedResidualMetric { + /** @brief Evaluate serial solver coordinates without interpreting individual physical fields. */ + [[nodiscard]] MetricEvaluation operator()(std::span residual) const; + }; + + /** @brief A modified direction must report whether it is usable and any measured merit slope. */ + struct DirectionAssessment { + bool usable{true}; + std::optional merit_slope; + std::string message; + }; + + /** @brief A numerical policy explicitly accepts its work or explains why this candidate must be rejected. */ + struct TrialAssessment { + bool usable{true}; + std::string message; + }; + + struct SolveDiagnostics; + + /** @brief Ordinary Newton supplies every strategy phase; custom policies may replace individual phases explicitly. */ + struct StandardNewton { + /** @brief Ordinary Newton adds no optional stellar preparation computations or inputs. */ + static auto computations(const auto&, const auto&, auto, const auto&) { return std::pair{std::tuple{}, std::tuple{}}; } + /** @brief Ordinary Newton needs no modification of the evaluated initial state. */ + TrialAssessment initialize(auto&, auto&, auto&, const auto&) { return {}; } + /** @brief Physical preparation supplies the default sampled admissibility checks. */ + TrialAssessment admissible(auto&, const auto&) { return {}; } + /** @brief Select ordinary Newton's checked linear equation on the bound problem. */ + linear::Report solve(auto& problem, const auto&, const auto& prepared, const auto& rhs, auto& correction, const auto& options) { + return problem.solve_linear(prepared, rhs, correction, options); + } + /** @brief Retain the checked Newton direction and ordinary sufficient-decrease criterion. */ + DirectionAssessment direction(auto&, const auto&, const auto&, const auto&, auto&, const auto&) { return {}; } + /** @brief Retain the original evaluated trial without an auxiliary correction. */ + TrialAssessment trial(auto&, const auto&, const auto&, const auto&, double, auto&, auto&, const auto&, double) { return {}; } + /** @brief Ordinary Newton has no additional strategy work to report. */ + void contribute_diagnostics(SolveDiagnostics&) const {} + }; + + /** @brief Select numerical policies while keeping equilibrium physics in the model. */ + template + struct SolverOptions { + Newton method{}; + LinearSolver linear_solver{}; + Backtracking line_search{}; + Metric metric{}; + Strategies strategies{}; + }; + + /** @brief Skip only existing backtracking rungs when an inverted trial gives a finite boundary estimate. */ + [[nodiscard]] double next_backtracking_step(double rejected_step, double accepted_determinant, + double rejected_determinant, bool inverted, const Backtracking& options) noexcept; +} diff --git a/src/include/serif/solver/nonlinear/qualification.hpp b/src/include/serif/solver/nonlinear/qualification.hpp new file mode 100644 index 00000000..624b4e54 --- /dev/null +++ b/src/include/serif/solver/nonlinear/qualification.hpp @@ -0,0 +1,55 @@ +#pragma once + +#include +#include +#include +#include + +#include "serif/solver/nonlinear/result.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::solver::detail { + /** @brief Preserve numerical preparation failures while propagating invalid physical configurations. */ + [[nodiscard]] inline std::optional preparation_rejection(const auto& problem, const auto& evaluation) { + if (evaluation) { + return std::nullopt; + } + utils::errors::serif_ensure(problem.recoverable(evaluation.error()), + "Equilibrium preparation has an invalid configuration: {}", problem.failure_message(evaluation.error())); + return problem.failure_message(evaluation.error()); + } + + /** @brief Explain why a candidate cannot provide a usable nonlinear metric. */ + struct NewtonFailure { + SolveStatus status{SolveStatus::initial_state_rejected}; + std::string message; + }; + + /** @brief Qualify one already evaluated candidate through the same physical and metric gates everywhere. */ + [[nodiscard]] inline std::expected qualify_candidate(auto& problem, const auto& state, + const auto& evaluation, auto& metric_function, auto& strategy, TrialAssessment assessment = {}) { + if (const auto rejection = preparation_rejection(problem, evaluation)) { + return std::unexpected(NewtonFailure{.message = *rejection}); + } + if (assessment.usable) { + assessment = strategy.admissible(problem, state); + } + if (!assessment.usable) { + return std::unexpected(NewtonFailure{.message = std::move(assessment.message)}); + } + const MetricEvaluation metric = metric_function((*evaluation)->residual()); + if (!metric.valid()) { + return std::unexpected(NewtonFailure{.status = SolveStatus::invalid_metric, + .message = "The nonlinear metric is not finite and nonnegative."}); + } + return metric; + } + + /** @brief Every custom backend must preserve the correction dimensions and finite coefficients. */ + [[nodiscard]] inline bool finite_correction(std::span correction, const int expected_size) { + using namespace utils::errors; + serif_ensure(static_cast(correction.size()) == expected_size, + "The linear backend changed the Newton correction size from {} to {}. Return one correction for every solved coordinate.", expected_size, static_cast(correction.size())); + return std::ranges::all_of(correction, [](double value) { return std::isfinite(value); }); + } +} diff --git a/src/include/serif/solver/nonlinear/result.hpp b/src/include/serif/solver/nonlinear/result.hpp new file mode 100644 index 00000000..2fc7f8c1 --- /dev/null +++ b/src/include/serif/solver/nonlinear/result.hpp @@ -0,0 +1,181 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "serif/geometry/deformation/geometry_step.hpp" +#include "serif/solver/nonlinear/options.hpp" +#include "serif/solver/linear/linear.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::solver { + /** @brief Numerical termination; invalid configuration and user callback errors remain exceptions. */ + enum class SolveStatus { converged, initial_state_rejected, invalid_metric, preparation_failure, + linear_solve_failure, globalization_failure, maximum_iterations, stopped }; + + /** @brief Bounded work totals; accepted physical snapshots are stored separately from diagnostics. */ + struct SolveDiagnostics { + int attempted_iterations{0}; + int accepted_iterations{0}; + int line_search_trials{0}; + int rejected_trials{0}; + int geometry_limited_iterations{0}; + int total_linear_iterations{0}; + std::optional initial_residual_norm; + std::optional final_residual_norm; + std::optional convergence_threshold; + std::optional last_linear_solve; + std::optional last_geometry_preflight; + double last_accepted_step{0.0}; + }; + + namespace detail { + /** @brief Only built-in stateless numerical policies currently have a complete restart contract. */ + template + inline constexpr bool checkpointable_options = + std::same_as().strategies)>, StandardNewton> && + std::same_as().metric)>, NormalizedResidualMetric> && + std::same_as().linear_solver)>, linear::FGMRES>; + } + + /** + * @brief Numerical continuation data at one accepted Newton boundary. + * @note Correction values are normalized solver coordinates. Reuse the same physical problem, + * normalization, metric and numerical options on restart; the total iteration budget may increase. + * Runtime operators and caches are rebuilt, never stored here. + */ + struct NewtonProgress { + MetricEvaluation accepted_metric; + std::vector correction; + double next_step{1.0}; + SolveDiagnostics diagnostics; + }; + + /** @brief Reject incomplete or inconsistent imported history before preparing any runtime operators. */ + [[nodiscard]] inline std::optional validate_newton_progress(const NewtonProgress& progress, const std::size_t state_size) { + const SolveDiagnostics& diagnostics = progress.diagnostics; + const auto nonnegative = [](const std::optional& value) { + return value && std::isfinite(*value) && *value >= 0.0; + }; + if (!progress.accepted_metric.valid() || !std::isfinite(progress.next_step) || progress.next_step <= 0.0 || + progress.correction.size() != state_size || !std::ranges::all_of(progress.correction, [](double value) { return std::isfinite(value); })) { + return "The accepted metric, correction size/values or next trial step is invalid."; + } + if (!nonnegative(diagnostics.initial_residual_norm) || !nonnegative(diagnostics.final_residual_norm) || + !nonnegative(diagnostics.convergence_threshold) || *diagnostics.final_residual_norm != progress.accepted_metric.residual_norm) { + return "The original convergence reference or accepted residual is missing or inconsistent."; + } + if (diagnostics.accepted_iterations < 0 || diagnostics.attempted_iterations != diagnostics.accepted_iterations || + diagnostics.line_search_trials < diagnostics.accepted_iterations || diagnostics.rejected_trials < 0 || + diagnostics.line_search_trials - diagnostics.accepted_iterations != diagnostics.rejected_trials || + diagnostics.geometry_limited_iterations < 0 || diagnostics.geometry_limited_iterations > diagnostics.attempted_iterations || + diagnostics.total_linear_iterations < 0 || + !std::isfinite(diagnostics.last_accepted_step) || diagnostics.last_accepted_step < 0.0) { + return "The work counters do not describe an accepted Newton boundary."; + } + return std::nullopt; + } + + /** @brief Own the accepted physical star and the small amount of Newton history needed to continue it. */ + template + struct NewtonCheckpoint { + State state; + NewtonProgress progress; + }; + + /** @brief Events expose immutable accepted/trial snapshots; callbacks never decide acceptance. */ + enum class EventKind { initial, before_iteration, trial_rejected, before_accept, accepted, finished, resumed }; + + /** @brief A stop request is remembered until the next initial, resumed or accepted Newton boundary. */ + enum class SolveControl { continue_solve, stop }; + + /** @brief Small event record suitable for logging without retaining every iterate. */ + struct SolveEvent { + EventKind kind; + int iteration{0}; + int trial{0}; + double step{0.0}; + std::optional metric; + std::string message; + const SolveDiagnostics& diagnostics; + }; + + /** @brief Disable observation without introducing runtime callback allocation. */ + struct NoObserver { + /** @brief Accept the same immutable event and snapshot arguments as a user observer. */ + void operator()(const SolveEvent&, const auto&, const auto&) const noexcept {} + }; + + /** @brief Marker for algebraic tests and custom problems without a stellar description. */ + struct NoModel {}; + + /** + * @brief Own the last accepted physical state, matching frozen system, and descriptive model. + * @note A failed initial evaluation has no accepted state. A later failure retains the last + * accepted star, which remains usable after the solver and its dependency cache are destroyed. + */ + template + class SolveResult { + public: + /** @brief Retain the physical definition before any numerical attempt is made. */ + explicit SolveResult(std::shared_ptr model = {}) : m_model(std::move(model)) {} + + /** @brief True only after the complete nonlinear convergence test succeeds. */ + [[nodiscard]] bool converged() const noexcept { return m_status == SolveStatus::converged; } + /** @brief Distinguish a usable accepted checkpoint from a rejected initial guess. */ + [[nodiscard]] bool has_state() const noexcept { return m_state.has_value(); } + /** @brief Read the last accepted physical coefficients; reject access when none were evaluated. */ + [[nodiscard]] const State& state() const { + utils::errors::serif_ensure(has_state(), + "This solve has no accepted stellar state. Inspect status() and message() before requesting its state."); + return *m_state; + } + /** @brief Retain or inspect the exact residual and Jacobian belonging to state(). */ + [[nodiscard]] const std::shared_ptr& system() const noexcept { return m_system; } + /** @brief Read the physical definition owned by a stellar result. */ + [[nodiscard]] const Model& model() const requires (!std::same_as) { return *m_model; } + /** @brief Inspect work counts and final convergence measurements. */ + [[nodiscard]] const SolveDiagnostics& diagnostics() const noexcept { return m_diagnostics; } + /** @brief Inspect the numerical reason for termination. */ + [[nodiscard]] SolveStatus status() const noexcept { return m_status; } + /** @brief Read the diagnostic associated with numerical termination. */ + [[nodiscard]] const std::string& message() const noexcept { return m_message; } + + /** @brief Whether this accepted boundary has supported, complete numerical restart data. */ + [[nodiscard]] bool has_checkpoint() const noexcept { return m_progress.has_value(); } + /** @brief Copy a portable accepted boundary without retaining runtime operators or caches. */ + [[nodiscard]] NewtonCheckpoint checkpoint() const { + utils::errors::serif_ensure(has_state() && has_checkpoint(), + "This result has no supported Newton checkpoint. Checkpointing requires an accepted state with StandardNewton, NormalizedResidualMetric and FGMRES; custom policies need an explicit history restoration implementation."); + return {.state = *m_state, .progress = *m_progress}; + } + /** @brief Record history only after all accepted-step correction and step-length updates are complete. */ + void checkpoint(NewtonProgress progress) { m_progress = std::move(progress); } + + /** @brief Publish a complete state/system pair only after candidate preparation succeeds. */ + void accept(State state, std::shared_ptr system) { + m_state = std::move(state); + m_system = std::move(system); + } + /** @brief Finalize bounded diagnostics without altering the accepted checkpoint. */ + void finish(SolveStatus status, std::string message, SolveDiagnostics diagnostics) { + m_status = status; + m_message = std::move(message); + m_diagnostics = std::move(diagnostics); + } + + private: + std::shared_ptr m_model; + std::optional m_state; + std::shared_ptr m_system; + std::optional m_progress; + SolveStatus m_status{SolveStatus::initial_state_rejected}; + std::string m_message; + SolveDiagnostics m_diagnostics; + }; +} diff --git a/src/include/serif/solver/nonlinear/steps.hpp b/src/include/serif/solver/nonlinear/steps.hpp new file mode 100644 index 00000000..8a0f4f32 --- /dev/null +++ b/src/include/serif/solver/nonlinear/steps.hpp @@ -0,0 +1,188 @@ +#pragma once + +#include + +#include "serif/solver/nonlinear/qualification.hpp" + +namespace serif::solver::detail { + /** @brief Validate controls and restored history before any physical preparation. */ + void validate_start(const auto& options, const std::optional& restart, std::size_t size) { + using namespace utils::errors; + + options.method.validate(); + options.line_search.validate(); + if constexpr (requires { options.linear_solver.validate(); }) { + options.linear_solver.validate(); + } + if (!restart) { + return; + } + serif_ensure(checkpointable_options>, + "Newton restart requires StandardNewton, NormalizedResidualMetric and FGMRES. Custom policies require restoration of their numerical history."); + if (const auto error = validate_newton_progress(*restart, size)) { + serif_ensure(false, "Invalid Newton checkpoint: {}", *error); + } + serif_ensure(restart->next_step <= options.line_search.initial_step && + *restart->diagnostics.convergence_threshold == options.method.convergence_threshold(*restart->diagnostics.initial_residual_norm), + "The restart changes Newton's original tolerance or trial-step bound. Restore the archived numerical options before continuing."); + } + + /** @brief Keep a freshly qualified metric beside the exact system from which it was measured. */ + template + struct InitialEvaluation { + std::shared_ptr system; + MetricEvaluation metric; + }; + + /** @brief Qualify the initial state and verify a restart against its original convergence reference. */ + template + [[nodiscard]] std::expected, NewtonFailure> evaluate_start( + Problem& problem, auto& state, auto& options, const std::optional& restart, SolveDiagnostics& diagnostics + ) { + using namespace utils::errors; + + auto evaluation = problem.evaluate(state); + TrialAssessment initialized; + if (evaluation && !restart) { + initialized = options.strategies.initialize(problem, state, evaluation, options); + } + const auto metric = qualify_candidate(problem, state, evaluation, options.metric, options.strategies, initialized); + if (!metric) { + return std::unexpected(metric.error()); + } + if (restart) { + const double allowance = 256.0 * std::numeric_limits::epsilon() * std::max(1.0, restart->accepted_metric.residual_norm); + serif_ensure(std::abs(metric->residual_norm - restart->accepted_metric.residual_norm) <= allowance && + std::abs(metric->merit - restart->accepted_metric.merit) <= allowance * std::max(1.0, restart->accepted_metric.residual_norm), + "The rebuilt star does not reproduce the checkpoint's accepted residual metric. Restore the same model, discretization, normalization, prescribed inputs and metric."); + } else { + diagnostics.initial_residual_norm = metric->residual_norm; + diagnostics.convergence_threshold = options.method.convergence_threshold(metric->residual_norm); + } + diagnostics.final_residual_norm = metric->residual_norm; + return InitialEvaluation{.system = std::move(*evaluation), .metric = *metric}; + } + + /** @brief Solve and qualify a normalized correction before exposing its physical direction to backtracking. */ + [[nodiscard]] std::expected compute_direction( + auto& problem, const auto& accepted, const auto& prepared, auto& options, SolveDiagnostics& diagnostics, + auto& rhs, auto& correction, std::span physical_direction + ) { + const int size = static_cast(physical_direction.size()); + utils::errors::serif_ensure( + std::ranges::size(rhs) == accepted.system()->residual().size(), + "The residual changed size during Newton. Preserve the equation layout throughout one solve."); + std::ranges::transform(accepted.system()->residual(), rhs.begin(), [](double value) { return -value; }); + const linear::Report report = options.strategies.solve(problem, accepted.state(), prepared, rhs, correction, options); + diagnostics.last_linear_solve = report; + diagnostics.total_linear_iterations += report.iterations; + if (!finite_correction(std::span(correction), size) || !report.converged()) { + return std::unexpected(NewtonFailure{.status = SolveStatus::linear_solve_failure, + .message = "The Newton linear solve did not produce a finite correction meeting its true-residual tolerance."}); + } + DirectionAssessment direction = options.strategies.direction(problem, accepted.state(), prepared, rhs, correction, options); + if (!direction.usable || !finite_correction(std::span(correction), size) || + (direction.merit_slope && (!std::isfinite(*direction.merit_slope) || *direction.merit_slope >= 0.0))) { + return std::unexpected(NewtonFailure{.status = SolveStatus::globalization_failure, + .message = direction.message.empty() ? "The numerical direction policy found no finite decreasing direction." : direction.message}); + } + accepted.system()->coordinates().denormalize_state(std::span(correction), physical_direction); + return direction; + } + + /** @brief Own a qualified trial; a terminal root deliberately has no prepared inverse. */ + template + struct QualifiedStep { + State state; + std::shared_ptr system; + MetricEvaluation metric; + std::optional preparation; + double step; + int trial; + }; + + /** + * @brief Find a safe, sufficiently decreasing candidate without changing the accepted state. + * @note A nonterminal candidate must prepare its inverse before publication. Preparation failure + * terminates this attempt; the caller still owns the previously accepted physical star. + */ + template + [[nodiscard]] auto find_acceptable_step( + Problem& problem, const SolveResult& accepted, MetricEvaluation accepted_metric, + const auto& correction, std::span physical_direction, const DirectionAssessment& direction, + double step, auto& options, SolveDiagnostics& diagnostics, auto& observation + ) { + using Prepared = typename decltype(problem.prepare(accepted.state()))::value_type; + using Step = QualifiedStep; + using Result = std::expected; + + const int iteration = diagnostics.attempted_iterations - 1; + if constexpr (requires { problem.geometry_limit(accepted.state(), physical_direction, step, + options.line_search.fraction_to_boundary, options.line_search.geometry_preflight); }) { + const auto limit = problem.geometry_limit(accepted.state(), physical_direction, step, + options.line_search.fraction_to_boundary, options.line_search.geometry_preflight); + if (!limit) { + return Result(std::unexpected(NewtonFailure{.status = SolveStatus::globalization_failure, + .message = problem.failure_message(limit.error())})); + } + diagnostics.last_geometry_preflight = *limit; + if (options.line_search.geometry_preflight) { + step = limit->step; + diagnostics.geometry_limited_iterations += limit->limited_by_geometry; + } + } + std::string rejection = "Backtracking exhausted its trial or step-length budget."; + for (int trial = 0; trial < options.line_search.maximum_trials && step >= options.line_search.minimum_step; ++trial) { + ++diagnostics.line_search_trials; + State candidate = accepted.state(); + for (std::size_t index = 0; index < physical_direction.size(); ++index) { + candidate.values()[index] += step * physical_direction[index]; + } + double required_merit = (1.0 - options.line_search.sufficient_decrease * step) * accepted_metric.merit; + if (direction.merit_slope) { + required_merit = accepted_metric.merit + options.line_search.sufficient_decrease * step * *direction.merit_slope; + } + auto evaluation = problem.evaluate(candidate); + TrialAssessment assessment; + if (evaluation) { + assessment = options.strategies.trial(problem, accepted.state(), accepted.system(), correction, step, + candidate, evaluation, options, required_merit); + } + const auto metric = qualify_candidate(problem, candidate, evaluation, options.metric, options.strategies, std::move(assessment)); + if (metric && metric->merit <= required_merit) { + std::optional preparation; + if (metric->residual_norm > *diagnostics.convergence_threshold) { + auto prepared = problem.prepare(candidate); + if (const auto failure = preparation_rejection(problem, prepared)) { + return Result(std::unexpected(NewtonFailure{.status = SolveStatus::preparation_failure, .message = *failure})); + } + preparation = std::move(*prepared); + } + return Result(Step{.state = std::move(candidate), .system = std::move(*evaluation), .metric = *metric, + .preparation = std::move(preparation), .step = step, .trial = trial}); + } + rejection = metric ? "The trial did not sufficiently reduce the complete residual merit." : metric.error().message; + ++diagnostics.rejected_trials; + observation.notify(EventKind::trial_rejected, candidate, evaluation ? *evaluation : nullptr, iteration, trial, step, + metric ? std::optional(*metric) : std::nullopt, rejection); + double following = step * options.line_search.contraction; + if constexpr (requires { problem.next_step(step, options.line_search, evaluation); }) { + if (!evaluation) { + following = problem.next_step(step, options.line_search, evaluation); + } + } + step = following; + } + return Result(std::unexpected(NewtonFailure{.status = SolveStatus::globalization_failure, .message = std::move(rejection)})); + } + + /** @brief Snapshot only fully committed numerical history; reporting has no influence on restart. */ + void save_progress(auto& result, auto& options, SolveDiagnostics& diagnostics, MetricEvaluation metric, + std::span correction, double next_step) { + if constexpr (checkpointable_options>) { + options.strategies.contribute_diagnostics(diagnostics); + result.checkpoint(NewtonProgress{.accepted_metric = metric, .correction = {correction.begin(), correction.end()}, + .next_step = next_step, .diagnostics = diagnostics}); + } + } +} diff --git a/src/lib/serif/solver/nonlinear/options.cpp b/src/lib/serif/solver/nonlinear/options.cpp new file mode 100644 index 00000000..716c2a8c --- /dev/null +++ b/src/lib/serif/solver/nonlinear/options.cpp @@ -0,0 +1,67 @@ +#include "serif/solver/nonlinear/options.hpp" + +#include +#include "serif/utils/error/errors.hpp" + +namespace serif::solver { + void Newton::validate() const { + using namespace utils::errors; + serif_ensure(std::isfinite(relative_tolerance) && relative_tolerance >= 0.0 && + std::isfinite(absolute_tolerance) && absolute_tolerance >= 0.0 && maximum_iterations >= 0, + "Newton tolerances must be finite and nonnegative, and its iteration budget must be nonnegative."); + linear_solve.validate(); + } + + double Newton::convergence_threshold(const double initial_norm) const { + using namespace utils::errors; + validate(); + const double relative = relative_tolerance * initial_norm; + serif_ensure(std::isfinite(initial_norm) && initial_norm >= 0.0 && std::isfinite(relative), + "Newton requires a finite initial norm and representable relative convergence threshold."); + return std::max(absolute_tolerance, relative); + } + + void Backtracking::validate() const { + using namespace utils::errors; + serif_ensure(std::isfinite(initial_step) && initial_step > 0.0 && initial_step <= 1.0 && + std::isfinite(contraction) && contraction > 0.0 && contraction < 1.0 && + std::isfinite(fraction_to_boundary) && fraction_to_boundary > 0.0 && fraction_to_boundary < 1.0 && + std::isfinite(sufficient_decrease) && sufficient_decrease > 0.0 && sufficient_decrease < 1.0 && + std::isfinite(minimum_step) && minimum_step > 0.0 && minimum_step <= initial_step && maximum_trials > 0, + "Backtracking needs 0 < minimum step <= initial step <= 1, contraction, boundary safety and sufficient decrease in (0,1), and a positive trial budget."); + } + + MetricEvaluation NormalizedResidualMetric::operator()(std::span residual) const { + double norm = 0.0; + for (const double value : residual) { + norm = std::hypot(norm, value); + } + return {.residual_norm = norm, .merit = 0.5 * norm * norm}; + } + + double next_backtracking_step(const double rejected_step, const double accepted_determinant, + const double rejected_determinant, const bool inverted, const Backtracking& options) noexcept { + const double contracted = rejected_step * options.contraction; + if (!inverted || !std::isfinite(accepted_determinant) || accepted_determinant <= 0.0 || + !std::isfinite(rejected_determinant) || rejected_determinant > 0.0) { + return contracted; + } + const double change = accepted_determinant - rejected_determinant; + if (!std::isfinite(change) || change <= 0.0) { + return contracted; + } + const double safe_step = options.fraction_to_boundary * rejected_step * accepted_determinant / change; + if (!std::isfinite(safe_step) || safe_step <= 0.0 || contracted <= safe_step) { + return contracted; + } + const double rung = std::ceil(std::log(safe_step / rejected_step) / std::log(options.contraction)); + double skipped = rejected_step * std::pow(options.contraction, rung); + if (!std::isfinite(skipped) || skipped <= 0.0 || skipped >= rejected_step) { + return contracted; + } + if (skipped > safe_step) { + skipped *= options.contraction; + } + return skipped; + } +}