/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2017-     Statoil ASA
//
//  ResInsight is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
//  FITNESS FOR A PARTICULAR PURPOSE.
//
//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
//  for more details.
//
/////////////////////////////////////////////////////////////////////////////////

#include "RigFractureTransmissibilityEquations.h"

#include "cvfMath.h"
#include "cvfVector2.h"

#include <cmath>

const double RigFractureTransmissibilityEquations::EPSILON = 1.0e-9;

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigFractureTransmissibilityEquations::centerToCenterFractureCellTrans( double conductivityCell1,
                                                                              double sideLengthParallellTransCell1,
                                                                              double sideLengthNormalTransCell1,
                                                                              double conductivityCell2,
                                                                              double sideLengthParallellTransCell2,
                                                                              double sideLengthNormalTransCell2,
                                                                              double cDarcyForRelevantUnit )
{
    double transCell1 = centerToEdgeFractureCellTrans( conductivityCell1,
                                                       sideLengthParallellTransCell1,
                                                       sideLengthNormalTransCell1,
                                                       cDarcyForRelevantUnit );
    double transCell2 = centerToEdgeFractureCellTrans( conductivityCell2,
                                                       sideLengthParallellTransCell2,
                                                       sideLengthNormalTransCell2,
                                                       cDarcyForRelevantUnit );

    double totalTrans = 1 / ( ( 1 / transCell1 ) + ( 1 / transCell2 ) );

    return totalTrans;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigFractureTransmissibilityEquations::fractureCellToWellRadialTrans( double fractureCellConductivity,
                                                                            double fractureCellSizeX,
                                                                            double fractureCellSizeZ,
                                                                            double wellRadius,
                                                                            double skinFactor,
                                                                            double cDarcyForRelevantUnit )
{
    double ro = 0.14 * cvf::Math::sqrt( pow( fractureCellSizeX, 2.0 ) + pow( fractureCellSizeZ, 2 ) );

    if ( ro < ( wellRadius * 1.01 ) )
    {
        ro = wellRadius * 1.01;
    }

    double Tc = 2 * cvf::PI_D * cDarcyForRelevantUnit * fractureCellConductivity / ( log( ro / wellRadius ) + skinFactor );

    CVF_TIGHT_ASSERT( Tc > 0 );
    return Tc;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigFractureTransmissibilityEquations::fractureCellToWellLinearTrans( double fractureConductivity,
                                                                            double fractureCellSizeX,
                                                                            double fractureCellSizeZ,
                                                                            double perforationLengthVertical,
                                                                            double perforationLengthHorizontal,
                                                                            double perforationEfficiency,
                                                                            double skinfactor,
                                                                            double cDarcyForRelevantUnit,
                                                                            double wellRadius )
{
    const double invalidTrans = 1.0e9;
    const double epsilon      = 1.0e-8;

    double TcPrefix = 8 * cDarcyForRelevantUnit * fractureConductivity;

    cvf::Vec2d wellOrientation  = cvf::Vec2d( perforationLengthHorizontal, perforationLengthVertical ).getNormalized();
    cvf::Vec2d wellRadialVector = wellOrientation.perpendicularVector() * wellRadius;

    double DzPerf = perforationLengthVertical * perforationEfficiency;
    double DxPerf = perforationLengthHorizontal * perforationEfficiency;

    double TcZ = 0.0;
    if ( DzPerf > epsilon )
    {
        double effectiveFlowLengthHorizontal = fractureCellSizeX - 4.0 * std::abs( wellRadialVector.x() );
        double denominatorZ                  = effectiveFlowLengthHorizontal + skinfactor * DzPerf / cvf::PI_D;
        if ( denominatorZ < epsilon )
        {
            return invalidTrans;
        }
        TcZ = TcPrefix * DzPerf / denominatorZ;
    }

    double TcX = 0.0;
    if ( DxPerf > epsilon )
    {
        double effectiveFlowLengthVertical = fractureCellSizeZ - 4.0 * std::abs( wellRadialVector.y() );
        double denominatorX                = effectiveFlowLengthVertical + skinfactor * DxPerf / cvf::PI_D;
        if ( denominatorX < epsilon )
        {
            return invalidTrans;
        }
        TcX = TcPrefix * DxPerf / denominatorX;
    }

    double Tc = cvf::Math::sqrt( pow( TcX, 2 ) + pow( TcZ, 2 ) );
    return Tc;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigFractureTransmissibilityEquations::matrixToFractureTrans( double perm,
                                                                    double NTG,
                                                                    double A,
                                                                    double cellSizeLength,
                                                                    double skinfactor,
                                                                    double fractureAreaWeightedlength,
                                                                    double cDarcy )
{
    double transmissibility;

    double slDivPi = 0.0;
    if ( cvf::Math::abs( skinfactor ) > EPSILON )
    {
        slDivPi = ( skinfactor * fractureAreaWeightedlength ) / cvf::PI_D;
    }

    transmissibility = 8 * cDarcy * ( perm * NTG ) * A / ( cellSizeLength + slDivPi );

    CVF_ASSERT( transmissibility == transmissibility );
    return transmissibility;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigFractureTransmissibilityEquations::effectiveInternalFractureToWellTransPDDHC( double sumScaledMatrixToFractureTrans,
                                                                                        double scaledMatrixToWellTrans )
{
    double divisor = sumScaledMatrixToFractureTrans - scaledMatrixToWellTrans;
    if ( cvf::Math::abs( divisor ) > EPSILON )
    {
        return ( sumScaledMatrixToFractureTrans * scaledMatrixToWellTrans ) / divisor;
    }
    return 0.0;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigFractureTransmissibilityEquations::effectiveMatrixToWellTransPDDHC( double sumOriginalMatrixToFractureTrans,
                                                                              double effectiveInternalFractureToWellTrans )
{
    double divisor = sumOriginalMatrixToFractureTrans + effectiveInternalFractureToWellTrans;
    if ( cvf::Math::abs( divisor ) > EPSILON )
    {
        return ( sumOriginalMatrixToFractureTrans * effectiveInternalFractureToWellTrans ) / divisor;
    }
    return 0.0;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigFractureTransmissibilityEquations::matrixPermeability( double permx, double permy, double NTG )
{
    double permxy = cvf::Math::sqrt( permx * permy );

    return permxy * NTG;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigFractureTransmissibilityEquations::centerToEdgeFractureCellTrans( double conductivity,
                                                                            double sideLengthParallellTrans,
                                                                            double sideLengthNormalTrans,
                                                                            double cDarcyForRelevantUnit )
{
    double transmissibility =
        cDarcyForRelevantUnit * conductivity * sideLengthNormalTrans / ( sideLengthParallellTrans / 2 );
    return transmissibility;
}