InitStateEquil.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:
/*
  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 3 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/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef OPM_INIT_STATE_EQUIL_HPP
#define OPM_INIT_STATE_EQUIL_HPP

#include <opm/models/utils/propertysystem.hh>

#include <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>

#include <opm/simulators/utils/ParallelCommunication.hpp>

#include <array>
#include <cstddef>
#include <memory>
#include <utility>
#include <vector>
#include <string>

namespace Opm {

class EclipseState;
class EquilRecord;
class NumericalAquifers;

namespace EQUIL {

template<class Scalar> struct CellCornerData {
    std::array<Scalar, 8> X;
    std::array<Scalar, 8> Y;
    std::array<Scalar, 8> Z;
    CellCornerData() = default;

    CellCornerData(const std::array<Scalar, 8>& x,
                   const std::array<Scalar, 8>& y,
                   const std::array<Scalar, 8>& z)
        : X(x), Y(y), Z(z)
    {}
};

template<class Scalar> class EquilReg;
namespace Miscibility { template<class Scalar> class RsFunction; }

namespace Details {
template <class Scalar, class RHS>
class RK4IVP
{
public:
    RK4IVP(const RHS& f,
           const std::array<Scalar,2>& span,
           const Scalar y0,
           const int N);

    Scalar operator()(const Scalar x) const;

private:
    int N_;
    std::array<Scalar,2> span_;
    std::vector<Scalar>  y_;
    std::vector<Scalar>  f_;

    Scalar stepsize() const;
};

namespace PhasePressODE {
template <class FluidSystem>
class Water
{
    using Scalar = typename FluidSystem::Scalar;
    using TabulatedFunction = Tabulated1DFunction<Scalar>;

public:
    Water(const TabulatedFunction& tempVdTable,
          const TabulatedFunction& saltVdTable,
          const int pvtRegionIdx,
          const Scalar normGrav);

    Scalar operator()(const Scalar depth,
                      const Scalar press) const;

private:
    const TabulatedFunction& tempVdTable_;
    const TabulatedFunction& saltVdTable_;
    const int pvtRegionIdx_;
    const Scalar g_;

    Scalar density(const Scalar depth,
                   const Scalar press) const;
};

template <class FluidSystem, class RS>
class Oil
{
    using Scalar = typename FluidSystem::Scalar;
    using TabulatedFunction = Tabulated1DFunction<Scalar>;

public:
    Oil(const TabulatedFunction& tempVdTable,
        const RS& rs,
        const int pvtRegionIdx,
        const Scalar normGrav);

    Scalar operator()(const Scalar depth,
                      const Scalar press) const;

private:
    const TabulatedFunction& tempVdTable_;
    const RS& rs_;
    const int pvtRegionIdx_;
    const Scalar g_;

    Scalar density(const Scalar depth,
                   const Scalar press) const;
};

template <class FluidSystem, class RV, class RVW>
class Gas
{
    using Scalar = typename FluidSystem::Scalar;
    using TabulatedFunction = Tabulated1DFunction<Scalar>;

public:
    Gas(const TabulatedFunction& tempVdTable,
        const RV& rv,
        const RVW& rvw,
        const int pvtRegionIdx,
        const Scalar normGrav);

    Scalar operator()(const Scalar depth,
                      const Scalar press) const;

private:
    const TabulatedFunction& tempVdTable_;
    const RV& rv_;
    const RVW& rvw_;
    const int pvtRegionIdx_;
    const Scalar g_;

    Scalar density(const Scalar depth,
                   const Scalar press) const;
};

} // namespace PhasePressODE

template <class FluidSystem, class Region>
class PressureTable
{
public:
    using Scalar = typename FluidSystem::Scalar;
    using VSpan = std::array<Scalar, 2>;

    explicit PressureTable(const Scalar gravity,
                           const int    samplePoints = 2000);

    PressureTable(const PressureTable& rhs);

    PressureTable(PressureTable&& rhs);

    PressureTable& operator=(const PressureTable& rhs);

    PressureTable& operator=(PressureTable&& rhs);

    void equilibrate(const Region& reg,
                     const VSpan&  span);

    bool oilActive() const;

    bool gasActive() const;

    bool waterActive() const;

    Scalar oil(const Scalar depth) const;

    Scalar gas(const Scalar depth) const;

    Scalar water(const Scalar depth) const;

private:
    template <class ODE>
    class PressureFunction
    {
    public:
        struct InitCond {
            Scalar depth;
            Scalar pressure;
        };

        explicit PressureFunction(const ODE&      ode,
                                  const InitCond& ic,
                                  const int       nsample,
                                  const VSpan&    span);

        PressureFunction(const PressureFunction& rhs);

        PressureFunction(PressureFunction&& rhs) = default;

        PressureFunction& operator=(const PressureFunction& rhs);

        PressureFunction& operator=(PressureFunction&& rhs);

        Scalar value(const Scalar depth) const;

    private:
        enum Direction : std::size_t { Up, Down, NumDir };

        using Distribution = Details::RK4IVP<Scalar,ODE>;
        using DistrPtr = std::unique_ptr<Distribution>;

        InitCond initial_;
        std::array<DistrPtr, Direction::NumDir> value_;
    };

    using OilPressODE = PhasePressODE::Oil<
        FluidSystem, typename Region::CalcDissolution
    >;

    using GasPressODE = PhasePressODE::Gas<
        FluidSystem, typename Region::CalcEvaporation, typename Region::CalcWaterEvaporation
    >;

    using WatPressODE = PhasePressODE::Water<FluidSystem>;

    using OPress = PressureFunction<OilPressODE>;
    using GPress = PressureFunction<GasPressODE>;
    using WPress = PressureFunction<WatPressODE>;

    using Strategy = void (PressureTable::*)
        (const Region&, const VSpan&);

    Scalar gravity_;
    int    nsample_;

    std::unique_ptr<OPress> oil_{};
    std::unique_ptr<GPress> gas_{};
    std::unique_ptr<WPress> wat_{};

    template <typename PressFunc>
    void checkPtr(const PressFunc*   phasePress,
                  const std::string& phaseName) const;

    Strategy selectEquilibrationStrategy(const Region& reg) const;

    void copyInPointers(const PressureTable& rhs);

    void equil_WOG(const Region& reg, const VSpan& span);
    void equil_GOW(const Region& reg, const VSpan& span);
    void equil_OWG(const Region& reg, const VSpan& span);

    void makeOilPressure(const typename OPress::InitCond& ic,
                         const Region&                    reg,
                         const VSpan&                     span);

    void makeGasPressure(const typename GPress::InitCond& ic,
                         const Region&                    reg,
                         const VSpan&                     span);

    void makeWatPressure(const typename WPress::InitCond& ic,
                         const Region&                    reg,
                         const VSpan&                     span);
};

// ===========================================================================

template<class Scalar>
struct PhaseQuantityValue {
    Scalar oil{0.0};
    Scalar gas{0.0};
    Scalar water{0.0};

    PhaseQuantityValue& axpy(const PhaseQuantityValue& rhs, const Scalar a)
    {
        this->oil   += a * rhs.oil;
        this->gas   += a * rhs.gas;
        this->water += a * rhs.water;

        return *this;
    }

    PhaseQuantityValue& operator/=(const Scalar x)
    {
        this->oil   /= x;
        this->gas   /= x;
        this->water /= x;

        return *this;
    }

    void reset()
    {
        this->oil = this->gas = this->water = 0.0;
    }
};

template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
class PhaseSaturations
{
public:
    using Scalar = typename FluidSystem::Scalar;
    struct Position {
        CellID cell;
        Scalar depth;
    };

    using PTable = PressureTable<FluidSystem, Region>;

    explicit PhaseSaturations(MaterialLawManager&        matLawMgr,
                              const std::vector<Scalar>& swatInit);

    PhaseSaturations(const PhaseSaturations& rhs);

    PhaseSaturations& operator=(const PhaseSaturations&) = delete;

    PhaseSaturations& operator=(PhaseSaturations&&) = delete;

    const PhaseQuantityValue<Scalar>&
    deriveSaturations(const Position& x,
                      const Region&   reg,
                      const PTable&   ptable);

    const PhaseQuantityValue<Scalar>& correctedPhasePressures() const
    {
        return this->press_;
    }

private:
    struct EvaluationPoint {
        const Position* position{nullptr};
        const Region*   region  {nullptr};
        const PTable*   ptable  {nullptr};
    };

    using FluidState = ::Opm::
        SimpleModularFluidState<Scalar, /*numPhases=*/3, /*numComponents=*/3,
                                FluidSystem,
                                /*storePressure=*/false,
                                /*storeTemperature=*/false,
                                /*storeComposition=*/false,
                                /*storeFugacity=*/false,
                                /*storeSaturation=*/true,
                                /*storeDensity=*/false,
                                /*storeViscosity=*/false,
                                /*storeEnthalpy=*/false>;

    using MaterialLaw = typename MaterialLawManager::MaterialLaw;

    using PhaseIdx = std::remove_cv_t<
        std::remove_reference_t<decltype(FluidSystem::oilPhaseIdx)>
    >;

    MaterialLawManager& matLawMgr_;

    const std::vector<Scalar>& swatInit_;

    PhaseQuantityValue<Scalar> sat_;

    PhaseQuantityValue<Scalar> press_;

    EvaluationPoint evalPt_;

    FluidState fluidState_;

    std::array<Scalar, FluidSystem::numPhases> matLawCapPress_;

    void setEvaluationPoint(const Position& x,
                            const Region&   reg,
                            const PTable&   ptable);

    void initializePhaseQuantities();

    void deriveOilSat();

    void deriveGasSat();

    void deriveWaterSat();

    void fixUnphysicalTransition();

    void accountForScaledSaturations();

    // --------------------------------------------------------------------
    // Note: Function 'applySwatInit' is non-const because the overload set
    // needs to mutate the 'matLawMgr_'.
    // --------------------------------------------------------------------

    std::pair<Scalar, bool> applySwatInit(const Scalar pcow);

    std::pair<Scalar, bool> applySwatInit(const Scalar pc, const Scalar sw);

    void computeMaterialLawCapPress();

    Scalar materialLawCapPressGasOil() const;

    Scalar materialLawCapPressOilWater() const;

    Scalar materialLawCapPressGasWater() const;

    bool isConstCapPress(const PhaseIdx phaseIdx) const;

    bool isOverlappingTransition() const;

    Scalar fromDepthTable(const Scalar   contactdepth,
                          const PhaseIdx phasePos,
                          const bool     isincr) const;

    Scalar invertCapPress(const Scalar   pc,
                          const PhaseIdx phasePos,
                          const bool     isincr) const;

    PhaseIdx oilPos() const
    {
        return FluidSystem::oilPhaseIdx;
    }

    PhaseIdx gasPos() const
    {
        return FluidSystem::gasPhaseIdx;
    }

    PhaseIdx waterPos() const
    {
        return FluidSystem::waterPhaseIdx;
    }
};

