struct PreconditionerTraits
{
- typedef std::shared_ptr PointerType;
+ typedef typename SmootherChooser::Type SmootherType;
+ typedef std::shared_ptr PointerType;
};
template
- typename PreconditionerTraits::PointerType
+ typename PreconditionerTraits::PointerType
makePreconditioner(O& opA, double relax, const C& comm, int iterations=1)
{
- typename Dune::Amg::SmootherTraits::Arguments args;
- typename Dune::Amg::ConstructionTraits
::Arguments cargs;
+ typedef typename SmootherChooser::Type SmootherType;
+ typedef typename PreconditionerTraits::PointerType PointerType;
+ typename Dune::Amg::SmootherTraits::Arguments args;
+ typename Dune::Amg::ConstructionTraits::Arguments cargs;
cargs.setMatrix(opA.getmat());
args.iterations=iterations;
args.relaxationFactor=relax;
cargs.setArgs(args);
cargs.setComm(comm);
- return std::shared_ptr(Dune::Amg::ConstructionTraits
::construct(cargs));
+ return PointerType(Dune::Amg::ConstructionTraits::construct(cargs));
}
template
@@ -323,19 +371,20 @@ namespace Opm
#endif
#if SMOOTHER_ILU
- typedef Dune::SeqILU0 Smoother;
+ typedef Dune::SeqILU0 SeqSmoother;
#else
- typedef Dune::SeqSOR Smoother;
+ typedef Dune::SeqSOR SeqSmoother;
#endif
+ typedef typename SmootherChooser::Type Smoother;
typedef Dune::Amg::CoarsenCriterion Criterion;
- typedef Dune::Amg::AMG Precond;
+ typedef Dune::Amg::AMG Precond;
// Construct preconditioner.
Criterion criterion;
- Precond::SmootherArgs smootherArgs;
+ typename Precond::SmootherArgs smootherArgs;
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity,
linsolver_smooth_steps);
- Precond precond(opA, criterion, smootherArgs);
+ Precond precond(opA, criterion, smootherArgs, comm);
// Construct linear solver.
Dune::CGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity);
@@ -360,6 +409,7 @@ namespace Opm
double linsolver_prolongate_factor, int linsolver_smooth_steps)
{
// Solve with AMG solver.
+ Dune::MatrixAdapter sOpA(opA.getmat());
#if FIRST_DIAGONAL
typedef Dune::Amg::FirstDiagonal CouplingMetric;
@@ -386,10 +436,10 @@ namespace Opm
Criterion criterion;
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity,
linsolver_smooth_steps);
- Precond precond(opA, criterion, smootherArgs);
+ Precond precond(sOpA, criterion, smootherArgs);
// Construct linear solver.
- Dune::GeneralizedPCGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity);
+ Dune::GeneralizedPCGSolver linsolve(sOpA, precond, tolerance, maxit, verbosity);
// Solve system.
Dune::InverseOperatorResult result;
@@ -409,6 +459,8 @@ namespace Opm
double linsolver_prolongate_factor)
{
// Solve with AMG solver.
+ typedef Dune::MatrixAdapter Operator;
+ Operator sOpA(opA.getmat());
#if FIRST_DIAGONAL
typedef Dune::Amg::FirstDiagonal CouplingMetric;
@@ -434,10 +486,10 @@ namespace Opm
parms.setNoPreSmoothSteps(smooth_steps);
parms.setNoPostSmoothSteps(smooth_steps);
parms.setProlongationDampingFactor(linsolver_prolongate_factor);
- Precond precond(opA, criterion, parms);
+ Precond precond(sOpA, criterion, parms);
// Construct linear solver.
- Dune::GeneralizedPCGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity);
+ Dune::GeneralizedPCGSolver linsolve(sOpA, precond, tolerance, maxit, verbosity);
// Solve system.
Dune::InverseOperatorResult result;
diff --git a/opm/core/linalg/LinearSolverIstl.hpp b/opm/core/linalg/LinearSolverIstl.hpp
index 4ab853cc7..4a948c087 100644
--- a/opm/core/linalg/LinearSolverIstl.hpp
+++ b/opm/core/linalg/LinearSolverIstl.hpp
@@ -24,12 +24,11 @@
#include
#include
#include
-
+#include
namespace Opm
{
-
/// Concrete class encapsulating some dune-istl linear solvers.
class LinearSolverIstl : public LinearSolverInterface
{
@@ -75,7 +74,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
- double* solution) const;
+ double* solution,
+ const boost::any& comm=boost::any()) const;
/// Set tolerance for the residual in dune istl linear solver.
/// \param[in] tol tolerance value
diff --git a/opm/core/linalg/LinearSolverUmfpack.cpp b/opm/core/linalg/LinearSolverUmfpack.cpp
index 8980d4b19..e7067b69f 100644
--- a/opm/core/linalg/LinearSolverUmfpack.cpp
+++ b/opm/core/linalg/LinearSolverUmfpack.cpp
@@ -46,7 +46,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
- double* solution) const
+ double* solution,
+ const boost::any&) const
{
CSRMatrix A = {
(size_t)size,
diff --git a/opm/core/linalg/LinearSolverUmfpack.hpp b/opm/core/linalg/LinearSolverUmfpack.hpp
index f2562a156..7126bc4a6 100644
--- a/opm/core/linalg/LinearSolverUmfpack.hpp
+++ b/opm/core/linalg/LinearSolverUmfpack.hpp
@@ -55,7 +55,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
- double* solution) const;
+ double* solution,
+ const boost::any& add=boost::any()) const;
/// Set tolerance for the linear solver.
/// \param[in] tol tolerance value
diff --git a/opm/core/linalg/ParallelIstlInformation.hpp b/opm/core/linalg/ParallelIstlInformation.hpp
new file mode 100644
index 000000000..38795774c
--- /dev/null
+++ b/opm/core/linalg/ParallelIstlInformation.hpp
@@ -0,0 +1,160 @@
+/*
+ Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
+
+ 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 .
+*/
+#ifndef OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
+#define OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
+
+#if HAVE_MPI && HAVE_DUNE_ISTL
+
+#include "mpi.h"
+#include
+#include
+#include
+#include
+
+#include
+
+namespace Opm
+{
+
+/// \brief Class that encapsulates the parallelization information needed by the
+/// ISTL solvers.
+class ParallelISTLInformation
+{
+public:
+ /// \brief The type of the parallel index set used.
+ typedef Dune::OwnerOverlapCopyCommunication::ParallelIndexSet ParallelIndexSet;
+ /// \brief The type of the remote indices information used.
+ typedef Dune::OwnerOverlapCopyCommunication::RemoteIndices RemoteIndices;
+
+ /// \brief Constructs an empty parallel information object using MPI_COMM_WORLD
+ ParallelISTLInformation()
+ : indexSet_(new ParallelIndexSet),
+ remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
+ communicator_(MPI_COMM_WORLD)
+ {}
+ /// \brief Constructs an empty parallel information object using a communicator.
+ /// \param communicator The communicator to use.
+ ParallelISTLInformation(MPI_Comm communicator)
+ : indexSet_(new ParallelIndexSet),
+ remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
+ communicator_(communicator)
+ {}
+ /// \brief Constructs a parallel information object from the specified information.
+ /// \param indexSet The parallel index set to use.
+ /// \param remoteIndices The remote indices information to use.
+ /// \param communicator The communicator to use.
+ ParallelISTLInformation(const std::shared_ptr& indexSet,
+ const std::shared_ptr& remoteIndices,
+ MPI_Comm communicator)
+ : indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
+ {}
+ /// \brief Copy constructor.
+ ///
+ /// The information will be shared by the the two objects.
+ ParallelISTLInformation(const ParallelISTLInformation& other)
+ : indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
+ communicator_(other.communicator_)
+ {}
+ /// \brief Get a pointer to the underlying index set.
+ std::shared_ptr indexSet() const
+ {
+ return indexSet_;
+ }
+ /// \brief Get a pointer to the remote indices information.
+ std::shared_ptr remoteIndices() const
+ {
+ return remoteIndices_;
+ }
+ /// \brief Get the MPI communicator that we use.
+ MPI_Comm communicator() const
+ {
+ return communicator_;
+ }
+ /// \brief Copy the information stored to the specified objects.
+ /// \param[out] indexSet The object to store the index set in.
+ /// \param[out] remoteIndices The object to store the remote indices information in.
+ void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices) const
+ {
+ indexSet.beginResize();
+ IndexSetInserter inserter(indexSet);
+ std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
+ indexSet.endResize();
+ remoteIndices.rebuild();
+ }
+ /// \brief Communcate the dofs owned by us to the other process.
+ ///
+ /// Afterwards all associated dofs will contain the same data.
+ template
+ void copyOwnerToAll (const T& source, T& dest) const
+ {
+ typedef Dune::Combine,Dune::EnumItem,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet;
+ typedef Dune::Combine,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet;
+ OwnerOverlapSet sourceFlags;
+ AllSet destFlags;
+ Dune::Interface interface(communicator_);
+ if(!remoteIndices_->isSynced())
+ remoteIndices_->rebuild();
+ interface.build(*remoteIndices_,sourceFlags,destFlags);
+ Dune::BufferedCommunicator communicator;
+ communicator.template build(interface);
+ communicator.template forward >(source,dest);
+ communicator.free();
+ }
+private:
+ /** \brief gather/scatter callback for communcation */
+ template
+ struct CopyGatherScatter
+ {
+ typedef typename Dune::CommPolicy::IndexedType V;
+
+ static V gather(const T& a, std::size_t i)
+ {
+ return a[i];
+ }
+
+ static void scatter(T& a, V v, std::size_t i)
+ {
+ a[i] = v;
+ }
+ };
+ template
+ class IndexSetInserter
+ {
+ public:
+ typedef T ParallelIndexSet;
+
+ IndexSetInserter(ParallelIndexSet& indexSet)
+ : indexSet_(&indexSet)
+ {}
+ void operator()(const typename ParallelIndexSet::IndexPair& pair)
+ {
+ indexSet_->add(pair.global(), pair.local());
+ }
+
+ private:
+ ParallelIndexSet* indexSet_;
+ };
+
+ std::shared_ptr indexSet_;
+ std::shared_ptr remoteIndices_;
+ MPI_Comm communicator_;
+};
+} // end namespace Opm
+#endif
+#endif
diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
index fb6f8f94f..cf340e7ec 100644
--- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
+++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
@@ -471,8 +471,7 @@ namespace Opm
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
-
- // Initialize saturation scaling parameter
+ // Initialize saturation scaling parameters
template
template
void SaturationPropsFromDeck::initEPS(const EclipseGridParser& deck,
@@ -1003,8 +1002,9 @@ namespace Opm
bool hasENPTVD = newParserDeck->hasKeyword("ENPTVD");
bool hasENKRVD = newParserDeck->hasKeyword("ENKRVD");
int itab = 0;
- std::vector > > table_dummy;
- std::vector > >& table = table_dummy;
+ std::vector > param_col;
+ std::vector > depth_col;
+ std::vector col_names;
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
@@ -1070,22 +1070,14 @@ namespace Opm
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
- Opm::EnptvdTable enptvd(newParserDeck->getKeyword("ENPTVD"));
- table.resize(1); // only one region supported so far
- for (unsigned i = 0; i < table.size(); ++i) {
- table[i].resize(9);
-
- for (int k = 0; k < enptvd.numRows(); ++k) {
- table[i][0] = enptvd.getDepthColumn();
- table[i][1] = enptvd.getSwcoColumn();
- table[i][2] = enptvd.getSwcritColumn();
- table[i][3] = enptvd.getSwmaxColumn();
- table[i][4] = enptvd.getSgcoColumn();
- table[i][5] = enptvd.getSgcritColumn();
- table[i][6] = enptvd.getSgmaxColumn();
- table[i][7] = enptvd.getSowcritColumn();
- table[i][8] = enptvd.getSogcritColumn();
- }
+ int num_tables = newParserDeck->getKeyword("ENPTVD")->size();
+ param_col.resize(num_tables);
+ depth_col.resize(num_tables);
+ col_names.resize(9);
+ for (int table_num=0; table_numgetKeyword("ENPTVD"), col_names, table_num);
+ depth_col[table_num] = enptvd.getColumn(0); // depth
+ param_col[table_num] = enptvd.getColumn(itab); // itab=[1-8]: swl swcr swu sgl sgcr sgu sowcr sogcr
}
}
} else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
@@ -1141,21 +1133,14 @@ namespace Opm
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
- Opm::EnkrvdTable enkrvd(newParserDeck->getKeyword("ENKRVD"));
- table.resize(1); // only one region supported so far
- for (unsigned i = 0; i < table.size(); ++i) {
- table[i].resize(8);
-
- for (int k = 0; k < enkrvd.numRows(); ++k) {
- table[i][0] = enkrvd.getDepthColumn();
- table[i][1] = enkrvd.getKrwmaxColumn();
- table[i][2] = enkrvd.getKrgmaxColumn();
- table[i][3] = enkrvd.getKromaxColumn();
- table[i][4] = enkrvd.getKrwcritColumn();
- table[i][5] = enkrvd.getKrgcritColumn();
- table[i][6] = enkrvd.getKrocritgColumn();
- table[i][7] = enkrvd.getKrocritwColumn();
- }
+ int num_tables = newParserDeck->getKeyword("ENKRVD")->size();
+ param_col.resize(num_tables);
+ depth_col.resize(num_tables);
+ col_names.resize(8);
+ for (int table_num=0; table_numgetKeyword("ENKRVD"), col_names, table_num);
+ depth_col[table_num] = enkrvd.getColumn(0); // depth
+ param_col[table_num] = enkrvd.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg
}
}
}
@@ -1172,25 +1157,15 @@ namespace Opm
scaleparam[c] = val[deck_pos];
}
} else {
- // TODO for new parser
- /*
- std::cout << "--- Scaling parameter '" << keyword << "' assigned via ";
- if (keyword[0] == 'S')
- newParserDeck.getENPTVD().write(std::cout);
- else
- newParserDeck.getENKRVD().write(std::cout);
- */
const int dim = dimensions;
for (int cell = 0; cell < number_of_cells; ++cell) {
int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
- if (table[itab][jtab][0] != -1.0) {
- std::vector& depth = table[0][jtab];
- std::vector& val = table[itab][jtab];
+ if (param_col[jtab][0] >= 0.0) {
double zc = UgGridHelpers
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
dim-1);
- if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval
- scaleparam[cell] = linearInterpolation(depth, val, zc);
+ if (zc >= depth_col[jtab].front() && zc <= depth_col[jtab].back()) { //don't want extrap outside depth interval
+ scaleparam[cell] = linearInterpolation(depth_col[jtab], param_col[jtab], zc);
}
}
}
diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp
new file mode 100644
index 000000000..d89dd8c3c
--- /dev/null
+++ b/opm/core/simulator/EquilibrationHelpers.hpp
@@ -0,0 +1,829 @@
+/*
+ Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+ 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 .
+*/
+
+#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
+#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+
+/*
+---- synopsis of EquilibrationHelpers.hpp ----
+
+namespace Opm
+{
+ namespace Equil {
+
+ template
+ class DensityCalculator;
+
+ template <>
+ class DensityCalculator< BlackoilPropertiesInterface >;
+
+ namespace Miscibility {
+ class RsFunction;
+ class NoMixing;
+ class RsVD;
+ class RsSatAtContact;
+ }
+
+ struct EquilRecord;
+
+ template
+ class EquilReg;
+
+
+ struct PcEq;
+
+ inline double satFromPc(const BlackoilPropertiesInterface& props,
+ const int phase,
+ const int cell,
+ const double target_pc,
+ const bool increasing = false);
+ struct PcEqSum
+ inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
+ const int phase1,
+ const int phase2,
+ const int cell,
+ const double target_pc);
+ } // namespace Equil
+} // namespace Opm
+
+---- end of synopsis of EquilibrationHelpers.hpp ----
+*/
+
+
+namespace Opm
+{
+ /**
+ * Types and routines that collectively implement a basic
+ * ECLIPSE-style equilibration-based initialisation scheme.
+ *
+ * This namespace is intentionally nested to avoid name clashes
+ * with other parts of OPM.
+ */
+ namespace Equil {
+
+
+ template
+ class DensityCalculator;
+
+ /**
+ * Facility for calculating phase densities based on the
+ * BlackoilPropertiesInterface.
+ *
+ * Implements the crucial operator()(p,svol)
+ * function that is expected by class EquilReg.
+ */
+ template <>
+ class DensityCalculator< BlackoilPropertiesInterface > {
+ public:
+ /**
+ * Constructor.
+ *
+ * \param[in] props Implementation of the
+ * BlackoilPropertiesInterface.
+ *
+ * \param[in] c Single cell used as a representative cell
+ * in a PVT region.
+ */
+ DensityCalculator(const BlackoilPropertiesInterface& props,
+ const int c)
+ : props_(props)
+ , c_(1, c)
+ {
+ }
+
+ /**
+ * Compute phase densities of all phases at phase point
+ * given by (pressure, surface volume) tuple.
+ *
+ * \param[in] p Fluid pressure.
+ *
+ * \param[in] z Surface volumes of all phases.
+ *
+ * \return Phase densities at phase point.
+ */
+ std::vector
+ operator()(const double p,
+ const std::vector& z) const
+ {
+ const int np = props_.numPhases();
+ std::vector A(np * np, 0);
+
+ assert (z.size() == std::vector::size_type(np));
+
+ double* dAdp = 0;
+ props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
+
+ std::vector rho(np, 0.0);
+ props_.density(1, &A[0], &rho[0]);
+
+ return rho;
+ }
+
+ private:
+ const BlackoilPropertiesInterface& props_;
+ const std::vector c_;
+ };
+
+
+ /**
+ * Types and routines relating to phase mixing in
+ * equilibration calculations.
+ */
+ namespace Miscibility {
+
+ /**
+ * Base class for phase mixing functions.
+ */
+ class RsFunction
+ {
+ public:
+ /**
+ * Function call operator.
+ *
+ * \param[in] depth Depth at which to calculate RS
+ * value.
+ *
+ * \param[in] press Pressure at which to calculate RS
+ * value.
+ *
+ * \return Dissolved gas-oil ratio (RS) at depth @c
+ * depth and pressure @c press.
+ */
+ virtual double operator()(const double depth,
+ const double press,
+ const double sat = 0.0) const = 0;
+ };
+
+
+ /**
+ * Type that implements "no phase mixing" policy.
+ */
+ class NoMixing : public RsFunction {
+ public:
+ /**
+ * Function call.
+ *
+ * \param[in] depth Depth at which to calculate RS
+ * value.
+ *
+ * \param[in] press Pressure at which to calculate RS
+ * value.
+ *
+ * \return Dissolved gas-oil ratio (RS) at depth @c
+ * depth and pressure @c press. In "no mixing
+ * policy", this is identically zero.
+ */
+ double
+ operator()(const double /* depth */,
+ const double /* press */,
+ const double sat = 0.0) const
+ {
+ return 0.0;
+ }
+ };
+
+
+ /**
+ * Type that implements "dissolved gas-oil ratio"
+ * tabulated as a function of depth policy. Data
+ * typically taken from keyword 'RSVD'.
+ */
+ class RsVD : public RsFunction {
+ public:
+ /**
+ * Constructor.
+ *
+ * \param[in] props property object
+ * \param[in] cell any cell in the pvt region
+ * \param[in] depth Depth nodes.
+ * \param[in] rs Dissolved gas-oil ratio at @c depth.
+ */
+ RsVD(const BlackoilPropertiesInterface& props,
+ const int cell,
+ const std::vector& depth,
+ const std::vector& rs)
+ : props_(props)
+ , cell_(cell)
+ , depth_(depth)
+ , rs_(rs)
+ {
+ auto pu = props_.phaseUsage();
+ std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+ z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
+ z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+ }
+
+ /**
+ * Function call.
+ *
+ * \param[in] depth Depth at which to calculate RS
+ * value.
+ *
+ * \param[in] press Pressure at which to calculate RS
+ * value.
+ *
+ * \return Dissolved gas-oil ratio (RS) at depth @c
+ * depth and pressure @c press.
+ */
+ double
+ operator()(const double depth,
+ const double press,
+ const double sat_gas = 0.0) const
+ {
+ if (sat_gas > 0.0) {
+ return satRs(press);
+ } else {
+ return std::min(satRs(press), linearInterpolation(depth_, rs_, depth));
+ }
+ }
+
+ private:
+ const BlackoilPropertiesInterface& props_;
+ const int cell_;
+ std::vector depth_; /**< Depth nodes */
+ std::vector rs_; /**< Dissolved gas-oil ratio */
+ double z_[BlackoilPhases::MaxNumPhases];
+ mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+ double satRs(const double press) const
+ {
+ props_.matrix(1, &press, z_, &cell_, A_, 0);
+ // Rs/Bo is in the gas row and oil column of A_.
+ // 1/Bo is in the oil row and column.
+ // Recall also that it is stored in column-major order.
+ const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+ const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+ const int np = props_.numPhases();
+ return A_[np*opos + gpos] / A_[np*opos + opos];
+ }
+ };
+
+
+ /**
+ * Type that implements "vaporized oil-gas ratio"
+ * tabulated as a function of depth policy. Data
+ * typically taken from keyword 'RVVD'.
+ */
+ class RvVD : public RsFunction {
+ public:
+ /**
+ * Constructor.
+ *
+ * \param[in] props property object
+ * \param[in] cell any cell in the pvt region
+ * \param[in] depth Depth nodes.
+ * \param[in] rv Dissolved gas-oil ratio at @c depth.
+ */
+ RvVD(const BlackoilPropertiesInterface& props,
+ const int cell,
+ const std::vector& depth,
+ const std::vector& rv)
+ : props_(props)
+ , cell_(cell)
+ , depth_(depth)
+ , rv_(rv)
+ {
+ auto pu = props_.phaseUsage();
+ std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+ z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+ z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
+ }
+
+ /**
+ * Function call.
+ *
+ * \param[in] depth Depth at which to calculate RV
+ * value.
+ *
+ * \param[in] press Pressure at which to calculate RV
+ * value.
+ *
+ * \return Vaporized oil-gas ratio (RV) at depth @c
+ * depth and pressure @c press.
+ */
+ double
+ operator()(const double depth,
+ const double press,
+ const double sat_oil = 0.0 ) const
+ {
+ if (sat_oil > 0.0) {
+ return satRv(press);
+ } else {
+ return std::min(satRv(press), linearInterpolation(depth_, rv_, depth));
+ }
+ }
+
+ private:
+ const BlackoilPropertiesInterface& props_;
+ const int cell_;
+ std::vector depth_; /**< Depth nodes */
+ std::vector rv_; /**< Vaporized oil-gas ratio */
+ double z_[BlackoilPhases::MaxNumPhases];
+ mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+ double satRv(const double press) const
+ {
+ props_.matrix(1, &press, z_, &cell_, A_, 0);
+ // Rv/Bg is in the oil row and gas column of A_.
+ // 1/Bg is in the gas row and column.
+ // Recall also that it is stored in column-major order.
+ const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+ const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+ const int np = props_.numPhases();
+ return A_[np*gpos + opos] / A_[np*gpos + gpos];
+ }
+ };
+
+
+ /**
+ * Class that implements "dissolved gas-oil ratio" (Rs)
+ * as function of depth and pressure as follows:
+ *
+ * 1. The Rs at the gas-oil contact is equal to the
+ * saturated Rs value, Rs_sat_contact.
+ *
+ * 2. The Rs elsewhere is equal to Rs_sat_contact, but
+ * constrained to the saturated value as given by the
+ * local pressure.
+ *
+ * This should yield Rs-values that are constant below the
+ * contact, and decreasing above the contact.
+ */
+ class RsSatAtContact : public RsFunction {
+ public:
+ /**
+ * Constructor.
+ *
+ * \param[in] props property object
+ * \param[in] cell any cell in the pvt region
+ * \param[in] p_contact oil pressure at the contact
+ */
+ RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
+ : props_(props), cell_(cell)
+ {
+ auto pu = props_.phaseUsage();
+ std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+ z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
+ z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+ rs_sat_contact_ = satRs(p_contact);
+ }
+
+ /**
+ * Function call.
+ *
+ * \param[in] depth Depth at which to calculate RS
+ * value.
+ *
+ * \param[in] press Pressure at which to calculate RS
+ * value.
+ *
+ * \return Dissolved gas-oil ratio (RS) at depth @c
+ * depth and pressure @c press.
+ */
+ double
+ operator()(const double /* depth */,
+ const double press,
+ const double sat_gas = 0.0) const
+ {
+ if (sat_gas > 0.0) {
+ return satRs(press);
+ } else {
+ return std::min(satRs(press), rs_sat_contact_);
+ }
+ }
+
+ private:
+ const BlackoilPropertiesInterface& props_;
+ const int cell_;
+ double z_[BlackoilPhases::MaxNumPhases];
+ double rs_sat_contact_;
+ mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+ double satRs(const double press) const
+ {
+ props_.matrix(1, &press, z_, &cell_, A_, 0);
+ // Rs/Bo is in the gas row and oil column of A_.
+ // 1/Bo is in the oil row and column.
+ // Recall also that it is stored in column-major order.
+ const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+ const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+ const int np = props_.numPhases();
+ return A_[np*opos + gpos] / A_[np*opos + opos];
+ }
+ };
+
+
+ /**
+ * Class that implements "vaporized oil-gas ratio" (Rv)
+ * as function of depth and pressure as follows:
+ *
+ * 1. The Rv at the gas-oil contact is equal to the
+ * saturated Rv value, Rv_sat_contact.
+ *
+ * 2. The Rv elsewhere is equal to Rv_sat_contact, but
+ * constrained to the saturated value as given by the
+ * local pressure.
+ *
+ * This should yield Rv-values that are constant below the
+ * contact, and decreasing above the contact.
+ */
+ class RvSatAtContact : public RsFunction {
+ public:
+ /**
+ * Constructor.
+ *
+ * \param[in] props property object
+ * \param[in] cell any cell in the pvt region
+ * \param[in] p_contact oil pressure at the contact
+ */
+ RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
+ : props_(props), cell_(cell)
+ {
+ auto pu = props_.phaseUsage();
+ std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+ z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+ z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
+ rv_sat_contact_ = satRv(p_contact);
+ }
+
+ /**
+ * Function call.
+ *
+ * \param[in] depth Depth at which to calculate RV
+ * value.
+ *
+ * \param[in] press Pressure at which to calculate RV
+ * value.
+ *
+ * \return Dissolved oil-gas ratio (RV) at depth @c
+ * depth and pressure @c press.
+ */
+ double
+ operator()(const double /*depth*/,
+ const double press,
+ const double sat_oil = 0.0) const
+ {
+ if (sat_oil > 0.0) {
+ return satRv(press);
+ } else {
+ return std::min(satRv(press), rv_sat_contact_);
+ }
+ }
+
+ private:
+ const BlackoilPropertiesInterface& props_;
+ const int cell_;
+ double z_[BlackoilPhases::MaxNumPhases];
+ double rv_sat_contact_;
+ mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+ double satRv(const double press) const
+ {
+ props_.matrix(1, &press, z_, &cell_, A_, 0);
+ // Rv/Bg is in the oil row and gas column of A_.
+ // 1/Bg is in the gas row and column.
+ // Recall also that it is stored in column-major order.
+ const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+ const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+ const int np = props_.numPhases();
+ return A_[np*gpos + opos] / A_[np*gpos + gpos];
+ }
+ };
+
+ } // namespace Miscibility
+
+
+
+ /**
+ * Equilibration record.
+ *
+ * Layout and contents inspired by first six items of
+ * ECLIPSE's 'EQUIL' records. This is the minimum amount of
+ * input data needed to define phase pressures in an
+ * equilibration region.
+ *
+ * Data consists of three pairs of depth and pressure values:
+ * 1. main
+ * - @c depth Main datum depth.
+ * - @c press Pressure at datum depth.
+ *
+ * 2. woc
+ * - @c depth Depth of water-oil contact
+ * - @c press water-oil capillary pressure at water-oil contact.
+ * Capillary pressure defined as "P_oil - P_water".
+ *
+ * 3. goc
+ * - @c depth Depth of gas-oil contact
+ * - @c press Gas-oil capillary pressure at gas-oil contact.
+ * Capillary pressure defined as "P_gas - P_oil".
+ */
+ struct EquilRecord {
+ struct {
+ double depth;
+ double press;
+ } main, woc, goc;
+ int live_oil_table_index;
+ int wet_gas_table_index;
+ int N;
+ };
+
+ /**
+ * Aggregate information base of an equilibration region.
+ *
+ * Provides inquiry methods for retrieving depths of contacs
+ * and pressure values as well as a means of calculating fluid
+ * densities, dissolved gas-oil ratio and vapourised oil-gas
+ * ratios.
+ *
+ * \tparam DensCalc Type that provides access to a phase
+ * density calculation facility. Must implement an operator()
+ * declared as
+ *
+ * std::vector
+ * operator()(const double press,
+ * const std::vector& svol )
+ *
+ * that calculates the phase densities of all phases in @c
+ * svol at fluid pressure @c press.
+ */
+ template
+ class EquilReg {
+ public:
+ /**
+ * Constructor.
+ *
+ * \param[in] rec Equilibration data of current region.
+ * \param[in] density Density calculator of current region.
+ * \param[in] rs Calculator of dissolved gas-oil ratio.
+ * \param[in] rv Calculator of vapourised oil-gas ratio.
+ * \param[in] pu Summary of current active phases.
+ */
+ EquilReg(const EquilRecord& rec,
+ const DensCalc& density,
+ std::shared_ptr rs,
+ std::shared_ptr rv,
+ const PhaseUsage& pu)
+ : rec_ (rec)
+ , density_(density)
+ , rs_ (rs)
+ , rv_ (rv)
+ , pu_ (pu)
+ {
+ }
+
+ /**
+ * Type of density calculator.
+ */
+ typedef DensCalc CalcDensity;
+
+ /**
+ * Type of dissolved gas-oil ratio calculator.
+ */
+ typedef Miscibility::RsFunction CalcDissolution;
+
+ /**
+ * Type of vapourised oil-gas ratio calculator.
+ */
+ typedef Miscibility::RsFunction CalcEvaporation;
+
+ /**
+ * Datum depth in current region
+ */
+ double datum() const { return this->rec_.main.depth; }
+
+ /**
+ * Pressure at datum depth in current region.
+ */
+ double pressure() const { return this->rec_.main.press; }
+
+ /**
+ * Depth of water-oil contact.
+ */
+ double zwoc() const { return this->rec_.woc .depth; }
+
+ /**
+ * water-oil capillary pressure at water-oil contact.
+ *
+ * \return P_o - P_w at WOC.
+ */
+ double pcow_woc() const { return this->rec_.woc .press; }
+
+ /**
+ * Depth of gas-oil contact.
+ */
+ double zgoc() const { return this->rec_.goc .depth; }
+
+ /**
+ * Gas-oil capillary pressure at gas-oil contact.
+ *
+ * \return P_g - P_o at GOC.
+ */
+ double pcgo_goc() const { return this->rec_.goc .press; }
+
+ /**
+ * Retrieve phase density calculator of current region.
+ */
+ const CalcDensity&
+ densityCalculator() const { return this->density_; }
+
+ /**
+ * Retrieve dissolved gas-oil ratio calculator of current
+ * region.
+ */
+ const CalcDissolution&
+ dissolutionCalculator() const { return *this->rs_; }
+
+ /**
+ * Retrieve vapourised oil-gas ratio calculator of current
+ * region.
+ */
+ const CalcEvaporation&
+ evaporationCalculator() const { return *this->rv_; }
+
+ /**
+ * Retrieve active fluid phase summary.
+ */
+ const PhaseUsage&
+ phaseUsage() const { return this->pu_; }
+
+ private:
+ EquilRecord rec_; /**< Equilibration data */
+ DensCalc density_; /**< Density calculator */
+ std::shared_ptr rs_; /**< RS calculator */
+ std::shared_ptr rv_; /**< RV calculator */
+ PhaseUsage pu_; /**< Active phase summary */
+ };
+
+
+
+ /// Functor for inverting capillary pressure function.
+ /// Function represented is
+ /// f(s) = pc(s) - target_pc
+ struct PcEq
+ {
+ PcEq(const BlackoilPropertiesInterface& props,
+ const int phase,
+ const int cell,
+ const double target_pc)
+ : props_(props),
+ phase_(phase),
+ cell_(cell),
+ target_pc_(target_pc)
+ {
+ std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
+ std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
+ }
+ double operator()(double s) const
+ {
+ s_[phase_] = s;
+ props_.capPress(1, s_, &cell_, pc_, 0);
+ return pc_[phase_] - target_pc_;
+ }
+ private:
+ const BlackoilPropertiesInterface& props_;
+ const int phase_;
+ const int cell_;
+ const int target_pc_;
+ mutable double s_[BlackoilPhases::MaxNumPhases];
+ mutable double pc_[BlackoilPhases::MaxNumPhases];
+ };
+
+
+
+ /// Compute saturation of some phase corresponding to a given
+ /// capillary pressure.
+ inline double satFromPc(const BlackoilPropertiesInterface& props,
+ const int phase,
+ const int cell,
+ const double target_pc,
+ const bool increasing = false)
+ {
+ // Find minimum and maximum saturations.
+ double sminarr[BlackoilPhases::MaxNumPhases];
+ double smaxarr[BlackoilPhases::MaxNumPhases];
+ props.satRange(1, &cell, sminarr, smaxarr);
+ const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
+ const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
+
+ // Create the equation f(s) = pc(s) - target_pc
+ const PcEq f(props, phase, cell, target_pc);
+ const double f0 = f(s0);
+ const double f1 = f(s1);
+ if (f0 <= 0.0) {
+ return s0;
+ } else if (f1 > 0.0) {
+ return s1;
+ } else {
+ const int max_iter = 30;
+ const double tol = 1e-6;
+ int iter_used = -1;
+ typedef RegulaFalsi ScalarSolver;
+ const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used);
+ return sol;
+ }
+ }
+
+
+ /// Functor for inverting a sum of capillary pressure functions.
+ /// Function represented is
+ /// f(s) = pc1(s) + pc2(1 - s) - target_pc
+ struct PcEqSum
+ {
+ PcEqSum(const BlackoilPropertiesInterface& props,
+ const int phase1,
+ const int phase2,
+ const int cell,
+ const double target_pc)
+ : props_(props),
+ phase1_(phase1),
+ phase2_(phase2),
+ cell_(cell),
+ target_pc_(target_pc)
+ {
+ std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
+ std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
+ }
+ double operator()(double s) const
+ {
+ s_[phase1_] = s;
+ s_[phase2_] = 1.0 - s;
+ props_.capPress(1, s_, &cell_, pc_, 0);
+ return pc_[phase1_] + pc_[phase2_] - target_pc_;
+ }
+ private:
+ const BlackoilPropertiesInterface& props_;
+ const int phase1_;
+ const int phase2_;
+ const int cell_;
+ const int target_pc_;
+ mutable double s_[BlackoilPhases::MaxNumPhases];
+ mutable double pc_[BlackoilPhases::MaxNumPhases];
+ };
+
+
+
+
+ /// Compute saturation of some phase corresponding to a given
+ /// capillary pressure, where the capillary pressure function
+ /// is given as a sum of two other functions.
+ inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
+ const int phase1,
+ const int phase2,
+ const int cell,
+ const double target_pc)
+ {
+ // Find minimum and maximum saturations.
+ double sminarr[BlackoilPhases::MaxNumPhases];
+ double smaxarr[BlackoilPhases::MaxNumPhases];
+ props.satRange(1, &cell, sminarr, smaxarr);
+ const double smin = sminarr[phase1];
+ const double smax = smaxarr[phase1];
+
+ // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
+ const PcEqSum f(props, phase1, phase2, cell, target_pc);
+ const double f0 = f(smin);
+ const double f1 = f(smax);
+ if (f0 <= 0.0) {
+ return smin;
+ } else if (f1 > 0.0) {
+ return smax;
+ } else {
+ const int max_iter = 30;
+ const double tol = 1e-6;
+ int iter_used = -1;
+ typedef RegulaFalsi ScalarSolver;
+ const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
+ return sol;
+ }
+ }
+
+ } // namespace Equil
+} // namespace Opm
+
+
+#endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp
new file mode 100644
index 000000000..c1f72b714
--- /dev/null
+++ b/opm/core/simulator/initStateEquil.hpp
@@ -0,0 +1,673 @@
+/*
+ Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+ 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 .
+*/
+
+#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
+#define OPM_INITSTATEEQUIL_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+/**
+ * \file
+ * Facilities for an ECLIPSE-style equilibration-based
+ * initialisation scheme (keyword 'EQUIL').
+ */
+struct UnstructuredGrid;
+
+namespace Opm
+{
+
+ /**
+ * Compute initial state by an equilibration procedure.
+ *
+ * The following state fields are modified:
+ * pressure(),
+ * saturation(),
+ * surfacevol(),
+ * gasoilratio(),
+ * rv().
+ *
+ * \param[in] grid Grid.
+ * \param[in] props Property object, pvt and capillary properties are used.
+ * \param[in] deck Simulation deck, used to obtain EQUIL and related data.
+ * \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
+ */
+ void initStateEquil(const UnstructuredGrid& grid,
+ const BlackoilPropertiesInterface& props,
+ const EclipseGridParser& deck,
+ const double gravity,
+ BlackoilState& state);
+
+
+ void initStateEquil(const UnstructuredGrid& grid,
+ const BlackoilPropertiesInterface& props,
+ const Opm::DeckConstPtr newParserDeck,
+ const double gravity,
+ BlackoilState& state);
+
+
+ /**
+ * Types and routines that collectively implement a basic
+ * ECLIPSE-style equilibration-based initialisation scheme.
+ *
+ * This namespace is intentionally nested to avoid name clashes
+ * with other parts of OPM.
+ */
+ namespace Equil {
+
+ /**
+ * Compute initial phase pressures by means of equilibration.
+ *
+ * This function uses the information contained in an
+ * equilibration record (i.e., depths and pressurs) as well as
+ * a density calculator and related data to vertically
+ * integrate the phase pressure ODE
+ * \f[
+ * \frac{\mathrm{d}p_{\alpha}}{\mathrm{d}z} =
+ * \rho_{\alpha}(z,p_{\alpha})\cdot g
+ * \f]
+ * in which \f$\rho_{\alpha}$ denotes the fluid density of
+ * fluid phase \f$\alpha\f$, \f$p_{\alpha}\f$ is the
+ * corresponding phase pressure, \f$z\f$ is the depth and
+ * \f$g\f$ is the acceleration due to gravity (assumed
+ * directed downwords, in the positive \f$z\f$ direction).
+ *
+ * \tparam Region Type of an equilibration region information
+ * base. Typically an instance of the EquilReg
+ * class template.
+ *
+ * \tparam CellRange Type of cell range that demarcates the
+ * cells pertaining to the current
+ * equilibration region. Must implement
+ * methods begin() and end() to bound the range
+ * as well as provide an inner type,
+ * const_iterator, to traverse the range.
+ *
+ * \param[in] G Grid.
+ * \param[in] reg Current equilibration region.
+ * \param[in] cells Range that spans the cells of the current
+ * equilibration region.
+ * \param[in] grav Acceleration of gravity.
+ *
+ * \return Phase pressures, one vector for each active phase,
+ * of pressure values in each cell in the current
+ * equilibration region.
+ */
+ template
+ std::vector< std::vector >
+ phasePressures(const UnstructuredGrid& G,
+ const Region& reg,
+ const CellRange& cells,
+ const double grav = unit::gravity);
+
+
+
+ /**
+ * Compute initial phase saturations by means of equilibration.
+ *
+ * \tparam Region Type of an equilibration region information
+ * base. Typically an instance of the EquilReg
+ * class template.
+ *
+ * \tparam CellRange Type of cell range that demarcates the
+ * cells pertaining to the current
+ * equilibration region. Must implement
+ * methods begin() and end() to bound the range
+ * as well as provide an inner type,
+ * const_iterator, to traverse the range.
+ *
+ * \param[in] reg Current equilibration region.
+ * \param[in] cells Range that spans the cells of the current
+ * equilibration region.
+ * \param[in] props Property object, needed for capillary functions.
+ * \param[in] phase_pressures Phase pressures, one vector for each active phase,
+ * of pressure values in each cell in the current
+ * equilibration region.
+ * \return Phase saturations, one vector for each phase, each containing
+ * one saturation value per cell in the region.
+ */
+ template
+ std::vector< std::vector >
+ phaseSaturations(const Region& reg,
+ const CellRange& cells,
+ const BlackoilPropertiesInterface& props,
+ const std::vector< std::vector >& phase_pressures);
+
+
+
+ /**
+ * Compute initial Rs values.
+ *
+ * \tparam CellRangeType Type of cell range that demarcates the
+ * cells pertaining to the current
+ * equilibration region. Must implement
+ * methods begin() and end() to bound the range
+ * as well as provide an inner type,
+ * const_iterator, to traverse the range.
+ *
+ * \param[in] grid Grid.
+ * \param[in] cells Range that spans the cells of the current
+ * equilibration region.
+ * \param[in] oil_pressure Oil pressure for each cell in range.
+ * \param[in] rs_func Rs as function of pressure and depth.
+ * \return Rs values, one for each cell in the 'cells' range.
+ */
+ template
+ std::vector computeRs(const UnstructuredGrid& grid,
+ const CellRangeType& cells,
+ const std::vector oil_pressure,
+ const Miscibility::RsFunction& rs_func,
+ const std::vector gas_saturation);
+
+ namespace DeckDependent {
+
+ inline
+ std::vector
+ getEquil(const EclipseGridParser& deck)
+ {
+ if (deck.hasField("EQUIL")) {
+ const EQUIL& eql = deck.getEQUIL();
+
+ typedef std::vector::size_type sz_t;
+ const sz_t nrec = eql.equil.size();
+
+ std::vector ret;
+ ret.reserve(nrec);
+ for (sz_t r = 0; r < nrec; ++r) {
+ const EquilLine& rec = eql.equil[r];
+
+ EquilRecord record =
+ {
+ { rec.datum_depth_ ,
+ rec.datum_depth_pressure_ }
+ ,
+ { rec.water_oil_contact_depth_ ,
+ rec.oil_water_cap_pressure_ }
+ ,
+ { rec.gas_oil_contact_depth_ ,
+ rec.gas_oil_cap_pressure_ }
+ ,
+ rec.live_oil_table_index_
+ ,
+ rec.wet_gas_table_index_
+ ,
+ rec.N_
+ };
+ if (record.N != 0) {
+ OPM_THROW(std::domain_error,
+ "kw EQUIL, item 9: Only N=0 supported.");
+ }
+ ret.push_back(record);
+ }
+
+ return ret;
+ }
+ else {
+ OPM_THROW(std::domain_error,
+ "Deck does not provide equilibration data.");
+ }
+ }
+
+ inline
+ std::vector
+ getEquil(const Opm::DeckConstPtr newParserDeck)
+ {
+ if (newParserDeck->hasKeyword("EQUIL")) {
+
+ Opm::EquilWrapper eql(newParserDeck->getKeyword("EQUIL"));
+
+ const int nrec = eql.numRegions();
+
+ std::vector ret;
+ ret.reserve(nrec);
+ for (int r = 0; r < nrec; ++r) {
+
+ EquilRecord record =
+ {
+ { eql.datumDepth(r) ,
+ eql.datumDepthPressure(r) }
+ ,
+ { eql.waterOilContactDepth(r) ,
+ eql.waterOilContactCapillaryPressure(r) }
+ ,
+ { eql.gasOilContactDepth(r) ,
+ eql.gasOilContactCapillaryPressure(r) }
+ ,
+ eql.liveOilInitProceedure(r)
+ ,
+ eql.wetGasInitProceedure(r)
+ ,
+ eql.initializationTargetAccuracy(r)
+ };
+ if (record.N != 0) {
+ OPM_THROW(std::domain_error,
+ "kw EQUIL, item 9: Only N=0 supported.");
+ }
+ ret.push_back(record);
+ }
+
+ return ret;
+ }
+ else {
+ OPM_THROW(std::domain_error,
+ "Deck does not provide equilibration data.");
+ }
+ }
+
+
+ inline
+ std::vector
+ equilnum(const EclipseGridParser& deck,
+ const UnstructuredGrid& G )
+ {
+ std::vector eqlnum;
+ if (deck.hasField("EQLNUM")) {
+ const std::vector& e = deck.getIntegerValue("EQLNUM");
+ eqlnum.reserve(e.size());
+ std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
+ std::bind2nd(std::minus(), 1));
+ }
+ else {
+ // No explicit equilibration region.
+ // All cells in region zero.
+ eqlnum.assign(G.number_of_cells, 0);
+ }
+
+ return eqlnum;
+ }
+
+
+ inline
+ std::vector
+ equilnum(const Opm::DeckConstPtr newParserDeck,
+ const UnstructuredGrid& G )
+ {
+ std::vector eqlnum;
+ if (newParserDeck->hasKeyword("EQLNUM")) {
+ const std::vector& e = newParserDeck->getKeyword("EQLNUM")->getIntData();
+ eqlnum.reserve(e.size());
+ std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
+ std::bind2nd(std::minus(), 1));
+ }
+ else {
+ // No explicit equilibration region.
+ // All cells in region zero.
+ eqlnum.assign(G.number_of_cells, 0);
+ }
+
+ return eqlnum;
+ }
+
+
+ template
+ class InitialStateComputer;
+
+ template <>
+ class InitialStateComputer {
+ public:
+ InitialStateComputer(const BlackoilPropertiesInterface& props,
+ const EclipseGridParser& deck ,
+ const UnstructuredGrid& G ,
+ const double grav = unit::gravity)
+ : pp_(props.numPhases(),
+ std::vector(G.number_of_cells)),
+ sat_(props.numPhases(),
+ std::vector(G.number_of_cells)),
+ rs_(G.number_of_cells),
+ rv_(G.number_of_cells)
+ {
+ // Get the equilibration records.
+ const std::vector rec = getEquil(deck);
+
+ // Create (inverse) region mapping.
+ const RegionMapping<> eqlmap(equilnum(deck, G));
+
+ // Create Rs functions.
+ rs_func_.reserve(rec.size());
+ if (deck.hasField("DISGAS")) {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ const int cell = *(eqlmap.cells(i).begin());
+ if (rec[i].live_oil_table_index > 0) {
+ if (deck.hasField("RSVD")) {
+ // TODO When this kw is actually parsed, also check for proper number of available tables
+ // For now, just use dummy ...
+ std::vector depth; depth.push_back(0.0); depth.push_back(100.0);
+ std::vector rs; rs.push_back(0.0); rs.push_back(100.0);
+ rs_func_.push_back(std::make_shared(props, cell, depth, rs));
+ } else {
+ OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
+ }
+ } else {
+ if (rec[i].goc.depth != rec[i].main.depth) {
+ OPM_THROW(std::runtime_error,
+ "Cannot initialise: when no explicit RSVD table is given, \n"
+ "datum depth must be at the gas-oil-contact. "
+ "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
+ }
+ const double p_contact = rec[i].main.press;
+ rs_func_.push_back(std::make_shared(props, cell, p_contact));
+ }
+ }
+ } else {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ rs_func_.push_back(std::make_shared());
+ }
+ }
+
+ rv_func_.reserve(rec.size());
+ if (deck.hasField("VAPOIL")) {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ const int cell = *(eqlmap.cells(i).begin());
+ if (rec[i].wet_gas_table_index > 0) {
+ if (deck.hasField("RVVD")) {
+ // TODO When this kw is actually parsed, also check for proper number of available tables
+ // For now, just use dummy ...
+ std::vector depth; depth.push_back(0.0); depth.push_back(100.0);
+ std::vector rv; rv.push_back(0.0); rv.push_back(0.0001);
+ rv_func_.push_back(std::make_shared(props, cell, depth, rv));
+ } else {
+ OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
+ }
+ } else {
+ if (rec[i].goc.depth != rec[i].main.depth) {
+ OPM_THROW(std::runtime_error,
+ "Cannot initialise: when no explicit RVVD table is given, \n"
+ "datum depth must be at the gas-oil-contact. "
+ "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
+ }
+ const double p_contact = rec[i].main.press + rec[i].goc.press;
+ rv_func_.push_back(std::make_shared(props, cell, p_contact));
+ }
+ }
+ } else {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ rv_func_.push_back(std::make_shared());
+ }
+ }
+
+ // Compute pressures, saturations, rs and rv factors.
+ calcPressSatRsRv(eqlmap, rec, props, G, grav);
+
+ // Modify oil pressure in no-oil regions so that the pressures of present phases can
+ // be recovered from the oil pressure and capillary relations.
+ }
+
+ typedef std::vector Vec;
+ typedef std::vector PVec; // One per phase.
+
+ const PVec& press() const { return pp_; }
+ const PVec& saturation() const { return sat_; }
+ const Vec& rs() const { return rs_; }
+ const Vec& rv() const { return rv_; }
+
+ private:
+ typedef DensityCalculator RhoCalc;
+ typedef EquilReg EqReg;
+
+ std::vector< std::shared_ptr > rs_func_;
+ std::vector< std::shared_ptr > rv_func_;
+
+ PVec pp_;
+ PVec sat_;
+ Vec rs_;
+ Vec rv_;
+
+ template
+ void
+ calcPressSatRsRv(const RMap& reg ,
+ const std::vector< EquilRecord >& rec ,
+ const Opm::BlackoilPropertiesInterface& props,
+ const UnstructuredGrid& G ,
+ const double grav)
+ {
+ typedef Miscibility::NoMixing NoMix;
+
+ for (typename RMap::RegionId
+ r = 0, nr = reg.numRegions();
+ r < nr; ++r)
+ {
+ const typename RMap::CellRange cells = reg.cells(r);
+
+ const int repcell = *cells.begin();
+ const RhoCalc calc(props, repcell);
+ const EqReg eqreg(rec[r], calc,
+ rs_func_[r], rv_func_[r],
+ props.phaseUsage());
+
+ PVec press = phasePressures(G, eqreg, cells, grav);
+
+ const PVec sat = phaseSaturations(eqreg, cells, props, press);
+
+ const int np = props.numPhases();
+ for (int p = 0; p < np; ++p) {
+ copyFromRegion(press[p], cells, pp_[p]);
+ copyFromRegion(sat[p], cells, sat_[p]);
+ }
+ if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
+ && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
+ const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+ const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+ const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
+ const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
+ copyFromRegion(rs, cells, rs_);
+ copyFromRegion(rv, cells, rv_);
+ }
+ }
+ }
+
+ template
+ void copyFromRegion(const Vec& source,
+ const CellRangeType& cells,
+ Vec& destination)
+ {
+ auto s = source.begin();
+ auto c = cells.begin();
+ const auto e = cells.end();
+ for (; c != e; ++c, ++s) {
+ destination[*c] = *s;
+ }
+ }
+
+ };
+
+
+ template <>
+ class InitialStateComputer {
+ public:
+ InitialStateComputer(const BlackoilPropertiesInterface& props,
+ const Opm::DeckConstPtr newParserDeck,
+ const UnstructuredGrid& G ,
+ const double grav = unit::gravity)
+ : pp_(props.numPhases(),
+ std::vector(G.number_of_cells)),
+ sat_(props.numPhases(),
+ std::vector(G.number_of_cells)),
+ rs_(G.number_of_cells),
+ rv_(G.number_of_cells)
+ {
+ // Get the equilibration records.
+ const std::vector rec = getEquil(newParserDeck);
+
+ // Create (inverse) region mapping.
+ const RegionMapping<> eqlmap(equilnum(newParserDeck, G));
+
+ // Create Rs functions.
+ rs_func_.reserve(rec.size());
+ if (newParserDeck->hasKeyword("DISGAS")) {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ const int cell = *(eqlmap.cells(i).begin());
+ if (rec[i].live_oil_table_index > 0) {
+ if (newParserDeck->hasKeyword("RSVD") && rec[i].live_oil_table_index <= newParserDeck->getKeyword("RSVD")->size()) {
+ Opm::SimpleTable rsvd(newParserDeck->getKeyword("RSVD"),std::vector{"vd", "rs"},rec[i].live_oil_table_index-1);
+ std::vector vd(rsvd.getColumn("vd"));
+ std::vector rs(rsvd.getColumn("rs"));
+ rs_func_.push_back(std::make_shared(props, cell, vd, rs));
+ } else {
+ OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
+ }
+ } else {
+ if (rec[i].goc.depth != rec[i].main.depth) {
+ OPM_THROW(std::runtime_error,
+ "Cannot initialise: when no explicit RSVD table is given, \n"
+ "datum depth must be at the gas-oil-contact. "
+ "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
+ }
+ const double p_contact = rec[i].main.press;
+ rs_func_.push_back(std::make_shared(props, cell, p_contact));
+ }
+ }
+ } else {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ rs_func_.push_back(std::make_shared());
+ }
+ }
+
+ rv_func_.reserve(rec.size());
+ if (newParserDeck->hasKeyword("VAPOIL")) {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ const int cell = *(eqlmap.cells(i).begin());
+ if (rec[i].wet_gas_table_index > 0) {
+ if (newParserDeck->hasKeyword("RVVD") && rec[i].wet_gas_table_index <= newParserDeck->getKeyword("RVVD")->size()) {
+ Opm::SimpleTable rvvd(newParserDeck->getKeyword("RVVD"),std::vector{"vd", "rv"},rec[i].wet_gas_table_index-1);
+ std::vector vd(rvvd.getColumn("vd"));
+ std::vector rv(rvvd.getColumn("rv"));
+ rv_func_.push_back(std::make_shared(props, cell, vd, rv));
+ } else {
+ OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
+ }
+ } else {
+ if (rec[i].goc.depth != rec[i].main.depth) {
+ OPM_THROW(std::runtime_error,
+ "Cannot initialise: when no explicit RVVD table is given, \n"
+ "datum depth must be at the gas-oil-contact. "
+ "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
+ }
+ const double p_contact = rec[i].main.press + rec[i].goc.press;
+ rv_func_.push_back(std::make_shared(props, cell, p_contact));
+ }
+ }
+ } else {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ rv_func_.push_back(std::make_shared());
+ }
+ }
+
+ // Compute pressures, saturations, rs and rv factors.
+ calcPressSatRsRv(eqlmap, rec, props, G, grav);
+
+ // Modify oil pressure in no-oil regions so that the pressures of present phases can
+ // be recovered from the oil pressure and capillary relations.
+ }
+
+ typedef std::vector Vec;
+ typedef std::vector PVec; // One per phase.
+
+ const PVec& press() const { return pp_; }
+ const PVec& saturation() const { return sat_; }
+ const Vec& rs() const { return rs_; }
+ const Vec& rv() const { return rv_; }
+
+ private:
+ typedef DensityCalculator RhoCalc;
+ typedef EquilReg EqReg;
+
+ std::vector< std::shared_ptr > rs_func_;
+ std::vector< std::shared_ptr > rv_func_;
+
+ PVec pp_;
+ PVec sat_;
+ Vec rs_;
+ Vec rv_;
+
+ template
+ void
+ calcPressSatRsRv(const RMap& reg ,
+ const std::vector< EquilRecord >& rec ,
+ const Opm::BlackoilPropertiesInterface& props,
+ const UnstructuredGrid& G ,
+ const double grav)
+ {
+ typedef Miscibility::NoMixing NoMix;
+
+ for (typename RMap::RegionId
+ r = 0, nr = reg.numRegions();
+ r < nr; ++r)
+ {
+ const typename RMap::CellRange cells = reg.cells(r);
+
+ const int repcell = *cells.begin();
+ const RhoCalc calc(props, repcell);
+ const EqReg eqreg(rec[r], calc,
+ rs_func_[r], rv_func_[r],
+ props.phaseUsage());
+
+ PVec press = phasePressures(G, eqreg, cells, grav);
+
+ const PVec sat = phaseSaturations(eqreg, cells, props, press);
+
+ const int np = props.numPhases();
+ for (int p = 0; p < np; ++p) {
+ copyFromRegion(press[p], cells, pp_[p]);
+ copyFromRegion(sat[p], cells, sat_[p]);
+ }
+ if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
+ && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
+ const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+ const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+ const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
+ const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
+ copyFromRegion(rs, cells, rs_);
+ copyFromRegion(rv, cells, rv_);
+ }
+ }
+ }
+
+ template
+ void copyFromRegion(const Vec& source,
+ const CellRangeType& cells,
+ Vec& destination)
+ {
+ auto s = source.begin();
+ auto c = cells.begin();
+ const auto e = cells.end();
+ for (; c != e; ++c, ++s) {
+ destination[*c] = *s;
+ }
+ }
+
+ };
+ } // namespace DeckDependent
+ } // namespace Equil
+} // namespace Opm
+
+#include
+
+#endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED
diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp
new file mode 100644
index 000000000..d227d54d1
--- /dev/null
+++ b/opm/core/simulator/initStateEquil_impl.hpp
@@ -0,0 +1,781 @@
+/*
+ Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+ 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 .
+*/
+
+#ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
+#define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+namespace Opm
+{
+ namespace Details {
+ template
+ class RK4IVP {
+ public:
+ RK4IVP(const RHS& f ,
+ const std::array& span,
+ const double y0 ,
+ const int N )
+ : N_(N)
+ , span_(span)
+ {
+ const double h = stepsize();
+ const double h2 = h / 2;
+ const double h6 = h / 6;
+
+ y_.reserve(N + 1);
+ f_.reserve(N + 1);
+
+ y_.push_back(y0);
+ f_.push_back(f(span_[0], y0));
+
+ for (int i = 0; i < N; ++i) {
+ const double x = span_[0] + i*h;
+ const double y = y_.back();
+
+ const double k1 = f_[i];
+ const double k2 = f(x + h2, y + h2*k1);
+ const double k3 = f(x + h2, y + h2*k2);
+ const double k4 = f(x + h , y + h*k3);
+
+ y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
+ f_.push_back(f(x + h, y_.back()));
+ }
+
+ assert (y_.size() == std::vector::size_type(N + 1));
+ }
+
+ double
+ operator()(const double x) const
+ {
+ // Dense output (O(h**3)) according to Shampine
+ // (Hermite interpolation)
+ const double h = stepsize();
+ int i = (x - span_[0]) / h;
+ const double t = (x - (span_[0] + i*h)) / h;
+
+ // Crude handling of evaluation point outside "span_";
+ if (i < 0) { i = 0; }
+ if (N_ <= i) { i = N_ - 1; }
+
+ const double y0 = y_[i], y1 = y_[i + 1];
+ const double f0 = f_[i], f1 = f_[i + 1];
+
+ double u = (1 - 2*t) * (y1 - y0);
+ u += h * ((t - 1)*f0 + t*f1);
+ u *= t * (t - 1);
+ u += (1 - t)*y0 + t*y1;
+
+ return u;
+ }
+
+ private:
+ int N_;
+ std::array span_;
+ std::vector y_;
+ std::vector f_;
+
+ double
+ stepsize() const { return (span_[1] - span_[0]) / N_; }
+ };
+
+ namespace PhasePressODE {
+ template
+ class Water {
+ public:
+ Water(const Density& rho,
+ const int np,
+ const int ix,
+ const double norm_grav)
+ : rho_(rho)
+ , svol_(np, 0)
+ , ix_(ix)
+ , g_(norm_grav)
+ {
+ svol_[ix_] = 1.0;
+ }
+
+ double
+ operator()(const double /* depth */,
+ const double press) const
+ {
+ return this->density(press) * g_;
+ }
+
+ private:
+ const Density& rho_;
+ std::vector svol_;
+ const int ix_;
+ const double g_;
+
+ double
+ density(const double press) const
+ {
+ const std::vector& rho = rho_(press, svol_);
+
+ return rho[ix_];
+ }
+ };
+
+ template
+ class Oil {
+ public:
+ Oil(const Density& rho,
+ const RS& rs,
+ const int np,
+ const int oix,
+ const int gix,
+ const double norm_grav)
+ : rho_(rho)
+ , rs_(rs)
+ , svol_(np, 0)
+ , oix_(oix)
+ , gix_(gix)
+ , g_(norm_grav)
+ {
+ svol_[oix_] = 1.0;
+ }
+
+ double
+ operator()(const double depth,
+ const double press) const
+ {
+ return this->density(depth, press) * g_;
+ }
+
+ private:
+ const Density& rho_;
+ const RS& rs_;
+ mutable std::vector svol_;
+ const int oix_;
+ const int gix_;
+ const double g_;
+
+ double
+ density(const double depth,
+ const double press) const
+ {
+ if (gix_ >= 0) {
+ svol_[gix_] = rs_(depth, press);
+ }
+
+ const std::vector& rho = rho_(press, svol_);
+ return rho[oix_];
+ }
+ };
+
+ template
+ class Gas {
+ public:
+ Gas(const Density& rho,
+ const RV& rv,
+ const int np,
+ const int gix,
+ const int oix,
+ const double norm_grav)
+ : rho_(rho)
+ , rv_(rv)
+ , svol_(np, 0)
+ , gix_(gix)
+ , oix_(oix)
+ , g_(norm_grav)
+ {
+ svol_[gix_] = 1.0;
+ }
+
+ double
+ operator()(const double depth,
+ const double press) const
+ {
+ return this->density(depth, press) * g_;
+ }
+
+ private:
+ const Density& rho_;
+ const RV& rv_;
+ mutable std::vector svol_;
+ const int gix_;
+ const int oix_;
+ const double g_;
+
+ double
+ density(const double depth,
+ const double press) const
+ {
+ if (oix_ >= 0) {
+ svol_[oix_] = rv_(depth, press);
+ }
+
+ const std::vector& rho = rho_(press, svol_);
+ return rho[gix_];
+ }
+ };
+ } // namespace PhasePressODE
+
+ namespace PhaseUsed {
+ inline bool
+ water(const PhaseUsage& pu)
+ {
+ return bool(pu.phase_used[ Opm::BlackoilPhases::Aqua ]);
+ }
+
+ inline bool
+ oil(const PhaseUsage& pu)
+ {
+ return bool(pu.phase_used[ Opm::BlackoilPhases::Liquid ]);
+ }
+
+ inline bool
+ gas(const PhaseUsage& pu)
+ {
+ return bool(pu.phase_used[ Opm::BlackoilPhases::Vapour ]);
+ }
+ } // namespace PhaseUsed
+
+ namespace PhaseIndex {
+ inline int
+ water(const PhaseUsage& pu)
+ {
+ int i = -1;
+ if (PhaseUsed::water(pu)) {
+ i = pu.phase_pos[ Opm::BlackoilPhases::Aqua ];
+ }
+
+ return i;
+ }
+
+ inline int
+ oil(const PhaseUsage& pu)
+ {
+ int i = -1;
+ if (PhaseUsed::oil(pu)) {
+ i = pu.phase_pos[ Opm::BlackoilPhases::Liquid ];
+ }
+
+ return i;
+ }
+
+ inline int
+ gas(const PhaseUsage& pu)
+ {
+ int i = -1;
+ if (PhaseUsed::gas(pu)) {
+ i = pu.phase_pos[ Opm::BlackoilPhases::Vapour ];
+ }
+
+ return i;
+ }
+ } // namespace PhaseIndex
+
+ namespace PhasePressure {
+ template
+ void
+ assign(const UnstructuredGrid& G ,
+ const std::array& f ,
+ const double split,
+ const CellRange& cells,
+ std::vector& p )
+ {
+ const int nd = G.dimensions;
+
+ enum { up = 0, down = 1 };
+
+ std::vector::size_type c = 0;
+ for (typename CellRange::const_iterator
+ ci = cells.begin(), ce = cells.end();
+ ci != ce; ++ci, ++c)
+ {
+ assert (c < p.size());
+
+ const double z = G.cell_centroids[(*ci)*nd + (nd - 1)];
+ p[c] = (z < split) ? f[up](z) : f[down](z);
+ }
+ }
+
+ template
+ void
+ water(const UnstructuredGrid& G ,
+ const Region& reg ,
+ const std::array& span ,
+ const double grav ,
+ const double po_woc,
+ const CellRange& cells ,
+ std::vector