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

#include "TwoPhaseFlowWithPrhoMaterialProperties.h"
#include <utility>
#include "BaseLib/Logging.h"
#include "MaterialLib/Fluid/FluidProperty.h"
#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
#include "MaterialLib/PorousMedium/Storage/Storage.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/PropertyVector.h"
#include "NumLib/NewtonRaphson.h"
#include "ParameterLib/Parameter.h"
#include "ParameterLib/SpatialPosition.h"

using MaterialLib::PhysicalConstant::MolarMass::H2;
using MaterialLib::PhysicalConstant::IdealGasConstant;
using MaterialLib::PhysicalConstant::HenryConstant::HenryConstantH2;
namespace ProcessLib
{
namespace TwoPhaseFlowWithPrho
{
TwoPhaseFlowWithPrhoMaterialProperties::TwoPhaseFlowWithPrhoMaterialProperties(
    boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
        liquid_density,
    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
        viscosity,
    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
        gas_density,
    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
        gas_viscosity,
    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>&&
        intrinsic_permeability_models,
    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
        porosity_models,
    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
        storage_models,
    std::vector<std::unique_ptr<
        MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
        capillary_pressure_models,
    std::vector<
        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
        relative_permeability_models)
    : _liquid_density(std::move(liquid_density)),
      _viscosity(std::move(viscosity)),
      _gas_density(std::move(gas_density)),
      _gas_viscosity(std::move(gas_viscosity)),
      _material_ids(material_ids),
      _intrinsic_permeability_models(std::move(intrinsic_permeability_models)),
      _porosity_models(std::move(porosity_models)),
      _storage_models(std::move(storage_models)),
      _capillary_pressure_models(std::move(capillary_pressure_models)),
      _relative_permeability_models(std::move(relative_permeability_models))
{
    DBUG("Create material properties for Two-Phase flow with P-RHO model.");
}

int TwoPhaseFlowWithPrhoMaterialProperties::getMaterialID(
    const std::size_t element_id)
{
    if (!_material_ids)
    {
        return 0;
    }

    assert(element_id < _material_ids->size());
    return (*_material_ids)[element_id];
}

double TwoPhaseFlowWithPrhoMaterialProperties::getLiquidDensity(
    const double p, const double T) const
{
    ArrayType vars;
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
    return _liquid_density->getValue(vars);
}

double TwoPhaseFlowWithPrhoMaterialProperties::getGasDensity(
    const double p, const double T) const
{
    ArrayType vars;
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
    return _gas_density->getValue(vars);
}

double TwoPhaseFlowWithPrhoMaterialProperties::getLiquidViscosity(
    const double p, const double T) const
{
    ArrayType vars;
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
    return _viscosity->getValue(vars);
}

double TwoPhaseFlowWithPrhoMaterialProperties::getGasViscosity(
    const double p, const double T) const
{
    ArrayType vars;
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
    return _gas_viscosity->getValue(vars);
}

Eigen::MatrixXd TwoPhaseFlowWithPrhoMaterialProperties::getPermeability(
    const int material_id, const double t,
    const ParameterLib::SpatialPosition& pos, const int /*dim*/) const
{
    return _intrinsic_permeability_models[material_id]->getValue(t, pos, 0.0,
                                                                 0.0);
}

double TwoPhaseFlowWithPrhoMaterialProperties::getPorosity(
    const int material_id, const double t,
    const ParameterLib::SpatialPosition& pos, const double /*p*/,
    const double T, const double porosity_variable) const
{
    return _porosity_models[material_id]->getValue(t, pos, porosity_variable,
                                                   T);
}

double TwoPhaseFlowWithPrhoMaterialProperties::getNonwetRelativePermeability(
    const double /*t*/, const ParameterLib::SpatialPosition& /*pos*/,
    const double /*p*/, const double /*T*/, const double saturation) const
{
    return _relative_permeability_models[0]->getValue(saturation);
}

double TwoPhaseFlowWithPrhoMaterialProperties::getWetRelativePermeability(
    const double /*t*/, const ParameterLib::SpatialPosition& /*pos*/,
    const double /*p*/, const double /*T*/, const double saturation) const
{
    return _relative_permeability_models[1]->getValue(saturation);
}

double TwoPhaseFlowWithPrhoMaterialProperties::getCapillaryPressure(
    const int material_id, const double /*t*/,
    const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
    const double /*T*/, const double saturation) const
{
    return _capillary_pressure_models[material_id]->getCapillaryPressure(
        saturation);
}

double TwoPhaseFlowWithPrhoMaterialProperties::getCapillaryPressureDerivative(
    const int material_id, const double /*t*/,
    const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
    const double /*T*/, const double saturation) const
{
    return _capillary_pressure_models[material_id]->getdPcdS(saturation);
}

bool TwoPhaseFlowWithPrhoMaterialProperties::computeConstitutiveRelation(
    double const t,
    ParameterLib::SpatialPosition const& x,
    int const material_id,
    double const pg,
    double const X,
    double const T,
    double& Sw,
    double& X_m,
    double& dsw_dpg,
    double& dsw_dX,
    double& dxm_dpg,
    double& dxm_dX)
{
    {  // Local Newton solver
        using LocalJacobianMatrix =
            Eigen::Matrix<double, 2, 2, Eigen::RowMajor>;
        using LocalResidualVector = Eigen::Matrix<double, 2, 1>;
        using LocalUnknownVector = Eigen::Matrix<double, 2, 1>;
        LocalJacobianMatrix J_loc;

        Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(2);
        auto const update_residual = [&](LocalResidualVector& residual) {
            calculateResidual(material_id, pg, X, T, Sw, X_m, residual);
        };

        auto const update_jacobian = [&](LocalJacobianMatrix& jacobian) {
            calculateJacobian(material_id, t, x, pg, X, T, jacobian, Sw,
                              X_m);  // for solution dependent Jacobians
        };

        auto const update_solution = [&](LocalUnknownVector const& increment) {
            // increment solution vectors
            Sw += increment[0];
            X_m += increment[1];
        };

        // TODO Make the following choice of maximum iterations and convergence
        // criteria available from the input file configuration. See Ehlers
        // material model implementation for the example.
        const int maximum_iterations(20);
        const double residuum_tolerance(1.e-14);
        const double increment_tolerance(0);

        auto newton_solver = NumLib::NewtonRaphson<
            decltype(linear_solver), LocalJacobianMatrix,
            decltype(update_jacobian), LocalResidualVector,
            decltype(update_residual), decltype(update_solution)>(
            linear_solver, update_jacobian, update_residual, update_solution,
            {maximum_iterations, residuum_tolerance, increment_tolerance});

        auto const success_iterations = newton_solver.solve(J_loc);

        if (!success_iterations)
        {
            return false;
        }
    }
    dsw_dpg = calculatedSwdP(pg, Sw, X_m, T, material_id);
    dsw_dX = calculatedSwdX(pg, X, Sw, X_m, T, material_id);
    dxm_dpg = calculatedXmdP(pg, Sw, X_m, dsw_dpg, material_id);
    dxm_dX = calculatedXmdX(pg, Sw, X_m, dsw_dX, material_id);
    return true;
}
void TwoPhaseFlowWithPrhoMaterialProperties::calculateResidual(
    const int material_id, double const pl, double const X, double const T,
    double const Sw, double const rho_h2_wet, ResidualVector& res)
{
    const double pg =
        pl + _capillary_pressure_models[material_id]->getCapillaryPressure(Sw);
    const double rho_h2_nonwet = pg * H2 / IdealGasConstant / T;

    // calculating residual
    res(0) = calculateEquilibiumRhoWetLight(pg, Sw, rho_h2_wet);
    res(1) = calculateSaturation(pl, X, Sw, rho_h2_wet, rho_h2_nonwet, T);
}

void TwoPhaseFlowWithPrhoMaterialProperties::calculateJacobian(
    const int material_id, double const /*t*/,
    ParameterLib::SpatialPosition const& /*x*/, double const pl,
    double const /*X*/, double const T, JacobianMatrix& Jac, double const Sw,
    double const rho_h2_wet)
{
    const double pg =
        pl + _capillary_pressure_models[material_id]->getCapillaryPressure(Sw);
    const double rho_h2_nonwet = pg * H2 / IdealGasConstant / T;
    double const rho_equili_h2_wet = pg * HenryConstantH2 * H2;
    double const dPC_dSw =
        _capillary_pressure_models[material_id]->getdPcdS(Sw);
    double const drhoh2wet_dpg = HenryConstantH2 * H2;
    Jac.setZero();
    if ((1 - Sw) < (rho_equili_h2_wet - rho_h2_wet))
    {
        Jac(0, 0) = -1;
    }
    else
    {
        Jac(0, 0) = drhoh2wet_dpg * dPC_dSw;
        Jac(0, 1) = -1;
    }

    Jac(1, 0) = rho_h2_nonwet - rho_h2_wet;
    Jac(1, 1) = -Sw;
}

/** Complementary condition 1
* for calculating molar fraction of light component in the liquid phase
*/
double TwoPhaseFlowWithPrhoMaterialProperties::calculateEquilibiumRhoWetLight(
    double const pg, double const Sw, double const rho_wet_h2)
{
    double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
    return std::min(1 - Sw, rho_equilibrium_wet_h2 - rho_wet_h2);
}

/** Complementary condition 2
* for calculating the saturation
*/
double TwoPhaseFlowWithPrhoMaterialProperties::calculateSaturation(
    double /*PL*/, double X, double Sw, double rho_wet_h2, double rho_nonwet_h2,
    double /*T*/)
{
    return X - (Sw * rho_wet_h2 + (1 - Sw) * rho_nonwet_h2);
}

/**
* Calculate the derivatives using the analytical way
*/
double TwoPhaseFlowWithPrhoMaterialProperties::calculatedSwdP(
    double pl, double S, double rho_wet_h2, double const T,
    int current_material_id) const
{
    const double pg =
        pl +
        _capillary_pressure_models[current_material_id]->getCapillaryPressure(
            S);
    double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
    if ((1 - S) < (rho_equilibrium_wet_h2 - rho_wet_h2))
    {
        return 0.0;
    }
    double const drhoh2wet_dpg = HenryConstantH2 * H2;
    double const drhoh2nonwet_dpg = H2 / IdealGasConstant / T;
    double const alpha =
        ((drhoh2nonwet_dpg - drhoh2wet_dpg) * (1 - S) + drhoh2wet_dpg);
    double const beta = (drhoh2nonwet_dpg - drhoh2wet_dpg) *
                        pg;  // NOTE here should be PG^h, but we ignore vapor
    double const dPC_dSw =
        _capillary_pressure_models[current_material_id]->getdPcdS(S);
    return alpha / (beta - alpha * dPC_dSw);
}
/**
* Calculate the derivatives using the analytical way
*/
double TwoPhaseFlowWithPrhoMaterialProperties::calculatedSwdX(
    double const pl, const double /*X*/, const double S,
    const double rho_wet_h2, double const T, int current_material_id) const
{
    const double pg =
        pl +
        _capillary_pressure_models[current_material_id]->getCapillaryPressure(
            S);
    double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
    if ((1 - S) < (rho_equilibrium_wet_h2 - rho_wet_h2))
    {
        return 0.0;
    }
    double const drhoh2wet_dpg = HenryConstantH2 * H2;
    double const drhoh2nonwet_dpg = H2 / IdealGasConstant / T;
    double const alpha =
        ((drhoh2nonwet_dpg - drhoh2wet_dpg) * (1 - S) + drhoh2wet_dpg);
    double const beta = (drhoh2nonwet_dpg - drhoh2wet_dpg) *
                        pg;  // NOTE here should be PG^h, but we ignore vapor
    double const dPC_dSw =
        _capillary_pressure_models[current_material_id]->getdPcdS(S);
    return -1 / (beta - alpha * dPC_dSw);
}
/**
* Calculate the derivatives using the analytical way
*/
double TwoPhaseFlowWithPrhoMaterialProperties::calculatedXmdX(
    double pl, double Sw, double rho_wet_h2, double dSwdX,
    int current_material_id) const
{
    const double pg =
        pl +
        _capillary_pressure_models[current_material_id]->getCapillaryPressure(
            Sw);
    double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
    double const dPC_dSw =
        _capillary_pressure_models[current_material_id]->getdPcdS(Sw);
    if ((1 - Sw) < (rho_equilibrium_wet_h2 - rho_wet_h2))
    {
        return 1.0;
    }
    return HenryConstantH2 * H2 * dPC_dSw * dSwdX;
}
/**
* Calculate the derivatives using the analytical way
*/
double TwoPhaseFlowWithPrhoMaterialProperties::calculatedXmdP(
    double pl, double Sw, double rho_wet_h2, double dSwdP,
    int current_material_id) const
{
    const double pg =
        pl +
        _capillary_pressure_models[current_material_id]->getCapillaryPressure(
            Sw);
    double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
    double const dPC_dSw =
        _capillary_pressure_models[current_material_id]->getdPcdS(Sw);
    if ((1 - Sw) < (rho_equilibrium_wet_h2 - rho_wet_h2))
    {
        return 0.0;
    }
    return HenryConstantH2 * H2 * (1 + dPC_dSw * dSwdP);
}
}  // namespace TwoPhaseFlowWithPrho
}  // namespace ProcessLib