/**
 * \file
 * \copyright
 * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 * Created on February 17, 2021, 4:22 PM
 */

#include <gtest/gtest.h>

#include <Eigen/Core>
#include <cmath>
#include <functional>
#include <limits>
#include <memory>
#include <vector>

#include "BaseLib/ConfigTree.h"
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Properties/ThermalConductivity/CreateSaturationWeightedThermalConductivity.h"
#include "MaterialLib/MPL/Properties/ThermalConductivity/SaturationWeightedThermalConductivity.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
#include "ParameterLib/ConstantParameter.h"
#include "Tests/TestTools.h"

namespace MPL = MaterialPropertyLib;

std::unique_ptr<MaterialPropertyLib::Property>
createTestSaturationWeightedThermalConductivityProperty(
    const char xml[], int const geometry_dimension,
    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
    std::function<std::unique_ptr<MaterialPropertyLib::Property>(
        int const geometry_dimension,
        BaseLib::ConfigTree const& config,
        std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
            parameters)>
        createProperty)
{
    auto ptree = Tests::readXml(xml);
    BaseLib::ConfigTree conf(std::move(ptree), "", BaseLib::ConfigTree::onerror,
                             BaseLib::ConfigTree::onwarning);
    auto const& sub_config = conf.getConfigSubtree("property");
    // Parsing the property name:
    auto const property_name =
        sub_config.getConfigParameter<std::string>("name");

    return createProperty(geometry_dimension, sub_config, parameters);
}

TEST(MaterialPropertyLib,
     SaturationWeightedThermalConductivity1Darithmetic_squareroot)
{
    const char xml[] =
        "<property>"
        "   <name>thermal_conductivity</name>"
        "   <type>SaturationWeightedThermalConductivity</type>"
        "   <mean_type>arithmetic_squareroot</mean_type>"
        "   <dry_thermal_conductivity>k_T_dry</dry_thermal_conductivity>"
        "   <wet_thermal_conductivity>k_T_wet</wet_thermal_conductivity>"
        "</property>";

    double const k_T_dry = 0.1;
    double const k_T_wet = 0.3;
    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> parameters;
    parameters.push_back(
        std::make_unique<ParameterLib::ConstantParameter<double>>("k_T_dry",
                                                                  k_T_dry));
    parameters.push_back(
        std::make_unique<ParameterLib::ConstantParameter<double>>("k_T_wet",
                                                                  k_T_wet));

    const int dimemsion = 1;
    std::unique_ptr<MPL::Property> const k_T_ptr =
        createTestSaturationWeightedThermalConductivityProperty(
            xml, dimemsion, parameters,
            MPL::createSaturationWeightedThermalConductivity);

    MPL::Property const& k_T_property = *k_T_ptr;

    MPL::VariableArray variable_array;
    ParameterLib::SpatialPosition const pos;
    double const t = std::numeric_limits<double>::quiet_NaN();
    double const dt = std::numeric_limits<double>::quiet_NaN();

    // For S <= 0.0
    {
        double const S_less_than_zero[] = {-1, 0.0};
        for (int i = 0; i < 2; i++)
        {
            variable_array.liquid_saturation = S_less_than_zero[i];

            double const computed_k_T =
                k_T_property.template value<double>(variable_array, pos, t, dt);
            ASSERT_LE(std::fabs(k_T_dry - computed_k_T), 1e-10)
                << "for expected thermal conductivity " << k_T_dry
                << " and for computed thermal conductivity." << computed_k_T;

            double const computed_dk_T_dS =
                k_T_property.template dValue<double>(
                    variable_array,
                    MaterialPropertyLib::Variable::liquid_saturation, pos, t,
                    dt);

            ASSERT_LE(std::fabs(0.0 - computed_dk_T_dS), 1e-10)
                << "for expected derivative of thermal conductivity with "
                   "respect to saturation "
                << 0.0
                << " and for computed derivative of thermal conductivity with "
                   "respect to saturation."
                << computed_dk_T_dS;
        }
    }
    // For S > 1.0
    {
        variable_array.liquid_saturation = 2.0;

        double const computed_k_T =
            k_T_property.template value<double>(variable_array, pos, t, dt);
        ASSERT_LE(std::fabs(k_T_wet - computed_k_T), 1e-10)
            << "for expected thermal conductivity " << k_T_wet
            << " and for computed thermal conductivity." << computed_k_T;

        double const computed_dk_T_dS = k_T_property.template dValue<double>(
            variable_array, MaterialPropertyLib::Variable::liquid_saturation,
            pos, t, dt);

        ASSERT_LE(std::fabs(0.0 - computed_dk_T_dS), 1e-10)
            << "for expected derivative of thermal conductivity with respect "
               "to saturation "
            << 0.0
            << " and for computed derivative of thermal conductivity with "
               "respect to saturation."
            << computed_dk_T_dS;
    }

    // For S in (0, 1)
    std::vector<double> const S_L = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
    std::vector<double> const expected_k_T = {
        1.632456e-01, 1.894427e-01, 2.095445e-01, 2.264911e-01,
        2.414214e-01, 2.549193e-01, 2.673320e-01, 2.788854e-01};

    double const perturbation = 1.e-6;
    for (std::size_t i = 0; i < S_L.size(); i++)
    {
        variable_array.liquid_saturation = S_L[i];

        double const k_T_i =
            k_T_property.template value<double>(variable_array, pos, t, dt);
        ASSERT_LE(std::fabs(expected_k_T[i] - k_T_i), 1e-7)
            << "for expected thermal conductivity " << expected_k_T[i]
            << " and for computed thermal conductivity." << k_T_i;

        variable_array.liquid_saturation = S_L[i] - perturbation;
        double const k_T_i_0 =
            k_T_property.template value<double>(variable_array, pos, t, dt);

        variable_array.liquid_saturation = S_L[i] + perturbation;

        double const k_T_i_1 =
            k_T_property.template value<double>(variable_array, pos, t, dt);

        double const approximated_dk_T_dS =
            0.5 * (k_T_i_1 - k_T_i_0) / perturbation;
        double const analytic_dk_T_dS = k_T_property.template dValue<double>(
            variable_array, MaterialPropertyLib::Variable::liquid_saturation,
            pos, t, dt);

        ASSERT_LE(std::fabs(analytic_dk_T_dS - approximated_dk_T_dS), 1e-5)
            << "for expected derivative of thermal conductivity with respect "
               "to saturation "
            << analytic_dk_T_dS
            << " and for computed derivative of thermal conductivity with "
               "respect to saturation."
            << approximated_dk_T_dS;
    }
}

