OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge from upstream and corrected conflictes.

Browse Source
This commit is contained in:
Halvor Møll Nilsen 2012-10-30 13:38:55 +01:00
commit ae05d3d7e8
45 changed files with 2408 additions and 582 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
opm/core/fluid/BlackoilPropertiesBasic.cpp
View File

@ -26,20 +26,20 @@ namespace Opm
{
BlackoilPropertiesBasic::BlackoilPropertiesBasic(const parameter::ParameterGroup& param,
const int dim,
const int num_cells)
const int dim,
const int num_cells)
{
double poro = param.getDefault("porosity", 1.0);
using namespace Opm::unit;
using namespace Opm::prefix;
double perm = param.getDefault("permeability", 100.0)*milli*darcy;
double poro = param.getDefault("porosity", 1.0);
using namespace Opm::unit;
using namespace Opm::prefix;
double perm = param.getDefault("permeability", 100.0)*milli*darcy;
rock_.init(dim, num_cells, poro, perm);
pvt_.init(param);
pvt_.init(param);
satprops_.init(param);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
}
BlackoilPropertiesBasic::~BlackoilPropertiesBasic()
@ -90,11 +90,11 @@ namespace Opm
/// \param[out] dmudp If non-null: array of nP viscosity derivative values,
/// array must be valid before calling.
void BlackoilPropertiesBasic::viscosity(const int n,
const double* p,
const double* z,
const int* /*cells*/,
double* mu,
double* dmudp) const
const double* p,
const double* z,
const int* /*cells*/,
double* mu,
double* dmudp) const
{
if (dmudp) {
THROW("BlackoilPropertiesBasic::viscosity() -- derivatives of viscosity not yet implemented.");
@ -114,16 +114,16 @@ namespace Opm
/// array must be valid before calling. The matrices are output
/// in Fortran order.
void BlackoilPropertiesBasic::matrix(const int n,
const double* /*p*/,
const double* /*z*/,
const int* /*cells*/,
double* A,
double* dAdp) const
const double* /*p*/,
const double* /*z*/,
const int* /*cells*/,
double* A,
double* dAdp) const
{
const int np = numPhases();
ASSERT(np <= 2);
double B[2]; // Must be enough since component classes do not handle more than 2.
pvt_.B(1, 0, 0, B);
const int np = numPhases();
ASSERT(np <= 2);
double B[2]; // Must be enough since component classes do not handle more than 2.
pvt_.B(1, 0, 0, B);
// Compute A matrix
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
@ -152,8 +152,8 @@ namespace Opm
/// of a call to the method matrix().
/// \param[out] rho Array of nP density values, array must be valid before calling.
void BlackoilPropertiesBasic::density(const int n,
const double* A,
double* rho) const
const double* A,
double* rho) const
{
const int np = numPhases();
const double* sdens = pvt_.surfaceDensities();
@ -186,10 +186,10 @@ namespace Opm
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void BlackoilPropertiesBasic::relperm(const int n,
const double* s,
const int* /*cells*/,
double* kr,
double* dkrds) const
const double* s,
const int* /*cells*/,
double* kr,
double* dkrds) const
{
satprops_.relperm(n, s, kr, dkrds);
}
@ -205,10 +205,10 @@ namespace Opm
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void BlackoilPropertiesBasic::capPress(const int n,
const double* s,
const int* /*cells*/,
double* pc,
double* dpcds) const
const double* s,
const int* /*cells*/,
double* pc,
double* dpcds) const
{
satprops_.capPress(n, s, pc, dpcds);
}
@ -226,7 +226,7 @@ namespace Opm
double* smin,
double* smax) const
{
satprops_.satRange(n, smin, smax);
satprops_.satRange(n, smin, smax);
}

20
opm/core/fluid/BlackoilPropertiesBasic.hpp
View File

@ -35,16 +35,16 @@ namespace Opm
{
public:
/// Construct from parameters.
/// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP
/// porosity (1.0) Porosity
/// permeability (100.0) Permeability in mD
/// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP
/// porosity (1.0) Porosity
/// permeability (100.0) Permeability in mD
BlackoilPropertiesBasic(const parameter::ParameterGroup& param,
const int dim,
const int num_cells);
const int dim,
const int num_cells);
/// Destructor.
virtual ~BlackoilPropertiesBasic();
@ -151,7 +151,7 @@ namespace Opm
double* dpcds) const;
/// Obtain the range of allowable saturation values.
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points.

80
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@ -19,6 +19,7 @@
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
namespace Opm
{
@ -26,15 +27,19 @@ namespace Opm
const UnstructuredGrid& grid,
bool init_rock)
{
if (init_rock){
if (init_rock){
rock_.init(deck, grid);
}
pvt_.init(deck);
satprops_.init(deck, grid);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
pvt_.init(deck, 200);
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, 200);
if (pvt_.numPhases() != satprops_->numPhases()) {
THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
@ -45,13 +50,56 @@ namespace Opm
if(init_rock){
rock_.init(deck, grid);
}
const int samples = param.getDefault("dead_tab_size", 1025);
pvt_.init(deck, samples);
satprops_.init(deck, grid, param);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
const int pvt_samples = param.getDefault("pvt_tab_size", 200);
pvt_.init(deck, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 200);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
if (sat_samples > 1) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
THROW("Unknown threephase_model: " << threephase_model);
}
} else {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
THROW("Unknown threephase_model: " << threephase_model);
}
}
if (pvt_.numPhases() != satprops_->numPhases()) {
THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
}
@ -256,7 +304,7 @@ namespace Opm
double* kr,
double* dkrds) const
{
satprops_.relperm(n, s, cells, kr, dkrds);
satprops_->relperm(n, s, cells, kr, dkrds);
}
@ -275,7 +323,7 @@ namespace Opm
double* pc,
double* dpcds) const
{
satprops_.capPress(n, s, cells, pc, dpcds);
satprops_->capPress(n, s, cells, pc, dpcds);
}
@ -291,7 +339,7 @@ namespace Opm
double* smin,
double* smax) const
{
satprops_.satRange(n, cells, smin, smax);
satprops_->satRange(n, cells, smin, smax);
}

20
opm/core/fluid/BlackoilPropertiesFromDeck.hpp
View File

@ -27,6 +27,7 @@
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <boost/scoped_ptr.hpp>
struct UnstructuredGrid;
@ -40,7 +41,7 @@ namespace Opm
public:
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
@ -48,12 +49,15 @@ namespace Opm
/// Initialize from deck, grid and parameters.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
/// \param[in] param Parameters. Accepted parameters include:
/// dead_tab_size (1025) number of uniform sample points for dead-oil pvt tables.
/// tab_size_kr (200) number of uniform sample points for saturation tables.
/// pvt_tab_size (200) number of uniform sample points for dead-oil pvt tables.
/// sat_tab_size (200) number of uniform sample points for saturation tables.
/// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2").
/// For both size parameters, a 0 or negative value indicates that no spline fitting is to
/// be done, and the input fluid data used directly for linear interpolation.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
@ -164,9 +168,9 @@ namespace Opm
double* dpcds) const;
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
@ -179,7 +183,7 @@ namespace Opm
private:
RockFromDeck rock_;
BlackoilPvtProperties pvt_;
SaturationPropsFromDeck satprops_;
boost::scoped_ptr<SaturationPropsInterface> satprops_;
mutable std::vector<double> B_;
mutable std::vector<double> dB_;
mutable std::vector<double> R_;

6
opm/core/fluid/BlackoilPropertiesInterface.hpp
View File

@ -138,9 +138,9 @@ namespace Opm
double* dpcds) const = 0;
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.

70
opm/core/fluid/IncompPropertiesBasic.cpp
View File

@ -28,22 +28,22 @@ namespace Opm
{
IncompPropertiesBasic::IncompPropertiesBasic(const parameter::ParameterGroup& param,
const int dim,
const int num_cells)
const int dim,
const int num_cells)
{
double poro = param.getDefault("porosity", 1.0);
using namespace Opm::unit;
using namespace Opm::prefix;
double perm = param.getDefault("permeability", 100.0)*milli*darcy;
double poro = param.getDefault("porosity", 1.0);
using namespace Opm::unit;
using namespace Opm::prefix;
double perm = param.getDefault("permeability", 100.0)*milli*darcy;
rock_.init(dim, num_cells, poro, perm);
pvt_.init(param);
pvt_.init(param);
satprops_.init(param);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]);
}
IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases,
@ -56,14 +56,14 @@ namespace Opm
const int num_cells)
{
rock_.init(dim, num_cells, por, perm);
pvt_.init(num_phases, rho, mu);
pvt_.init(num_phases, rho, mu);
satprops_.init(num_phases, relpermfunc);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]);
}
IncompPropertiesBasic::~IncompPropertiesBasic()
@ -109,7 +109,7 @@ namespace Opm
/// \return Array of P viscosity values.
const double* IncompPropertiesBasic::viscosity() const
{
return &viscosity_[0];
return &viscosity_[0];
}
/// \return Array of P density values.
@ -117,7 +117,7 @@ namespace Opm
{
// No difference between reservoir and surface densities
// modelled by this class.
return pvt_.surfaceDensities();
return pvt_.surfaceDensities();
}
/// \return Array of P density values.
@ -125,7 +125,7 @@ namespace Opm
{
// No difference between reservoir and surface densities
// modelled by this class.
return pvt_.surfaceDensities();
return pvt_.surfaceDensities();
}
/// \param[in] n Number of data points.
@ -138,10 +138,10 @@ namespace Opm
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesBasic::relperm(const int n,
const double* s,
const int* /*cells*/,
double* kr,
double* dkrds) const
const double* s,
const int* /*cells*/,
double* kr,
double* dkrds) const
{
satprops_.relperm(n, s, kr, dkrds);
}
@ -157,10 +157,10 @@ namespace Opm
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesBasic::capPress(const int n,
const double* s,
const int* /*cells*/,
double* pc,
double* dpcds) const
const double* s,
const int* /*cells*/,
double* pc,
double* dpcds) const
{
satprops_.capPress(n, s, pc, dpcds);
}
@ -174,11 +174,11 @@ namespace Opm
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void IncompPropertiesBasic::satRange(const int n,
const int* /*cells*/,
double* smin,
double* smax) const
const int* /*cells*/,
double* smin,
double* smax) const
{
satprops_.satRange(n, smin, smax);
satprops_.satRange(n, smin, smax);
}
} // namespace Opm

34
opm/core/fluid/IncompPropertiesBasic.hpp
View File

@ -42,29 +42,29 @@ namespace Opm
{
public:
/// Construct from parameters.
/// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP
/// porosity (1.0) Porosity
/// permeability (100.0) Permeability in mD
/// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP
/// porosity (1.0) Porosity
/// permeability (100.0) Permeability in mD
IncompPropertiesBasic(const parameter::ParameterGroup& param,
const int dim,
const int num_cells);
const int dim,
const int num_cells);
/// Construct from arguments a basic two phase fluid.
IncompPropertiesBasic(const int num_phases,
const SaturationPropsBasic::RelPermFunc& relpermfunc,
const std::vector<double>& rho,
const std::vector<double>& mu,
const std::vector<double>& mu,
const double porosity,
const double permeability,
const int dim,
const int num_cells);
const int num_cells);
/// Destructor.
/// Destructor.
virtual ~IncompPropertiesBasic();
// ---- Rock interface ----
@ -132,9 +132,9 @@ namespace Opm
double* dpcds) const;
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
@ -145,9 +145,9 @@ namespace Opm
double* smax) const;
private:
RockBasic rock_;
PvtPropertiesBasic pvt_;
PvtPropertiesBasic pvt_;
SaturationPropsBasic satprops_;
std::vector<double> viscosity_;
std::vector<double> viscosity_;
};

44
opm/core/fluid/IncompPropertiesFromDeck.cpp
View File

