diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 6dbb6612..5f4004e7 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -169,6 +169,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_blackoilstate.cpp
tests/test_wellsmanager.cpp
tests/test_wellcontrols.cpp
+ tests/test_equil.cpp
)
# originally generated with the command:
@@ -323,6 +324,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/simulator/WellState.hpp
opm/core/simulator/initState.hpp
opm/core/simulator/initState_impl.hpp
+ opm/core/simulator/initStateEquil.hpp
+ opm/core/simulator/initStateEquil_impl.hpp
opm/core/tof/DGBasis.hpp
opm/core/tof/TofReorder.hpp
opm/core/tof/TofDiscGalReorder.hpp
diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp
new file mode 100644
index 00000000..f9b5a65a
--- /dev/null
+++ b/opm/core/simulator/initStateEquil.hpp
@@ -0,0 +1,710 @@
+/*
+ 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
+
+/**
+ * \file
+ * Facilities for an ECLIPSE-style equilibration-based
+ * initialisation scheme (keyword 'EQUIL').
+ */
+struct UnstructuredGrid;
+
+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 {
+ /**
+ * Type that implements "no phase mixing" policy.
+ */
+ struct NoMixing {
+ /**
+ * 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
+ {
+ 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:
+ /**
+ * Constructor.
+ *
+ * \param[in] depth Depth nodes.
+ * \param[in] rs Dissolved gas-oil ratio at @c depth.
+ */
+ RsVD(const std::vector& depth,
+ const std::vector& rs)
+ : depth_(depth)
+ , rs_(rs)
+ {
+ }
+
+ /**
+ * 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
+ {
+ return linearInterpolation(depth_, rs_, depth);
+ }
+
+ private:
+ std::vector depth_; /**< Depth nodes */
+ std::vector rs_; /**< Dissolved gas-oil ratio */
+ };
+ } // namespace miscibility
+
+ /**
+ * Forward and reverse mappings between cells and
+ * regions/partitions (e.g., the ECLIPSE-style 'SATNUM',
+ * 'PVTNUM', or 'EQUILNUM' arrays).
+ *
+ * \tparam Region Type of a forward region mapping. Expected
+ * to provide indexed access through
+ * operator[]() as well as inner types
+ * 'value_type', 'size_type', and
+ * 'const_iterator'.
+ */
+ template < class Region = std::vector >
+ class RegionMapping {
+ public:
+ /**
+ * Constructor.
+ *
+ * \param[in] reg Forward region mapping, restricted to
+ * active cells only.
+ */
+ explicit
+ RegionMapping(const Region& reg)
+ : reg_(reg)
+ {
+ rev_.init(reg_);
+ }
+
+ /**
+ * Type of forward (cell-to-region) mapping result.
+ * Expected to be an integer.
+ */
+ typedef typename Region::value_type RegionId;
+
+ /**
+ * Type of reverse (region-to-cell) mapping (element)
+ * result.
+ */
+ typedef typename Region::size_type CellId;
+
+ /**
+ * Type of reverse region-to-cell range bounds and
+ * iterators.
+ */
+ typedef typename std::vector::const_iterator CellIter;
+
+ /**
+ * Range of cells. Result from reverse (region-to-cell)
+ * mapping.
+ */
+ class CellRange {
+ public:
+ /**
+ * Constructor.
+ *
+ * \param[in] b Beginning of range.
+ * \param[in] e One past end of range.
+ */
+ CellRange(const CellIter b,
+ const CellIter e)
+ : b_(b), e_(e)
+ {}
+
+ /**
+ * Read-only iterator on cell ranges.
+ */
+ typedef CellIter const_iterator;
+
+ /**
+ * Beginning of cell range.
+ */
+ const_iterator begin() const { return b_; }
+
+ /**
+ * One past end of cell range.
+ */
+ const_iterator end() const { return e_; }
+
+ private:
+ const_iterator b_;
+ const_iterator e_;
+ };
+
+ /**
+ * Number of declared regions in cell-to-region mapping.
+ */
+ RegionId
+ numRegions() const { return RegionId(rev_.p.size()) - 1; }
+
+ /**
+ * Compute region number of given active cell.
+ *
+ * \param[in] c Active cell
+ * \return Region to which @c c belongs.
+ */
+ RegionId
+ region(const CellId c) const { return reg_[c]; }
+
+ /**
+ * Extract active cells in particular region.
+ *
+ * \param[in] r Region number
+ * \returns Range of active cells in region @c r.
+ */
+ CellRange
+ cells(const RegionId r) const {
+ const RegionId i = r - rev_.low;
+ return CellRange(rev_.c.begin() + rev_.p[i + 0],
+ rev_.c.begin() + rev_.p[i + 1]);
+ }
+
+ private:
+ /**
+ * Copy of forward region mapping (cell-to-region).
+ */
+ Region reg_;
+
+ /**
+ * Reverse mapping (region-to-cell).
+ */
+ struct {
+ typedef typename std::vector::size_type Pos;
+ std::vector p; /**< Region start pointers */
+ std::vector c; /**< Region cells */
+ RegionId low; /**< Smallest region number */
+
+ /**
+ * Compute reverse mapping. Standard linear insertion
+ * sort algorithm.
+ */
+ void
+ init(const Region& reg)
+ {
+ typedef typename Region::const_iterator CI;
+ const std::pair
+ m = std::minmax_element(reg.begin(), reg.end());
+
+ low = *m.first;
+
+ const typename Region::size_type
+ n = *m.second - low + 1;
+
+ p.resize(n + 1); std::fill(p.begin(), p.end(), Pos(0));
+ for (CellId i = 0, nc = reg.size(); i < nc; ++i) {
+ p[ reg[i] - low + 1 ] += 1;
+ }
+
+ for (typename std::vector::size_type
+ i = 1, sz = p.size(); i < sz; ++i) {
+ p[0] += p[i];
+ p[i] = p[0] - p[i];
+ }
+
+ assert (p[0] ==
+ static_cast(reg.size()));
+
+ c.resize(reg.size());
+ for (CellId i = 0, nc = reg.size(); i < nc; ++i) {
+ c[ p[ reg[i] - low + 1 ] ++ ] = i;
+ }
+
+ p[0] = 0;
+ }
+ } rev_; /**< Reverse mapping instance */
+ };
+
+ /**
+ * 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;
+ };
+
+ /**
+ * 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.
+ *
+ * \tparam RS Type that provides access to a calculator for
+ * (initial) dissolved gas-oil ratios as a function of depth
+ * and (oil) pressure. Must implement an operator() declared
+ * as
+ *
+ * double
+ * operator()(const double depth,
+ * const double press)
+ *
+ * that calculates the dissolved gas-oil ratio at depth @c
+ * depth and (oil) pressure @c press.
+ *
+ * \tparam RV Type that provides access to a calculator for
+ * (initial) vapourised oil-gas ratios as a function of depth
+ * and (gas) pressure. Must implement an operator() declared
+ * as
+ *
+ * double
+ * operator()(const double depth,
+ * const double press)
+ *
+ * that calculates the vapourised oil-gas ratio at depth @c
+ * depth and (gas) 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,
+ const RS& rs,
+ const RV& 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 RS CalcDissolution;
+
+ /**
+ * Type of vapourised oil-gas ratio calculator.
+ */
+ typedef RV 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 */
+ RS rs_; /**< RS calculator */
+ RV rv_; /**< RV calculator */
+ PhaseUsage pu_; /**< Active phase summary */
+ };
+
+ /**
+ * 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);
+
+ 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_ }
+ };
+
+ 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")) {
+ eqlnum = deck.getIntegerValue("EQLNUM");
+ }
+ else {
+ // No explicit equilibration region.
+ // All cells in region zero.
+ eqlnum.assign(G.number_of_cells, 0);
+ }
+
+ return eqlnum;
+ }
+
+ template
+ class PhasePressureComputer;
+
+ template <>
+ class PhasePressureComputer {
+ public:
+ PhasePressureComputer(const BlackoilPropertiesInterface& props,
+ const EclipseGridParser& deck ,
+ const UnstructuredGrid& G )
+ : pp_(props.numPhases(),
+ std::vector(G.number_of_cells))
+ {
+ const std::vector rec = getEquil(deck);
+ const RegionMapping<> eqlmap(equilnum(deck, G));
+
+ calcII(eqlmap, rec, props, G);
+ }
+
+ typedef std::vector PVal;
+ typedef std::vector PPress;
+
+ const PPress& press() const { return pp_; }
+
+ private:
+ typedef DensityCalculator RhoCalc;
+ typedef EquilReg EqReg;
+
+ PPress pp_;
+
+ template
+ void
+ calcII(const RMap& reg ,
+ const std::vector< EquilRecord >& rec ,
+ const Opm::BlackoilPropertiesInterface& props,
+ const UnstructuredGrid& G )
+ {
+ 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, NoMix(), NoMix(),
+ props.phaseUsage());
+
+ const PPress& res = phasePressures(G, eqreg, cells);
+
+ for (int p = 0, np = props.numPhases(); p < np; ++p) {
+ PVal& d = pp_[p];
+ PVal::const_iterator s = res[p].begin();
+ for (typename RMap::CellRange::const_iterator
+ c = cells.begin(),
+ e = cells.end();
+ c != e; ++c, ++s)
+ {
+ d[*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 00000000..ffc566be
--- /dev/null
+++ b/opm/core/simulator/initStateEquil_impl.hpp
@@ -0,0 +1,572 @@
+/*
+ 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 std::binary_function {
+ 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& press )
+ {
+ using PhasePressODE::Water;
+ typedef Water ODE;
+
+ const PhaseUsage& pu = reg.phaseUsage();
+
+ const int wix = PhaseIndex::water(pu);
+ ODE drho(reg.densityCalculator(), pu.num_phases, wix, grav);
+
+ const double z0 = reg.zwoc();
+ const double p0 = po_woc - reg.pcow_woc(); // Pcow = Po - Pw
+
+ std::array up = {{ z0, span[0] }};
+ std::array down = {{ z0, span[1] }};
+
+ typedef Details::RK4IVP WPress;
+ std::array wpress = {
+ {
+ WPress(drho, up , p0, 100)
+ ,
+ WPress(drho, down, p0, 100)
+ }
+ };
+
+ assign(G, wpress, z0, cells, press);
+ }
+
+ template
+ void
+ oil(const UnstructuredGrid& G ,
+ const Region& reg ,
+ const std::array& span ,
+ const double grav ,
+ const CellRange& cells ,
+ std::vector& press ,
+ double& po_woc,
+ double& po_goc)
+ {
+ using PhasePressODE::Oil;
+ typedef Oil ODE;
+
+ const PhaseUsage& pu = reg.phaseUsage();
+
+ const int oix = PhaseIndex::oil(pu);
+ const int gix = PhaseIndex::gas(pu);
+ ODE drho(reg.densityCalculator(),
+ reg.dissolutionCalculator(),
+ pu.num_phases, oix, gix, grav);
+
+ const double z0 = reg.datum();
+ const double p0 = reg.pressure();
+
+ std::array up = {{ z0, span[0] }};
+ std::array down = {{ z0, span[1] }};
+
+ typedef Details::RK4IVP OPress;
+ std::array opress = {
+ {
+ OPress(drho, up , p0, 100)
+ ,
+ OPress(drho, down, p0, 100)
+ }
+ };
+
+ assign(G, opress, z0, cells, press);
+
+ po_woc = -std::numeric_limits::max();
+ po_goc = -std::numeric_limits::max();
+
+ const double woc = reg.zwoc();
+ // Compute Oil pressure at WOC
+ if (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum
+ else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum
+ else { po_woc = p0; } // WOC *at* datum
+
+ const double goc = reg.zgoc();
+ // Compute Oil pressure at GOC
+ if (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum
+ else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum
+ else { po_goc = p0; } // GOC *at* datum
+ }
+
+ template
+ void
+ gas(const UnstructuredGrid& G ,
+ const Region& reg ,
+ const std::array& span ,
+ const double grav ,
+ const double po_goc,
+ const CellRange& cells ,
+ std::vector& press )
+ {
+ using PhasePressODE::Gas;
+ typedef Gas ODE;
+
+ const PhaseUsage& pu = reg.phaseUsage();
+
+ const int gix = PhaseIndex::gas(pu);
+ const int oix = PhaseIndex::oil(pu);
+
+ ODE drho(reg.densityCalculator(),
+ reg.evaporationCalculator(),
+ pu.num_phases, gix, oix, grav);
+
+ const double z0 = reg.zgoc();
+ const double p0 = po_goc + reg.pcgo_goc(); // Pcog = Pg - Po
+
+ std::array up = {{ z0, span[0] }};
+ std::array down = {{ z0, span[1] }};
+
+ typedef Details::RK4IVP GPress;
+ std::array gpress = {
+ {
+ GPress(drho, up , p0, 100)
+ ,
+ GPress(drho, down, p0, 100)
+ }
+ };
+
+ assign(G, gpress, z0, cells, press);
+ }
+ } // namespace PhasePressure
+
+ template
+ void
+ equilibrateOWG(const UnstructuredGrid& G,
+ const Region& reg,
+ const double grav,
+ const std::array& span,
+ const CellRange& cells,
+ std::vector< std::vector >& press)
+ {
+ const PhaseUsage& pu = reg.phaseUsage();
+
+ double po_woc = -1, po_goc = -1;
+ if (PhaseUsed::oil(pu)) {
+ const int oix = PhaseIndex::oil(pu);
+ PhasePressure::oil(G, reg, span, grav, cells,
+ press[ oix ], po_woc, po_goc);
+ }
+
+ if (PhaseUsed::water(pu)) {
+ const int wix = PhaseIndex::water(pu);
+ PhasePressure::water(G, reg, span, grav, po_woc,
+ cells, press[ wix ]);
+ }
+
+ if (PhaseUsed::gas(pu)) {
+ const int gix = PhaseIndex::gas(pu);
+ PhasePressure::gas(G, reg, span, grav, po_goc,
+ cells, press[ gix ]);
+ }
+ }
+ } // namespace Details
+
+ namespace equil {
+ template
+ std::vector< std::vector >
+ phasePressures(const UnstructuredGrid& G,
+ const Region& reg,
+ const CellRange& cells,
+ const double grav)
+ {
+ std::array span =
+ {{ std::numeric_limits::max() ,
+ -std::numeric_limits::max() }}; // Symm. about 0.
+
+ int ncell = 0;
+ {
+ // This code is only supported in three space dimensions
+ assert (G.dimensions == 3);
+
+ const int nd = G.dimensions;
+
+ // Define short-name aliases to reduce visual clutter.
+ const double* const nc = & G.node_coordinates[0];
+
+ const int* const cfp = & G.cell_facepos[0];
+ const int* const cf = & G.cell_faces[0];
+
+ const int* const fnp = & G.face_nodepos[0];
+ const int* const fn = & G.face_nodes[0];
+
+ // Define vertical span as
+ //
+ // [minimum(node depth(cells)), maximum(node depth(cells))]
+ //
+ // Note: We use a sledgehammer approach--looping all
+ // the nodes of all the faces of all the 'cells'--to
+ // compute those bounds. This necessarily entails
+ // visiting some nodes (and faces) multiple times.
+ //
+ // Note: The implementation of 'RK4IVP<>' implicitly
+ // imposes the requirement that cell centroids are all
+ // within this vertical span. That requirement is not
+ // checked.
+ for (typename CellRange::const_iterator
+ ci = cells.begin(), ce = cells.end();
+ ci != ce; ++ci, ++ncell)
+ {
+ for (const int
+ *fi = & cf[ cfp[*ci + 0] ],
+ *fe = & cf[ cfp[*ci + 1] ];
+ fi != fe; ++fi)
+ {
+ for (const int
+ *i = & fn[ fnp[*fi + 0] ],
+ *e = & fn[ fnp[*fi + 1] ];
+ i != e; ++i)
+ {
+ const double z = nc[(*i)*nd + (nd - 1)];
+
+ if (z < span[0]) { span[0] = z; }
+ if (z > span[1]) { span[1] = z; }
+ }
+ }
+ }
+ }
+
+ const int np = reg.phaseUsage().num_phases;
+
+ typedef std::vector pval;
+ std::vector press(np, pval(ncell, 0.0));
+
+ const double z0 = reg.datum();
+ const double zwoc = reg.zwoc ();
+ const double zgoc = reg.zgoc ();
+
+ if (! ((zgoc > z0) || (z0 > zwoc))) {
+ // Datum depth in oil zone (zgoc <= z0 <= zwoc)
+ Details::equilibrateOWG(G, reg, grav, span, cells, press);
+ }
+
+ return press;
+ }
+ } // namespace equil
+} // namespace Opm
+
+#endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp
new file mode 100644
index 00000000..27b46b55
--- /dev/null
+++ b/tests/test_equil.cpp
@@ -0,0 +1,315 @@
+/*
+ Copyright 2014 SINTEF ICT, Applied Mathematics.
+*/
+
+#include "config.h"
+
+/* --- Boost.Test boilerplate --- */
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE // Suppress own messages when throw()ing
+
+#define BOOST_TEST_MODULE UnitsTest
+#include
+#include
+
+/* --- our own headers --- */
+
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+BOOST_AUTO_TEST_SUITE ()
+
+BOOST_AUTO_TEST_CASE (PhasePressure)
+{
+ typedef std::vector PVal;
+ typedef std::vector PPress;
+
+ std::shared_ptr
+ G(create_grid_cart3d(10, 1, 10), destroy_grid);
+
+ Opm::parameter::ParameterGroup param;
+ {
+ using Opm::unit::kilogram;
+ using Opm::unit::meter;
+ using Opm::unit::cubic;
+
+ std::stringstream dens; dens << 700*kilogram/cubic(meter);
+ param.insertParameter("rho2", dens.str());
+ }
+
+ typedef Opm::BlackoilPropertiesBasic Props;
+ Props props(param, G->dimensions, G->number_of_cells);
+
+ typedef Opm::equil::DensityCalculator RhoCalc;
+ RhoCalc calc(props, 0);
+
+ Opm::equil::EquilRecord record =
+ {
+ { 0 , 1e5 } , // Datum depth, pressure
+ { 5 , 0 } , // Zwoc , Pcow_woc
+ { 0 , 0 } // Zgoc , Pcgo_goc
+ };
+
+ Opm::equil::EquilReg
+ region(record, calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage());
+
+ std::vector cells(G->number_of_cells);
+ std::iota(cells.begin(), cells.end(), 0);
+
+ const double grav = 10;
+ const PPress ppress = Opm::equil::phasePressures(*G, region, cells, grav);
+
+ const int first = 0, last = G->number_of_cells - 1;
+ const double reltol = 1.0e-8;
+ BOOST_CHECK_CLOSE(ppress[0][first] , 90e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[0][last ] , 180e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
+}
+
+BOOST_AUTO_TEST_CASE (CellSubset)
+{
+ typedef std::vector PVal;
+ typedef std::vector PPress;
+
+ std::shared_ptr
+ G(create_grid_cart3d(10, 1, 10), destroy_grid);
+
+ Opm::parameter::ParameterGroup param;
+ {
+ using Opm::unit::kilogram;
+ using Opm::unit::meter;
+ using Opm::unit::cubic;
+
+ std::stringstream dens; dens << 700*kilogram/cubic(meter);
+ param.insertParameter("rho2", dens.str());
+ }
+
+ typedef Opm::BlackoilPropertiesBasic Props;
+ Props props(param, G->dimensions, G->number_of_cells);
+
+ typedef Opm::equil::DensityCalculator RhoCalc;
+ RhoCalc calc(props, 0);
+
+ Opm::equil::EquilRecord record[] =
+ {
+ {
+ { 0 , 1e5 } , // Datum depth, pressure
+ { 2.5 , -0.075e5 } , // Zwoc , Pcow_woc
+ { 0 , 0 } // Zgoc , Pcgo_goc
+ }
+ ,
+ {
+ { 5 , 1.35e5 } , // Datum depth, pressure
+ { 7.5 , -0.225e5 } , // Zwoc , Pcow_woc
+ { 5 , 0 } // Zgoc , Pcgo_goc
+ }
+ };
+
+ Opm::equil::EquilReg region[] =
+ {
+ Opm::equil::EquilReg(record[0], calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage())
+ ,
+ Opm::equil::EquilReg(record[0], calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage())
+ ,
+ Opm::equil::EquilReg(record[1], calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage())
+ ,
+ Opm::equil::EquilReg(record[1], calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage())
+ };
+
+ const int cdim[] = { 2, 1, 2 };
+ int ncoarse = cdim[0];
+ for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; }
+
+ std::vector< std::vector > cells(ncoarse);
+ for (int c = 0; c < G->number_of_cells; ++c) {
+ int ci = c;
+ const int i = ci % G->cartdims[0]; ci /= G->cartdims[0];
+ const int j = ci % G->cartdims[1];
+ const int k = ci / G->cartdims[1];
+
+ const int ic = (i / (G->cartdims[0] / cdim[0]));
+ const int jc = (j / (G->cartdims[1] / cdim[1]));
+ const int kc = (k / (G->cartdims[2] / cdim[2]));
+ const int ix = ic + cdim[0]*(jc + cdim[1]*kc);
+
+ assert ((0 <= ix) && (ix < ncoarse));
+ cells[ix].push_back(c);
+ }
+
+ PPress ppress(2, PVal(G->number_of_cells, 0));
+ for (std::vector< std::vector >::const_iterator
+ r = cells.begin(), e = cells.end();
+ r != e; ++r)
+ {
+ const int rno = int(r - cells.begin());
+ const double grav = 10;
+ const PPress p =
+ Opm::equil::phasePressures(*G, region[rno], *r, grav);
+
+ PVal::size_type i = 0;
+ for (std::vector::const_iterator
+ c = r->begin(), ce = r->end();
+ c != ce; ++c, ++i)
+ {
+ assert (i < p[0].size());
+
+ ppress[0][*c] = p[0][i];
+ ppress[1][*c] = p[1][i];
+ }
+ }
+
+ const int first = 0, last = G->number_of_cells - 1;
+ const double reltol = 1.0e-8;
+ BOOST_CHECK_CLOSE(ppress[0][first] , 105e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
+}
+
+BOOST_AUTO_TEST_CASE (RegMapping)
+{
+ typedef std::vector PVal;
+ typedef std::vector PPress;
+
+ std::shared_ptr
+ G(create_grid_cart3d(10, 1, 10), destroy_grid);
+
+ Opm::parameter::ParameterGroup param;
+ {
+ using Opm::unit::kilogram;
+ using Opm::unit::meter;
+ using Opm::unit::cubic;
+
+ std::stringstream dens; dens << 700*kilogram/cubic(meter);
+ param.insertParameter("rho2", dens.str());
+ }
+
+ typedef Opm::BlackoilPropertiesBasic Props;
+ Props props(param, G->dimensions, G->number_of_cells);
+
+ typedef Opm::equil::DensityCalculator RhoCalc;
+ RhoCalc calc(props, 0);
+
+ Opm::equil::EquilRecord record[] =
+ {
+ {
+ { 0 , 1e5 } , // Datum depth, pressure
+ { 2.5 , -0.075e5 } , // Zwoc , Pcow_woc
+ { 0 , 0 } // Zgoc , Pcgo_goc
+ }
+ ,
+ {
+ { 5 , 1.35e5 } , // Datum depth, pressure
+ { 7.5 , -0.225e5 } , // Zwoc , Pcow_woc
+ { 5 , 0 } // Zgoc , Pcgo_goc
+ }
+ };
+
+ Opm::equil::EquilReg region[] =
+ {
+ Opm::equil::EquilReg(record[0], calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage())
+ ,
+ Opm::equil::EquilReg(record[0], calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage())
+ ,
+ Opm::equil::EquilReg(record[1], calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage())
+ ,
+ Opm::equil::EquilReg(record[1], calc,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage())
+ };
+
+ std::vector eqlnum(G->number_of_cells);
+ {
+ std::vector cells(G->number_of_cells);
+ std::iota(cells.begin(), cells.end(), 0);
+
+ const int cdim[] = { 2, 1, 2 };
+ int ncoarse = cdim[0];
+ for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; }
+
+ partition_unif_idx(G->dimensions, G->number_of_cells,
+ G->cartdims, cdim,
+ &cells[0], &eqlnum[0]);
+ }
+ Opm::equil::RegionMapping<> eqlmap(eqlnum);
+
+ PPress ppress(2, PVal(G->number_of_cells, 0));
+ for (int r = 0, e = eqlmap.numRegions(); r != e; ++r)
+ {
+ const Opm::equil::RegionMapping<>::CellRange&
+ rng = eqlmap.cells(r);
+
+ const int rno = r;
+ const double grav = 10;
+ const PPress p =
+ Opm::equil::phasePressures(*G, region[rno], rng, grav);
+
+ PVal::size_type i = 0;
+ for (Opm::equil::RegionMapping<>::CellRange::const_iterator
+ c = rng.begin(), ce = rng.end();
+ c != ce; ++c, ++i)
+ {
+ assert (i < p[0].size());
+
+ ppress[0][*c] = p[0][i];
+ ppress[1][*c] = p[1][i];
+ }
+ }
+
+ const int first = 0, last = G->number_of_cells - 1;
+ const double reltol = 1.0e-8;
+ BOOST_CHECK_CLOSE(ppress[0][first] , 105e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
+ BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
+}
+
+BOOST_AUTO_TEST_SUITE_END()