// ===========================================================================

template <typename CellRange, class Scalar>
void verticalExtent(const CellRange&      cells,
                    const std::vector<std::pair<Scalar, Scalar>>& cellZMinMax,
                    const Parallel::Communication& comm,
                    std::array<Scalar,2>& span);

template <class Scalar, class Element>
std::pair<Scalar,Scalar> cellZMinMax(const Element& element);

} // namespace Details

namespace DeckDependent {

template<class FluidSystem,
         class Grid,
         class GridView,
         class ElementMapper,
         class CartesianIndexMapper>
class InitialStateComputer
{
    using Element = typename GridView::template Codim<0>::Entity;
    using Scalar = typename FluidSystem::Scalar;
public:
    template<class MaterialLawManager>
    InitialStateComputer(MaterialLawManager& materialLawManager,
                         const EclipseState& eclipseState,
                         const Grid& grid,
                         const GridView& gridView,
                         const CartesianIndexMapper& cartMapper,
                         const Scalar grav,
                         const int num_pressure_points = 2000,
                         const bool applySwatInit = true);

    using Vec = std::vector<Scalar>;
    using PVec = std::vector<Vec>; // One per phase.

    const Vec& temperature() const { return temperature_; }
    const Vec& saltConcentration() const { return saltConcentration_; }
    const Vec& saltSaturation() const { return saltSaturation_; }
    const PVec& press() const { return pp_; }
    const PVec& saturation() const { return sat_; }
    const Vec& rs() const { return rs_; }
    const Vec& rv() const { return rv_; }
    const Vec& rvw() const { return rvw_; }

private:
    template <class RMap>
    void updateInitialTemperature_(const EclipseState& eclState, const RMap& reg);