@ -27,15 +27,15 @@ namespace Opm
{
IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
const UnstructuredGrid& grid)
{
rock_.init(deck, grid);
pvt_.init(deck);
satprops_.init(deck, grid);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
pvt_.init(deck);
satprops_.init(deck, grid, 200);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
}
IncompPropertiesFromDeck::~IncompPropertiesFromDeck()
@ -81,19 +81,19 @@ namespace Opm
/// \return Array of P viscosity values.
const double* IncompPropertiesFromDeck::viscosity() const
{
return pvt_.viscosity();
return pvt_.viscosity();
}
/// \return Array of P density values.
const double* IncompPropertiesFromDeck::density() const
{
return pvt_.reservoirDensities();
return pvt_.reservoirDensities();
}
/// \return Array of P density values.
const double* IncompPropertiesFromDeck::surfaceDensity() const
{
return pvt_.surfaceDensities();
return pvt_.surfaceDensities();
}
/// \param[in] n Number of data points.
@ -106,10 +106,10 @@ namespace Opm
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesFromDeck::relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const
const double* s,
const int* cells,
double* kr,
double* dkrds) const
{
satprops_.relperm(n, s, cells, kr, dkrds);
}
@ -125,10 +125,10 @@ namespace Opm
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesFromDeck::capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const
const double* s,
const int* cells,
double* pc,
double* dpcds) const
{
satprops_.capPress(n, s, cells, pc, dpcds);
}
@ -142,11 +142,11 @@ namespace Opm
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void IncompPropertiesFromDeck::satRange(const int n,
const int* cells,
double* smin,
double* smax) const
const int* cells,
double* smin,
double* smax) const
{
satprops_.satRange(n, cells, smin, smax);
satprops_.satRange(n, cells, smin, smax);
}
} // namespace Opm

16
opm/core/fluid/IncompPropertiesFromDeck.hpp
View File

@ -47,13 +47,13 @@ namespace Opm
public:
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \param grid Grid to which property object applies, needed for the
/// \param grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
IncompPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
const UnstructuredGrid& grid);
/// Destructor.
/// Destructor.
virtual ~IncompPropertiesFromDeck();
// ---- Rock interface ----
@ -121,9 +121,9 @@ namespace Opm
double* dpcds) const;
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
@ -134,8 +134,8 @@ namespace Opm
double* smax) const;
private:
RockFromDeck rock_;
PvtPropertiesIncompFromDeck pvt_;
SaturationPropsFromDeck satprops_;
PvtPropertiesIncompFromDeck pvt_;
SaturationPropsFromDeck<SatFuncStone2Uniform> satprops_;
};

8
opm/core/fluid/IncompPropertiesInterface.hpp
View File

@ -62,7 +62,7 @@ namespace Opm
/// \return Array of P viscosity values.
virtual const double* viscosity() const = 0;
/// Densities of fluid phases at surface conditions.
/// Densities of fluid phases at reservoir conditions.
/// \return Array of P density values.
virtual const double* density() const = 0;
@ -109,9 +109,9 @@ namespace Opm
double* pc,
double* dpcds) const = 0;
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.

112
opm/core/fluid/PvtPropertiesBasic.cpp
View File

@ -34,41 +34,41 @@ namespace Opm
void PvtPropertiesBasic::init(const parameter::ParameterGroup& param)
{
int num_phases = param.getDefault("num_phases", 2);
if (num_phases > 3 || num_phases < 1) {
THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases);
}
density_.resize(num_phases);
viscosity_.resize(num_phases);
// We currently do not allow the user to set B.
formation_volume_factor_.clear();
formation_volume_factor_.resize(num_phases, 1.0);
int num_phases = param.getDefault("num_phases", 2);
if (num_phases > 3 || num_phases < 1) {
THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases);
}
density_.resize(num_phases);
viscosity_.resize(num_phases);
// We currently do not allow the user to set B.
formation_volume_factor_.clear();
formation_volume_factor_.resize(num_phases, 1.0);
// Setting mu and rho from parameters
using namespace Opm::prefix;
using namespace Opm::unit;
const double kgpm3 = kilogram/cubic(meter);
const double cP = centi*Poise;
std::string rname[3] = { "rho1", "rho2", "rho3" };
double rdefault[3] = { 1.0e3, 1.0e3, 1.0e3 };
std::string vname[3] = { "mu1", "mu2", "mu3" };
double vdefault[3] = { 1.0, 1.0, 1.0 };
for (int phase = 0; phase < num_phases; ++phase) {
density_[phase] = kgpm3*param.getDefault(rname[phase], rdefault[phase]);
viscosity_[phase] = cP*param.getDefault(vname[phase], vdefault[phase]);
}
// Setting mu and rho from parameters
using namespace Opm::prefix;
using namespace Opm::unit;
const double kgpm3 = kilogram/cubic(meter);
const double cP = centi*Poise;
std::string rname[3] = { "rho1", "rho2", "rho3" };
double rdefault[3] = { 1.0e3, 1.0e3, 1.0e3 };
std::string vname[3] = { "mu1", "mu2", "mu3" };
double vdefault[3] = { 1.0, 1.0, 1.0 };
for (int phase = 0; phase < num_phases; ++phase) {
density_[phase] = kgpm3*param.getDefault(rname[phase], rdefault[phase]);
viscosity_[phase] = cP*param.getDefault(vname[phase], vdefault[phase]);
}
}
void PvtPropertiesBasic::init(const int num_phases,
const std::vector<double>& rho,
const std::vector<double>& visc)
{
if (num_phases > 3 || num_phases < 1) {
THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases);
}
// We currently do not allow the user to set B.
formation_volume_factor_.clear();
formation_volume_factor_.resize(num_phases, 1.0);
if (num_phases > 3 || num_phases < 1) {
THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases);
}
// We currently do not allow the user to set B.
formation_volume_factor_.clear();
formation_volume_factor_.resize(num_phases, 1.0);
density_ = rho;
viscosity_ = visc;
}
@ -87,69 +87,69 @@ namespace Opm
void PvtPropertiesBasic::mu(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_mu) const
const double* /*p*/,
const double* /*z*/,
double* output_mu) const
{
const int np = numPhases();
const int np = numPhases();
for (int phase = 0; phase < np; ++phase) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[np*i + phase] = viscosity_[phase];
}
}
}
}
void PvtPropertiesBasic::B(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_B) const
const double* /*p*/,
const double* /*z*/,
double* output_B) const
{
const int np = numPhases();
const int np = numPhases();
for (int phase = 0; phase < np; ++phase) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[np*i + phase] = formation_volume_factor_[phase];
}
}
}
}
void PvtPropertiesBasic::dBdp(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
const double* /*p*/,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
const int np = numPhases();
const int np = numPhases();
for (int phase = 0; phase < np; ++phase) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[np*i + phase] = formation_volume_factor_[phase];
output_dBdp[np*i + phase] = 0.0;
}
}
}
}
void PvtPropertiesBasic::R(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
const int np = numPhases();
std::fill(output_R, output_R + n*np, 0.0);
const int np = numPhases();
std::fill(output_R, output_R + n*np, 0.0);
}
void PvtPropertiesBasic::dRdp(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
const int np = numPhases();
std::fill(output_R, output_R + n*np, 0.0);
std::fill(output_dRdp, output_dRdp + n*np, 0.0);
const int np = numPhases();
std::fill(output_R, output_R + n*np, 0.0);
std::fill(output_dRdp, output_dRdp + n*np, 0.0);
}
} // namespace Opm

18
opm/core/fluid/PvtPropertiesBasic.hpp
View File

@ -38,11 +38,11 @@ namespace Opm
PvtPropertiesBasic();
/// Initialize from parameters.
/// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1, 2 or 3.
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP
void init(const parameter::ParameterGroup& param);
/// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1, 2 or 3.
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP
void init(const parameter::ParameterGroup& param);
/// Initialize from arguments.
/// Basic multi phase fluid pvt properties.
@ -55,7 +55,7 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of size numPhases().
const double* surfaceDensities() const;
const double* surfaceDensities() const;
/// Viscosity as a function of p and z.
void mu(const int n,
@ -90,9 +90,9 @@ namespace Opm
double* output_dRdp) const;
private:
std::vector<double> density_;
std::vector<double> viscosity_;
std::vector<double> formation_volume_factor_;
std::vector<double> density_;
std::vector<double> viscosity_;
std::vector<double> formation_volume_factor_;
};
}

70
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
View File

@ -38,54 +38,54 @@ namespace Opm
{
typedef std::vector<std::vector<std::vector<double> > > table_t;
// If we need multiple regions, this class and the SinglePvt* classes must change.
int region_number = 0;
int region_number = 0;
PhaseUsage phase_usage = phaseUsageFromDeck(deck);
if (phase_usage.phase_used[PhaseUsage::Vapour] ||
!phase_usage.phase_used[PhaseUsage::Aqua] ||
!phase_usage.phase_used[PhaseUsage::Liquid]) {
THROW("PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n");
}
if (phase_usage.phase_used[PhaseUsage::Vapour] ||
!phase_usage.phase_used[PhaseUsage::Aqua] ||
!phase_usage.phase_used[PhaseUsage::Liquid]) {
THROW("PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n");
}
// Surface densities. Accounting for different orders in eclipse and our code.
if (deck.hasField("DENSITY")) {
const std::vector<double>& d = deck.getDENSITY().densities_[region_number];
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water];
surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
} else {
THROW("Input is missing DENSITY\n");
}
// Surface densities. Accounting for different orders in eclipse and our code.
if (deck.hasField("DENSITY")) {
const std::vector<double>& d = deck.getDENSITY().densities_[region_number];
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water];
surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
} else {
THROW("Input is missing DENSITY\n");
}
// Make reservoir densities the same as surface densities initially.
// We will modify them with formation volume factors if found.
reservoir_density_ = surface_density_;
// Water viscosity.
if (deck.hasField("PVTW")) {
const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number];
if (pvtw[2] != 0.0 || pvtw[4] != 0.0) {
MESSAGE("Compressibility effects in PVTW are ignored.");
}
if (deck.hasField("PVTW")) {
const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number];
if (pvtw[2] != 0.0 || pvtw[4] != 0.0) {
MESSAGE("Compressibility effects in PVTW are ignored.");
}
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3];
} else {
// Eclipse 100 default.
// viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise;
THROW("Input is missing PVTW\n");
}
viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3];
} else {
// Eclipse 100 default.
// viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise;
THROW("Input is missing PVTW\n");
}
// Oil viscosity.
if (deck.hasField("PVCDO")) {
const std::vector<double>& pvcdo = deck.getPVCDO().pvcdo_[region_number];
if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) {
MESSAGE("Compressibility effects in PVCDO are ignored.");
}
if (deck.hasField("PVCDO")) {
const std::vector<double>& pvcdo = deck.getPVCDO().pvcdo_[region_number];
if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) {
MESSAGE("Compressibility effects in PVCDO are ignored.");
}
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3];
} else {
THROW("Input is missing PVCDO\n");
}
viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3];
} else {
THROW("Input is missing PVCDO\n");
}
}
const double* PvtPropertiesIncompFromDeck::surfaceDensities() const

12
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp
View File

@ -39,14 +39,14 @@ namespace Opm
PvtPropertiesIncompFromDeck();
/// Initialize from deck.
void init(const EclipseGridParser& deck);
void init(const EclipseGridParser& deck);
/// Number of active phases.
int numPhases() const;
/// Densities of stock components at surface conditions.
/// \return Array of size numPhases().
const double* surfaceDensities() const;
const double* surfaceDensities() const;
/// Densities of stock components at reservoir conditions.
/// Note: a reasonable question to ask is why there can be
@ -58,15 +58,15 @@ namespace Opm
/// reporting and using data given in terms of surface values,
/// we need to handle this difference.
/// \return Array of size numPhases().
const double* reservoirDensities() const;
const double* reservoirDensities() const;
/// Viscosities.
const double* viscosity() const;
private:
std::tr1::array<double, 2> surface_density_;
std::tr1::array<double, 2> reservoir_density_;
std::tr1::array<double, 2> viscosity_;
std::tr1::array<double, 2> surface_density_;
std::tr1::array<double, 2> reservoir_density_;
std::tr1::array<double, 2> viscosity_;
};
}

8
opm/core/fluid/RockBasic.hpp
View File

@ -35,9 +35,9 @@ namespace Opm
/// Initialize with homogenous porosity and permeability.
void init(const int dimensions,
const int num_cells,
const double poro,
const double perm);
const int num_cells,
const double poro,
const double perm);
/// \return D, the number of spatial dimensions.
int numDimensions() const
@ -66,7 +66,7 @@ namespace Opm
}
private:
int dimensions_;
int dimensions_;
std::vector<double> porosity_;
std::vector<double> permeability_;
};

