MeanField/libmeanfield/interface/eos/polytropic.cppm

module;
#include <cmath>
#include <format>
#include <stdexcept>
export module mean_field:eos.polytrope;
export import :eos.evaluation;

export namespace mean_field::eos {
    class Polytrope final {
    public:
        struct Parameters final {
            double n;
            double K;
        };

        using Relations = RelationCatalog<
            PressureFromDensity,
            PressureFromSpecificEnthalpy,
            SpecificEnthalpyFromDensity,
            SpecificEnthalpyFromPressure,
            DensityFromSpecificEnthalpy>;

        explicit Polytrope(const Parameters parameters)
            : Polytrope(
                  parameters.n,
                  parameters.K
              ) {
        }

        Polytrope(
            const double polytropic_index,
            const double polytropic_constant
        )
            : m_polytropic_index(polytropic_index),
              m_polytropic_constant(polytropic_constant),
              m_enthalpy_scale((polytropic_index + 1.0) * polytropic_constant) {
            if (!std::isfinite(polytropic_index) || polytropic_index < 1.0) {
                throw std::invalid_argument(
                    std::format(
                        "The differentiable polytropic closure requires a "
                        "finite polytropic index greater than or equal to one. "
                        "Instead a value of {} has been provided",
                        polytropic_index
                    )
                );
            }

            if (!std::isfinite(polytropic_constant) || polytropic_constant <= 0.0) {
                throw std::invalid_argument(
                    std::format(
                        "The polytropic constant must be finite and positive. "
                        "Instead a value of {} has been provided",
                        polytropic_constant
                    )
                );
            }
        };

        [[nodiscard]] double polytropic_index() const noexcept {
            return m_polytropic_index;
        }

        [[nodiscard]] double polytropic_constant() const noexcept {
            return m_polytropic_constant;
        }

        [[nodiscard]] double enthalpy_scale() const noexcept {
            return m_enthalpy_scale;
        }

        [[nodiscard]] dimensions::PressureValue evaluate(
            PressureFromDensity,
            const dimensions::DensityValue density
        ) const {
            validate_nonnegativity(density.value(), "density");
            if (density.value() == 0.0) {
                return dimensions::PressureValue{0.0};
            }

            return dimensions::PressureValue{
                m_polytropic_constant * std::pow(density.value(), 1.0 + 1.0 / m_polytropic_index)
            };
        }

        [[nodiscard]] dimensions::SpecificEnthalpyValue evaluate(
            SpecificEnthalpyFromDensity,
            const dimensions::DensityValue density
        ) const {
            validate_nonnegativity(density.value(), "density");
            if (density.value() == 0.0) {
                return dimensions::SpecificEnthalpyValue{0.0};
            }

            return dimensions::SpecificEnthalpyValue{
                m_enthalpy_scale * std::pow(density.value(), 1.0 / m_polytropic_index)
            };
        }

        [[nodiscard]] dimensions::DensityValue evaluate(
            DensityFromSpecificEnthalpy,
            const dimensions::SpecificEnthalpyValue specificEnthalpy
        ) const {
            validate_finite(specificEnthalpy.value(), "specific enthalpy");

            if (specificEnthalpy.value() <= 0.0) {
                return dimensions::DensityValue{0.0};
            }

            return dimensions::DensityValue{std::pow(specificEnthalpy.value() / m_enthalpy_scale, m_polytropic_index)};
        }

        [[nodiscard]] dimensions::PressureValue evaluate(
            PressureFromSpecificEnthalpy,
            const dimensions::SpecificEnthalpyValue specificEnthalpy
        ) const {
            const dimensions::DensityValue density = evaluate(DensityFromSpecificEnthalpy{}, specificEnthalpy);

            if (specificEnthalpy.value() <= 0.0) {
                return dimensions::PressureValue{0.0};
            }

            return dimensions::PressureValue{density.value() * specificEnthalpy.value() / (m_polytropic_index + 1.0)};
        }

