module;

#include <concepts>
#include <type_traits>
#include <utility>

export module mean_field:solver.newton_numerics;

export import :solver.enthalpy_admissibility;
export import :solver.linear_backend;

export namespace mean_field::solver::nonlinear {
    struct SampledEnthalpy final { };
    struct PolynomialEnthalpy final {
        EnthalpyAdmissibilityOptions options{};
    };

    struct MonolithicGravity final { };

    // These tolerances control an auxiliary equation solve, not the physical
    // model's acceptance threshold. The complete residual remains authoritative.
    struct GravitySolveOptions final {
        LinearSolveControl control{.relativeTolerance = 1.0e-10, .absoluteTolerance = 1.0e-14,
                                   .maximumIterations = 300};
        linear::FGMRESOptions backend{.restartLength = 80};

        void Validate() const { control.Validate(); backend.Validate(); }
    };

    struct GravityCurvatureCorrection final {
        GravitySolveOptions gravity{};
    };

    struct EliminatedGravity final {
        GravitySolveOptions gravity{};
    };

    template <typename Candidate>
    concept EnthalpyNumerics = std::same_as<Candidate, SampledEnthalpy> ||
                              std::same_as<Candidate, PolynomialEnthalpy>;

    template <typename Candidate>
    concept GravityNumerics = std::same_as<Candidate, MonolithicGravity> ||
                              std::same_as<Candidate, GravityCurvatureCorrection> ||
                              std::same_as<Candidate, EliminatedGravity>;

    // Types select algorithms; tolerances and the mesh-dependent work sizes
    // remain runtime data. Adding a strategy cannot silently change a model's
    // equation set or impose an additional specification such as rho_c.
    template <EnthalpyNumerics Enthalpy = SampledEnthalpy, GravityNumerics Gravity = MonolithicGravity>
    struct NewtonNumerics final {
        using EnthalpyType = Enthalpy;
        using GravityType = Gravity;
        [[no_unique_address]] Enthalpy enthalpy{};
        [[no_unique_address]] Gravity gravity{};
    };

    template <typename Enthalpy, typename Gravity>
    NewtonNumerics(Enthalpy, Gravity) -> NewtonNumerics<Enthalpy, Gravity>;

    template <typename> struct IsNewtonNumerics : std::false_type { };
    template <typename Enthalpy, typename Gravity>
    struct IsNewtonNumerics<NewtonNumerics<Enthalpy, Gravity>> : std::true_type { };
    template <typename Candidate>
    concept NewtonNumericalConfiguration = IsNewtonNumerics<std::remove_cvref_t<Candidate>>::value;

    struct NumericalStrategyDiagnostics final {
        int polynomialChecks{0};
        int directionRepairs{0};
        double lastRepairFraction{0.0};
        double lastMeritDirectionalDerivative{0.0};
        int gravitySolves{0};
        int gravityIterations{0};
        int gravityFailures{0};
        double gravitySolveSeconds{0.0};
        double lastGravityTrueResidual{0.0};
        double lastGravityCorrectionNorm{0.0};
        double lastGravityPreparedDefect{0.0};
        double initialGravityResidual{0.0};
        double completedInitialGravityResidual{0.0};
        double lastLiftedRelativeTrueResidual{0.0};
    };
}