    template <class RMap>
    void updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg);

    template <class RMap>
    void updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg);

    void updateCellProps_(const GridView& gridView,
                          const NumericalAquifers& aquifer);

    void applyNumericalAquifers_(const GridView& gridView,
                                 const NumericalAquifers& aquifer,
                                 const bool co2store_or_h2store);

    template<class RMap>
    void setRegionPvtIdx(const EclipseState& eclState, const RMap& reg);

    template <class RMap, class MaterialLawManager, class Comm>
    void calcPressSatRsRv(const RMap& reg,
                          const std::vector<EquilRecord>& rec,
                          MaterialLawManager& materialLawManager,
                          const GridView& gridView,
                          const Comm& comm,
                          const Scalar grav);

    template <class CellRange, class EquilibrationMethod>
    void cellLoop(const CellRange&      cells,
                  EquilibrationMethod&& eqmethod);

    template <class CellRange, class PressTable, class PhaseSat>
    void equilibrateCellCentres(const CellRange&        cells,
                                const EquilReg<Scalar>& eqreg,
                                const PressTable&       ptable,
                                PhaseSat&               psat);

    template <class CellRange, class PressTable, class PhaseSat>
    void equilibrateHorizontal(const CellRange&        cells,
                               const EquilReg<Scalar>& eqreg,
                               const int               acc,
                               const PressTable&       ptable,
                               PhaseSat&               psat);

     template<class CellRange, class PressTable, class PhaseSat>
     void equilibrateTiltedFaultBlock(const CellRange& cells, 
                            const EquilReg<Scalar>& eqreg,
                            const GridView& gridView, const int numLevels,
                            const PressTable& ptable, PhaseSat& psat);

     template<class CellRange, class PressTable, class PhaseSat>
     void equilibrateTiltedFaultBlockSimple(const CellRange& cells,
                           const EquilReg<Scalar>& eqreg,
                           const GridView& gridView, const int numLevels,
                           const PressTable& ptable, PhaseSat& psat);

    std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rsFunc_;
    std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rvFunc_;
    std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rvwFunc_;
    using TabulatedFunction = Tabulated1DFunction<Scalar>;
    std::vector<TabulatedFunction> tempVdTable_;
    std::vector<TabulatedFunction> saltVdTable_;
    std::vector<TabulatedFunction> saltpVdTable_;
    std::vector<int> regionPvtIdx_;
    Vec temperature_;
    Vec saltConcentration_;
    Vec saltSaturation_;
    PVec pp_;
    PVec sat_;
    Vec rs_;
    Vec rv_;
    Vec rvw_;
    const CartesianIndexMapper& cartesianIndexMapper_;
    Vec swatInit_;
    Vec cellCenterDepth_;
    std::vector<std::pair<Scalar,Scalar>> cellCenterXY_;
    std::vector<std::pair<Scalar,Scalar>> cellZSpan_;
    std::vector<std::pair<Scalar,Scalar>> cellZMinMax_;
    std::vector<CellCornerData<Scalar>> cellCorners_;
    int num_pressure_points_;
};

} // namespace DeckDependent
} // namespace EQUIL
} // namespace Opm

#endif // OPM_INIT_STATE_EQUIL_HPP