MeanField/tests/physics/equation_of_state_runtime_view.cpp

#include <array>
#include <concepts>
#include <expected>
#include <limits>
#include <memory>
#include <span>
#include <string_view>
#include <type_traits>

#include <catch2/catch_test_macros.hpp>

import mean_field;
import test_helpers;

namespace {
    namespace eos = mean_field::eos;

    class LinearPressureEquationOfState final {
    public:
        using Relations = eos::RelationCatalog<eos::PressureFromDensity>;

        [[nodiscard]] constexpr eos::PressureValue evaluate(
            eos::PressureFromDensity,
            const eos::DensityValue density
        ) const noexcept {
            return eos::PressureValue{2.0 * density.value() + 0.5};
        }
    };

    struct DensityAlias final : eos::ThermodynamicQuantity {
        static constexpr std::string_view identifier = "density";
    };

    class AmbiguouslyIdentifiedEquationOfState final {
    public:
        using Relations = eos::RelationCatalog<eos::Relation<DensityAlias, eos::quantity::Density>>;

        [[nodiscard]] constexpr eos::QuantityValue<DensityAlias> evaluate(
            eos::Relation<
                DensityAlias,
                eos::quantity::Density>,
            const eos::DensityValue density
        ) const noexcept {
            return eos::QuantityValue<DensityAlias>{density.value()};
        }
    };

    [[nodiscard]] std::expected<
        eos::PressureValue,
        eos::EvaluationError>
    pressureAtDensity(
        const eos::EquationOfStateView equationOfState,
        const eos::DensityValue density
    ) {
        return equationOfState.tryEvaluate<eos::quantity::Pressure>(density);
    }

    [[nodiscard]] const eos::RuntimeRelationDescriptor *findRelation(
        const eos::EquationOfStateView equationOfState,
        const eos::ThermodynamicQuantityId output,
        const eos::ThermodynamicQuantityId input
    ) {
        for (const eos::RuntimeRelationDescriptor &relation : equationOfState.relations()) {
            if (relation.outputQuantity == output && relation.inputQuantities.size() == 1 &&
                relation.inputQuantities[0] == input) {
                return std::addressof(relation);
            }
        }

        return nullptr;
    }
} // namespace

TEST_CASE(
    "Runtime EOS View Generates The Polytropic Relation Catalog",
    tags::equation_of_state_runtime_contract
) {
    STATIC_CHECK(eos::RuntimeEquationOfStateModel<eos::Polytrope>);
    STATIC_CHECK(eos::RuntimeEquationOfStateModel<LinearPressureEquationOfState>);
    STATIC_CHECK(eos::EquationOfStateModel<AmbiguouslyIdentifiedEquationOfState>);
    STATIC_CHECK_FALSE(eos::RuntimeEquationOfStateModel<AmbiguouslyIdentifiedEquationOfState>);
    STATIC_CHECK(std::is_trivially_copyable_v<eos::EquationOfStateView>);
    STATIC_CHECK_FALSE(std::constructible_from<eos::EquationOfStateView, eos::Polytrope &&>);

    const eos::Polytrope equationOfState(3.0, 0.25);
    const eos::Polytrope secondEquationOfState(1.5, 0.73);

    const eos::EquationOfStateView view{equationOfState};
    const eos::EquationOfStateView secondView{secondEquationOfState};

    REQUIRE(view.relations().size() == eos::Polytrope::Relations::size);
    CHECK(view.relations().data() == secondView.relations().data());

    CHECK(eos::thermodynamicQuantityId<eos::quantity::Density>.name() == "density");
    CHECK(eos::thermodynamicQuantityId<eos::quantity::Pressure>.name() == "pressure");
    CHECK(eos::thermodynamicQuantityId<eos::quantity::SpecificEnthalpy>.name() == "specific_enthalpy");

    const eos::RuntimeRelationDescriptor *pressureFromDensity = findRelation(
        view, eos::thermodynamicQuantityId<eos::quantity::Pressure>,
        eos::thermodynamicQuantityId<eos::quantity::Density>
    );

    REQUIRE(pressureFromDensity != nullptr);
    CHECK(pressureFromDensity->hasPartialDerivative(0));

    const eos::RuntimeRelationDescriptor *specificEnthalpyFromPressure = findRelation(
        view, eos::thermodynamicQuantityId<eos::quantity::SpecificEnthalpy>,
        eos::thermodynamicQuantityId<eos::quantity::Pressure>
    );

    REQUIRE(specificEnthalpyFromPressure != nullptr);
    CHECK_FALSE(specificEnthalpyFromPressure->hasPartialDerivative(0));

    const eos::RuntimeRelationDescriptor *pressureFromSpecificEnthalpy = findRelation(
        view, eos::thermodynamicQuantityId<eos::quantity::Pressure>,
        eos::thermodynamicQuantityId<eos::quantity::SpecificEnthalpy>
    );

    REQUIRE(pressureFromSpecificEnthalpy != nullptr);
    CHECK(pressureFromSpecificEnthalpy->hasPartialDerivative(0));

    const eos::RuntimeRelationDescriptor *specificEnthalpyFromDensity = findRelation(
        view, eos::thermodynamicQuantityId<eos::quantity::SpecificEnthalpy>,
        eos::thermodynamicQuantityId<eos::quantity::Density>
    );

    REQUIRE(specificEnthalpyFromDensity != nullptr);
    CHECK_FALSE(specificEnthalpyFromDensity->hasPartialDerivative(0));

    const eos::RuntimeRelationDescriptor *densityFromSpecificEnthalpy = findRelation(
        view, eos::thermodynamicQuantityId<eos::quantity::Density>,
        eos::thermodynamicQuantityId<eos::quantity::SpecificEnthalpy>
    );

    REQUIRE(densityFromSpecificEnthalpy != nullptr);
    CHECK(densityFromSpecificEnthalpy->hasPartialDerivative(0));
}