        [[nodiscard]] dimensions::SpecificEnthalpyValue evaluate(
            SpecificEnthalpyFromPressure,
            const dimensions::PressureValue pressure
        ) const {
            validate_nonnegativity(pressure.value(), "pressure");
            if (pressure.value() == 0.0) {
                return dimensions::SpecificEnthalpyValue{0.0};
            }

            const double indexPlusOne = m_polytropic_index + 1.0;

            return dimensions::SpecificEnthalpyValue{
                indexPlusOne * std::pow(m_polytropic_constant, m_polytropic_index / indexPlusOne) *
                std::pow(pressure.value(), 1.0 / indexPlusOne)
            };
        }

        [[nodiscard]] PartialDerivative<
            dimensions::quantity::Density,
            dimensions::quantity::SpecificEnthalpy>
        partialDerivative(
            DensityFromSpecificEnthalpy,
            WithRespectTo<dimensions::quantity::SpecificEnthalpy>,
            const dimensions::SpecificEnthalpyValue specificEnthalpy
        ) const {
            validate_finite(specificEnthalpy.value(), "specific enthalpy");
            if (specificEnthalpy.value() < 0.0) {
                return PartialDerivative<dimensions::quantity::Density, dimensions::quantity::SpecificEnthalpy>{0.0};
            }

            if (specificEnthalpy.value() == 0.0) {
                return PartialDerivative<dimensions::quantity::Density, dimensions::quantity::SpecificEnthalpy>{
                    m_polytropic_index == 1.0 ? 1.0 / m_enthalpy_scale : 0.0
                };
            }

            return PartialDerivative<dimensions::quantity::Density, dimensions::quantity::SpecificEnthalpy>{
                m_polytropic_index / m_enthalpy_scale *
                std::pow(specificEnthalpy.value() / m_enthalpy_scale, m_polytropic_index - 1.0)
            };
        }

        [[nodiscard]] PartialDerivative<
            dimensions::quantity::Pressure,
            dimensions::quantity::SpecificEnthalpy>
        partialDerivative(
            PressureFromSpecificEnthalpy,
            WithRespectTo<dimensions::quantity::SpecificEnthalpy>,
            const dimensions::SpecificEnthalpyValue specificEnthalpy
        ) const {
            const dimensions::DensityValue density = evaluate(DensityFromSpecificEnthalpy{}, specificEnthalpy);

            return PartialDerivative<dimensions::quantity::Pressure, dimensions::quantity::SpecificEnthalpy>{
                density.value()
            };
        }

        [[nodiscard]] PartialDerivative<
            dimensions::quantity::Pressure,
            dimensions::quantity::Density>
        partialDerivative(
            PressureFromDensity,
            WithRespectTo<dimensions::quantity::Density>,
            const dimensions::DensityValue density
        ) const {
            validate_nonnegativity(density.value(), "density");
            if (density.value() == 0.0) {
                return PartialDerivative<dimensions::quantity::Pressure, dimensions::quantity::Density>{0.0};
            }

            return PartialDerivative<dimensions::quantity::Pressure, dimensions::quantity::Density>{
                m_polytropic_constant * (1.0 + 1.0 / m_polytropic_index) *
                std::pow(density.value(), 1.0 / m_polytropic_index)
            };
        }

    private:
        static void validate_finite(
            const double value,
            const char *quantity
        ) {
            if (!std::isfinite(value)) {
                throw EvaluationError(
                    EvaluationErrorCode::nonfinite_input, std::format(
                                                              "The {} must be finite. Instead a value of {} has been "
                                                              "provided",
                                                              quantity, value
                                                          )
                );
            }
        }

        static void validate_nonnegativity(
            const double value,
            const char *quantity
        ) {
            validate_finite(value, quantity);
            if (value < 0.0) {
                throw EvaluationError(
                    EvaluationErrorCode::outside_domain, std::format(
                                                             "The {} must be non-negative. Instead a value of {} "
                                                             "has been "
                                                             "provided",
                                                             quantity, value
                                                         )
                );
            }
        }
        double m_polytropic_index;
        double m_polytropic_constant;
        double m_enthalpy_scale;
    };
} // namespace mean_field::eos