6
opm/core/fluid/RockCompressibility.cpp
View File

@ -69,8 +69,8 @@ namespace Opm
const double cpnorm = rock_comp_*(pressure - pref_);
return (1.0 + cpnorm + 0.5*cpnorm*cpnorm);
} else {
// return Opm::linearInterpolation(p_, poromult_, pressure);
return Opm::linearInterpolationExtrap(p_, poromult_, pressure);
// return Opm::linearInterpolation(p_, poromult_, pressure);
return Opm::linearInterpolationExtrap(p_, poromult_, pressure);
}
}
@ -81,7 +81,7 @@ namespace Opm
} else {
//const double poromult = Opm::linearInterpolation(p_, poromult_, pressure);
//const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure);
const double poromult = Opm::linearInterpolationExtrap(p_, poromult_, pressure);
const double poromult = Opm::linearInterpolationExtrap(p_, poromult_, pressure);
const double dporomultdp = Opm::linearInterpolationDerivativeExtrap(p_, poromult_, pressure);
return dporomultdp/poromult;

2
opm/core/fluid/RockFromDeck.cpp
View File

@ -51,7 +51,7 @@ namespace Opm
/// Initialize from deck and cell mapping.
/// \param deck Deck input parser
/// \param grid grid to which property object applies, needed for the
/// \param grid grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
void RockFromDeck::init(const EclipseGridParser& deck,

2
opm/core/fluid/RockFromDeck.hpp
View File

@ -37,7 +37,7 @@ namespace Opm
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \param grid Grid to which property object applies, needed for the
/// \param grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
void init(const EclipseGridParser& deck,

204
opm/core/fluid/SaturationPropsBasic.cpp
View File

@ -29,64 +29,64 @@ namespace Opm
namespace {
struct KrFunConstant
{
double kr(double)
{
return 1.0;
}
double dkrds(double)
{
return 0.0;
}
};
struct KrFunConstant
{
double kr(double)
{
return 1.0;
}
double dkrds(double)
{
return 0.0;
}
};
struct KrFunLinear
{
double kr(double s)
{
return s;
}
double dkrds(double)
{
return 1.0;
}
};
struct KrFunLinear
{
double kr(double s)
{
return s;
}
double dkrds(double)
{
return 1.0;
}
};
struct KrFunQuadratic
{
double kr(double s)
{
return s*s;
}
double dkrds(double s)
{
return 2.0*s;
}
};
struct KrFunQuadratic
{
double kr(double s)
{
return s*s;
}
double dkrds(double s)
{
return 2.0*s;
}
};
template <class Fun>
static inline void evalAllKrDeriv(const int n, const int np,
const double* s, double* kr, double* dkrds, Fun fun)
{
if (dkrds == 0) {
template <class Fun>
static inline void evalAllKrDeriv(const int n, const int np,
const double* s, double* kr, double* dkrds, Fun fun)
{
if (dkrds == 0) {
// #pragma omp parallel for
for (int i = 0; i < n*np; ++i) {
kr[i] = fun.kr(s[i]);
}
return;
}
for (int i = 0; i < n*np; ++i) {
kr[i] = fun.kr(s[i]);
}
return;
}
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0);
for (int phase = 0; phase < np; ++phase) {
kr[i*np + phase] = fun.kr(s[i*np + phase]);
// Only diagonal elements in derivative.
dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]);
}
}
}
for (int i = 0; i < n; ++i) {
std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0);
for (int phase = 0; phase < np; ++phase) {
kr[i*np + phase] = fun.kr(s[i*np + phase]);
// Only diagonal elements in derivative.
dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]);
}
}
}
} // anon namespace
@ -109,25 +109,25 @@ namespace Opm
/// Initialize from parameters.
void SaturationPropsBasic::init(const parameter::ParameterGroup& param)
{
int num_phases = param.getDefault("num_phases", 2);
if (num_phases > 2 || num_phases < 1) {
THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases);
}
int num_phases = param.getDefault("num_phases", 2);
if (num_phases > 2 || num_phases < 1) {
THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases);
}
num_phases_ = num_phases;
//std::string rpf = param.getDefault("relperm_func", std::string("Unset"));
std::string rpf = param.getDefault("relperm_func", std::string("Linear"));
if (rpf == "Constant") {
relperm_func_ = Constant;
if(num_phases!=1){
THROW("Constant relperm with more than one phase???");
}
} else if (rpf == "Linear") {
relperm_func_ = Linear;
} else if (rpf == "Quadratic") {
relperm_func_ = Quadratic;
} else {
THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf);
}
//std::string rpf = param.getDefault("relperm_func", std::string("Unset"));
std::string rpf = param.getDefault("relperm_func", std::string("Linear"));
if (rpf == "Constant") {
relperm_func_ = Constant;
if(num_phases!=1){
THROW("Constant relperm with more than one phase???");
}
} else if (rpf == "Linear") {
relperm_func_ = Linear;
} else if (rpf == "Quadratic") {
relperm_func_ = Quadratic;
} else {
THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf);
}
}
@ -136,7 +136,7 @@ namespace Opm
/// \return P, the number of phases.
int SaturationPropsBasic::numPhases() const
{
return num_phases_;
return num_phases_;
}
@ -152,29 +152,29 @@ namespace Opm
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void SaturationPropsBasic::relperm(const int n,
const double* s,
double* kr,
double* dkrds) const
const double* s,
double* kr,
double* dkrds) const
{
switch (relperm_func_) {
case Constant:
{
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant());
break;
}
case Linear:
{
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear());
break;
}
case Quadratic:
{
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic());
break;
}
default:
THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_);
}
switch (relperm_func_) {
case Constant:
{
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant());
break;
}
case Linear:
{
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear());
break;
}
case Quadratic:
{
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic());
break;
}
default:
THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_);
}
}
@ -190,13 +190,13 @@ namespace Opm
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void SaturationPropsBasic::capPress(const int n,
const double* /*s*/,
double* pc,
double* dpcds) const
const double* /*s*/,
double* pc,
double* dpcds) const
{
std::fill(pc, pc + num_phases_*n, 0.0);
std::fill(pc, pc + num_phases_*n, 0.0);
if (dpcds) {
std::fill(dpcds, dpcds + num_phases_*num_phases_*n, 0.0);
std::fill(dpcds, dpcds + num_phases_*num_phases_*n, 0.0);
}
}
@ -207,11 +207,11 @@ namespace Opm
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void SaturationPropsBasic::satRange(const int n,
double* smin,
double* smax) const
double* smin,
double* smax) const
{
std::fill(smin, smin + num_phases_*n, 0.0);
std::fill(smax, smax + num_phases_*n, 1.0);
std::fill(smin, smin + num_phases_*n, 0.0);
std::fill(smax, smax + num_phases_*n, 1.0);
}

22
opm/core/fluid/SaturationPropsBasic.hpp
View File

@ -40,16 +40,16 @@ namespace Opm
SaturationPropsBasic();
/// Initialize from parameters.
/// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
/// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
void init(const parameter::ParameterGroup& param);
enum RelPermFunc { Constant, Linear, Quadratic };
enum RelPermFunc { Constant, Linear, Quadratic };
/// Initialize from arguments a basic Saturation property.
void init(const int num_phases,
const RelPermFunc& relperm_func)
const RelPermFunc& relperm_func)
{
num_phases_ = num_phases;
relperm_func_ = relperm_func;
@ -86,18 +86,18 @@ namespace Opm
double* pc,
double* dpcds) const;
/// Obtain the range of allowable saturation values.
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void satRange(const int n,
double* smin,
double* smax) const;
void satRange(const int n,
double* smin,
double* smax) const;
private:
int num_phases_;
RelPermFunc relperm_func_;
int num_phases_;
RelPermFunc relperm_func_;
};

42
opm/core/fluid/SaturationPropsFromDeck.hpp
View File

@ -19,12 +19,14 @@
#ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
#define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
#include <opm/core/fluid/SaturationPropsInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SatFuncStone2.hpp>
#include <opm/core/fluid/SatFuncSimple.hpp>
#include <opm/core/fluid/SatFuncGwseg.hpp>
#include <vector>
struct UnstructuredGrid;
@ -32,23 +34,31 @@ struct UnstructuredGrid;
namespace Opm
{
class SaturationPropsFromDeck : public BlackoilPhases
/// Interface to saturation functions from deck.
/// Possible values for template argument (for now):
/// SatFuncSetStone2Nonuniform,
/// SatFuncSetStone2Uniform.
/// SatFuncSetSimpleNonuniform,
/// SatFuncSetSimpleUniform.
template <class SatFuncSet>
class SaturationPropsFromDeck : public SaturationPropsInterface
{
public:
/// Default constructor.
SaturationPropsFromDeck();
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \param grid Grid to which property object applies, needed for the
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
void init(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
void init(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param);
const int samples);
/// \return P, the number of phases.
int numPhases() const;
@ -83,22 +93,23 @@ namespace Opm
double* pc,
double* dpcds) const;
/// Obtain the range of allowable saturation values.
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void satRange(const int n,
void satRange(const int n,
const int* cells,
double* smin,
double* smax) const;
double* smin,
double* smax) const;
private:
PhaseUsage phase_usage_;
typedef SatFuncSimple satfunc_t;
std::vector<satfunc_t> satfuncset_;
std::vector<SatFuncSet> satfuncset_;
std::vector<int> cell_to_func_; // = SATNUM - 1
const satfunc_t& funcForCell(const int cell) const;
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
};
@ -106,6 +117,7 @@ namespace Opm
} // namespace Opm
#include <opm/core/fluid/SaturationPropsFromDeck_impl.hpp>
#endif // OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED

221
opm/core/fluid/SaturationPropsFromDeck_impl.hpp Normal file
View File

