Loading...
Searching...
No Matches
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> 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 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);
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>> cellZSpan_;
std::vector<std::pair<Scalar,Scalar>> cellZMinMax_;
int num_pressure_points_;
};
} // namespace DeckDependent
} // namespace EQUIL
} // namespace Opm
#endif // OPM_INIT_STATE_EQUIL_HPP
Opm::EQUIL::Details::PhasePressODE::Gas
Definition InitStateEquil.hpp:138
Opm::EQUIL::Details::PhasePressODE::Oil
Definition InitStateEquil.hpp:113
Opm::EQUIL::Details::PhasePressODE::Water
Definition InitStateEquil.hpp:88
Opm::EQUIL::Details::PhaseSaturations::deriveSaturations
const PhaseQuantityValue< Scalar > & deriveSaturations(const Position &x, const Region &reg, const PTable &ptable)
Calculate phase saturations at particular point of the simulation model geometry.
Definition InitStateEquil_impl.hpp:596
Opm::EQUIL::Details::PhaseSaturations::PhaseSaturations
PhaseSaturations(MaterialLawManager &matLawMgr, const std::vector< Scalar > &swatInit)
Constructor.
Definition InitStateEquil_impl.hpp:572
Opm::EQUIL::Details::PhaseSaturations::PTable
PressureTable< FluidSystem, Region > PTable
Convenience type alias.
Definition InitStateEquil.hpp:398
Opm::EQUIL::Details::PhaseSaturations::operator=
PhaseSaturations & operator=(const PhaseSaturations &)=delete
Disabled assignment operator.
Opm::EQUIL::Details::PhaseSaturations::operator=
PhaseSaturations & operator=(PhaseSaturations &&)=delete
Disabled move-assignment operator.
Opm::EQUIL::Details::PhaseSaturations::correctedPhasePressures
const PhaseQuantityValue< Scalar > & correctedPhasePressures() const
Retrieve saturation-corrected phase pressures.
Definition InitStateEquil.hpp:442
Opm::EQUIL::Details::PressureTable
Definition InitStateEquil.hpp:167
Opm::EQUIL::Details::PressureTable::operator=
PressureTable & operator=(const PressureTable &rhs)
Assignment operator.
Definition InitStateEquil_impl.hpp:1017
Opm::EQUIL::Details::PressureTable::water
Scalar water(const Scalar depth) const
Evaluate water phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1097
Opm::EQUIL::Details::PressureTable::gas
Scalar gas(const Scalar depth) const
Evaluate gas phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1086
Opm::EQUIL::Details::PressureTable::waterActive
bool waterActive() const
Predicate for whether or not water is an active phase.
Definition InitStateEquil_impl.hpp:1068
Opm::EQUIL::Details::PressureTable::gasActive
bool gasActive() const
Predicate for whether or not gas is an active phase.
Definition InitStateEquil_impl.hpp:1061
Opm::EQUIL::Details::PressureTable::oil
Scalar oil(const Scalar depth) const
Evaluate oil phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1076
Opm::EQUIL::Details::PressureTable::oilActive
bool oilActive() const
Predicate for whether or not oil is an active phase.
Definition InitStateEquil_impl.hpp:1054
Opm::EQUIL::Details::PressureTable::PressureTable
PressureTable(const Scalar gravity, const int samplePoints=2000)
Constructor.
Definition InitStateEquil_impl.hpp:987
Opm::EQUIL::Details::RK4IVP
Definition InitStateEquil.hpp:67
Opm::EQUIL::EquilReg
Aggregate information base of an equilibration region.
Definition EquilibrationHelpers.hpp:627
Opm::EQUIL::Miscibility::RsFunction
Base class for phase mixing functions.
Definition EquilibrationHelpers.hpp:101
Opm::EQUIL::Miscibility
Types and routines relating to phase mixing in equilibration calculations.
Definition EquilibrationHelpers.hpp:94
Opm::EQUIL
Types and routines that collectively implement a basic ECLIPSE-style equilibration-based initialisati...
Definition EquilibrationHelpers.cpp:28
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
propertysystem.hh
The Opm property system, traits with inheritance.
Opm::EQUIL::Details::PhaseQuantityValue
Simple set of per-phase (named by primary component) quantities.
Definition InitStateEquil.hpp:338
Opm::EQUIL::Details::PhaseSaturations::Position
Evaluation point within a model geometry.
Definition InitStateEquil.hpp:392
Opm::EQUIL::Details::PressureTable::PressureFunction::InitCond
Definition InitStateEquil.hpp:255