TEST_CASE(
    "Runtime EOS View Matches Typed Polytropic Evaluation",
    tags::equation_of_state_runtime_compatibility
) {
    const eos::Polytrope equationOfState(3.0, 0.25);
    const eos::EquationOfStateView view{equationOfState};

    const eos::DensityValue density{0.7};
    const eos::SpecificEnthalpyValue specificEnthalpy{0.9};
    const eos::PressureValue pressure{0.04};

    const auto runtimePressureFromDensity          = view.tryEvaluate<eos::quantity::Pressure>(density);
    const auto runtimePressureFromSpecificEnthalpy = view.tryEvaluate<eos::quantity::Pressure>(specificEnthalpy);
    const auto runtimeSpecificEnthalpyFromDensity  = view.tryEvaluate<eos::quantity::SpecificEnthalpy>(density);
    const auto runtimeSpecificEnthalpyFromPressure = view.tryEvaluate<eos::quantity::SpecificEnthalpy>(pressure);
    const auto runtimeDensityFromSpecificEnthalpy  = view.tryEvaluate<eos::quantity::Density>(specificEnthalpy);

    REQUIRE(runtimePressureFromDensity.has_value());
    REQUIRE(runtimePressureFromSpecificEnthalpy.has_value());
    REQUIRE(runtimeSpecificEnthalpyFromDensity.has_value());
    REQUIRE(runtimeSpecificEnthalpyFromPressure.has_value());
    REQUIRE(runtimeDensityFromSpecificEnthalpy.has_value());

    CHECK(
        runtimePressureFromDensity->value() == eos::evaluate<eos::quantity::Pressure>(equationOfState, density).value()
    );
    CHECK(
        runtimePressureFromSpecificEnthalpy->value() ==
        eos::evaluate<eos::quantity::Pressure>(equationOfState, specificEnthalpy).value()
    );
    CHECK(
        runtimeSpecificEnthalpyFromDensity->value() ==
        eos::evaluate<eos::quantity::SpecificEnthalpy>(equationOfState, density).value()
    );
    CHECK(
        runtimeSpecificEnthalpyFromPressure->value() ==
        eos::evaluate<eos::quantity::SpecificEnthalpy>(equationOfState, pressure).value()
    );
    CHECK(
        runtimeDensityFromSpecificEnthalpy->value() ==
        eos::evaluate<eos::quantity::Density>(equationOfState, specificEnthalpy).value()
    );

    const std::array runtimeDensityInput{
        eos::RuntimeQuantityValue{eos::thermodynamicQuantityId<eos::quantity::Density>, density.value()}
    };

    const auto erasedPressureFromDensity = view.tryEvaluate(
        eos::thermodynamicQuantityId<eos::quantity::Pressure>,
        std::span<const eos::RuntimeQuantityValue>{runtimeDensityInput}
    );

    REQUIRE(erasedPressureFromDensity.has_value());
    CHECK(erasedPressureFromDensity->quantity == eos::thermodynamicQuantityId<eos::quantity::Pressure>);
    CHECK(erasedPressureFromDensity->value == runtimePressureFromDensity->value());

    const auto runtimePressureDerivative =
        view.tryPartialDerivative<eos::quantity::Pressure, eos::quantity::SpecificEnthalpy>(specificEnthalpy);

    const auto runtimeDensityDerivative =
        view.tryPartialDerivative<eos::quantity::Density, eos::quantity::SpecificEnthalpy>(specificEnthalpy);

    const auto runtimePressureDensityDerivative =
        view.tryPartialDerivative<eos::quantity::Pressure, eos::quantity::Density>(density);

    REQUIRE(runtimePressureDerivative.has_value());
    REQUIRE(runtimeDensityDerivative.has_value());
    REQUIRE(runtimePressureDensityDerivative.has_value());

    CHECK(
        runtimePressureDerivative->value() ==
        eos::partialDerivative<eos::quantity::Pressure, eos::quantity::SpecificEnthalpy>(
            equationOfState, specificEnthalpy
        )
            .value()
    );
    CHECK(
        runtimeDensityDerivative->value() ==
        eos::partialDerivative<eos::quantity::Density, eos::quantity::SpecificEnthalpy>(
            equationOfState, specificEnthalpy
        )
            .value()
    );
    CHECK(
        runtimePressureDensityDerivative->value() ==
        eos::partialDerivative<eos::quantity::Pressure, eos::quantity::Density>(equationOfState, density).value()
    );

    const auto erasedPressureDensityDerivative = view.tryPartialDerivative(
        eos::thermodynamicQuantityId<eos::quantity::Pressure>, eos::thermodynamicQuantityId<eos::quantity::Density>,
        std::span<const eos::RuntimeQuantityValue>{runtimeDensityInput}
    );

    REQUIRE(erasedPressureDensityDerivative.has_value());
    CHECK(*erasedPressureDensityDerivative == runtimePressureDensityDerivative->value());
}

