From f6b9dbc57b48b5eafb0523d192cf96de445e99f7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 17 Jan 2012 09:30:43 +0100 Subject: [PATCH 1/7] Added operators += and *= for SmallVec. --- opm/core/fluid/blackoil/BlackoilDefs.hpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) 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) { From e3dc0dfa0e6d0b0824d6e03fb277ac48f88ef1af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 17 Jan 2012 09:31:05 +0100 Subject: [PATCH 2/7] Minor typo fix. --- opm/core/fluid/RockFromDeck.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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. From c065b8633a570b0f36e47ed94bce50cbc26c758f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 17 Jan 2012 09:31:25 +0100 Subject: [PATCH 3/7] Added RockBasic class. --- opm/core/fluid/RockBasic.hpp | 79 ++++++++++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) create mode 100644 opm/core/fluid/RockBasic.hpp 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 From 125cb33c368a95dfb07af69088817187f7dbb976 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 17 Jan 2012 10:25:49 +0100 Subject: [PATCH 4/7] Added SaturationPropsBasic class. --- opm/core/fluid/SaturationPropsBasic.cpp | 196 ++++++++++++++++++++++++ opm/core/fluid/SaturationPropsBasic.hpp | 86 +++++++++++ 2 files changed, 282 insertions(+) create mode 100644 opm/core/fluid/SaturationPropsBasic.cpp create mode 100644 opm/core/fluid/SaturationPropsBasic.hpp 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 From bf7b0ce5ecd647a6cda7594d509f718ce084d39a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 17 Jan 2012 11:17:40 +0100 Subject: [PATCH 5/7] Added class PvtPropertiesBasic. --- opm/core/fluid/PvtPropertiesBasic.cpp | 141 ++++++++++++++++++++++++++ opm/core/fluid/PvtPropertiesBasic.hpp | 96 ++++++++++++++++++ 2 files changed, 237 insertions(+) create mode 100644 opm/core/fluid/PvtPropertiesBasic.cpp create mode 100644 opm/core/fluid/PvtPropertiesBasic.hpp 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 From 2d07dd206281c375ac3f77587b0fa6886980df18 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 17 Jan 2012 12:19:50 +0100 Subject: [PATCH 6/7] Added class BlackoilPropertiesBasic. --- opm/core/fluid/BlackoilPropertiesBasic.cpp | 261 +++++++++++++++++++++ opm/core/fluid/BlackoilPropertiesBasic.hpp | 152 ++++++++++++ 2 files changed, 413 insertions(+) create mode 100644 opm/core/fluid/BlackoilPropertiesBasic.cpp create mode 100644 opm/core/fluid/BlackoilPropertiesBasic.hpp diff --git a/opm/core/fluid/BlackoilPropertiesBasic.cpp b/opm/core/fluid/BlackoilPropertiesBasic.cpp new file mode 100644 index 000000000..cd7a95f4b --- /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..7a81a8ada --- /dev/null +++ b/opm/core/fluid/BlackoilPropertiesBasic.hpp @@ -0,0 +1,152 @@ +/* + 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. + 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 From 5bc77babb45db5dd40197caa53f6dd0bb0b5713f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 17 Jan 2012 14:16:28 +0100 Subject: [PATCH 7/7] Documented parameters, and made unit behaviour according to doc (perm in mD etc.) --- opm/core/fluid/BlackoilPropertiesBasic.cpp | 2 +- opm/core/fluid/BlackoilPropertiesBasic.hpp | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/opm/core/fluid/BlackoilPropertiesBasic.cpp b/opm/core/fluid/BlackoilPropertiesBasic.cpp index cd7a95f4b..7a70d2c6d 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.cpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.cpp @@ -32,7 +32,7 @@ namespace Opm double poro = param.getDefault("porosity", 1.0); using namespace Dune::unit; using namespace Dune::prefix; - double perm = param.getDefault("permeability", 100*milli*darcy); + double perm = param.getDefault("permeability", 100)*milli*darcy; rock_.init(dim, num_cells, poro, perm); pvt_.init(param); satprops_.init(param); diff --git a/opm/core/fluid/BlackoilPropertiesBasic.hpp b/opm/core/fluid/BlackoilPropertiesBasic.hpp index 7a81a8ada..b23d8d8e4 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.hpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.hpp @@ -35,6 +35,13 @@ 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 BlackoilPropertiesBasic(const Dune::parameter::ParameterGroup& param, const int dim, const int num_cells);