@ -0,0 +1,221 @@
/*
Copyright 2012 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
#define OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/grid.h>
namespace Opm
{
// ----------- Methods of SaturationPropsFromDeck ---------
/// Default constructor.
template <class SatFuncSet>
SaturationPropsFromDeck<SatFuncSet>::SaturationPropsFromDeck()
{
}
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::init(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const int samples)
{
phase_usage_ = phaseUsageFromDeck(deck);
// Extract input data.
// Oil phase should be active.
if (!phase_usage_.phase_used[Liquid]) {
THROW("SaturationPropsFromDeck::init() -- oil phase must be active.");
}
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
const int num_cells = grid.number_of_cells;
cell_to_func_.resize(num_cells);
const int* gc = grid.global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell];
cell_to_func_[cell] = satnum[deck_pos] - 1;
}
}
// Find number of tables, check for consistency.
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
num_tables = swof_table.size();
if (num_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
}
if (phase_usage_.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
int num_sgof_tables = sgof_table.size();
if (num_sgof_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
}
if (num_tables == Uninitialized) {
num_tables = num_sgof_tables;
} else if (num_tables != num_sgof_tables) {
THROW("Inconsistent number of tables in SWOF and SGOF.");
}
}
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_, samples);
}
}
/// \return P, the number of phases.
template <class SatFuncSet>
int SaturationPropsFromDeck<SatFuncSet>::numPhases() const
{
return phase_usage_.num_phases;
}
/// Relative permeability.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] kr Array of nP relperm values, array must be valid before calling.
/// \param[out] dkrds If non-null: array of nP^2 relperm derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
if (dkrds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
}
}
}
/// Capillary pressure.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] pc Array of nP capillary pressure values, array must be valid before calling.
/// \param[out] dpcds If non-null: array of nP^2 derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
if (dpcds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i);
}
}
}
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::satRange(const int n,
const int* cells,
double* smin,
double* smax) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
for (int i = 0; i < n; ++i) {
for (int p = 0; p < np; ++p) {
smin[np*i + p] = funcForCell(cells[i]).smin_[p];
smax[np*i + p] = funcForCell(cells[i]).smax_[p];
}
}
}
// Map the cell number to the correct function set.
template <class SatFuncSet>
const typename SaturationPropsFromDeck<SatFuncSet>::Funcs&
SaturationPropsFromDeck<SatFuncSet>::funcForCell(const int cell) const
{
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
} // namespace Opm
#endif // OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED

85
opm/core/fluid/SaturationPropsInterface.hpp Normal file
View File

@ -0,0 +1,85 @@
/*
Copyright 2012 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED
#define OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
namespace Opm
{
class SaturationPropsInterface : public BlackoilPhases
{
public:
/// Virtual destructor.
virtual ~SaturationPropsInterface() {};
/// \return P, the number of phases.
virtual int numPhases() const = 0;
/// Relative permeability.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[out] kr Array of nP relperm values, array must be valid before calling.
/// \param[out] dkrds If non-null: array of nP^2 relperm derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
virtual void relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const = 0;
/// Capillary pressure.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[out] pc Array of nP capillary pressure values, array must be valid before calling.
/// \param[out] dpcds If non-null: array of nP^2 derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
virtual void capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const = 0;
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
virtual void satRange(const int n,
const int* cells,
double* smin,
double* smax) const = 0;
};
} // namespace Opm
#endif // OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED

14
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp
View File

@ -47,7 +47,13 @@ namespace Opm
BlackoilPvtProperties();
/// Initialize from deck.
void init(const EclipseGridParser& deck, const int samples = 16);
/// \param deck An input deck.
/// \param samples If greater than zero, indicates the number of
/// uniform samples to be taken from monotone spline
/// curves interpolating the fluid data.
/// Otherwise, interpolate linearly in the original
/// data without fitting a spline.
void init(const EclipseGridParser& deck, const int samples);
/// Number of active phases.
int numPhases() const;
@ -64,7 +70,7 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of size numPhases().
const double* surfaceDensities() const;
const double* surfaceDensities() const;
/// Viscosity as a function of p and z.
void mu(const int n,
@ -105,11 +111,11 @@ namespace Opm
PhaseUsage phase_usage_;
int region_number_;
int region_number_;
std::vector<std::tr1::shared_ptr<SinglePvtInterface> > props_;
double densities_[MaxNumPhases];
double densities_[MaxNumPhases];
mutable std::vector<double> data1_;
mutable std::vector<double> data2_;
};

4
opm/core/linalg/LinearSolverIstl.cpp
View File

@ -26,6 +26,10 @@
#include <opm/core/utility/have_boost_redef.hpp>
// Silence compatibility warning from DUNE headers since we don't use
// the deprecated member anyway (in this compilation unit)
#define DUNE_COMMON_FIELDVECTOR_SIZE_IS_METHOD 1
// TODO: clean up includes.
#include <dune/common/deprecated.hh>
#include <dune/istl/bvector.hh>

13
opm/core/newwells.h
View File

@ -235,6 +235,19 @@ void
destroy_wells(struct Wells *W);
/**
* Create a deep-copy (i.e., clone) of an existing Wells object, including its
* controls.
*
* @param[in] W Existing Wells object.
* @return Complete clone of the input object. Dispose of resources using
* function destroy_wells() when no longer needed. Returns @c NULL in case of
* allocation failure.
*/
struct Wells *
clone_wells(const struct Wells *W);
#ifdef __cplusplus
}
#endif

81
opm/core/pressure/CompressibleTpfa.cpp
View File

@ -36,6 +36,7 @@
#include <cmath>
#include <iostream>
#include <iomanip>
#include <numeric>
namespace Opm
{
@ -122,7 +123,7 @@ namespace Opm
WellState& well_state)
{
const int nc = grid_.number_of_cells;
const int nw = wells_->number_of_wells;
const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
// Set up dynamic data.
computePerSolveDynamicData(dt, state, well_state);
@ -446,30 +447,54 @@ namespace Opm
/// Compute per-iteration dynamic properties for wells.
void CompressibleTpfa::computeWellDynamicData(const double /*dt*/,
const BlackoilState& /*state*/,
const WellState& /*well_state*/)
const WellState& well_state)
{
// These are the variables that get computed by this function:
//
// std::vector<double> wellperf_A_;
// std::vector<double> wellperf_phasemob_;
const int np = props_.numPhases();
const int nw = wells_->number_of_wells;
const int nperf = wells_->well_connpos[nw];
const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
const int nperf = (wells_ != 0) ? wells_->well_connpos[nw] : 0;
wellperf_A_.resize(nperf*np*np);
wellperf_phasemob_.resize(nperf*np);
// The A matrix is set equal to the perforation grid cells'
// matrix, for both injectors and producers.
// matrix for producers, computed from bhp and injection
// component fractions from
// The mobilities are set equal to the perforation grid cells'
// mobilities, for both injectors and producers.
// mobilities for producers.
std::vector<double> mu(np);
for (int w = 0; w < nw; ++w) {
bool producer = (wells_->type[w] == PRODUCER);
const double* comp_frac = &wells_->comp_frac[np*w];
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
const int c = wells_->well_cells[j];
const double* cA = &cell_A_[np*np*c];
double* wpA = &wellperf_A_[np*np*j];
std::copy(cA, cA + np*np, wpA);
const double* cM = &cell_phasemob_[np*c];
double* wpM = &wellperf_phasemob_[np*j];
std::copy(cM, cM + np, wpM);
if (producer) {
const double* cA = &cell_A_[np*np*c];
std::copy(cA, cA + np*np, wpA);
const double* cM = &cell_phasemob_[np*c];
std::copy(cM, cM + np, wpM);
} else {
const double bhp = well_state.bhp()[w];
double perf_p = bhp;
for (int phase = 0; phase < np; ++phase) {
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
}
// Hack warning: comp_frac is used as a component
// surface-volume variable in calls to matrix() and
// viscosity(), but as a saturation in the call to
// relperm(). This is probably ok as long as injectors
// only inject pure fluids.
props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL);
props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL);
ASSERT(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6);
props_.relperm (1, comp_frac, &c, wpM , NULL);
for (int phase = 0; phase < np; ++phase) {
wpM[phase] /= mu[phase];
}
}
}
}
}
@ -487,9 +512,9 @@ namespace Opm
const double* z = &state.surfacevol()[0];
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = &wellperf_gpot_[0];
completion_data.A = &wellperf_A_[0];
completion_data.phasemob = &wellperf_phasemob_[0];
completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0;
completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
cfs_tpfa_res_wells wells_tmp;
wells_tmp.W = const_cast<Wells*>(wells_);
wells_tmp.data = &completion_data;
@ -574,9 +599,9 @@ namespace Opm
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = const_cast<double*>(&wellperf_gpot_[0]);
completion_data.A = const_cast<double*>(&wellperf_A_[0]);
completion_data.phasemob = const_cast<double*>(&wellperf_phasemob_[0]);
completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast<double*>(&wellperf_gpot_[0]) : 0;
completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
cfs_tpfa_res_wells wells_tmp;
wells_tmp.W = const_cast<Wells*>(wells_);
wells_tmp.data = &completion_data;
@ -584,6 +609,9 @@ namespace Opm
forces.wells = &wells_tmp;
forces.src = NULL;
double* wpress = ! well_state.bhp ().empty() ? & well_state.bhp ()[0] : 0;
double* wflux = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0;
cfs_tpfa_res_flux(gg,
&forces,
props_.numPhases(),
@ -592,9 +620,9 @@ namespace Opm
&face_phasemob_[0],
&face_gravcap_[0],
&state.pressure()[0],
&well_state.bhp()[0],
wpress,
&state.faceflux()[0],
&well_state.perfRates()[0]);
wflux);
cfs_tpfa_res_fpress(gg,
props_.numPhases(),
&htrans_[0],
@ -604,6 +632,23 @@ namespace Opm
&state.pressure()[0],
&state.faceflux()[0],
&state.facepressure()[0]);
// Compute well perforation pressures (not done by the C code).
if (wells_ != 0) {
const int nw = wells_->number_of_wells;
const int np = props_.numPhases();
for (int w = 0; w < nw; ++w) {
const double* comp_frac = &wells_->comp_frac[np*w];
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
const double bhp = well_state.bhp()[w];
double perf_p = bhp;
for (int phase = 0; phase < np; ++phase) {
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
}
well_state.perfPress()[j] = perf_p;
}
}
}
}
} // namespace Opm

115
opm/core/pressure/tpfa/cfs_tpfa_residual.c
View File

