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

feat(solver/nonlinear): Implements a damped newton line-search with backtracking

Browse Source 
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)
This commit is contained in:
Emily Boudreaux
2026-11-09 13:27:14 -05:00
parent f520fe5f03
commit 9da779417b
7 changed files with 823 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

119
src/include/serif/solver/nonlinear/newton.hpp Normal file
View File

@@ -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 <typename Problem, typename State, typename Options>
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<linear::Report>;
};
/**
* @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 <typename Problem, typename State, typename Options, typename Observer>
[[nodiscard]] auto run_newton(Problem& problem, State initial_state, Options options, Observer observer,
std::optional<NewtonProgress> restart) {
using Vector = typename Problem::Vector;
using Prepared = typename decltype(problem.prepare(initial_state))::value_type;
using Result = SolveResult<State, typename Problem::System, typename Problem::Model>;
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<int>(result.state().values().size());
Vector correction(size), physical_direction(size);
Vector rhs(static_cast<int>(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> 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 Problem, typename State, typename Options = SolverOptions<>, typename Observer = NoObserver>
requires detail::NewtonStrategyFor<Problem, State, Options>
[[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 Problem, typename State, typename Options = SolverOptions<>, typename Observer = NoObserver>
requires detail::NewtonStrategyFor<Problem, State, Options>
[[nodiscard]] auto resume_newton(Problem& problem, NewtonCheckpoint<State> checkpoint, Options options = {}, Observer observer = {}) {
return detail::run_newton(problem, std::move(checkpoint.state), std::move(options), std::move(observer), std::move(checkpoint.progress));
}
}

107
src/include/serif/solver/nonlinear/observers.hpp Normal file
View File

@@ -0,0 +1,107 @@
#pragma once
#include <chrono>
#include <functional>
#include <type_traits>
#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 <typename Observer = NoObserver>
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<double>(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 <typename Observer, typename Strategies>
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<MetricEvaluation> 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<std::invoke_result_t<Observer&, SolveEvent, decltype(state), decltype(system)>, 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};
};
}
}

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

@@ -0,0 +1,106 @@
#pragma once
#include <cmath>
#include <span>
#include <optional>
#include <string>
#include <tuple>
#include <utility>
#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<const double> residual) const;
};
/** @brief A modified direction must report whether it is usable and any measured merit slope. */
struct DirectionAssessment {
bool usable{true};
std::optional<double> 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 <typename Metric = NormalizedResidualMetric, typename Strategies = StandardNewton, typename LinearSolver = linear::FGMRES>
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;
}

55
src/include/serif/solver/nonlinear/qualification.hpp Normal file
View File

@@ -0,0 +1,55 @@
#pragma once
#include <algorithm>
#include <cmath>
#include <expected>
#include <span>
#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<std::string> preparation_rejection(const auto& problem, const auto& evaluation) {
if (evaluation) {
return std::nullopt;
}
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT>(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<MetricEvaluation, NewtonFailure> 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<const double> correction, const int expected_size) {
using namespace utils::errors;
serif_ensure<SERiFErrorCode::INVALID_DOF_COUNT>(static_cast<int>(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<int>(correction.size()));
return std::ranges::all_of(correction, [](double value) { return std::isfinite(value); });
}
}

181
src/include/serif/solver/nonlinear/result.hpp Normal file
View File

