#include <print>
#include <utility>

#include "serif/eos/models/polytropic.hpp"
#include "serif/geometry/surface/surface.hpp"
#include "serif/utils/dependencies/computation_cache.hpp"

/**
 * @brief Prescribe one pressure and evaluate its residual over several already sampled boundary points.
 *
 * The EOS catalog and available state select the quantity enforcing the condition. This example binds density, as the
 * stellar model does; other models may bind another state. IsobaricSurface itself only knows the prescribed pressure.
 *
 * The ordinary dependency cache owns preparation. The EOS response is computed once for this pressure-only relation,
 * then reused across all boundary samples and subsequent changes to their density. No EOS evaluation occurs at graph construction.
 * Sampling a finite-element field on the boundary is deliberately a separate numerical operation.
 */
int main() {
    namespace surface = serif::geometry::surface;
    namespace dependencies = serif::utils::dependencies;
    namespace quantities = serif::dimensions;
    using Density = serif::discretization::blocks::density::mass;
    using EOS = serif::eos::models::Polytrope;
    using State = dependencies::Bindings<dependencies::Binding<quantities::MassDensity, Density>>;

    const EOS equation_of_state{3.0, 0.25};
    const surface::IsobaricSurface boundary{quantities::PressureValue{0.03125}};
    const auto context = dependencies::make_dependency_context(
        dependencies::borrow<surface::keys::EquationOfState>(equation_of_state, dependencies::track_by([](const EOS& model) {
            return std::pair{model.polytropic_index(), model.polytropic_constant()};
        })),
        surface::pressure_condition(boundary),
        surface::boundary_values<Density>({0.05, 0.1, 0.2})
    );

    auto graph = surface::make_pressure_surface_graph<EOS, State>();
    auto cache = dependencies::make_computation_cache<surface::SurfaceFailure>(std::move(graph), context);
    const auto preparation = cache.try_prepare(context);
    if (!preparation) {
        std::println(stderr, "Surface preparation failed: {}", preparation.error().error.message);
        return 1;
    }

    for (const auto& point : cache.get<surface::keys::Linearization>()) {
        std::println("Boundary density residual: {}; derivative: {}", point.residual, point.field_derivatives.get<Density>());
    }
    return 0;
}