@ -142,7 +142,7 @@ impl_allocate(struct UnstructuredGrid *G ,
nnu = G->number_of_cells;
nwperf = 0;
if (wells != NULL) { assert (wells->W != NULL);
if ((wells != NULL) && (wells->W != NULL)) {
nnu += wells->W->number_of_wells;
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
}
@ -185,13 +185,15 @@ construct_matrix(struct UnstructuredGrid *G ,
/* ---------------------------------------------------------------------- */
{
int f, c1, c2, w, i, nc, nnu;
int wells_present;
size_t nnz;
struct CSRMatrix *A;
nc = nnu = G->number_of_cells;
if (wells != NULL) {
assert (wells->W != NULL);
wells_present = (wells != NULL) && (wells->W != NULL);
if (wells_present) {
nnu += wells->W->number_of_wells;
}
@ -214,7 +216,7 @@ construct_matrix(struct UnstructuredGrid *G ,
}
}
if (wells != NULL) {
if (wells_present) {
/* Well <-> cell connections */
struct Wells *W = wells->W;
@ -252,7 +254,7 @@ construct_matrix(struct UnstructuredGrid *G ,
}
}
if (wells != NULL) {
if (wells_present) {
/* Fill well <-> cell connections */
struct Wells *W = wells->W;
@ -741,6 +743,29 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt,
}
/* ---------------------------------------------------------------------- */
static void
welleq_coeff_shut(int np, struct cfs_tpfa_res_data *h,
double *res, double *w2c, double *w2w)
/* ---------------------------------------------------------------------- */
{
int p;
double fwi;
const double *dpflux_w;
/* Sum reservoir phase flux derivatives set by
* compute_darcyflux_and_deriv(). */
dpflux_w = h->pimpl->flux_work + (1 * np);
for (p = 0, fwi = 0.0; p < np; p++) {
fwi += dpflux_w[ p ];
}
*res = 0.0;
*w2c = 0.0;
*w2w = fwi;
}
/* ---------------------------------------------------------------------- */
static void
welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h,
@ -795,7 +820,34 @@ welleq_coeff_resv(int np, struct cfs_tpfa_res_data *h,
/* ---------------------------------------------------------------------- */
static void
assemble_completion_to_well(int w, int c, int nc, int np,
welleq_coeff_surfrate(int i, int np, struct cfs_tpfa_res_data *h,
struct WellControls *ctrl,
double *res, double *w2c, double *w2w)
/* ---------------------------------------------------------------------- */
{
int p;
double f;
const double *pflux, *dpflux_w, *dpflux_c, *distr;
pflux = h->pimpl->compflux_p + (i * (1 * np));
dpflux_w = h->pimpl->compflux_deriv_p + (i * (2 * np));
dpflux_c = dpflux_w + (1 * (1 * np));
distr = ctrl->distr + (ctrl->current * (1 * np));
*res = *w2c = *w2w = 0.0;
for (p = 0; p < np; p++) {
f = distr[ p ];
*res += f * pflux [ p ];
*w2w += f * dpflux_w[ p ];
*w2c += f * dpflux_c[ p ];
}
}
/* ---------------------------------------------------------------------- */
static void
assemble_completion_to_well(int i, int w, int c, int nc, int np,
double pw, double dt,
struct cfs_tpfa_res_wells *wells,
struct cfs_tpfa_res_data *h )
@ -815,23 +867,25 @@ assemble_completion_to_well(int w, int c, int nc, int np,
ctrl = W->ctrls[ w ];
if (ctrl->current < 0) {
assert (0); /* Shut wells currently not supported */
/* Interpreting a negative current control index to mean a shut well */
welleq_coeff_shut(np, h, &res, &w2c, &w2w);
}
else {
switch (ctrl->type[ ctrl->current ]) {
case BHP :
welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
h, &res, &w2c, &w2w);
break;
switch (ctrl->type[ ctrl->current ]) {
case BHP :
welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
h, &res, &w2c, &w2w);
break;
case RESERVOIR_RATE:
assert (W->number_of_phases == np);
welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w);
break;
case RESERVOIR_RATE:
assert (W->number_of_phases == np);
welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w);
break;
case SURFACE_RATE:
assert (0); /* Surface rate currently not supported */
break;
case SURFACE_RATE:
welleq_coeff_surfrate(i, np, h, ctrl, &res, &w2c, &w2w);
break;
}
}
/* Assemble completion contributions */
@ -854,7 +908,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
struct cfs_tpfa_res_data *h )
{
int w, i, c, np, np2, nc;
int is_neumann;
int is_neumann, is_open;
double pw, dp;
double *WI, *gpot, *pmobp;
const double *Ac, *dAc;
@ -876,6 +930,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
for (w = i = 0; w < W->number_of_wells; w++) {
pw = wpress[ w ];
is_open = W->ctrls[w]->current >= 0;
for (; i < W->well_connpos[w + 1];
i++, gpot += np, pmobp += np) {
@ -888,14 +943,16 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
init_completion_contrib(i, np, Ac, dAc, h->pimpl);
assemble_completion_to_cell(c, nc + w, np, dt, h);
if (is_open) {
assemble_completion_to_cell(c, nc + w, np, dt, h);
}
/* Prepare for RESV controls */
compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
h->pimpl->flux_work,
h->pimpl->flux_work + np);
assemble_completion_to_well(w, c, nc, np, pw, dt, wells, h);
assemble_completion_to_well(i, w, c, nc, np, pw, dt, wells, h);
}
ctrl = W->ctrls[ w ];
@ -1127,8 +1184,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
nf = G->number_of_faces;
nwperf = 0;
if (wells != NULL) {
assert (wells->W != NULL);
if ((wells != NULL) && (wells->W != NULL)) {
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
}
@ -1194,7 +1250,9 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
assemble_cell_contrib(G, c, h);
}
if ((forces != NULL) && (forces->wells != NULL)) {
if ((forces != NULL) &&
(forces->wells != NULL) &&
(forces->wells->W != NULL)) {
compute_well_compflux_and_deriv(forces->wells, cq->nphases,
cpress, wpress, h->pimpl);
@ -1297,8 +1355,9 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
{
compute_flux(G, np, trans, pmobf, gravcap_f, cpress, fflux);
if ((forces != NULL) && (forces->wells != NULL) &&
(wpress != NULL) && (wflux != NULL)) {
if ((forces != NULL) &&
(forces->wells != NULL) &&
(forces->wells->W != NULL) && (wpress != NULL) && (wflux != NULL)) {
compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux);
}

21
opm/core/simulator/WellState.hpp
View File

@ -37,12 +37,20 @@ namespace Opm
if (wells) {
const int nw = wells->number_of_wells;
bhp_.resize(nw);
// Initialize bhp to be pressure in first perforation cell.
// Initialize bhp to be target pressure
// if bhp-controlled well, otherwise set
// to pressure in first perforation cell.
for (int w = 0; w < nw; ++w) {
const int cell = wells->well_cells[wells->well_connpos[w]];
bhp_[w] = state.pressure()[cell];
const WellControls* ctrl = wells->ctrls[w];
if (ctrl->type[ctrl->current] == BHP) {
bhp_[w] = ctrl->target[ctrl->current];
} else {
const int cell = wells->well_cells[wells->well_connpos[w]];
bhp_[w] = state.pressure()[cell];
}
}
perfrates_.resize(wells->well_connpos[nw]);
perfrates_.resize(wells->well_connpos[nw], 0.0);
perfpress_.resize(wells->well_connpos[nw], -1e100);
}
}
@ -54,9 +62,14 @@ namespace Opm
std::vector<double>& perfRates() { return perfrates_; }
const std::vector<double>& perfRates() const { return perfrates_; }
/// One pressure per well connection.
std::vector<double>& perfPress() { return perfpress_; }
const std::vector<double>& perfPress() const { return perfpress_; }
private:
std::vector<double> bhp_;
std::vector<double> perfrates_;
std::vector<double> perfpress_;
};
} // namespace Opm

41
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
View File

@ -152,8 +152,8 @@ namespace Opm
B_cell = 1.0/tm.A_[np*np*cell + 0];
double src_flux = -tm.source_[cell];
bool src_is_inflow = src_flux < 0.0;
influx = src_is_inflow ? B_cell*src_flux : 0.0;
outflux = !src_is_inflow ? B_cell*src_flux : 0.0;
influx = src_is_inflow ? B_cell* src_flux : 0.0;
outflux = !src_is_inflow ? src_flux : 0.0;
comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell];
dtpv = tm.dt_/tm.porevolume0_[cell];
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
@ -349,7 +349,7 @@ namespace Opm
gf[1] = gravflux[pos];
}
s0 = tm.saturation_[cell];
dtpv = tm.dt_/tm.porevolume0_[cell];
dtpv = tm.dt_/tm.porevolume_[cell];
}
double operator()(double s) const
@ -380,8 +380,8 @@ namespace Opm
{
double sat[2] = { s, 1.0 - s };
props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[0];
mob[1] /= visc_[1];
mob[0] /= visc_[2*cell + 0];
mob[1] /= visc_[2*cell + 1];
}
@ -407,7 +407,7 @@ namespace Opm
{
// Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f)
// + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ]
// But b_w,i * rho_w,S = rho_w,i, which we conmpute with a call to props_.density().
// But b_w,i * rho_w,S = rho_w,i, which we compute with a call to props_.density().
// We assume that we already have stored T_ij in trans_.
// We also assume that the A_ matrices are updated from an earlier call to solve().
const int nc = grid_.number_of_cells;
@ -437,16 +437,15 @@ namespace Opm
void TransportModelCompressibleTwophase::solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
const double* gravflux)
const int pos,
const double* gravflux)
{
const int cell = cells[pos];
GravityResidual res(*this, cells, pos, gravflux);
if (std::fabs(res(saturation_[cell])) > tol_) {
int iters_used;
saturation_[cell] = RootFinder::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used);
}
saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]);
mobility(saturation_[cell], cell, &mob_[2*cell]);
}
@ -506,8 +505,6 @@ namespace Opm
void TransportModelCompressibleTwophase::solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure,
const double* porevolume0,
const double dt,
std::vector<double>& saturation,
std::vector<double>& surfacevol)
@ -522,18 +519,22 @@ namespace Opm
cells[c] = c;
}
mob_.resize(2*nc);
std::vector<double> boths;
Opm::toBothSat(saturation, boths);
props_.relperm(cells.size(), &boths[0], &cells[0], &mob_[0], 0);
props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
for (int c = 0; c < nc; ++c) {
mob_[2*c + 0] /= visc_[2*c + 0];
mob_[2*c + 1] /= visc_[2*c + 1];
// props_.relperm(cells.size(), &saturation[0], &cells[0], &mob_[0], 0);
// props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
// for (int c = 0; c < nc; ++c) {
// mob_[2*c + 0] /= visc_[2*c + 0];
// mob_[2*c + 1] /= visc_[2*c + 1];
// }
const int np = props_.numPhases();
for (int cell = 0; cell < nc; ++cell) {
mobility(saturation_[cell], cell, &mob_[np*cell]);
}
// Set up other variables.
porevolume0_ = porevolume0;
dt_ = dt;
toWaterSat(saturation, saturation_);

3
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
View File

@ -74,13 +74,10 @@ namespace Opm
/// vertical stack, that do not interact with other columns (for
/// gravity segregation.
/// \param[in] columns Vector of cell-columns.
/// \param[in] porevolume0 Array of pore volumes at start of timestep.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
/// \param[in, out] surfacevol Surface volume densities for each phase.
void solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure,
const double* porevolume0,
const double dt,
std::vector<double>& saturation,
std::vector<double>& surfacevol);

5
opm/core/transport/reorder/TransportModelInterface.cpp
View File