TEST_CASE(
    "Runtime EOS View Reports Unsupported And Invalid Requests",
    tags::equation_of_state_runtime_contract
) {
    const eos::Polytrope equationOfState(3.0, 0.25);
    const eos::EquationOfStateView view{equationOfState};

    constexpr eos::ThermodynamicQuantityId temperature{"temperature"};

    const std::array densityInput{eos::RuntimeQuantityValue{eos::thermodynamicQuantityId<eos::quantity::Density>, 0.7}};

    const std::array pressureInput{
        eos::RuntimeQuantityValue{eos::thermodynamicQuantityId<eos::quantity::Pressure>, 0.04}
    };

    const std::array<eos::RuntimeQuantityValue, 0> noInputs{};

    const auto unsupportedOutput =
        view.tryEvaluate(temperature, std::span<const eos::RuntimeQuantityValue>{densityInput});
    REQUIRE_FALSE(unsupportedOutput.has_value());
    CHECK(unsupportedOutput.error().code() == eos::EvaluationErrorCode::unsupported_relation);

    const auto wrongInputCount = view.tryEvaluate(
        eos::thermodynamicQuantityId<eos::quantity::Pressure>, std::span<const eos::RuntimeQuantityValue>{noInputs}
    );
    REQUIRE_FALSE(wrongInputCount.has_value());
    CHECK(wrongInputCount.error().code() == eos::EvaluationErrorCode::wrong_input_count);

    const auto wrongInputQuantity = view.tryEvaluate(
        eos::thermodynamicQuantityId<eos::quantity::Density>, std::span<const eos::RuntimeQuantityValue>{pressureInput}
    );
    REQUIRE_FALSE(wrongInputQuantity.has_value());
    CHECK(wrongInputQuantity.error().code() == eos::EvaluationErrorCode::wrong_input_quantity);

    const auto unsupportedDerivative = view.tryPartialDerivative(
        eos::thermodynamicQuantityId<eos::quantity::SpecificEnthalpy>,
        eos::thermodynamicQuantityId<eos::quantity::Pressure>, std::span<const eos::RuntimeQuantityValue>{pressureInput}
    );
    REQUIRE_FALSE(unsupportedDerivative.has_value());
    CHECK(unsupportedDerivative.error().code() == eos::EvaluationErrorCode::unsupported_derivative);

    const auto invalidDensity = view.tryEvaluate<eos::quantity::Pressure>(eos::DensityValue{-0.1});
    REQUIRE_FALSE(invalidDensity.has_value());
    CHECK(invalidDensity.error().code() == eos::EvaluationErrorCode::outside_domain);

    const auto nonfiniteDensity =
        view.tryEvaluate<eos::quantity::Pressure>(eos::DensityValue{std::numeric_limits<double>::quiet_NaN()});
    REQUIRE_FALSE(nonfiniteDensity.has_value());
    CHECK(nonfiniteDensity.error().code() == eos::EvaluationErrorCode::nonfinite_input);
}

TEST_CASE(
    "One Runtime EOS Function Accepts Heterogeneous Concrete Models",
    tags::equation_of_state_runtime_compatibility
) {
    const eos::Polytrope polytrope(3.0, 0.25);
    const LinearPressureEquationOfState linearEquationOfState;

    const std::array views{eos::EquationOfStateView{polytrope}, eos::EquationOfStateView{linearEquationOfState}};

    const eos::DensityValue density{0.7};

    const auto polytropicPressure = pressureAtDensity(views[0], density);
    const auto linearPressure     = pressureAtDensity(views[1], density);

    REQUIRE(polytropicPressure.has_value());
    REQUIRE(linearPressure.has_value());

    CHECK(polytropicPressure->value() == eos::evaluate<eos::quantity::Pressure>(polytrope, density).value());
    CHECK(linearPressure->value() == 1.9);
}

TEST_CASE(
    "Runtime EOS View Remains Valid When Stable Ownership Moves",
    tags::equation_of_state_runtime_contract
) {
    auto owner = std::make_unique<const eos::Polytrope>(3.0, 0.25);
    const eos::EquationOfStateView view{*owner};

    auto movedOwner     = std::move(owner);

    const auto pressure = view.tryEvaluate<eos::quantity::Pressure>(eos::DensityValue{0.7});

    REQUIRE(movedOwner != nullptr);
    REQUIRE(pressure.has_value());
    CHECK(pressure->value() == eos::evaluate<eos::quantity::Pressure>(*movedOwner, eos::DensityValue{0.7}).value());
}