EclDefaultMaterial.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  Copyright (C) 2008-2013 by Andreas Lauser

  This file is part of the Open Porous Media project (OPM).

  OPM 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 2 of the License, or
  (at your option) any later version.

  OPM 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 for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECL_DEFAULT_MATERIAL_HPP
#define OPM_ECL_DEFAULT_MATERIAL_HPP

#include "EclDefaultMaterialParams.hpp"

#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/MathToolbox.hpp>

#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>

#include <algorithm>

namespace Opm {

template <class TraitsT,
          class GasOilMaterialLawT,
          class OilWaterMaterialLawT,
          class ParamsT = EclDefaultMaterialParams<TraitsT,
                                                   typename GasOilMaterialLawT::Params,
                                                   typename OilWaterMaterialLawT::Params> >
class EclDefaultMaterial : public TraitsT
{
public:
    typedef GasOilMaterialLawT GasOilMaterialLaw;
    typedef OilWaterMaterialLawT OilWaterMaterialLaw;

    // some safety checks
    static_assert(TraitsT::numPhases == 3,
                  "The number of phases considered by this capillary pressure "
                  "law is always three!");
    static_assert(GasOilMaterialLaw::numPhases == 2,
                  "The number of phases considered by the gas-oil capillary "
                  "pressure law must be two!");
    static_assert(OilWaterMaterialLaw::numPhases == 2,
                  "The number of phases considered by the oil-water capillary "
                  "pressure law must be two!");
    static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
                               typename OilWaterMaterialLaw::Scalar>::value,
                  "The two two-phase capillary pressure laws must use the same "
                  "type of floating point values.");

    static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
                  "The gas-oil material law must implement the two-phase saturation "
                  "only API to for the default Ecl capillary pressure law!");
    static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
                  "The oil-water material law must implement the two-phase saturation "
                  "only API to for the default Ecl capillary pressure law!");

    typedef TraitsT Traits;
    typedef ParamsT Params;
    typedef typename Traits::Scalar Scalar;

    static const int numPhases = 3;
    static const int waterPhaseIdx = Traits::wettingPhaseIdx;
    static const int oilPhaseIdx = Traits::nonWettingPhaseIdx;
    static const int gasPhaseIdx = Traits::gasPhaseIdx;

    static const bool implementsTwoPhaseApi = false;

    static const bool implementsTwoPhaseSatApi = false;

    static const bool isSaturationDependent = true;

    static const bool isPressureDependent = false;

    static const bool isTemperatureDependent = false;

    static const bool isCompositionDependent = false;

    template <class ContainerT, class FluidState>
    static void capillaryPressures(ContainerT &values,
                                   const Params &params,
                                   const FluidState &state)
    {
        typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
        values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
        values[oilPhaseIdx] = 0;
        values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);

        Valgrind::CheckDefined(values[gasPhaseIdx]);
        Valgrind::CheckDefined(values[oilPhaseIdx]);
        Valgrind::CheckDefined(values[waterPhaseIdx]);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation pcgn(const Params &params,
                           const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const auto& Sw = 1.0 - FsToolbox::template toLhs<Evaluation>(fs.saturation(gasPhaseIdx));
        return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation pcnw(const Params &params,
                           const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(waterPhaseIdx));
        return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
    }

    template <class ContainerT, class FluidState>
    static void saturations(ContainerT& /*values*/,
                            const Params& /*params*/,
                            const FluidState& /*fluidState*/)
    {
        OPM_THROW(std::logic_error, "Not implemented: saturations()");
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sg(const Params& /*params*/,
                         const FluidState& /*fluidState*/)
    {
        OPM_THROW(std::logic_error, "Not implemented: Sg()");
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sn(const Params& /*params*/,
                         const FluidState& /*fluidState*/)
    {
        OPM_THROW(std::logic_error, "Not implemented: Sn()");
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sw(const Params& /*params*/,
                         const FluidState& /*fluidState*/)
    {
        OPM_THROW(std::logic_error, "Not implemented: Sw()");
    }

    template <class ContainerT, class FluidState>
    static void relativePermeabilities(ContainerT &values,
                                       const Params &params,
                                       const FluidState &fluidState)
    {
        typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;

        values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
        values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
        values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krg(const Params &params,
                          const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const Evaluation& Sw = 1 - FsToolbox::template toLhs<Evaluation>(fluidState.saturation(gasPhaseIdx));
        return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krw(const Params &params,
                          const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const Evaluation& Sw = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(waterPhaseIdx));
        return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krn(const Params &params,
                          const FluidState &fluidState)
    {
        typedef MathToolbox<Evaluation> Toolbox;
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        Scalar Swco = params.Swl();

        Evaluation Sw =
            Toolbox::max(Evaluation(Swco),
                         FsToolbox::template toLhs<Evaluation>(fluidState.saturation(waterPhaseIdx)));
        Evaluation Sg = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(gasPhaseIdx));

        Evaluation Sw_ow = Sg + Sw;
        Evaluation So_go = 1.0 - Sw_ow;
        const Evaluation& kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
        const Evaluation& kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);

        // avoid the division by zero: chose a regularized kro which is used if Sw - Swco
        // < epsilon/2 and interpolate between the oridinary and the regularized kro between
        // epsilon and epsilon/2
        const Scalar epsilon = 1e-5;
        if (Toolbox::value(Sw_ow) - Swco < epsilon) {
            Evaluation kro2 = (kro_ow + kro_go)/2;;
            if (Toolbox::value(Sw_ow) - Swco > epsilon/2) {
                Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
                Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
                return kro2*alpha + kro1*(1 - alpha);
            }

            return kro2;
        }
        else
            return (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
    }

    template <class FluidState>
    static void updateHysteresis(Params &params, const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
        Scalar So = FsToolbox::value(fluidState.saturation(oilPhaseIdx));
        Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx));

        if (params.inconsistentHysteresisUpdate()) {
            Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
            // NOTE: the saturations which are passed to update the hysteresis curves are
            // inconsistent with the ones used to calculate the relative permabilities. We do
            // it like this anyway because (a) the saturation functions of opm-core do it
            // this way (b) the simulations seem to converge better (which is not too much
            // surprising actually, because the time step does not start on a kink in the
            // solution) and (c) the Eclipse 100 simulator may do the same.
            //
            // Though be aware that from a physical perspective this is definitively
            // incorrect!
            params.oilWaterParams().update(/*pcSw=*/1 - So, /*krwSw=*/1 - So, /*krn_Sw=*/1 - So);
            params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg);
        }
        else {
            Scalar Swco = params.Swl();
            Sw = std::min(Scalar(1.0), std::max(Scalar(0.0), Sw));
            Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));

            Scalar Sw_ow = Sg + std::max(Swco, Sw);
            Scalar So_go = 1 + Sw_ow;

            params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
            params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);
        }
    }
};
} // namespace Opm

#endif