@ -20,6 +20,7 @@
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/grid.h>
#include <opm/core/utility/StopWatch.hpp>
#include <vector>
#include <cassert>
@ -31,7 +32,11 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g
std::vector<int> sequence(grid.number_of_cells);
std::vector<int> components(grid.number_of_cells + 1);
int ncomponents;
time::StopWatch clock;
clock.start();
compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
clock.stop();
std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl;
// Invoke appropriate solve method for each interdependent component.
for (int comp = 0; comp < ncomponents; ++comp) {

122
opm/core/transport/reorder/TransportModelTracerTof.cpp Normal file
View File

@ -0,0 +1,122 @@
/*
Copyright 2012 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 <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <numeric>
#include <cmath>
namespace Opm
{
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid)
: grid_(grid)
{
}
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values.
void TransportModelTracerTof::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
#ifndef NDEBUG
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
THROW("Sources do not sum to zero: " << cum_src);
}
#endif
tof.resize(grid_.number_of_cells);
std::fill(tof.begin(), tof.end(), 0.0);
tof_ = &tof[0];
reorderAndTransport(grid_, darcyflux);
}
void TransportModelTracerTof::solveSingleCell(const int cell)
{
// Compute flux terms.
// Sources have zero tof, and therefore do not contribute
// to upwind_term. Sinks on the other hand, must be added
// to the downwind_flux (note sign change resulting from
// different sign conventions: pos. source is injection,
// pos. flux is outflow).
double upwind_term = 0.0;
double downwind_flux = std::max(-source_[cell], 0.0);
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
int f = grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == grid_.face_cells[2*f]) {
flux = darcyflux_[f];
other = grid_.face_cells[2*f+1];
} else {
flux =-darcyflux_[f];
other = grid_.face_cells[2*f];
}
// Add flux to upwind_term or downwind_flux, if interior.
if (other != -1) {
if (flux < 0.0) {
upwind_term += flux*tof_[other];
} else {
downwind_flux += flux;
}
}
}
// Compute tof.
tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux;
}
void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
}
}
} // namespace Opm

74
opm/core/transport/reorder/TransportModelTracerTof.hpp Normal file
View File

@ -0,0 +1,74 @@
/*
Copyright 2012 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <vector>
#include <map>
#include <ostream>
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
/// Implements a first-order finite volume solver for
/// (single-phase) time-of-flight using reordering.
/// The equation solved is:
/// v \cdot \grad\tau = \phi
/// where v is the fluid velocity, \tau is time-of-flight and
/// \phi is the porosity. This is a boundary value problem, where
/// \tau is specified to be zero on all inflow boundaries.
class TransportModelTracerTof : public TransportModelInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTof(const UnstructuredGrid& grid);
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values.
void solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof);
private:
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
private:
const UnstructuredGrid& grid_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
double* tof_;
};
} // namespace Opm
#endif // OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED

657
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp Normal file
View File

@ -0,0 +1,657 @@
/*
Copyright 2012 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 <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/VelocityInterpolation.hpp>
#include <opm/core/linalg/blas_lapack.h>
#include <algorithm>
#include <cmath>
#include <numeric>
namespace Opm
{
// --------------- Helpers for TransportModelTracerTofDiscGal ---------------
/// A class providing discontinuous Galerkin basis functions.
struct DGBasis
{
static int numBasisFunc(const int dimensions,
const int degree)
{
switch (dimensions) {
case 1:
return degree + 1;
case 2:
return (degree + 2)*(degree + 1)/2;
case 3:
return (degree + 3)*(degree + 2)*(degree + 1)/6;
default:
THROW("Dimensions must be 1, 2 or 3.");
}
}
/// Evaluate all nonzero basis functions at x,
/// writing to f_x. The array f_x must have
/// size numBasisFunc(grid.dimensions, degree).
///
/// The basis functions are the following
/// Degree 0: 1.
/// Degree 1: x - xc, y - yc, z - zc etc.
/// Further degrees await development.
static void eval(const UnstructuredGrid& grid,
const int cell,
const int degree,
const double* x,
double* f_x)
{
const int dim = grid.dimensions;
const double* cc = grid.cell_centroids + dim*cell;
// Note intentional fallthrough in this switch statement!
switch (degree) {
case 1:
for (int ix = 0; ix < dim; ++ix) {
f_x[1 + ix] = x[ix] - cc[ix];
}
case 0:
f_x[0] = 1;
break;
default:
THROW("Maximum degree is 1 for now.");
}
}
/// Evaluate gradients of all nonzero basis functions at x,
/// writing to grad_f_x. The array grad_f_x must have size
/// numBasisFunc(grid.dimensions, degree) * grid.dimensions.
/// The <grid.dimensions> components of the first basis function
/// gradient come before the components of the second etc.
static void evalGrad(const UnstructuredGrid& grid,
const int /*cell*/,
const int degree,
const double* /*x*/,
double* grad_f_x)
{
const int dim = grid.dimensions;
const int num_basis = numBasisFunc(dim, degree);
std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
if (degree > 1) {
THROW("Maximum degree is 1 for now.");
} else if (degree == 1) {
for (int ix = 0; ix < dim; ++ix) {
grad_f_x[dim*(ix + 1) + ix] = 1.0;
}
}
}
};
static void cross(const double* a, const double* b, double* res)
{
res[0] = a[1]*b[2] - a[2]*b[1];
res[1] = a[2]*b[0] - a[0]*b[2];
res[2] = a[0]*b[1] - a[1]*b[0];
}
static double triangleArea3d(const double* p0,
const double* p1,
const double* p2)
{
double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
double cr[3];
cross(a, b, cr);
return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]);
}
/// Calculates the determinant of a 3 x 3 matrix, represented as
/// three three-dimensional arrays.
static double determinantOf(const double* a0,
const double* a1,
const double* a2)
{
return
a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) -
a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) +
a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]);
}
/// Computes the volume of a tetrahedron consisting of 4 vertices
/// with 3-dimensional coordinates
static double tetVolume(const double* p0,
const double* p1,
const double* p2,
const double* p3)
{
double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] };
return std::fabs(determinantOf(a, b, c) / 6.0);
}
/// A class providing numerical quadrature for cells.
/// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by cell volume,
/// so weights always sum to cell volume.
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = cell volume, x_0 = cell centroid
/// Degree 2 method:
/// Based on subdivision of each cell face into triangles
/// with the face centroid as a common vertex, and then
/// subdividing the cell into tetrahedra with the cell
/// centroid as a common vertex. Then apply the tetrahedron
/// rule with the following 4 nodes (uniform weights):
/// a = 0.138196601125010515179541316563436
/// x_i has all barycentric coordinates = a, except for
/// the i'th coordinate which is = 1 - 3a.
/// This rule is from http://nines.cs.kuleuven.be/ecf,
/// it is the second degree 2 4-point rule for tets,
/// referenced to Stroud(1971).
/// The tetrahedra are numbered T_{i,j}, and are given by the
/// cell centroid C, the face centroid FC_i, and two nodes
/// of face i: FN_{i,j}, FN_{i,j+1}.
class CellQuadrature
{
public:
CellQuadrature(const UnstructuredGrid& grid,
const int cell,
const int degree)
: grid_(grid), cell_(cell), degree_(degree)
{
if (degree > 2) {
THROW("CellQuadrature exact for polynomial degrees > 1 not implemented.");
}
if (degree == 2) {
// Prepare subdivision.
}
}
int numQuadPts() const
{
if (degree_ < 2) {
return 1;
}
// Degree 2 case.
int sumnodes = 0;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
}
return 4*sumnodes;
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
if (degree_ < 2) {
std::copy(cc, cc + dim, coord);
return;
}
// Degree 2 case.
int tetindex = index / 4;
const int subindex = index % 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
if (nfn <= tetindex) {
// Our tet is not associated with this face.
tetindex -= nfn;
continue;
}
const double* fc = grid_.face_centroids + dim*face;
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
const int node0 = fnodes[tetindex];
const int node1 = fnodes[(tetindex + 1) % nfn];
const double* n0c = nc + dim*node0;
const double* n1c = nc + dim*node1;
const double a = 0.138196601125010515179541316563436;
// Barycentric coordinates of our point in the tet.
double baryc[4] = { a, a, a, a };
baryc[subindex] = 1.0 - 3.0*a;
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd];
}
return;
}
THROW("Should never reach this point.");
}
double quadPtWeight(const int index) const
{
if (degree_ < 2) {
return grid_.cell_volumes[cell_];
}
// Degree 2 case.
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
int tetindex = index / 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
if (nfn <= tetindex) {
// Our tet is not associated with this face.
tetindex -= nfn;
continue;
}
const double* fc = grid_.face_centroids + dim*face;
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
const int node0 = fnodes[tetindex];
const int node1 = fnodes[(tetindex + 1) % nfn];
const double* n0c = nc + dim*node0;
const double* n1c = nc + dim*node1;
return 0.25*tetVolume(cc, fc, n0c, n1c);
}
THROW("Should never reach this point.");
}
private:
const UnstructuredGrid& grid_;
const int cell_;
const int degree_;
};
/// A class providing numerical quadrature for faces.
/// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by face area,
/// so weights always sum to face area.
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = face area, x_0 = face centroid
/// Degree 2 method:
/// Based on subdivision of the face into triangles,
/// with the centroid as a common vertex, and the triangle
/// edge midpoint rule.
/// Triangle i consists of the centroid C, nodes N_i and N_{i+1}.
/// Its area is A_i.
/// n = 2 * nn (nn = num nodes in face)
/// For i = 0..(nn-1):
/// w_i = 1/3 A_i.
/// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i
/// x_i = (N_i + N_{i+1})/2
/// x_{nn+i} = (C + N_i)/2
/// All N and A indices are interpreted cyclic, modulus nn.
class FaceQuadrature
{
public:
FaceQuadrature(const UnstructuredGrid& grid,
const int face,
const int degree)
: grid_(grid), face_(face), degree_(degree)
{
if (grid_.dimensions != 3) {
THROW("FaceQuadrature only implemented for 3D case.");
}
if (degree_ > 2) {
THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented.");
}
}
int numQuadPts() const
{
if (degree_ < 2) {
return 1;
}
// Degree 2 case.
return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
if (degree_ < 2) {
std::copy(fc, fc + dim, coord);
return;
}
// Degree 2 case.
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
const int node1 = fnodes[(index + 1)%nn];
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]);
}
} else {
// Interiour edge midpoint.
// Recall that index is now in [nn, 2*nn).
const int node = fnodes[index - nn];
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]);
}
}
}
double quadPtWeight(const int index) const
{
if (degree_ < 2) {
return grid_.face_areas[face_];
}
// Degree 2 case.
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
const int node1 = fnodes[(index + 1)%nn];
const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
return area / 3.0;
} else {
// Interiour edge midpoint.
// Recall that index is now in [nn, 2*nn).
const int node0 = fnodes[(index - 1) % nn];
const int node1 = fnodes[index - nn];
const int node2 = fnodes[(index + 1) % nn];
const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc);
return (area0 + area1) / 3.0;
}
}
private:
const UnstructuredGrid& grid_;
const int face_;
const int degree_;
};
// --------------- Methods of TransportModelTracerTofDiscGal ---------------
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi)
: grid_(grid),
coord_(grid.dimensions),
velocity_(grid.dimensions)
{
if (use_cvi) {
velocity_interpolation_.reset(new VelocityInterpolationECVI(grid));
} else {
velocity_interpolation_.reset(new VelocityInterpolationConstant(grid));
}
}
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] degree Polynomial degree of DG basis functions used.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
const int degree,
std::vector<double>& tof_coeff)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
#ifndef NDEBUG
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
THROW("Sources do not sum to zero: " << cum_src);
}
#endif
degree_ = degree;
const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_);
tof_coeff.resize(num_basis*grid_.number_of_cells);
std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0);
tof_coeff_ = &tof_coeff[0];
rhs_.resize(num_basis);
jac_.resize(num_basis*num_basis);
orig_jac_.resize(num_basis*num_basis);
basis_.resize(num_basis);
basis_nb_.resize(num_basis);
grad_basis_.resize(num_basis*grid_.dimensions);
velocity_interpolation_->setupFluxes(darcyflux);
reorderAndTransport(grid_, darcyflux);
}
void TransportModelTracerTofDiscGal::solveSingleCell(const int cell)
{
// Residual:
// For each cell K, basis function b_j (spanning V_h),
// writing the solution u_h|K = \sum_i c_i b_i
// Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx
// + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds
// - \int_K \phi b_j
// This is linear in c_i, so we do not need any nonlinear iterations.
// We assemble the jacobian and the right-hand side. The residual is
// equal to Res = Jac*c - rhs, and we compute rhs directly.
const int dim = grid_.dimensions;
const int num_basis = DGBasis::numBasisFunc(dim, degree_);
std::fill(rhs_.begin(), rhs_.end(), 0.0);
std::fill(jac_.begin(), jac_.end(), 0.0);
// Compute cell residual contribution.
// Note: Assumes that \int_K b_j = 0 for all j > 0
rhs_[0] += porevolume_[cell];
// Compute upstream residual contribution.
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
int upstream_cell = -1;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
upstream_cell = grid_.face_cells[2*face+1];
} else {
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
if (upstream_cell < 0) {
// This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
continue;
}
if (flux >= 0.0) {
// This is an outflow boundary.
continue;
}
// Do quadrature over the face to compute
// \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds
// (where u_h^{ext} is the upstream unknown (tof)).
const double normal_velocity = flux / grid_.face_areas[face];
FaceQuadrature quad(grid_, face, degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]);
const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(),
tof_coeff_ + num_basis*upstream_cell, 0.0);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j];
}
}
}
// Compute cell jacobian contribution. We use Fortran ordering
// for jac_, i.e. rows cycling fastest.
{
CellQuadrature quad(grid_, cell, 2*degree_ - 1);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// b_i (v \cdot \grad b_j)
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]);
velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
for (int dd = 0; dd < dim; ++dd) {
jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd];
}
}
}
}
}
// Compute downstream jacobian contribution from faces.
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux <= 0.0) {
// This is an inflow boundary.
continue;
}
// Do quadrature over the face to compute
// \int_{\partial K} b_i (v(x) \cdot n) b_j ds
const double normal_velocity = flux / grid_.face_areas[face];
FaceQuadrature quad(grid_, face, 2*degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// u^ext flux B (B = {b_j})
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
jac_[j*num_basis + i] += w * basis_[i] * normal_velocity * basis_[j];
}
}
}
}
// Compute downstream jacobian contribution from sink terms.
// Contribution from inflow sources would be
// similar to the contribution from upstream faces, but
// it is zero since we let all external inflow be associated
// with a zero tof.
if (source_[cell] < 0.0) {
// A sink.
const double flux = -source_[cell]; // Sign convention for flux: outflux > 0.
const double flux_density = flux / grid_.cell_volumes[cell];
// Do quadrature over the cell to compute
// \int_{K} b_i flux b_j dx
CellQuadrature quad(grid_, cell, 2*degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j];
}
}
}
}
// Solve linear equation.
MAT_SIZE_T n = num_basis;
MAT_SIZE_T nrhs = 1;
MAT_SIZE_T lda = num_basis;
std::vector<MAT_SIZE_T> piv(num_basis);
MAT_SIZE_T ldb = num_basis;
MAT_SIZE_T info = 0;
orig_jac_ = jac_;
dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info);
if (info != 0) {
// Print the local matrix and rhs.
std::cerr << "Failed solving single-cell system Ax = b in cell " << cell
<< " with A = \n";
for (int row = 0; row < n; ++row) {
for (int col = 0; col < n; ++col) {
std::cerr << " " << orig_jac_[row + n*col];
}
std::cerr << '\n';
}
std::cerr << "and b = \n";
for (int row = 0; row < n; ++row) {
std::cerr << " " << rhs_[row] << '\n';
}
THROW("Lapack error: " << info << " encountered in cell " << cell);
}
// The solution ends up in rhs_, so we must copy it.
std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell);
}
void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
}
}
} // namespace Opm

105
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp Normal file
View File

