diff --git a/opm/core/fluid/BlackoilPropertiesBasic.cpp b/opm/core/fluid/BlackoilPropertiesBasic.cpp
new file mode 100644
index 000000000..7a70d2c6d
--- /dev/null
+++ b/opm/core/fluid/BlackoilPropertiesBasic.cpp
@@ -0,0 +1,261 @@
+/*
+ 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 .
+*/
+
+#include
+#include
+#include
+#include
+
+namespace Opm
+{
+
+ BlackoilPropertiesBasic::BlackoilPropertiesBasic(const Dune::parameter::ParameterGroup& param,
+ const int dim,
+ const int num_cells)
+ {
+ double poro = param.getDefault("porosity", 1.0);
+ using namespace Dune::unit;
+ using namespace Dune::prefix;
+ double perm = param.getDefault("permeability", 100)*milli*darcy;
+ rock_.init(dim, num_cells, poro, perm);
+ pvt_.init(param);
+ satprops_.init(param);
+ }
+
+ BlackoilPropertiesBasic::~BlackoilPropertiesBasic()
+ {
+ }
+
+
+ /// \return D, the number of spatial dimensions.
+ int BlackoilPropertiesBasic::numDimensions() const
+ {
+ return rock_.numDimensions();
+ }
+
+ /// \return N, the number of cells.
+ int BlackoilPropertiesBasic::numCells() const
+ {
+ return rock_.numCells();
+ }
+
+ /// \return Array of N porosity values.
+ const double* BlackoilPropertiesBasic::porosity() const
+ {
+ return rock_.porosity();
+ }
+
+ /// \return Array of ND^2 permeability values.
+ /// The D^2 permeability values for a cell are organized as a matrix,
+ /// which is symmetric (so ordering does not matter).
+ const double* BlackoilPropertiesBasic::permeability() const
+ {
+ return rock_.permeability();
+ }
+
+
+ // ---- Fluid interface ----
+
+ /// \return P, the number of phases (also the number of components).
+ int BlackoilPropertiesBasic::numPhases() const
+ {
+ return pvt_.numPhases();
+ }
+
+ /// \param[in] n Number of data points.
+ /// \param[in] p Array of n pressure values.
+ /// \param[in] z Array of nP surface volume values.
+ /// \param[in] cells Array of n cell indices to be associated with the p and z values.
+ /// \param[out] mu Array of nP viscosity values, array must be valid before calling.
+ /// \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
+ {
+ if (dmudp) {
+ THROW("BlackoilPropertiesBasic::viscosity() -- derivatives of viscosity not yet implemented.");
+ } else {
+ pvt_.mu(n, p, z, mu);
+ }
+ }
+
+ /// \param[in] n Number of data points.
+ /// \param[in] p Array of n pressure values.
+ /// \param[in] z Array of nP surface volume values.
+ /// \param[in] cells Array of n cell indices to be associated with the p and z values.
+ /// \param[out] A Array of nP^2 values, array must be valid before calling.
+ /// The P^2 values for a cell give the matrix A = RB^{-1} which
+ /// relates z to u by z = Au. The matrices are output in Fortran order.
+ /// \param[out] dAdp If non-null: array of nP^2 matrix derivative values,
+ /// 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 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) {
+ double* m = A + i*np*np;
+ std::fill(m, m + np*np, 0.0);
+ // Diagonal entries only.
+ for (int phase = 0; phase < np; ++phase) {
+ m[phase + phase*np] = 1.0/B[phase];
+ }
+ }
+
+ // Derivative of A matrix.
+ if (dAdp) {
+#pragma omp parallel for
+ for (int i = 0; i < n; ++i) {
+ double* m = dAdp + i*np*np;
+ std::fill(m, m + np*np, 0.0);
+ }
+ }
+#if 0
+ // This is copied from BlackoilPropertiesFromDeck.
+ const int np = numPhases();
+ B_.resize(n*np);
+ R_.resize(n*np);
+ if (dAdp) {
+ dB_.resize(n*np);
+ dR_.resize(n*np);
+ pvt_.dBdp(n, p, z, &B_[0], &dB_[0]);
+ pvt_.dRdp(n, p, z, &R_[0], &dR_[0]);
+ } else {
+ pvt_.B(n, p, z, &B_[0]);
+ pvt_.R(n, p, z, &R_[0]);
+ }
+ const int* phase_pos = pvt_.phasePosition();
+ bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
+ pvt_.phaseUsed()[BlackoilPhases::Vapour];
+ const int o = phase_pos[BlackoilPhases::Liquid];
+ const int g = phase_pos[BlackoilPhases::Vapour];
+
+ // Compute A matrix
+#pragma omp parallel for
+ for (int i = 0; i < n; ++i) {
+ double* m = A + i*np*np;
+ std::fill(m, m + np*np, 0.0);
+ // Diagonal entries.
+ for (int phase = 0; phase < np; ++phase) {
+ m[phase + phase*np] = 1.0/B_[i*np + phase];
+ }
+ // Off-diagonal entries.
+ if (oil_and_gas) {
+ m[o + g*np] = R_[i*np + g]/B_[i*np + g];
+ m[g + o*np] = R_[i*np + o]/B_[i*np + o];
+ }
+ }
+
+ // Derivative of A matrix.
+ if (dAdp) {
+#pragma omp parallel for
+ for (int i = 0; i < n; ++i) {
+ double* m = dAdp + i*np*np;
+ std::fill(m, m + np*np, 0.0);
+ // Diagonal entries.
+ for (int phase = 0; phase < np; ++phase) {
+ m[phase + phase*np] = -dB_[i*np + phase]/B_[i*np + phase]*B_[i*np + phase];
+ }
+ // Off-diagonal entries.
+ if (oil_and_gas) {
+ m[o + g*np] = m[g + g*np]*R_[i*np + g] + dR_[i*np + g]/B_[i*np + g];
+ m[g + o*np] = m[o + o*np]*R_[i*np + o] + dR_[i*np + o]/B_[i*np + o];
+ }
+ }
+ }
+#endif
+ }
+
+ /// \param[in] n Number of data points.
+ /// \param[in] A Array of nP^2 values, where the P^2 values for a cell give the
+ /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
+ /// are assumed to be in Fortran order, and are typically the result
+ /// 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 int np = numPhases();
+ const double* sdens = pvt_.surfaceDensities();
+#pragma omp parallel for
+ for (int i = 0; i < n; ++i) {
+ for (int phase = 0; phase < np; ++phase) {
+ rho[np*i + phase] = 0.0;
+ for (int comp = 0; comp < np; ++comp) {
+ rho[np*i + phase] += A[n*np*np + np*phase + comp]*sdens[comp];
+ }
+ }
+ }
+ }
+
+ /// \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 ...)
+ void BlackoilPropertiesBasic::relperm(const int n,
+ const double* s,
+ const int* /*cells*/,
+ double* kr,
+ double* dkrds) const
+ {
+ satprops_.relperm(n, s, kr, dkrds);
+ }
+
+
+ /// \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 ...)
+ void BlackoilPropertiesBasic::capPress(const int n,
+ const double* s,
+ const int* /*cells*/,
+ double* pc,
+ double* dpcds) const
+ {
+ satprops_.relperm(n, s, pc, dpcds);
+ }
+
+
+
+} // namespace Opm
+
diff --git a/opm/core/fluid/BlackoilPropertiesBasic.hpp b/opm/core/fluid/BlackoilPropertiesBasic.hpp
new file mode 100644
index 000000000..b23d8d8e4
--- /dev/null
+++ b/opm/core/fluid/BlackoilPropertiesBasic.hpp
@@ -0,0 +1,159 @@
+/*
+ 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 .
+*/
+
+#ifndef OPM_BLACKOILPROPERTIESBASIC_HEADER_INCLUDED
+#define OPM_BLACKOILPROPERTIESBASIC_HEADER_INCLUDED
+
+
+#include
+#include
+#include
+#include
+
+namespace Opm
+{
+
+ /// Concrete class implementing the blackoil property interface,
+ /// reading all necessary input from parameters.
+ class BlackoilPropertiesBasic : public BlackoilPropertiesInterface
+ {
+ 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
+ BlackoilPropertiesBasic(const Dune::parameter::ParameterGroup& param,
+ const int dim,
+ const int num_cells);
+
+ /// Destructor.
+ virtual ~BlackoilPropertiesBasic();
+
+
+ // ---- Rock interface ----
+
+ /// \return D, the number of spatial dimensions.
+ virtual int numDimensions() const;
+
+ /// \return N, the number of cells.
+ virtual int numCells() const;
+
+ /// \return Array of N porosity values.
+ virtual const double* porosity() const;
+
+ /// \return Array of ND^2 permeability values.
+ /// The D^2 permeability values for a cell are organized as a matrix,
+ /// which is symmetric (so ordering does not matter).
+ virtual const double* permeability() const;
+
+
+ // ---- Fluid interface ----
+
+ /// \return P, the number of phases (also the number of components).
+ virtual int numPhases() const;
+
+ /// \param[in] n Number of data points.
+ /// \param[in] p Array of n pressure values.
+ /// \param[in] z Array of nP surface volume values.
+ /// \param[in] cells Array of n cell indices to be associated with the p and z values.
+ /// \param[out] mu Array of nP viscosity values, array must be valid before calling.
+ /// \param[out] dmudp If non-null: array of nP viscosity derivative values,
+ /// array must be valid before calling.
+ virtual void viscosity(const int n,
+ const double* p,
+ const double* z,
+ const int* cells,
+ double* mu,
+ double* dmudp) const;
+
+ /// \param[in] n Number of data points.
+ /// \param[in] p Array of n pressure values.
+ /// \param[in] z Array of nP surface volume values.
+ /// \param[in] cells Array of n cell indices to be associated with the p and z values.
+ /// \param[out] A Array of nP^2 values, array must be valid before calling.
+ /// The P^2 values for a cell give the matrix A = RB^{-1} which
+ /// relates z to u by z = Au. The matrices are output in Fortran order.
+ /// \param[out] dAdp If non-null: array of nP^2 matrix derivative values,
+ /// array must be valid before calling. The matrices are output
+ /// in Fortran order.
+ virtual void matrix(const int n,
+ const double* p,
+ const double* z,
+ const int* cells,
+ double* A,
+ double* dAdp) const;
+
+
+ /// \param[in] n Number of data points.
+ /// \param[in] A Array of nP^2 values, where the P^2 values for a cell give the
+ /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
+ /// are assumed to be in Fortran order, and are typically the result
+ /// of a call to the method matrix().
+ /// \param[out] rho Array of nP density values, array must be valid before calling.
+ virtual void density(const int n,
+ const double* A,
+ double* rho) const;
+
+ /// \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 ...)
+ virtual void relperm(const int n,
+ const double* s,
+ const int* cells,
+ double* kr,
+ double* dkrds) const;
+
+
+ /// \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 ...)
+ virtual void capPress(const int n,
+ const double* s,
+ const int* cells,
+ double* pc,
+ double* dpcds) const;
+
+ private:
+ RockBasic rock_;
+ PvtPropertiesBasic pvt_;
+ SaturationPropsBasic satprops_;
+ };
+
+
+
+} // namespace Opm
+
+
+#endif // OPM_BLACKOILPROPERTIESBASIC_HEADER_INCLUDED
diff --git a/opm/core/fluid/PvtPropertiesBasic.cpp b/opm/core/fluid/PvtPropertiesBasic.cpp
new file mode 100644
index 000000000..5cd0f727b
--- /dev/null
+++ b/opm/core/fluid/PvtPropertiesBasic.cpp
@@ -0,0 +1,141 @@
+/*
+ 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 .
+*/
+
+
+
+#include
+#include
+#include
+
+
+namespace Opm
+{
+
+ PvtPropertiesBasic::PvtPropertiesBasic()
+ {
+ }
+
+
+ void PvtPropertiesBasic::init(const Dune::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);
+
+ // Setting mu and rho from parameters
+ using namespace Dune::prefix;
+ using namespace Dune::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]);
+ }
+ }
+
+ const double* PvtPropertiesBasic::surfaceDensities() const
+ {
+ return &density_[0];
+ }
+
+
+ int PvtPropertiesBasic::numPhases() const
+ {
+ return density_.size();
+ }
+
+
+
+ void PvtPropertiesBasic::mu(const int n,
+ const double* /*p*/,
+ const double* /*z*/,
+ double* output_mu) const
+ {
+ 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 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 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 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 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
diff --git a/opm/core/fluid/PvtPropertiesBasic.hpp b/opm/core/fluid/PvtPropertiesBasic.hpp
new file mode 100644
index 000000000..9314417db
--- /dev/null
+++ b/opm/core/fluid/PvtPropertiesBasic.hpp
@@ -0,0 +1,96 @@
+/*
+ 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 .
+*/
+
+#ifndef OPM_PVTPROPERTIESBASIC_HEADER_INCLUDED
+#define OPM_PVTPROPERTIESBASIC_HEADER_INCLUDED
+
+#include
+
+namespace Opm
+{
+
+ /// Class collecting simple pvt properties for 1-3 phases.
+ /// All phases are incompressible and have constant viscosities.
+ /// For all the methods, the following apply: p and z are unused.
+ /// Output arrays shall be of size n*numPhases(), and must be valid
+ /// before calling the method.
+ /// NOTE: This class is intentionally similar to BlackoilPvtProperties.
+ class PvtPropertiesBasic
+ {
+ public:
+ /// Default constructor.
+ 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 Dune::parameter::ParameterGroup& param);
+
+ /// Number of active phases.
+ int numPhases() const;
+
+ /// Densities of stock components at surface conditions.
+ /// \return Array of size numPhases().
+ const double* surfaceDensities() const;
+
+ /// Viscosity as a function of p and z.
+ void mu(const int n,
+ const double* p,
+ const double* z,
+ double* output_mu) const;
+
+ /// Formation volume factor as a function of p and z.
+ void B(const int n,
+ const double* p,
+ const double* z,
+ double* output_B) const;
+
+ /// Formation volume factor and p-derivative as functions of p and z.
+ void dBdp(const int n,
+ const double* p,
+ const double* z,
+ double* output_B,
+ double* output_dBdp) const;
+
+ /// Solution factor as a function of p and z.
+ void R(const int n,
+ const double* p,
+ const double* z,
+ double* output_R) const;
+
+ /// Solution factor and p-derivative as functions of p and z.
+ void dRdp(const int n,
+ const double* p,
+ const double* z,
+ double* output_R,
+ double* output_dRdp) const;
+
+ private:
+ std::vector density_;
+ std::vector viscosity_;
+ std::vector formation_volume_factor_;
+ };
+
+}
+
+
+
+#endif // OPM_PVTPROPERTIESBASIC_HEADER_INCLUDED
diff --git a/opm/core/fluid/RockBasic.hpp b/opm/core/fluid/RockBasic.hpp
new file mode 100644
index 000000000..965637807
--- /dev/null
+++ b/opm/core/fluid/RockBasic.hpp
@@ -0,0 +1,79 @@
+/*
+ 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 .
+*/
+
+#ifndef OPM_ROCKBASIC_HEADER_INCLUDED
+#define OPM_ROCKBASIC_HEADER_INCLUDED
+
+
+#include
+
+
+namespace Opm
+{
+
+ class RockBasic
+ {
+ public:
+ /// Default constructor.
+ RockBasic();
+
+ /// Initialize with homogenous porosity and permeability.
+ void init(const int dimensions,
+ const int num_cells,
+ const double poro,
+ const double perm);
+
+ /// \return D, the number of spatial dimensions.
+ int numDimensions() const
+ {
+ return dimensions_;
+ }
+
+ /// \return N, the number of cells.
+ int numCells() const
+ {
+ return porosity_.size();
+ }
+
+ /// \return Array of N porosity values.
+ const double* porosity() const
+ {
+ return &porosity_[0];
+ }
+
+ /// \return Array of ND^2 permeability values.
+ /// The D^2 permeability values for a cell are organized as a matrix,
+ /// which is symmetric (so ordering does not matter).
+ const double* permeability() const
+ {
+ return &permeability_[0];
+ }
+
+ private:
+ int dimensions_;
+ std::vector porosity_;
+ std::vector permeability_;
+ };
+
+
+
+} // namespace Opm
+
+
+#endif // OPM_ROCKBASIC_HEADER_INCLUDED
diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp
index 9cd227038..2fa0829b3 100644
--- a/opm/core/fluid/RockFromDeck.cpp
+++ b/opm/core/fluid/RockFromDeck.cpp
@@ -39,7 +39,7 @@ namespace Opm
- // ---- ReadFromDeck methods ----
+ // ---- RockFromDeck methods ----
/// Default constructor.
diff --git a/opm/core/fluid/SaturationPropsBasic.cpp b/opm/core/fluid/SaturationPropsBasic.cpp
new file mode 100644
index 000000000..077576b04
--- /dev/null
+++ b/opm/core/fluid/SaturationPropsBasic.cpp
@@ -0,0 +1,196 @@
+/*
+ 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 .
+*/
+
+#include
+#include
+#include
+
+namespace Opm
+{
+
+
+ // ---------- Helper functions ----------
+
+ namespace {
+
+ 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 KrFunQuadratic
+ {
+ double kr(double s)
+ {
+ return s*s;
+ }
+ double dkrds(double s)
+ {
+ return 2.0*s;
+ }
+ };
+
+
+ template
+ 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;
+ }
+#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]);
+ }
+ }
+ }
+
+
+ } // anon namespace
+
+
+
+ // ---------- Class methods ----------
+
+
+
+ /// Default constructor.
+ SaturationPropsBasic::SaturationPropsBasic()
+ {
+ }
+
+
+
+
+ /// Initialize from parameters.
+ void SaturationPropsBasic::init(const Dune::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);
+ }
+ num_phases_ = num_phases;
+ std::string rpf = param.getDefault("relperm_func", std::string("Linear"));
+ if (rpf == "Constant") {
+ relperm_func_ = Constant;
+ } else if (rpf == "Linear") {
+ relperm_func_ = Linear;
+ } else if (rpf == "Quadratic") {
+ relperm_func_ = Quadratic;
+ } else {
+ THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf);
+ }
+ }
+
+
+
+
+ /// 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 ...)
+ void SaturationPropsBasic::relperm(const int n,
+ 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_);
+ }
+ }
+
+
+
+
+ /// 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 ...)
+ void SaturationPropsBasic::capPress(const int n,
+ const double* /*s*/,
+ double* pc,
+ double* dpcds) const
+ {
+ std::fill(pc, pc + num_phases_*n, 0.0);
+ if (dpcds) {
+ std::fill(dpcds, dpcds + num_phases_*num_phases_*n, 0.0);
+ }
+ }
+
+
+
+
+
+
+} // namespace Opm
+
+
diff --git a/opm/core/fluid/SaturationPropsBasic.hpp b/opm/core/fluid/SaturationPropsBasic.hpp
new file mode 100644
index 000000000..3a8e97fa4
--- /dev/null
+++ b/opm/core/fluid/SaturationPropsBasic.hpp
@@ -0,0 +1,86 @@
+/*
+ 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 .
+*/
+
+#ifndef OPM_SATURATIONPROPSBASIC_HEADER_INCLUDED
+#define OPM_SATURATIONPROPSBASIC_HEADER_INCLUDED
+
+#include
+
+namespace Opm
+{
+
+
+ /// Class encapsulating basic saturation function behaviour,
+ /// by which we mean constant, linear or quadratic relative
+ /// permeability functions for a maximum of two phases,
+ /// and zero capillary pressure.
+ class SaturationPropsBasic
+ {
+ public:
+ /// Default constructor.
+ 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".
+ void init(const Dune::parameter::ParameterGroup& param);
+
+ /// 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 ...)
+ void relperm(const int n,
+ const double* s,
+ double* kr,
+ double* dkrds) const;
+
+ /// 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 ...)
+ void capPress(const int n,
+ const double* s,
+ double* pc,
+ double* dpcds) const;
+
+ private:
+ enum RelPermFunc { Constant, Linear, Quadratic };
+ int num_phases_;
+ RelPermFunc relperm_func_;
+ };
+
+
+
+} // namespace Opm
+
+
+
+
+#endif // OPM_SATURATIONPROPSBASIC_HEADER_INCLUDED
diff --git a/opm/core/fluid/blackoil/BlackoilDefs.hpp b/opm/core/fluid/blackoil/BlackoilDefs.hpp
index 08e44cb83..c60d66d0d 100644
--- a/opm/core/fluid/blackoil/BlackoilDefs.hpp
+++ b/opm/core/fluid/blackoil/BlackoilDefs.hpp
@@ -54,6 +54,18 @@ namespace Opm
{ return data_[i]; }
T& operator[](int i)
{ return data_[i]; }
+ void operator += (const SmallVec& sm)
+ {
+ for (int i = 0; i < N; ++i) {
+ data_[i] += sm.data_[i];
+ }
+ }
+ void operator *= (const T& scalar)
+ {
+ for (int i = 0; i < N; ++i) {
+ data_[i] *= scalar;
+ }
+ }
template
void assign(const U& elem)
{