TEST(MaterialPropertyLib,
     SaturationWeightedThermalConductivity3Darithmetic_squareroot)
{
    const char xml[] =
        "<property>"
        "   <name>thermal_conductivity</name>"
        "   <type>SaturationWeightedThermalConductivity</type>"
        "   <mean_type>arithmetic_squareroot</mean_type>"
        "   <dry_thermal_conductivity>k_T_dry</dry_thermal_conductivity>"
        "   <wet_thermal_conductivity>k_T_wet</wet_thermal_conductivity>"
        "</property>";

    std::vector<double> k_T_dry{0.1, 0.1, 0.1};
    std::vector<double> k_T_wet{0.3, 0.3, 0.3};
    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> parameters;
    parameters.push_back(
        std::make_unique<ParameterLib::ConstantParameter<double>>("k_T_dry",
                                                                  k_T_dry));
    parameters.push_back(
        std::make_unique<ParameterLib::ConstantParameter<double>>("k_T_wet",
                                                                  k_T_wet));

    const int dimemsion = 3;
    std::unique_ptr<MPL::Property> const k_T_ptr =
        createTestSaturationWeightedThermalConductivityProperty(
            xml, dimemsion, parameters,
            MPL::createSaturationWeightedThermalConductivity);

    MPL::Property const& k_T_property = *k_T_ptr;

    MPL::VariableArray variable_array;
    ParameterLib::SpatialPosition const pos;
    double const t = std::numeric_limits<double>::quiet_NaN();
    double const dt = std::numeric_limits<double>::quiet_NaN();

    // For S <= 0.0
    {
        double const S_less_than_zero[] = {-1, 0.0};
        for (int i = 0; i < 2; i++)
        {
            variable_array.liquid_saturation = S_less_than_zero[i];

            auto const computed_k_T = MPL::formEigenTensor<3>(
                k_T_property.value(variable_array, pos, t, dt));

            for (std::size_t k = 0; k < k_T_dry.size(); k++)
            {
                ASSERT_LE(std::fabs(k_T_dry[k] - computed_k_T(k, k)), 1e-10)
                    << "for expected thermal conductivity " << k_T_dry[k]
                    << " and for computed thermal conductivity."
                    << computed_k_T(k, k);
            }

            auto const computed_dk_T_dS =
                MPL::formEigenTensor<3>(k_T_property.dValue(
                    variable_array,
                    MaterialPropertyLib::Variable::liquid_saturation, pos, t,
                    dt));

            for (std::size_t k = 0; k < k_T_dry.size(); k++)
            {
                ASSERT_LE(std::fabs(0.0 - computed_dk_T_dS(k, k)), 1e-10)
                    << "for expected derivative of thermal conductivity with "
                       "respect to saturation "
                    << 0.0
                    << " and for computed derivative of thermal conductivity "
                       "with respect to saturation."
                    << computed_dk_T_dS(k, k);
            }
        }
    }
    // For S > 1.0
    {
        variable_array.liquid_saturation = 2.0;

        auto const computed_k_T = MPL::formEigenTensor<3>(
            k_T_property.value(variable_array, pos, t, dt));

        for (std::size_t k = 0; k < k_T_dry.size(); k++)
        {
            ASSERT_LE(std::fabs(k_T_wet[k] - computed_k_T(k, k)), 1e-10)
                << "for expected thermal conductivity " << k_T_wet[k]
                << " and for computed thermal conductivity."
                << computed_k_T(k, k);
        }

        auto const computed_dk_T_dS =
            MPL::formEigenTensor<3>(k_T_property.dValue(
                variable_array,
                MaterialPropertyLib::Variable::liquid_saturation, pos, t, dt));

        for (std::size_t k = 0; k < k_T_dry.size(); k++)
        {
            ASSERT_LE(std::fabs(0.0 - computed_dk_T_dS(k, k)), 1e-10)
                << "for expected derivative of thermal conductivity with "
                   "respect to saturation "
                << 0.0
                << " and for computed derivative of thermal conductivity with "
                   "respect to saturation."
                << computed_dk_T_dS(k, k);
        }
    }

    // For S in (0, 1)
    std::vector<double> const S_L = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
    std::vector<double> const expected_k_T = {
        1.632456e-01, 1.894427e-01, 2.095445e-01, 2.264911e-01,
        2.414214e-01, 2.549193e-01, 2.673320e-01, 2.788854e-01};

    double const perturbation = 1.e-6;
    for (std::size_t i = 0; i < S_L.size(); i++)
    {
        variable_array.liquid_saturation = S_L[i];

        auto const k_T_i = MPL::formEigenTensor<3>(
            k_T_property.value(variable_array, pos, t, dt));

        for (std::size_t k = 0; k < k_T_dry.size(); k++)
        {
            ASSERT_LE(std::fabs(expected_k_T[i] - k_T_i(k, k)), 1e-7)
                << "for expected thermal conductivity " << expected_k_T[i]
                << " and for computed thermal conductivity." << k_T_i(k, k);
        }

        variable_array.liquid_saturation = S_L[i] - perturbation;
        auto const k_T_i_0 = MPL::formEigenTensor<3>(
            k_T_property.value(variable_array, pos, t, dt));

        variable_array.liquid_saturation = S_L[i] + perturbation;

        auto const k_T_i_1 = MPL::formEigenTensor<3>(
            k_T_property.value(variable_array, pos, t, dt));

        auto const analytic_dk_T_dS =
            MPL::formEigenTensor<3>(k_T_property.dValue(
                variable_array,
                MaterialPropertyLib::Variable::liquid_saturation, pos, t, dt));

        for (std::size_t k = 0; k < k_T_dry.size(); k++)
        {
            auto const approximated_dk_T_dS =
                0.5 * (k_T_i_1(k, k) - k_T_i_0(k, k)) / perturbation;
            ASSERT_LE(std::fabs(analytic_dk_T_dS(k, k) - approximated_dk_T_dS),
                      1e-5)
                << "for expected derivative of thermal conductivity with "
                   "respect to saturation "
                << analytic_dk_T_dS(k, k)
                << " and for computed derivative of thermal conductivity with "
                   "respect to saturation."
                << approximated_dk_T_dS;
        }
    }
}
TEST(MaterialPropertyLib,
     SaturationWeightedThermalConductivity1Darithmetic_linear)
{
    ParameterLib::ConstantParameter<double> const k_dry{"k_dry", {0.2}};
    ParameterLib::ConstantParameter<double> const k_wet{"k_wet", {1.5}};

    auto const k_model_eff = MPL::SaturationWeightedThermalConductivity<
        MaterialPropertyLib::MeanType::ARITHMETIC_LINEAR, 1>(
        "thermal_conductivity", k_dry, k_wet);

    ParameterLib::SpatialPosition const pos;
    double const t = std::numeric_limits<double>::quiet_NaN();
    double const dt = std::numeric_limits<double>::quiet_NaN();
    MPL::VariableArray vars;

    /// check derivative
    {
        double const S_0 = -0.1;
        double const S_max = 1.1;
        int const n_steps = 10000;
        double const eps = 1e-6;
        for (int i = 0; i <= n_steps; i++)
        {
            double const S_L = S_0 + i * (S_max - S_0) / n_steps;
            vars.liquid_saturation = S_L;

            auto const lambda =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            switch (i)
            {
                case 0:
                    ASSERT_LE(std::abs(lambda - 0.2), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
                case 5000:
                    ASSERT_LE(std::abs(lambda - 0.85), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
                case 10000:
                    ASSERT_LE(std::abs(lambda - 1.5), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
            }

            auto const dlambda = std::get<double>(k_model_eff.dValue(
                vars, MaterialPropertyLib::Variable::liquid_saturation, pos, t,
                dt));

            vars.liquid_saturation = S_L - eps;
            auto const lambda_minus =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            vars.liquid_saturation = S_L + eps;
            auto const lambda_plus =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            double const Dlambda = (lambda_plus - lambda_minus) / 2 / eps;

            ASSERT_LE(std::abs(dlambda - Dlambda), 1e-9)
                << "for conductivity " << lambda << " and saturation " << S_L;
        }
    }
}
TEST(MaterialPropertyLib,
     SaturationWeightedThermalConductivity3Darithmetic_linear)
{
    ParameterLib::ConstantParameter<double> const k_dry{"k_dry", {0.2}};
    ParameterLib::ConstantParameter<double> const k_wet{"k_wet", {1.5}};

    auto const k_model_eff = MPL::SaturationWeightedThermalConductivity<
        MaterialPropertyLib::MeanType::ARITHMETIC_LINEAR, 3>(
        "thermal_conductivity", k_dry, k_wet);

    ParameterLib::SpatialPosition const pos;
    double const t = std::numeric_limits<double>::quiet_NaN();
    double const dt = std::numeric_limits<double>::quiet_NaN();
    MPL::VariableArray vars;

    {
        double const S_0 = -0.1;
        double const S_max = 1.1;
        int const n_steps = 10000;
        double const eps = 1e-6;
        for (int i = 0; i <= n_steps; i++)
        {
            double const S_L = S_0 + i * (S_max - S_0) / n_steps;
            vars.liquid_saturation = S_L;

            auto const lambda =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            switch (i)
            {
                case 0:
                    ASSERT_LE(std::abs(lambda - 0.2), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
                case 5000:
                    ASSERT_LE(std::abs(lambda - 0.85), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
                case 10000:
                    ASSERT_LE(std::abs(lambda - 1.5), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
            }

            auto const dlambda = std::get<double>(k_model_eff.dValue(
                vars, MaterialPropertyLib::Variable::liquid_saturation, pos, t,
                dt));

            vars.liquid_saturation = S_L - eps;
            auto const lambda_minus =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            vars.liquid_saturation = S_L + eps;
            auto const lambda_plus =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            double const Dlambda = (lambda_plus - lambda_minus) / 2 / eps;

            ASSERT_LE(std::abs(dlambda - Dlambda), 1e-9)
                << "for conductivity " << lambda << " and saturation " << S_L;
        }
    }
}
TEST(MaterialPropertyLib, SaturationWeightedThermalConductivity1Dgeometric)
{
    ParameterLib::ConstantParameter<double> const k_dry{"k_dry", {0.2}};
    ParameterLib::ConstantParameter<double> const k_wet{"k_wet", {1.5}};

    auto const k_model_eff = MPL::SaturationWeightedThermalConductivity<
        MaterialPropertyLib::MeanType::GEOMETRIC, 1>(
        "thermal_conductivity", k_dry, k_wet);

    ParameterLib::SpatialPosition const pos;
    double const t = std::numeric_limits<double>::quiet_NaN();
    double const dt = std::numeric_limits<double>::quiet_NaN();
    MPL::VariableArray vars;

    {
        double const S_0 = -0.1;
        double const S_max = 1.1;
        int const n_steps = 10000;
        double const eps = 1e-6;
        for (int i = 0; i <= n_steps; i++)
        {
            double const S_L = S_0 + i * (S_max - S_0) / n_steps;
            vars.liquid_saturation = S_L;

            auto const lambda =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            switch (i)
            {
                case 0:
                    ASSERT_LE(std::abs(lambda - 0.2), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
                case 5000:
                    ASSERT_LE(std::abs(lambda - 0.5477225575051661), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
                case 10000:
                    ASSERT_LE(std::abs(lambda - 1.5), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
            }

            auto const dlambda = std::get<double>(k_model_eff.dValue(
                vars, MaterialPropertyLib::Variable::liquid_saturation, pos, t,
                dt));

            vars.liquid_saturation = S_L - eps;
            auto const lambda_minus =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            vars.liquid_saturation = S_L + eps;
            auto const lambda_plus =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            double const Dlambda = (lambda_plus - lambda_minus) / 2 / eps;

            ASSERT_LE(std::abs(dlambda - Dlambda), 1e-9)
                << "for conductivity " << dlambda << Dlambda << lambda
                << " and saturation " << S_L;
        }
    }
}
TEST(MaterialPropertyLib, SaturationWeightedThermalConductivity3Dgeometric)
{
    ParameterLib::ConstantParameter<double> const k_dry{"k_dry", {0.2}};
    ParameterLib::ConstantParameter<double> const k_wet{"k_wet", {1.5}};

    auto const k_model_eff = MPL::SaturationWeightedThermalConductivity<
        MaterialPropertyLib::MeanType::GEOMETRIC, 3>(
        "thermal_conductivity", k_dry, k_wet);

    ParameterLib::SpatialPosition const pos;
    double const t = std::numeric_limits<double>::quiet_NaN();
    double const dt = std::numeric_limits<double>::quiet_NaN();
    MPL::VariableArray vars;

    {
        double const S_0 = -0.1;
        double const S_max = 1.1;
        int const n_steps = 10000;
        double const eps = 1e-6;
        for (int i = 0; i <= n_steps; i++)
        {
            double const S_L = S_0 + i * (S_max - S_0) / n_steps;
            vars.liquid_saturation = S_L;

            auto const lambda =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            switch (i)
            {
                case 0:
                    ASSERT_LE(std::abs(lambda - 0.2), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
                case 5000:
                    ASSERT_LE(std::abs(lambda - 0.5477225575051661), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
                case 10000:
                    ASSERT_LE(std::abs(lambda - 1.5), 1e-9)
                        << "for conductivity " << lambda << " and saturation "
                        << S_L;
                    break;
            }

            auto const dlambda = std::get<double>(k_model_eff.dValue(
                vars, MaterialPropertyLib::Variable::liquid_saturation, pos, t,
                dt));

            vars.liquid_saturation = S_L - eps;
            auto const lambda_minus =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            vars.liquid_saturation = S_L + eps;
            auto const lambda_plus =
                std::get<double>(k_model_eff.value(vars, pos, t, dt));

            double const Dlambda = (lambda_plus - lambda_minus) / 2 / eps;

            ASSERT_LE(std::abs(dlambda - Dlambda), 1e-9)
                << "for conductivity " << dlambda << Dlambda << lambda
                << " and saturation " << S_L;
        }
    }
}