@ -0,0 +1,105 @@
/*
Copyright 2012 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <boost/shared_ptr.hpp>
#include <vector>
#include <map>
#include <ostream>
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
class VelocityInterpolationInterface;
/// Implements a discontinuous Galerkin solver for
/// (single-phase) time-of-flight using reordering.
/// The equation solved is:
/// v \cdot \grad\tau = \phi
/// where v is the fluid velocity, \tau is time-of-flight and
/// \phi is the porosity. This is a boundary value problem, where
/// \tau is specified to be zero on all inflow boundaries.
/// The user may specify the polynomial degree of the basis function space
/// used, but only degrees 0 and 1 are supported so far.
class TransportModelTracerTofDiscGal : public TransportModelInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi);
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] degree Polynomial degree of DG basis functions used.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
void solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
const int degree,
std::vector<double>& tof_coeff);
private:
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
private:
// Disable copying and assignment.
TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&);
TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&);
const UnstructuredGrid& grid_;
boost::shared_ptr<VelocityInterpolationInterface> velocity_interpolation_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
int degree_;
double* tof_coeff_;
std::vector<double> rhs_; // single-cell right-hand-side
std::vector<double> jac_; // single-cell jacobian
std::vector<double> orig_jac_; // single-cell jacobian (copy)
// Below: storage for quantities needed by solveSingleCell().
std::vector<double> coord_;
std::vector<double> basis_;
std::vector<double> basis_nb_;
std::vector<double> grad_basis_;
std::vector<double> velocity_;
};
} // namespace Opm
#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED

66
opm/core/utility/initState_impl.hpp
View File

@ -29,6 +29,7 @@
#include <opm/core/utility/Units.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <cmath>
namespace Opm
@ -512,15 +513,23 @@ namespace Opm
State& state)
{
const int num_phases = props.numPhases();
const PhaseUsage pu = phaseUsageFromDeck(deck);
if (num_phases != pu.num_phases) {
THROW("initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
state.init(grid, num_phases);
if (deck.hasField("EQUIL")) {
if (num_phases != 2) {
THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
THROW("initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
}
// Set saturations depending on oil-water contact.
const EQUIL& equil= deck.getEQUIL();
if (equil.equil.size() != 1) {
THROW("No region support yet.");
THROW("initStateFromDeck(): No region support yet.");
}
const double woc = equil.equil[0].water_oil_contact_depth_;
initWaterOilContact(grid, props, woc, WaterBelow, state);
@ -528,37 +537,64 @@ namespace Opm
const double datum_z = equil.equil[0].datum_depth_;
const double datum_p = equil.equil[0].datum_depth_pressure_;
initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state);
} else if (deck.hasField("SWAT") && deck.hasField("PRESSURE")) {
// Set saturations from SWAT, pressure from PRESSURE.
} else if (deck.hasField("PRESSURE")) {
// Set saturations from SWAT/SGAS, pressure from PRESSURE.
std::vector<double>& s = state.saturation();
std::vector<double>& p = state.pressure();
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
const int num_cells = grid.number_of_cells;
if (num_phases == 2) {
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[2*c] = sw_deck[c_deck];
s[2*c + 1] = 1.0 - s[2*c];
p[c] = p_deck[c_deck];
if (!pu.phase_used[BlackoilPhases::Aqua]) {
// oil-gas: we require SGAS
if (!deck.hasField("SGAS")) {
THROW("initStateFromDeck(): missing SGAS keyword in 2-phase init");
}
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[2*c + gpos] = sg_deck[c_deck];
s[2*c + opos] = 1.0 - sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
// water-oil or water-gas: we require SWAT
if (!deck.hasField("SWAT")) {
THROW("initStateFromDeck(): missing SWAT keyword in 2-phase init");
}
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int nwpos = (wpos + 1) % 2;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[2*c + wpos] = sw_deck[c_deck];
s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
p[c] = p_deck[c_deck];
}
}
} else if (num_phases == 3) {
if (!deck.hasField("SGAS")) {
THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found).");
const bool has_swat_sgas = deck.hasField("SWAT") && deck.hasField("SGAS");
if (!has_swat_sgas) {
THROW("initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
}
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[3*c] = sw_deck[c_deck];
s[3*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
s[3*c + 2] = sg_deck[c_deck];
s[3*c + wpos] = sw_deck[c_deck];
s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
s[3*c + gpos] = sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
}
} else {
THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE.");
THROW("initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS.");
}
// Finally, init face pressures.

6
opm/core/utility/miscUtilities.cpp
View File

@ -374,7 +374,10 @@ namespace Opm
if (perf_rate > 0.0) {
// perf_rate is a total inflow rate, we want a water rate.
if (wells->type[w] != INJECTOR) {
std::cout << "**** Warning: crossflow in well with index " << w << " ignored." << std::endl;
std::cout << "**** Warning: crossflow in well "
<< w << " perf " << perf - wells->well_connpos[w]
<< " ignored. Rate was "
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
perf_rate = 0.0;
} else {
ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
@ -550,6 +553,7 @@ namespace Opm
{
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
ASSERT(int(flow_rates_per_well_cell.size()) == wells.well_connpos[nw]);
phase_flow_per_well.resize(nw * np);
for (int wix = 0; wix < nw; ++wix) {
for (int phase = 0; phase < np; ++phase) {

124
opm/core/utility/miscUtilitiesBlackoil.cpp
View File

@ -21,7 +21,10 @@
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/grid.h>
#include <opm/core/newwells.h>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <functional>
@ -31,53 +34,74 @@
namespace Opm
{
/// @brief Computes injected and produced volumes of all phases.
/// @brief Computes injected and produced surface volumes of all phases.
/// Note 1: assumes that only the first phase is injected.
/// Note 2: assumes that transport has been done with an
/// implicit method, i.e. that the current state
/// gives the mobilities used for the preceding timestep.
/// @param[in] props fluid and rock properties.
/// @param[in] p pressure (one value per cell)
/// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all P phases)
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements.
/// Note 3: Gives surface volume values, not reservoir volumes
/// (as the incompressible version of the function does).
/// Also, assumes that transport_src is given in surface volumes
/// for injector terms!
/// @param[in] props fluid and rock properties.
/// @param[in] state state variables (pressure, sat, surfvol)
/// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements.
void computeInjectedProduced(const BlackoilPropertiesInterface& props,
const std::vector<double>& press,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& src,
const BlackoilState& state,
const std::vector<double>& transport_src,
const double dt,
double* injected,
double* produced)
{
const int num_cells = src.size();
const int np = s.size()/src.size();
if (int(s.size()) != num_cells*np) {
THROW("Sizes of s and src vectors do not match.");
const int num_cells = transport_src.size();
if (props.numCells() != num_cells) {
THROW("Size of transport_src vector does not match number of cells in props.");
}
const int np = props.numPhases();
if (int(state.saturation().size()) != num_cells*np) {
THROW("Sizes of state vectors do not match number of cells.");
}
const std::vector<double>& press = state.pressure();
const std::vector<double>& s = state.saturation();
const std::vector<double>& z = state.surfacevol();
std::fill(injected, injected + np, 0.0);
std::fill(produced, produced + np, 0.0);
std::vector<double> visc(np);
std::vector<double> mob(np);
std::vector<double> A(np*np);
std::vector<double> prod_resv_phase(np);
std::vector<double> prod_surfvol(np);
for (int c = 0; c < num_cells; ++c) {
if (src[c] > 0.0) {
injected[0] += src[c]*dt;
} else if (src[c] < 0.0) {
const double flux = -src[c]*dt;
if (transport_src[c] > 0.0) {
// Inflowing transport source is a surface volume flux
// for the first phase.
injected[0] += transport_src[c]*dt;
} else if (transport_src[c] < 0.0) {
// Outflowing transport source is a total reservoir
// volume flux.
const double flux = -transport_src[c]*dt;
const double* sat = &s[np*c];
props.relperm(1, sat, &c, &mob[0], 0);
props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0);
props.matrix(1, &press[c], &z[np*c], &c, &A[0], 0);
double totmob = 0.0;
for (int p = 0; p < np; ++p) {
mob[p] /= visc[p];
totmob += mob[p];
}
std::fill(prod_surfvol.begin(), prod_surfvol.end(), 0.0);
for (int p = 0; p < np; ++p) {
produced[p] += (mob[p]/totmob)*flux;
prod_resv_phase[p] = (mob[p]/totmob)*flux;
for (int q = 0; q < np; ++q) {
prod_surfvol[q] += prod_resv_phase[p]*A[q + np*p];
}
}
for (int p = 0; p < np; ++p) {
produced[p] += prod_surfvol[p];
}
}
}
@ -251,4 +275,58 @@ namespace Opm
}
}
/// Compute two-phase transport source terms from well terms.
/// Note: Unlike the incompressible version of this function,
/// this version computes surface volume injection rates,
/// production rates are still total reservoir volumes.
/// \param[in] props Fluid and rock properties.
/// \param[in] wells Wells data structure.
/// \param[in] well_state Well pressures and fluxes.
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
/// (+) positive inflow of first (water) phase (surface volume),
/// (-) negative total outflow of both phases (reservoir volume).
void computeTransportSource(const BlackoilPropertiesInterface& props,
const Wells* wells,
const WellState& well_state,
std::vector<double>& transport_src)
{
int nc = props.numCells();
transport_src.clear();
transport_src.resize(nc, 0.0);
// Well contributions.
if (wells) {
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases;
if (np != 2) {
THROW("computeTransportSource() requires a 2 phase case.");
}
std::vector<double> A(np*np);
for (int w = 0; w < nw; ++w) {
const double* comp_frac = wells->comp_frac + np*w;
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
const int perf_cell = wells->well_cells[perf];
double perf_rate = well_state.perfRates()[perf];
if (perf_rate > 0.0) {
// perf_rate is a total inflow reservoir rate, we want a surface water rate.
if (wells->type[w] != INJECTOR) {
std::cout << "**** Warning: crossflow in well "
<< w << " perf " << perf - wells->well_connpos[w]
<< " ignored. Reservoir rate was "
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
perf_rate = 0.0;
} else {
ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
perf_rate *= comp_frac[0]; // Water reservoir volume rate.
props.matrix(1, &well_state.perfPress()[perf], comp_frac, &perf_cell, &A[0], 0);
perf_rate *= A[0]; // Water surface volume rate.
}
}
transport_src[perf_cell] += perf_rate;
}
}
}
}
} // namespace Opm

51
opm/core/utility/miscUtilitiesBlackoil.hpp
View File

@ -22,36 +22,40 @@
#include <vector>
struct UnstructuredGrid;
struct Wells;
namespace Opm
{
class BlackoilPropertiesInterface;
class BlackoilState;
class WellState;
/// @brief Computes injected and produced volumes of all phases.
/// @brief Computes injected and produced surface volumes of all phases.
/// Note 1: assumes that only the first phase is injected.
/// Note 2: assumes that transport has been done with an
/// implicit method, i.e. that the current state
/// gives the mobilities used for the preceding timestep.
/// @param[in] props fluid and rock properties.
/// @param[in] p pressure (one value per cell)
/// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all P phases)
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements.
/// Note 3: Gives surface volume values, not reservoir volumes
/// (as the incompressible version of the function does).
/// Also, assumes that transport_src is given in surface volumes
/// for injector terms!
/// @param[in] props fluid and rock properties.
/// @param[in] state state variables (pressure, sat, surfvol)
/// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements.
void computeInjectedProduced(const BlackoilPropertiesInterface& props,
const std::vector<double>& p,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& src,
const BlackoilState& state,
const std::vector<double>& transport_src,
const double dt,
double* injected,
double* produced);
/// @brief Computes total mobility for a set of saturation values.
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
@ -66,6 +70,7 @@ namespace Opm
const std::vector<double>& s,
std::vector<double>& totmob);
/// @brief Computes total mobility and omega for a set of saturation values.
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
@ -131,6 +136,22 @@ namespace Opm
const double* saturation,
double* surfacevol);
/// Compute two-phase transport source terms from well terms.
/// Note: Unlike the incompressible version of this function,
/// this version computes surface volume injection rates,
/// production rates are still total reservoir volumes.
/// \param[in] props Fluid and rock properties.
/// \param[in] wells Wells data structure.
/// \param[in] well_state Well pressures and fluxes.
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
/// (+) positive inflow of first (water) phase (surface volume),
/// (-) negative total outflow of both phases (reservoir volume).
void computeTransportSource(const BlackoilPropertiesInterface& props,
const Wells* wells,
const WellState& well_state,
std::vector<double>& transport_src);
} // namespace Opm
#endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED

14
opm/core/utility/writeVtkData.hpp
View File

@ -26,25 +26,23 @@
#include <vector>
#include <tr1/array>
#include <iosfwd>
#include <opm/core/utility/DataMap.hpp>
struct UnstructuredGrid;
namespace Opm
{
/// Intended to map strings (giving the output field names) to data.
typedef std::map<std::string, const std::vector<double>*> DataMap;
/// Vtk output for cartesian grids.
void writeVtkData(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size,
const DataMap& data,
std::ostream& os);
const std::tr1::array<double, 3>& cell_size,
const DataMap& data,
std::ostream& os);
/// Vtk output for general grids.
void writeVtkData(const UnstructuredGrid& grid,
const DataMap& data,
std::ostream& os);
const DataMap& data,
std::ostream& os);
} // namespace Opm
#endif // OPM_WRITEVTKDATA_HEADER_INCLUDED

26
opm/core/wells/WellsGroup.cpp
View File

@ -572,6 +572,7 @@ namespace Opm
break;
}
const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase);
const double total_reinjected = - total_produced; // Production negative, injection positive
const double my_guide_rate = injectionGuideRate(true);
for (size_t i = 0; i < children_.size(); ++i) {
// Apply for all children.
@ -580,11 +581,11 @@ namespace Opm
const double children_guide_rate = children_[i]->injectionGuideRate(true);
#ifdef DIRTY_WELLCTRL_HACK
children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_,
false);
#else
children_[i]->applyInjGroupControl(InjectionSpecification::RATE,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_,
false);
#endif
}
@ -600,15 +601,15 @@ namespace Opm
if (phaseUsage().phase_used[BlackoilPhases::Vapour]) {
total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour);
}
const double my_guide_rate = injectionGuideRate(true);
const double total_reinjected = - total_produced; // Production negative, injection positive
const double my_guide_rate = injectionGuideRate(true);
for (size_t i = 0; i < children_.size(); ++i) {
// Apply for all children.
// Note, we do _not_ want to call the applyProdGroupControl in this object,
// as that would check if we're under group control, something we're not.
const double children_guide_rate = children_[i]->injectionGuideRate(true);
children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().voidage_replacment_fraction_,
false);
}
@ -760,16 +761,17 @@ namespace Opm
wells_->ctrls[self_index_]->current = ~ wells_->ctrls[self_index_]->current;
}
}
std::pair<WellNode*, double> WellNode::getWorstOffending(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase,
ProductionSpecification::ControlMode mode)
{
const int np = phaseUsage().num_phases;
const int index = self_index_*np;
return std::make_pair<WellNode*, double>(this, rateByMode(&well_reservoirrates_phase[index],
&well_surfacerates_phase[index],
mode));
return std::pair<WellNode*, double>(this,
rateByMode(&well_reservoirrates_phase[index],
&well_surfacerates_phase[index],
mode));
}
void WellNode::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,
@ -781,7 +783,7 @@ namespace Opm
&& (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) {
return;
}
if (!wells_->type[self_index_] == INJECTOR) {
if (wells_->type[self_index_] != INJECTOR) {
ASSERT(target == 0.0);
return;
}
@ -858,12 +860,12 @@ namespace Opm
std::cout << "Returning" << std::endl;
return;
}
if (!wells_->type[self_index_] == PRODUCER) {
if (wells_->type[self_index_] != PRODUCER) {
ASSERT(target == 0.0);
return;
}
// We're a producer, so we need to negate the input
double ntarget = target;
double ntarget = -target;
double distr[3] = { 0.0, 0.0, 0.0 };
const int* phase_pos = phaseUsage().phase_pos;

61
opm/core/wells/WellsManager.cpp
View File

@ -230,6 +230,14 @@ namespace Opm
/// Construct from existing wells object.
WellsManager::WellsManager(struct Wells* W)
: w_(clone_wells(W))
{
}
/// Construct wells from deck.
WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
@ -548,7 +556,7 @@ namespace Opm
cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
} else if (wci_line.injector_type_[0] == 'G') {
if (!pu.phase_used[BlackoilPhases::Vapour]) {
THROW("Water phase not used, yet found water-injecting well.");
THROW("Gas phase not used, yet found gas-injecting well.");
}
cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
}
@ -720,32 +728,31 @@ namespace Opm
#endif
if (deck.hasField("WELOPEN")) {
const WELOPEN& welopen = deck.getWELOPEN();
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
WelopenLine line = welopen.welopen[i];
std::string wellname = line.well_;
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
if (it == well_names_to_index.end()) {
THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
}
int index = it->second;
if (line.openshutflag_ == "SHUT") {
// We currently don't care if the well is open or not.
/// \TODO Should this perhaps be allowed? I.e. should it be if(well_shut) { shutwell(); } else { /* do nothing*/ }?
//ASSERT(w_->ctrls[index]->current < 0);
} else if (line.openshutflag_ == "OPEN") {
//ASSERT(w_->ctrls[index]->current >= 0);
} else {
THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
}
// We revert back to it's original control.
// Note that this is OK as ~~ = id.
w_->ctrls[index]->current = ~w_->ctrls[index]->current;
}
const WELOPEN& welopen = deck.getWELOPEN();
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
WelopenLine line = welopen.welopen[i];
std::string wellname = line.well_;
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
if (it == well_names_to_index.end()) {
THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
}
const int index = it->second;
if (line.openshutflag_ == "SHUT") {
int& cur_ctrl = w_->ctrls[index]->current;
if (cur_ctrl >= 0) {
cur_ctrl = ~cur_ctrl;
}
ASSERT(w_->ctrls[index]->current < 0);
} else if (line.openshutflag_ == "OPEN") {
int& cur_ctrl = w_->ctrls[index]->current;
if (cur_ctrl < 0) {
cur_ctrl = ~cur_ctrl;
}
ASSERT(w_->ctrls[index]->current >= 0);
} else {
THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
}
}
}
// Build the well_collection_ well group hierarchy.