@@ -0,0 +1,181 @@
#pragma once
#include <algorithm>
#include <cmath>
#include <memory>
#include <optional>
#include <string>
#include <utility>
#include <vector>
#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<double> initial_residual_norm;
std::optional<double> final_residual_norm;
std::optional<double> convergence_threshold;
std::optional<linear::Report> last_linear_solve;
std::optional<geometry::deformation::GeometryStepLimit> 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 <typename Options>
inline constexpr bool checkpointable_options =
std::same_as<std::remove_cvref_t<decltype(std::declval<Options>().strategies)>, StandardNewton> &&
std::same_as<std::remove_cvref_t<decltype(std::declval<Options>().metric)>, NormalizedResidualMetric> &&
std::same_as<std::remove_cvref_t<decltype(std::declval<Options>().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<double> correction;
double next_step{1.0};
SolveDiagnostics diagnostics;
};
/** @brief Reject incomplete or inconsistent imported history before preparing any runtime operators. */
[[nodiscard]] inline std::optional<std::string> validate_newton_progress(const NewtonProgress& progress, const std::size_t state_size) {
const SolveDiagnostics& diagnostics = progress.diagnostics;
const auto nonnegative = [](const std::optional<double>& 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 <typename State>
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<MetricEvaluation> 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 <typename State, typename System, typename Model = NoModel>
class SolveResult {
public:
/** @brief Retain the physical definition before any numerical attempt is made. */
explicit SolveResult(std::shared_ptr<const Model> 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<utils::errors::SERiFErrorCode::INVALID_ARGUMENT>(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<const System>& 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<Model, NoModel>) { 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<State> checkpoint() const {
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_ARGUMENT>(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<const System> 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<const Model> m_model;
std::optional<State> m_state;
std::shared_ptr<const System> m_system;
std::optional<NewtonProgress> m_progress;
SolveStatus m_status{SolveStatus::initial_state_rejected};
std::string m_message;
SolveDiagnostics m_diagnostics;
};
}

188
src/include/serif/solver/nonlinear/steps.hpp Normal file
View File

@@ -0,0 +1,188 @@
#pragma once
#include <limits>
#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<NewtonProgress>& 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<SERiFErrorCode::INVALID_ARGUMENT>(checkpointable_options<std::remove_cvref_t<decltype(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<SERiFErrorCode::INVALID_ARGUMENT>(false, "Invalid Newton checkpoint: {}", *error);
}
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(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 <typename System>
struct InitialEvaluation {
std::shared_ptr<const System> system;
MetricEvaluation metric;
};
/** @brief Qualify the initial state and verify a restart against its original convergence reference. */
template <typename Problem>
[[nodiscard]] std::expected<InitialEvaluation<typename Problem::System>, NewtonFailure> evaluate_start(
Problem& problem, auto& state, auto& options, const std::optional<NewtonProgress>& 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<double>::epsilon() * std::max(1.0, restart->accepted_metric.residual_norm);
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(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<typename Problem::System>{.system = std::move(*evaluation), .metric = *metric};
}
/** @brief Solve and qualify a normalized correction before exposing its physical direction to backtracking. */
[[nodiscard]] std::expected<DirectionAssessment, NewtonFailure> compute_direction(
auto& problem, const auto& accepted, const auto& prepared, auto& options, SolveDiagnostics& diagnostics,
auto& rhs, auto& correction, std::span<double> physical_direction
) {
const int size = static_cast<int>(physical_direction.size());
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_DOF_COUNT>(
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 <typename State, typename System, typename Prepared>
struct QualifiedStep {
State state;
std::shared_ptr<const System> system;
MetricEvaluation metric;
std::optional<Prepared> 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 <typename Problem, typename State, typename System, typename Model>
[[nodiscard]] auto find_acceptable_step(
Problem& problem, const SolveResult<State, System, Model>& accepted, MetricEvaluation accepted_metric,
const auto& correction, std::span<const double> 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<State, System, Prepared>;
using Result = std::expected<Step, NewtonFailure>;
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<Prepared> 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<const double> correction, double next_step) {
if constexpr (checkpointable_options<std::remove_cvref_t<decltype(options)>>) {
options.strategies.contribute_diagnostics(diagnostics);
result.checkpoint(NewtonProgress{.accepted_metric = metric, .correction = {correction.begin(), correction.end()},
.next_step = next_step, .diagnostics = diagnostics});
}
}
}

67
src/lib/serif/solver/nonlinear/options.cpp Normal file
View File

@@ -0,0 +1,67 @@
#include "serif/solver/nonlinear/options.hpp"
#include <algorithm>
#include "serif/utils/error/errors.hpp"
namespace serif::solver {
void Newton::validate() const {
using namespace utils::errors;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(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<SERiFErrorCode::INVALID_ARGUMENT>(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<SERiFErrorCode::INVALID_ARGUMENT>(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<const double> 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;
}
}