8
opm/core/wells/WellsManager.hpp
View File

@ -42,8 +42,12 @@ namespace Opm
public:
/// Default constructor -- no wells.
WellsManager();
/// Construct from mrst type output.
/// Wellmanger is not properly initialized
/// Construct from existing wells object.
/// WellsManager is not properly initialised in the sense that the logic to
/// manage control switching does not exist.
///
/// @param[in] W Existing wells object.
WellsManager(struct Wells* W);
/// Construct from input deck and grid.

205
tests/test_wells.cpp Normal file
View File

@ -0,0 +1,205 @@
/*
Copyright 2012 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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE WellsModuleTest
#include <boost/test/unit_test.hpp>
#include <opm/core/newwells.h>
#include <iostream>
#include <vector>
#include <boost/shared_ptr.hpp>
BOOST_AUTO_TEST_CASE(Construction)
{
const int nphases = 2;
const int nwells = 2;
const int nperfs = 2;
boost::shared_ptr<Wells> W(create_wells(nphases, nwells, nperfs),
destroy_wells);
if (W) {
int cells[] = { 0, 9 };
double WI = 1.0;
const double ifrac[] = { 1.0, 0.0 };
const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0],
&WI, "INJECTOR", W.get());
const double pfrac[] = { 0.0, 0.0 };
const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1],
&WI, "PRODUCER", W.get());
if (ok0 && ok1) {
BOOST_CHECK_EQUAL(W->number_of_phases, nphases);
BOOST_CHECK_EQUAL(W->number_of_wells , nwells );
BOOST_CHECK_EQUAL(W->well_connpos[0], 0);
BOOST_CHECK_EQUAL(W->well_connpos[1], 1);
BOOST_CHECK_EQUAL(W->well_connpos[W->number_of_wells], nperfs);
BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[0]], cells[0]);
BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[1]], cells[1]);
BOOST_CHECK_EQUAL(W->WI[W->well_connpos[0]], WI);
BOOST_CHECK_EQUAL(W->WI[W->well_connpos[1]], WI);
using std::string;
BOOST_CHECK_EQUAL(string(W->name[0]), string("INJECTOR"));
BOOST_CHECK_EQUAL(string(W->name[1]), string("PRODUCER"));
}
}
}
BOOST_AUTO_TEST_CASE(Controls)
{
const int nphases = 2;
const int nwells = 1;
const int nperfs = 2;
boost::shared_ptr<Wells> W(create_wells(nphases, nwells, nperfs),
destroy_wells);
if (W) {
int cells[] = { 0 , 9 };
double WI [] = { 1.0, 1.0 };
const double ifrac[] = { 1.0, 0.0 };
const bool ok = add_well(INJECTOR, 0.0, nperfs, &ifrac[0], &cells[0],
&WI[0], "INJECTOR", W.get());
if (ok) {
const double distr[] = { 1.0, 0.0 };
const bool ok1 = append_well_controls(BHP, 1, &distr[0],
0, W.get());
const bool ok2 = append_well_controls(SURFACE_RATE, 1,
&distr[0], 0, W.get());
if (ok1 && ok2) {
WellControls* ctrls = W->ctrls[0];
BOOST_CHECK_EQUAL(ctrls->num , 2);
BOOST_CHECK_EQUAL(ctrls->current, -1);
set_current_control(0, 0, W.get());
BOOST_CHECK_EQUAL(ctrls->current, 0);
set_current_control(0, 1, W.get());
BOOST_CHECK_EQUAL(ctrls->current, 1);
BOOST_CHECK_EQUAL(ctrls->type[0], BHP);
BOOST_CHECK_EQUAL(ctrls->type[1], SURFACE_RATE);
BOOST_CHECK_EQUAL(ctrls->target[0], 1.0);
BOOST_CHECK_EQUAL(ctrls->target[1], 1.0);
}
}
}
}
BOOST_AUTO_TEST_CASE(Copy)
{
const int nphases = 2;
const int nwells = 2;
const int nperfs = 2;
boost::shared_ptr<Wells> W1(create_wells(nphases, nwells, nperfs),
destroy_wells);
boost::shared_ptr<Wells> W2;
if (W1) {
int cells[] = { 0, 9 };
const double WI = 1.0;
const double ifrac[] = { 1.0, 0.0 };
const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0],
&WI, "INJECTOR", W1.get());
const double pfrac[] = { 0.0, 0.0 };
const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1],
&WI, "PRODUCER", W1.get());
bool ok = ok0 && ok1;
for (int w = 0; ok && (w < W1->number_of_wells); ++w) {
const double distr[] = { 1.0, 0.0 };
const bool okc1 = append_well_controls(BHP, 1, &distr[0],
w, W1.get());
const bool okc2 = append_well_controls(SURFACE_RATE, 1,
&distr[0], w,
W1.get());
ok = okc1 && okc2;
}
if (ok) {
W2.reset(clone_wells(W1.get()), destroy_wells);
}
}
if (W2) {
BOOST_CHECK_EQUAL(W2->number_of_phases, W1->number_of_phases);
BOOST_CHECK_EQUAL(W2->number_of_wells , W1->number_of_wells );
BOOST_CHECK_EQUAL(W2->well_connpos[0] , W1->well_connpos[0] );
for (int w = 0; w < W1->number_of_wells; ++w) {
using std::string;
BOOST_CHECK_EQUAL(string(W2->name[w]), string(W1->name[w]));
BOOST_CHECK_EQUAL( W2->type[w] , W1->type[w] );
BOOST_CHECK_EQUAL(W2->well_connpos[w + 1],
W1->well_connpos[w + 1]);
for (int j = W1->well_connpos[w];
j < W1->well_connpos[w + 1]; ++j) {
BOOST_CHECK_EQUAL(W2->well_cells[j], W1->well_cells[j]);
BOOST_CHECK_EQUAL(W2->WI [j], W1->WI [j]);
}
BOOST_CHECK(W1->ctrls[w] != 0);
BOOST_CHECK(W2->ctrls[w] != 0);
WellControls* c1 = W1->ctrls[w];
WellControls* c2 = W2->ctrls[w];
BOOST_CHECK_EQUAL(c2->num , c1->num );
BOOST_CHECK_EQUAL(c2->current, c1->current);
for (int c = 0; c < c1->num; ++c) {
BOOST_CHECK_EQUAL(c2->type [c], c1->type [c]);
BOOST_CHECK_EQUAL(c2->target[c], c1->target[c]);
for (int p = 0; p < W1->number_of_phases; ++p) {
BOOST_CHECK_EQUAL(c2->distr[c*W1->number_of_phases + p],
c1->distr[c*W1->number_of_phases + p]);
}
}
}
}
}