From 63a57ff03b571ac7a90be5654d81684eec64d68f Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 16 Apr 2012 17:49:02 +0200 Subject: [PATCH 01/51] Added Doxygen documentation to Units.hpp. --- opm/core/utility/Units.hpp | 67 ++++++++++++++++++++++++++++---------- 1 file changed, 50 insertions(+), 17 deletions(-) diff --git a/opm/core/utility/Units.hpp b/opm/core/utility/Units.hpp index 5e6b9f3d..a2939783 100644 --- a/opm/core/utility/Units.hpp +++ b/opm/core/utility/Units.hpp @@ -38,6 +38,7 @@ namespace Opm { namespace prefix + /// Conversion prefix for units. { const double micro = 1.0e-6; const double milli = 1.0e-3; @@ -49,68 +50,99 @@ namespace Opm } // namespace prefix namespace unit + /// Definition of various units. + /// All the units are defined in terms of international standard units (SI). + /// Example of use: We define a variable \c k which gives a permeability. We want to set \c k to \f$1\,mD\f$. + /// \code + /// using namespace Opm::unit + /// double k = 0.001*darcy; + /// \endcode + /// We can also use one of the prefixes defined in Opm::prefix + /// \code + /// using namespace Opm::unit + /// using namespace Opm::prefix + /// double k = 1.0*milli*darcy; + /// \endcode { - // Common powers + ///\name Common powers + /// @{ inline double square(double v) { return v * v; } inline double cubic (double v) { return v * v * v; } + /// @} // -------------------------------------------------------------- // Basic (fundamental) units and conversions // -------------------------------------------------------------- - // Length: + /// \name Length + /// @{ const double meter = 1; const double inch = 2.54 * prefix::centi*meter; const double feet = 12 * inch; + /// @} - // Time: + /// \name Time + /// @{ const double second = 1; const double minute = 60 * second; const double hour = 60 * minute; const double day = 24 * hour; const double year = 365 * day; + /// @} - // Volume + /// \name Volume + /// @{ const double stb = 0.158987294928 * cubic(meter); + /// @} - - // Mass: + /// \name Mass + /// @{ const double kilogram = 1; - // http://en.wikipedia.org/wiki/Pound_(mass)#Avoirdupois_pound const double pound = 0.45359237 * kilogram; + /// @} // -------------------------------------------------------------- // Standardised constants // -------------------------------------------------------------- + /// \name Standardised constant + /// @{ const double gravity = 9.80665 * meter/square(second); + /// @} // -------------------------------------------------------------- // Derived units and conversions // -------------------------------------------------------------- - // Force: + /// \name Force + /// @{ const double Newton = kilogram*meter / square(second); // == 1 const double lbf = pound * gravity; // Pound-force + /// @} - // Pressure: + /// \name Pressure + /// @{ const double Pascal = Newton / square(meter); // == 1 const double barsa = 100000 * Pascal; const double atm = 101325 * Pascal; const double psia = lbf / square(inch); + /// @} - // Viscosity: + /// \name Viscosity + /// @{ const double Pas = Pascal * second; // == 1 const double Poise = prefix::deci*Pas; + /// @} - // Permeability: - // - // A porous medium with a permeability of 1 darcy permits a - // flow (flux) of 1 cm³/s of a fluid with viscosity 1 cP (1 - // mPa·s) under a pressure gradient of 1 atm/cm acting across - // an area of 1 cm². - // + /// \name Permeability + /// @{ + /// + /// A porous medium with a permeability of 1 darcy permits a + /// flow (flux) of \f$1\,cm^2/s\f$ of a fluid with viscosity \f$1\,cP\f$ ( + /// \f$1\,mPa\cdot s\f$) under a pressure gradient of \f$1\,atm/cm\f$ acting across + /// an area of \f$1\,cm^2\f$. + /// namespace perm_details { const double p_grad = atm / (prefix::centi*meter); const double area = square(prefix::centi*meter); @@ -122,6 +154,7 @@ namespace Opm // == 9.869232667160130e-13 [m^2] } const double darcy = perm_details::darcy; + /// @} // Unit conversion support. // From 9d8cb130c8f884964933960b7c1e5aa12d93a867 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 16 Apr 2012 17:49:33 +0200 Subject: [PATCH 02/51] Update python script to generate documentation figures. --- generate_doc_figures.py | 42 ++++++++++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 9 deletions(-) diff --git a/generate_doc_figures.py b/generate_doc_figures.py index 9364c6f4..2829d816 100755 --- a/generate_doc_figures.py +++ b/generate_doc_figures.py @@ -17,14 +17,14 @@ dp.AmbientColor = [1, 0, 0] view = GetActiveView() view.Background = [1, 1, 1] camera = GetActiveCamera() -camera.SetPosition(20, 20, 110) -camera.SetViewUp(0.5,0.3,0.7) +camera.SetPosition(4, -6, 5) +camera.SetViewUp(-0.19, 0.4, 0.9) camera.SetViewAngle(30) -camera.SetFocalPoint(1,1,0.5) +camera.SetFocalPoint(1.5, 1.5, 1) Render() WriteImage("Figure/tutorial1.png") Hide(grid) -# remove("tutorial1.vtu") +remove("tutorial1.vtu") # tutorial 2 call("tutorials/tutorial2") @@ -33,12 +33,10 @@ grid.UpdatePipeline() Show(grid) dp = GetDisplayProperties(grid) dp.Representation = 'Surface' -data = grid.GetCellDataInformation() -pressure = data.GetArray('pressure') -pmin = pressure.GetRange()[0] -pmax = pressure.GetRange()[1] -dp.LookupTable = MakeBlueToRedLT(0.5*(pmin+pmax)-0.2*(pmax-pmin), 0.5*(pmin+pmax)-0.2*(pmax+pmin)) dp.ColorArrayName = 'pressure' +pres = grid.CellData.GetArray(0) +pres_lookuptable = GetLookupTableForArray( "pressure", 1, RGBPoints=[pres.GetRange()[0], 1, 0, 0, pres.GetRange()[1], 0, 0, 1] ) +dp.LookupTable = pres_lookuptable view = GetActiveView() view.Background = [1, 1, 1] camera = GetActiveCamera() @@ -48,4 +46,30 @@ camera.SetViewAngle(30) camera.SetFocalPoint(20, 20, 0.5) Render() WriteImage("Figure/tutorial2.png") +Hide(grid) # remove("tutorial2.vtu") + +# # tutorial 3 +call("tutorials/tutorial3") +cases = ["000", "005", "010", "015", "019"] +for case in cases: + grid = servermanager.sources.XMLUnstructuredGridReader(FileName="tutorial3-"+case+".vtu") + grid.UpdatePipeline() + Show(grid) + dp = GetDisplayProperties(grid) + dp.Representation = 'Surface' + dp.ColorArrayName = 'saturation' + sat = grid.CellData.GetArray(1) + sat_lookuptable = GetLookupTableForArray( "saturation", 1, RGBPoints=[0, 1, 0, 0, 1, 0, 0, 1]) + dp.LookupTable = sat_lookuptable + view.Background = [1, 1, 1] + camera = GetActiveCamera() + camera.SetPosition(100, 100, 550) + camera.SetViewUp(0, 1, 0) + camera.SetViewAngle(30) + camera.SetFocalPoint(100, 100, 5) + Render() + WriteImage("Figure/tutorial3-"+case+".png") +Hide(grid) +# remove("tutorial3.vtu") + From 6b37e011ae6b5917e684824629e3f56ec6d5dced Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 16 Apr 2012 17:50:29 +0200 Subject: [PATCH 03/51] Added Multiphase tutorial. --- tutorials/Makefile.am | 5 + tutorials/tutorial3.cpp | 283 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 288 insertions(+) create mode 100644 tutorials/tutorial3.cpp diff --git a/tutorials/Makefile.am b/tutorials/Makefile.am index de6e91e4..7bd7f06e 100644 --- a/tutorials/Makefile.am +++ b/tutorials/Makefile.am @@ -13,3 +13,8 @@ if UMFPACK noinst_PROGRAMS += tutorial2 tutorial2_SOURCES = tutorial2.cpp endif + +if UMFPACK +noinst_PROGRAMS += tutorial3 +tutorial3_SOURCES = tutorial3.cpp +endif diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp new file mode 100644 index 00000000..76d2cbdb --- /dev/null +++ b/tutorials/tutorial3.cpp @@ -0,0 +1,283 @@ +/* + 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 . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +/// \page tutorial3 Multiphase flow +/// Multiphase flow + +class TwophaseFluid +{ +public: + TwophaseFluid(const Opm::IncompPropertiesInterface& props) + : props_(props), + smin_(props.numCells()*props.numPhases()), + smax_(props.numCells()*props.numPhases()) + { + const int num_cells = props.numCells(); + std::vector cells(num_cells); + for (int c = 0; c < num_cells; ++c) { + cells[c] = c; + } + props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]); + } + + double density(int phase) const + { + return props_.density()[phase]; + } + + template + void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const + { + props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]); + const double* mu = props_.viscosity(); + mob[0] /= mu[0]; + mob[1] /= mu[1]; + // Recall that we use Fortran ordering for kr derivatives, + // therefore dmob[i*2 + j] is row j and column i of the + // matrix. + // Each row corresponds to a kr function, so which mu to + // divide by also depends on the row, j. + dmob[0*2 + 0] /= mu[0]; + dmob[0*2 + 1] /= mu[1]; + dmob[1*2 + 0] /= mu[0]; + dmob[1*2 + 1] /= mu[1]; + } + + template + void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const + { + double pcow[2]; + double dpcow[4]; + props_.capPress(1, &s[0], &c, pcow, dpcow); + pcap = pcow[0]; + dpcap = dpcow[0]; + } + + double s_min(int c) const + { + return smin_[2*c + 0]; + } + + double s_max(int c) const + { + return smax_[2*c + 0]; + } + +private: + const Opm::IncompPropertiesInterface& props_; + std::vector smin_; + std::vector smax_; +}; + +typedef Opm::SinglePointUpwindTwoPhase TransportModel; + +using namespace Opm::ImplicitTransportDefault; + +typedef NewtonVectorCollection< ::std::vector > NVecColl; +typedef JacobianSystem < struct CSRMatrix, NVecColl > JacSys; + +template +class MaxNorm { +public: + static double + norm(const Vector& v) { + return AccumulationNorm ::norm(v); + } +}; + +typedef Opm::ImplicitTransport TransportSolver; + + + +int main () +{ + + int dim = 3; + int nx = 20; + int ny = 20; + int nz = 1; + double dx = 10.; + double dy = 10.; + double dz = 10.; + using namespace Opm; + GridManager grid(nx, ny, nz, dx, dy, dz); + int num_cells = grid.c_grid()->number_of_cells; + int num_faces = grid.c_grid()->number_of_faces; + + int num_phases = 2; + using namespace unit; + using namespace prefix; + std::vector rho(2, 1000.); + std::vector mu(2, 1.*centi*Poise); + double porosity = 0.5; + double k = 10*milli*darcy; + + SaturationPropsBasic::RelPermFunc rel_perm_func; + rel_perm_func = SaturationPropsBasic::Linear; + + IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu, + porosity, k, dim, num_cells); + + + + const double *grav = 0; + std::vector omega; + + LinearSolverUmfpack linsolver; + IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver); + + + TwophaseFluid fluid(props); + + std::vector src(num_cells, 0.0); + src[0] = 1.; + src[num_cells-1] = -1.; + + + std::vector empty_wdp; + std::vector empty_well_bhp; + std::vector empty_well_flux; + + TransportSource* tsrc = create_transport_source(2, 2); + double ssrc[] = { 1.0, 0.0 }; + double ssink[] = { 0.0, 1.0 }; + double zdummy[] = { 0.0, 0.0 }; + for (int cell = 0; cell < num_cells; ++cell) { + if (src[cell] > 0.0) { + append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc); + } else if (src[cell] < 0.0) { + append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc); + } + } + + std::vector porevol; + computePorevolume(*grid.c_grid(), props, porevol); + + const bool guess_old_solution = true; + TransportModel model (fluid, *grid.c_grid(), porevol, grav, guess_old_solution); + TransportSolver tsolver(model); + + double dt = 0.1*day; + int num_time_steps = 20; + + TwophaseState state; + ImplicitTransportDetails::NRReport rpt; + ImplicitTransportDetails::NRControl ctrl; + std::vector totmob; + + std::vector allcells(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells[cell] = cell; + } + + FlowBCManager bcs; + + // Linear solver init. + using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver; + CSRMatrixUmfpackSolver linsolve; + + state.init(*grid.c_grid()); + // By default: initialise water saturation to minimum everywhere. + state.setWaterSat(allcells, props, TwophaseState::MinSat); + + + std::ostringstream vtkfilename; + + for (int i = 0; i < num_time_steps; ++i) { + + computeTotalMobility(props, allcells, state.saturation(), totmob); + psolver.solve(totmob, omega, src, empty_wdp, bcs.c_bcs(), + state.pressure(), state.faceflux(), empty_well_bhp, + empty_well_flux); + tsolver.solve(*grid.c_grid(), tsrc, dt, ctrl, state, linsolve, rpt); + + vtkfilename.str(""); + vtkfilename << "tutorial3-" << std::setw(3) << std::setfill('0') << i << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + Opm::writeVtkData(*grid.c_grid(), dm, vtkfile); + + } +} + + +/// \page tutorial3 +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +///
\image html tutorial3-000.png \image html tutorial3-005.png
\image html tutorial3-010.png \image html tutorial3-015.png
\image html tutorial3-019.png
+ +/// \page tutorial3 +/// \section sourcecode Complete source code. +/// \include tutorial3.cpp From e8ad818ce29c2d3d291e28e24f0ae5cb48b0ef0e Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 16 Apr 2012 17:53:10 +0200 Subject: [PATCH 04/51] Added setup functions which take arguments directly (do not use ParamterGroup) --- opm/core/fluid/IncompPropertiesBasic.cpp | 20 ++++++++++++++++++++ opm/core/fluid/IncompPropertiesBasic.hpp | 11 +++++++++++ opm/core/fluid/PvtPropertiesBasic.cpp | 14 ++++++++++++++ opm/core/fluid/PvtPropertiesBasic.hpp | 5 +++++ opm/core/fluid/SaturationPropsBasic.hpp | 13 ++++++++++++- 5 files changed, 62 insertions(+), 1 deletion(-) diff --git a/opm/core/fluid/IncompPropertiesBasic.cpp b/opm/core/fluid/IncompPropertiesBasic.cpp index 1b78481a..8201dacd 100644 --- a/opm/core/fluid/IncompPropertiesBasic.cpp +++ b/opm/core/fluid/IncompPropertiesBasic.cpp @@ -46,6 +46,26 @@ namespace Opm pvt_.mu(1, 0, 0, &viscosity_[0]); } + IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases, + const SaturationPropsBasic::RelPermFunc& relpermfunc, + const std::vector& rho, + const std::vector& mu, + const double porosity, + const double permeability, + const int dim, + const int num_cells) + { + rock_.init(dim, num_cells, porosity, permeability); + 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]); + } + IncompPropertiesBasic::~IncompPropertiesBasic() { } diff --git a/opm/core/fluid/IncompPropertiesBasic.hpp b/opm/core/fluid/IncompPropertiesBasic.hpp index e52044af..9784398f 100644 --- a/opm/core/fluid/IncompPropertiesBasic.hpp +++ b/opm/core/fluid/IncompPropertiesBasic.hpp @@ -53,6 +53,17 @@ namespace Opm 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& rho, + const std::vector& mu, + const double porosity, + const double permeability, + const int dim, + const int num_cells); + /// Destructor. virtual ~IncompPropertiesBasic(); diff --git a/opm/core/fluid/PvtPropertiesBasic.cpp b/opm/core/fluid/PvtPropertiesBasic.cpp index 56e20cfe..53b6cbfd 100644 --- a/opm/core/fluid/PvtPropertiesBasic.cpp +++ b/opm/core/fluid/PvtPropertiesBasic.cpp @@ -59,6 +59,20 @@ namespace Opm } } + void PvtPropertiesBasic::init(const int num_phases, + const std::vector rho, + const std::vector mu) + { + 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_ = mu; + } + const double* PvtPropertiesBasic::surfaceDensities() const { return &density_[0]; diff --git a/opm/core/fluid/PvtPropertiesBasic.hpp b/opm/core/fluid/PvtPropertiesBasic.hpp index d628344b..7bf234a7 100644 --- a/opm/core/fluid/PvtPropertiesBasic.hpp +++ b/opm/core/fluid/PvtPropertiesBasic.hpp @@ -44,6 +44,11 @@ namespace Opm /// mu1 [mu2, mu3] (1.0) Viscosity in cP void init(const parameter::ParameterGroup& param); + /// Initialize from argument Basic multi phase fluid pvt properties. + void init(const int num_phases, + const std::vector rho, + const std::vector mu); + /// Number of active phases. int numPhases() const; diff --git a/opm/core/fluid/SaturationPropsBasic.hpp b/opm/core/fluid/SaturationPropsBasic.hpp index a00db3f2..789e0a35 100644 --- a/opm/core/fluid/SaturationPropsBasic.hpp +++ b/opm/core/fluid/SaturationPropsBasic.hpp @@ -45,6 +45,16 @@ namespace Opm /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". void init(const parameter::ParameterGroup& param); + enum RelPermFunc { Constant, Linear, Quadratic }; + + /// Initialize from arguments a basic Saturation property. + void init(const int num_phases, + const RelPermFunc& relperm_func) + { + num_phases_ = num_phases; + relperm_func_ = relperm_func; + } + /// \return P, the number of phases. int numPhases() const; @@ -83,8 +93,9 @@ namespace Opm void satRange(const int n, double* smin, double* smax) const; + + private: - enum RelPermFunc { Constant, Linear, Quadratic }; int num_phases_; RelPermFunc relperm_func_; }; From f9675df45051c3374f8ef27e7e7a33dedf26afad Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 16 Apr 2012 17:54:42 +0200 Subject: [PATCH 05/51] Added Doxygen comments in tutorials. --- tutorials/tutorial1.cpp | 29 +++++++++++------------- tutorials/tutorial2.cpp | 49 +++++++++++++++++++++++------------------ 2 files changed, 40 insertions(+), 38 deletions(-) diff --git a/tutorials/tutorial1.cpp b/tutorials/tutorial1.cpp index b0bec824..567a594c 100644 --- a/tutorials/tutorial1.cpp +++ b/tutorials/tutorial1.cpp @@ -26,43 +26,40 @@ #include "config.h" #endif // HAVE_CONFIG_H -/// \page tutorial1 A simple carthesian grid -/// This tutorial explains how to construct a simple carthesian grid.\n\n -/// We construct a 2x2 two dimensional carthesian grid with 4 blocks of equal size. +/// \page tutorial1 A simple cartesian grid +/// This tutorial explains how to construct a simple cartesian grid. #include #include #include -#include -#include #include -#include #include #include // ----------------- Main program ----------------- -/// \page tutorial1 +/// \page tutorial1 /// \code int main() { /// \endcode /// \page tutorial1 - /// By setting nz = 1, we make the grid two dimensional + /// We set the number of blocks in each direction. /// \code int nx = 3; int ny = 3; - int nz = 1; + int nz = 2; /// \endcode - /// The size of each block is 1x1x1. We use standard units (SI) + /// The size of each block is 1x1x1. The default units are allways the + /// standard units (SI). But other units can easily be dealt with, see Opm::unit. /// \code double dx = 1.0; double dy = 1.0; double dz = 1.0; - /// \endcode - /// \page tutorial1 - /// One of the constructors of the class Opm::GridManager takes nx,ny,nz,dx,dy,dz - /// and construct the corresponding carthesian grid. + /// \endcode + /// \page tutorial1 + /// One of the constructors of the class Opm::GridManager takes nx,ny,nz,dx,dy,dz + /// and construct the corresponding cartesian grid. /// \code Opm::GridManager grid(nx, ny, nz, dx, dy, dz); /// \endcode @@ -87,7 +84,7 @@ int main() /// We read the the vtu output file in \a Paraview and obtain the following grid. /// \image html tutorial1.png -/// \page tutorial1 +/// \page tutorial1 /// \section sourcecode Source code. -/// \include tutorial1.cpp +/// \include tutorial1.cpp diff --git a/tutorials/tutorial2.cpp b/tutorials/tutorial2.cpp index 27f69030..ce3d1233 100644 --- a/tutorials/tutorial2.cpp +++ b/tutorials/tutorial2.cpp @@ -18,14 +18,14 @@ */ -/// \page tutorial2 Flow Solver for a single phase +/// \page tutorial2 Flow Solver for a single phase /// \details The flow equations consist of the mass conservation equation /// \f[\nabla\cdot u=q\f] and the Darcy law \f[u=-\frac{1}{\mu}K\nabla p.\f] Here, /// \f$u\f$ denotes the velocity and \f$p\f$ the pressure. The permeability tensor is -/// given by \f$K\f$ and \f$\mu\f$ denotes the viscosity. -/// -/// We solve the flow equations for a carthesian grid and we set the source term -/// \f$q\f$ be zero except at the left-lower and right-upper corner, where it is equal +/// given by \f$K\f$ and \f$\mu\f$ denotes the viscosity. +/// +/// We solve the flow equations for a cartesian grid and we set the source term +/// \f$q\f$ be zero except at the left-lower and right-upper corner, where it is equal /// with opposite sign (inflow equal to outflow). @@ -36,16 +36,14 @@ #include #include #include -#include -#include #include -#include #include #include #include #include #include #include +#include /// \page tutorial2 /// \section commentedcode Commented code: @@ -54,7 +52,7 @@ int main() { /// \endcode /// \page tutorial2 - /// We construct a carthesian grid + /// We construct a cartesian grid /// \code int dim = 3; int nx = 40; @@ -71,19 +69,25 @@ int main() int num_faces = grid.c_grid()->number_of_faces; /// \endcode /// \page tutorial2 - /// We define the viscosity (unit: cP). + /// \details + /// We define a fluid viscosity equal to \f$1\,cP\f$. We use + /// the namespaces Opm::unit + /// and Opm::prefix to deal with the units. /// \code - double mu = 1.0; + using namespace Opm::unit; + using namespace Opm::prefix; + double mu = 1.0*centi*Poise; /// \endcode /// \page tutorial2 - /// We define the permeability (unit: mD). + /// \details + /// We define a permeability equal to \f$100\,mD\f$. /// \code - double k = 100.0; + double k = 100.0*milli*darcy; /// \endcode /// \page tutorial2 /// \details /// We set up the permeability tensor and compute the mobility for each cell. - /// The permeability tensor is flattened in a vector. + /// The resulting permeability matrix is flattened in a vector. /// \code std::vector permeability(num_cells*dim*dim, 0.); std::vector mob(num_cells); @@ -96,7 +100,8 @@ int main() /// \endcode /// \page tutorial2 - /// We choose the UMFPACK linear solver for the pressure solver. + /// \details + /// We take UMFPACK as the linear solver for the pressure solver (This library has therefore to be installed.) /// \code Opm::LinearSolverUmfpack linsolver; /// \endcode @@ -115,7 +120,7 @@ int main() src[num_cells-1] = -100.; /// \endcode /// \page tutorial2 - /// \details We set up the boundary conditions. We do not modify them. + /// \details We set up the boundary conditions. We do not modify them. /// By default, we obtain no outflow boundary conditions. /// \code /// \code @@ -133,21 +138,21 @@ int main() std::vector faceflux(num_faces); std::vector well_bhp; std::vector well_flux; - /// \endcode + /// \endcode /// \page tutorial2 /// \details /// We declare the gravity term which is required by the pressure solver (see /// Opm::IncompTpfa.solve()). In the absence of gravity, an empty vector is required. /// \code - std::vector omega; + std::vector omega; /// \endcode - + /// \page tutorial2 /// \details /// We declare the wdp term which is required by the pressure solver (see /// Opm::IncompTpfa.solve()). In the absence of wells, an empty vector is required. /// \code - std::vector wdp; + std::vector wdp; /// \endcode /// \page tutorial2 @@ -170,10 +175,10 @@ int main() } /// \endcode /// \page tutorial2 -/// We read the the vtu output file in \a Paraview and obtain the following pressure +/// We read the vtu output file in \a Paraview and obtain the following pressure /// distribution. \image html tutorial2.png /// \page tutorial2 /// \section sourcecode Complete source code. -/// \include tutorial2.cpp +/// \include tutorial2.cpp From e1df8085c3011a6e190f700345976f698d149f04 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 16 Apr 2012 17:57:29 +0200 Subject: [PATCH 06/51] Let Doxyfile generate documentation of whole opm-core. --- Doxyfile | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Doxyfile b/Doxyfile index 0be5414c..ed1b6bb6 100644 --- a/Doxyfile +++ b/Doxyfile @@ -610,9 +610,7 @@ WARN_LOGFILE = # directories like "/usr/src/myproject". Separate the files or directories # with spaces. -# INPUT = tutorials/tutorial1.cpp opm/core/GridManager.hpp opm/core/utility/writeVtkData.hpp opm/core/utility/writeVtkData.cpp -INPUT = tutorials opm/core/linalg/LinearSolverUmfpack.hpp opm/core/pressure/IncompTpfa.hpp opm/core/pressure/FlowBCManager.hpp opm/core/utility/miscUtilities.hpp opm/core/GridManager.hpp opm/core/utility/writeVtkData.hpp opm/core/grid.h - +INPUT = opm/core tutorials examples # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is From 156e88a7e28969bcf0421463994130e256ddb08e Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 17 Apr 2012 17:24:19 +0200 Subject: [PATCH 07/51] Updated tutorials. --- tutorials/tutorial1.cpp | 3 +- tutorials/tutorial2.cpp | 4 +- tutorials/tutorial3.cpp | 387 ++++++++++++++++++++++++++++++++-------- 3 files changed, 316 insertions(+), 78 deletions(-) diff --git a/tutorials/tutorial1.cpp b/tutorials/tutorial1.cpp index 567a594c..080244fd 100644 --- a/tutorials/tutorial1.cpp +++ b/tutorials/tutorial1.cpp @@ -39,6 +39,7 @@ // ----------------- Main program ----------------- /// \page tutorial1 +/// \section commentedsource1 Commented source code: /// \code int main() { @@ -85,6 +86,6 @@ int main() /// \image html tutorial1.png /// \page tutorial1 -/// \section sourcecode Source code. +/// \section completecode1 Complete source code: /// \include tutorial1.cpp diff --git a/tutorials/tutorial2.cpp b/tutorials/tutorial2.cpp index ce3d1233..c34d9658 100644 --- a/tutorials/tutorial2.cpp +++ b/tutorials/tutorial2.cpp @@ -46,7 +46,7 @@ #include /// \page tutorial2 -/// \section commentedcode Commented code: +/// \section commentedcode2 Commented code: /// \code int main() { @@ -180,5 +180,5 @@ int main() /// \page tutorial2 -/// \section sourcecode Complete source code. +/// \section completecode2 Complete source code: /// \include tutorial2.cpp diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp index 76d2cbdb..85e8ac23 100644 --- a/tutorials/tutorial3.cpp +++ b/tutorials/tutorial3.cpp @@ -51,83 +51,129 @@ #include /// \page tutorial3 Multiphase flow -/// Multiphase flow +/// The Darcy law gives +/// \f[u_\alpha= -\frac1{\mu_\alpha} K_\alpha\nabla p_\alpha\f] +/// where \f$\mu_\alpha\f$ and \f$K_\alpha\f$ represent the viscosity +/// and the permeability tensor for each phase \f$\alpha\f$. In the two phase +/// case, we have either \f$\alpha=w\f$ or \f$\alpha=o\f$. +/// In this tutorial, we do not take into account capillary pressure so that +/// \f$p=p_w=p_o\f$ and gravity +/// effects. We denote by \f$K\f$ the absolute permeability tensor and each phase +/// permeability is defined through its relative permeability by the expression +/// \f[K_\alpha=k_{r\alpha}K.\f] +/// The phase mobility are defined as +/// \f[\lambda_\alpha=\frac{k_{r\alpha}}{\mu_\alpha}\f] +/// so that the Darcy law may be rewritten as +/// \f[u_\alpha= -\lambda_\alpha K\nabla p.\f] +/// The conservation of mass for each phase writes: +/// \f[\frac{\partial}{\partial t}(\phi\rho_\alpha s_\alpha)+\nabla\cdot(\rho_\alpha u_\alpha)=q_\alpha\f] +/// where \f$s_\alpha\f$ denotes the saturation of the phase \f$\alpha\f$ and \f$q_\alpha\f$ is a source term. Let +/// us consider a two phase flow with oil and water. We assume that the phases are incompressible. Since +/// \f$s_w+s_o=1\f$, we get +/// \f[ \nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o}.\f] +/// where we define +/// \f[u=u_w+u_o.\f] +/// Let the total mobility be equal to +/// \f[\lambda=\lambda_w+\lambda_o\f] +/// Then, we have +/// \f[u=-\lambda K\nabla p.\f] +/// The set of equations +/// \f[\nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o},\quad u=-\lambda K\nabla p.\f] +/// is referred to as the pressure equation. We introduce +/// the fractional flow \f$f_w\f$ +/// as +/// \f[f_w=\frac{\lambda_w}{\lambda_w+\lambda_o}\f] +/// and obtain +/// \f[\phi\frac{\partial s}{\partial t}+\nabla\cdot(f_w u)=\frac{q_w}{\rho_w}\f] +/// which is referred to as the transport equation. The pressure and +/// transport equation are coupled. In this tutorial, we implement a splitting scheme, +/// where, at each time step, we decouple the two equations. We solve first +/// the pressure equation and then update the water saturation by solving +/// the transport equation assuming that \f$u\f$ is constant in time in the time step +/// interval we are considering. +/// \page tutorial3 +/// \section commentedcode3 Commented code: +/// \page tutorial3 +/// \details +/// We define a class which computes mobility, capillary pressure and +/// the minimum and maximum saturation value for each cell. +/// \code class TwophaseFluid { public: - TwophaseFluid(const Opm::IncompPropertiesInterface& props) - : props_(props), - smin_(props.numCells()*props.numPhases()), - smax_(props.numCells()*props.numPhases()) - { - const int num_cells = props.numCells(); - std::vector cells(num_cells); - for (int c = 0; c < num_cells; ++c) { - cells[c] = c; - } - props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]); - } + /// \endcode + /// \page tutorial3 + /// \details Constructor operator. Takes in the fluid properties defined + /// \c props + /// \code + TwophaseFluid(const Opm::IncompPropertiesInterface& props); + /// \endcode - double density(int phase) const - { - return props_.density()[phase]; - } + /// \page tutorial3 + /// \details Density for each phase. + /// \code + double density(int phase) const; + /// \endcode - template - void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const - { - props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]); - const double* mu = props_.viscosity(); - mob[0] /= mu[0]; - mob[1] /= mu[1]; - // Recall that we use Fortran ordering for kr derivatives, - // therefore dmob[i*2 + j] is row j and column i of the - // matrix. - // Each row corresponds to a kr function, so which mu to - // divide by also depends on the row, j. - dmob[0*2 + 0] /= mu[0]; - dmob[0*2 + 1] /= mu[1]; - dmob[1*2 + 0] /= mu[0]; - dmob[1*2 + 1] /= mu[1]; - } + /// \page tutorial3 + /// \details Computes the mobility and its derivative with respect to saturation + /// for a given cell \c c and saturation \c sat. The template classes \c Sat, + /// \c Mob, \c DMob are typically arrays. By using templates, we do not have to + /// investigate how these array objects are implemented + /// (as long as they have an \c operator[] method). + /// \code + template + void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const; + /// \endcode - template - void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const - { - double pcow[2]; - double dpcow[4]; - props_.capPress(1, &s[0], &c, pcow, dpcow); - pcap = pcow[0]; - dpcap = dpcow[0]; - } + /// \page tutorial3 + /// \details Computes the capillary pressure and its derivative with respect + /// to saturation + /// for a given cell \c c and saturation \c sat. + /// \code + template + void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const; + /// \endcode - double s_min(int c) const - { - return smin_[2*c + 0]; - } + /// \page tutorial3 + /// \details Returns the minimum permitted saturation. + /// \code + double s_min(int c) const; + /// \endcode - double s_max(int c) const - { - return smax_[2*c + 0]; - } + /// \page tutorial3 + /// \details Returns the maximum permitted saturation + /// \code + double s_max(int c) const; + /// \endcode + /// \page tutorial3 + /// \details Private variables + /// \code private: const Opm::IncompPropertiesInterface& props_; std::vector smin_; std::vector smax_; }; +/// \endcode +/// \page tutorial3 +/// \details We set up the transport model. +/// \code typedef Opm::SinglePointUpwindTwoPhase TransportModel; +/// \endcode +/// \page tutorial3 +/// \details +/// The transport equation is nonlinear. We use an implicit transport solver +/// which implements a Newton-Raphson solver. +/// We define the format of the objects +/// which will be used by the solver. +/// \code using namespace Opm::ImplicitTransportDefault; - -typedef NewtonVectorCollection< ::std::vector > NVecColl; -typedef JacobianSystem < struct CSRMatrix, NVecColl > JacSys; +typedef NewtonVectorCollection< ::std::vector > NVecColl; +typedef JacobianSystem< struct CSRMatrix, NVecColl > JacSys; template class MaxNorm { @@ -137,7 +183,12 @@ public: return AccumulationNorm ::norm(v); } }; +/// \endcode +/// \page tutorial3 +/// \details +/// We set up the solver. +/// \code typedef Opm::ImplicitTransport TransportSolver; - +/// \endcode +/// \page tutorial3 +/// \details +/// Main function +/// \code int main () { - + /// \endcode + /// \page tutorial3 + /// \details + /// We define the grid. A cartesian grid with 1200 cells. + /// \code int dim = 3; int nx = 20; int ny = 20; @@ -162,41 +221,88 @@ int main () GridManager grid(nx, ny, nz, dx, dy, dz); int num_cells = grid.c_grid()->number_of_cells; int num_faces = grid.c_grid()->number_of_faces; + /// \endcode + /// \page tutorial3 + /// \details + /// We define the properties of the fluid.\n + /// Number of phases. + /// \code int num_phases = 2; using namespace unit; using namespace prefix; + /// \endcode + /// \page tutorial3 + /// \details density vector (one component per phase). + /// \code std::vector rho(2, 1000.); + /// \endcode + /// \page tutorial3 + /// \details viscosity vector (one component per phase). + /// \code std::vector mu(2, 1.*centi*Poise); + /// \endcode + /// \page tutorial3 + /// \details porosity and permeability of the rock. + /// \code double porosity = 0.5; double k = 10*milli*darcy; - + /// \endcode + + /// \page tutorial3 + /// \details We define the relative permeability function. We use a basic fluid + /// description and set this function to be linear. For more realistic fluid, the + /// saturation function is given by the data. + /// \code SaturationPropsBasic::RelPermFunc rel_perm_func; rel_perm_func = SaturationPropsBasic::Linear; + /// \endcode + /// \page tutorial3 + /// \details We construct a basic fluid with the properties we have defined above. + /// Each property is constant and hold for all cells. + /// \code IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu, porosity, k, dim, num_cells); + TwophaseFluid fluid(props); + /// \endcode - - + + /// \page tutorial3 + /// \details Gravity parameters. Here, we set zero gravity. + /// \code const double *grav = 0; std::vector omega; + /// \endcode + /// \page tutorial3 + /// \details We set up the pressure solver. + /// \code LinearSolverUmfpack linsolver; IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver); + /// \endcode + - TwophaseFluid fluid(props); - + /// \page tutorial3 + /// \details We set up the source term + /// \code std::vector src(num_cells, 0.0); src[0] = 1.; src[num_cells-1] = -1.; + /// \endcode - + /// \page tutorial3 + /// \details We set up the wells. Here, there are no well and we let them empty. + /// \code std::vector empty_wdp; std::vector empty_well_bhp; std::vector empty_well_flux; + /// \endcode + /// \page tutorial3 + /// \details We set up the source term for the transport solver. + /// \code TransportSource* tsrc = create_transport_source(2, 2); double ssrc[] = { 1.0, 0.0 }; double ssink[] = { 0.0, 1.0 }; @@ -208,48 +314,112 @@ int main () append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc); } } - + /// \endcode + + /// \page tutorial3 + /// \details We compute the pore volume + /// \code std::vector porevol; computePorevolume(*grid.c_grid(), props, porevol); + /// \endcode + /// \page tutorial3 + /// \details We set up the transport solver. + /// \code const bool guess_old_solution = true; TransportModel model (fluid, *grid.c_grid(), porevol, grav, guess_old_solution); TransportSolver tsolver(model); + /// \endcode + /// \page tutorial3 + /// \details Time integration parameters + /// \code double dt = 0.1*day; int num_time_steps = 20; + /// \endcode - TwophaseState state; + /// \page tutorial3 + /// \details Control paramaters for the implicit solver. + /// \code ImplicitTransportDetails::NRReport rpt; ImplicitTransportDetails::NRControl ctrl; - std::vector totmob; + /// \endcode + + + /// \page tutorial3 + /// \details We define a vector which contains all cell indexes. We use this + /// vector to set up parameters on the whole domains. std::vector allcells(num_cells); for (int cell = 0; cell < num_cells; ++cell) { allcells[cell] = cell; } - FlowBCManager bcs; - // Linear solver init. + /// \page tutorial3 + /// \details We set up the boundary conditions. Letting bcs empty is equivalent + /// to no flow boundary conditions. + /// \code + FlowBCManager bcs; + /// \endcode + + /// \page tutorial3 + /// \details + /// Linear solver init. + /// \code using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver; CSRMatrixUmfpackSolver linsolve; + /// \endcode + /// \page tutorial3 + /// \details + /// We initialise water saturation to minimum everywhere. + /// \code + TwophaseState state; state.init(*grid.c_grid()); - // By default: initialise water saturation to minimum everywhere. state.setWaterSat(allcells, props, TwophaseState::MinSat); + /// \endcode + /// \page tutorial3 + /// \details We introduce a vector which contains the total mobility + /// on all cells. + /// \code + std::vector totmob; + /// \endcode + /// \page tutorial3 + /// \details This string will contain the name of a VTK output vector. + /// \code std::ostringstream vtkfilename; + /// \endcode + + /// \page tutorial3 + /// \details Loop over the time steps. + /// \code for (int i = 0; i < num_time_steps; ++i) { - + /// \endcode + /// \page tutorial3 + /// \details Compute the total mobility. It is needed by the pressure solver + /// \code computeTotalMobility(props, allcells, state.saturation(), totmob); + /// \endcode + /// \page tutorial3 + /// \details Solve the pressure equation + /// \code psolver.solve(totmob, omega, src, empty_wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), empty_well_bhp, empty_well_flux); + /// \endcode + /// \page tutorial3 + /// \details Transport solver + /// \code tsolver.solve(*grid.c_grid(), tsrc, dt, ctrl, state, linsolve, rpt); + /// \endcode + /// \page tutorial3 + /// \details Write the output to file. + /// \code vtkfilename.str(""); vtkfilename << "tutorial3-" << std::setw(3) << std::setfill('0') << i << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); @@ -257,12 +427,78 @@ int main () dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); Opm::writeVtkData(*grid.c_grid(), dm, vtkfile); - } } +/// \endcode + +/// \page tutorial3 +/// \details Implementation of the TwophaseFluid class. +/// \code +TwophaseFluid::TwophaseFluid(const Opm::IncompPropertiesInterface& props) + : props_(props), + smin_(props.numCells()*props.numPhases()), + smax_(props.numCells()*props.numPhases()) +{ + const int num_cells = props.numCells(); + std::vector cells(num_cells); + for (int c = 0; c < num_cells; ++c) { + cells[c] = c; + } + props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]); + } + +double TwophaseFluid::density(int phase) const +{ + return props_.density()[phase]; +} + +template +void TwophaseFluid::mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const +{ + props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]); + const double* mu = props_.viscosity(); + mob[0] /= mu[0]; + mob[1] /= mu[1]; + /// \endcode + /// \page tutorial3 + /// \details We + /// recall that we use Fortran ordering for kr derivatives, + /// therefore dmob[i*2 + j] is row j and column i of the + /// matrix. + /// Each row corresponds to a kr function, so which mu to + /// divide by also depends on the row, j. + /// \code + dmob[0*2 + 0] /= mu[0]; + dmob[0*2 + 1] /= mu[1]; + dmob[1*2 + 0] /= mu[0]; + dmob[1*2 + 1] /= mu[1]; +} + +template +void TwophaseFluid::pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const +{ + pcap = 0.; + dpcap = 0.; +} + +double TwophaseFluid::s_min(int c) const +{ + return smin_[2*c + 0]; +} + +double TwophaseFluid::s_max(int c) const +{ + return smax_[2*c + 0]; +} +/// \endcode /// \page tutorial3 +/// \section results3 Results. /// /// /// @@ -279,5 +515,6 @@ int main () ///
\image html tutorial3-000.png
/// \page tutorial3 -/// \section sourcecode Complete source code. +/// \details +/// \section completecode3 Complete source code: /// \include tutorial3.cpp From 2c541fee291fc2fe44f2d08651c9ad928f0259e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 17 Apr 2012 19:08:45 +0200 Subject: [PATCH 08/51] Darcy definition: Fix misprint introduced in Doxygenisation. --- opm/core/utility/Units.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/Units.hpp b/opm/core/utility/Units.hpp index a2939783..3d4f4338 100644 --- a/opm/core/utility/Units.hpp +++ b/opm/core/utility/Units.hpp @@ -139,7 +139,7 @@ namespace Opm /// @{ /// /// A porous medium with a permeability of 1 darcy permits a - /// flow (flux) of \f$1\,cm^2/s\f$ of a fluid with viscosity \f$1\,cP\f$ ( + /// flow (flux) of \f$1\,cm^3/s\f$ of a fluid with viscosity \f$1\,cP\f$ ( /// \f$1\,mPa\cdot s\f$) under a pressure gradient of \f$1\,atm/cm\f$ acting across /// an area of \f$1\,cm^2\f$. /// From 80787e57f11f7718e39cb4eee260f6f6dfaa0585 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 17 Apr 2012 19:09:16 +0200 Subject: [PATCH 09/51] Darcy definition: Split long lines to fit in 80 columns. --- opm/core/utility/Units.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/opm/core/utility/Units.hpp b/opm/core/utility/Units.hpp index 3d4f4338..bbf6ce78 100644 --- a/opm/core/utility/Units.hpp +++ b/opm/core/utility/Units.hpp @@ -139,9 +139,10 @@ namespace Opm /// @{ /// /// A porous medium with a permeability of 1 darcy permits a - /// flow (flux) of \f$1\,cm^3/s\f$ of a fluid with viscosity \f$1\,cP\f$ ( - /// \f$1\,mPa\cdot s\f$) under a pressure gradient of \f$1\,atm/cm\f$ acting across - /// an area of \f$1\,cm^2\f$. + /// flow (flux) of \f$1\,cm^3/s\f$ of a fluid with viscosity + /// \f$1\,cP\f$ ( \f$1\,mPa\cdot s\f$) under a pressure + /// gradient of \f$1\,atm/cm\f$ acting across an area of + /// \f$1\,cm^2\f$. /// namespace perm_details { const double p_grad = atm / (prefix::centi*meter); From 41235ed3ad182415794eb9220d1f286c43e6dd07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 17 Apr 2012 19:10:09 +0200 Subject: [PATCH 10/51] Darcy definition: Remove a leading blank inside parentheses. --- opm/core/utility/Units.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/Units.hpp b/opm/core/utility/Units.hpp index bbf6ce78..cb48985c 100644 --- a/opm/core/utility/Units.hpp +++ b/opm/core/utility/Units.hpp @@ -140,7 +140,7 @@ namespace Opm /// /// A porous medium with a permeability of 1 darcy permits a /// flow (flux) of \f$1\,cm^3/s\f$ of a fluid with viscosity - /// \f$1\,cP\f$ ( \f$1\,mPa\cdot s\f$) under a pressure + /// \f$1\,cP\f$ (\f$1\,mPa\cdot s\f$) under a pressure /// gradient of \f$1\,atm/cm\f$ acting across an area of /// \f$1\,cm^2\f$. /// From 8dfcadc36c992b6446fc405f6b19e7d4e76d9eb9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 17 Apr 2012 19:27:29 +0200 Subject: [PATCH 11/51] Units.hpp: Split a few more long lines to fit in 80 columns. --- opm/core/utility/Units.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/core/utility/Units.hpp b/opm/core/utility/Units.hpp index cb48985c..a57ead73 100644 --- a/opm/core/utility/Units.hpp +++ b/opm/core/utility/Units.hpp @@ -51,8 +51,9 @@ namespace Opm namespace unit /// Definition of various units. - /// All the units are defined in terms of international standard units (SI). - /// Example of use: We define a variable \c k which gives a permeability. We want to set \c k to \f$1\,mD\f$. + /// All the units are defined in terms of international standard + /// units (SI). Example of use: We define a variable \c k which + /// gives a permeability. We want to set \c k to \f$1\,mD\f$. /// \code /// using namespace Opm::unit /// double k = 0.001*darcy; From 0e275bea4e48540f7816b6b43ec0c0c5b1bf875c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 17 Apr 2012 19:31:20 +0200 Subject: [PATCH 12/51] Assert copyright for years 2011 and 2012. --- opm/core/utility/Units.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/utility/Units.hpp b/opm/core/utility/Units.hpp index a57ead73..5530df4d 100644 --- a/opm/core/utility/Units.hpp +++ b/opm/core/utility/Units.hpp @@ -13,8 +13,8 @@ //=========================================================================== /* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. + Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010, 2011, 2012 Statoil ASA. This file is part of The Open Reservoir Simulator Project (OpenRS). From 5a57f3147406b3660e2c9324f721572bda076b9d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 17 Apr 2012 19:45:39 +0200 Subject: [PATCH 13/51] Volume: Use official definition of a standard barrel; 42 U.S. gallon. --- opm/core/utility/Units.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/utility/Units.hpp b/opm/core/utility/Units.hpp index 5530df4d..2e8e37b2 100644 --- a/opm/core/utility/Units.hpp +++ b/opm/core/utility/Units.hpp @@ -93,7 +93,8 @@ namespace Opm /// \name Volume /// @{ - const double stb = 0.158987294928 * cubic(meter); + const double gallon = 231 * cubic(inch); + const double stb = 42 * gallon; /// @} /// \name Mass From 9bb4b2c51de16315165466ae53f0ddb2beb857ca Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 18 Apr 2012 14:32:20 +0200 Subject: [PATCH 14/51] Corrected typos. --- tutorials/tutorial1.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/tutorial1.cpp b/tutorials/tutorial1.cpp index 080244fd..e7c0649e 100644 --- a/tutorials/tutorial1.cpp +++ b/tutorials/tutorial1.cpp @@ -47,7 +47,7 @@ int main() /// \page tutorial1 /// We set the number of blocks in each direction. /// \code - int nx = 3; + int nx = 4; int ny = 3; int nz = 2; /// \endcode @@ -82,7 +82,7 @@ int main() } /// \endcode /// \page tutorial1 -/// We read the the vtu output file in \a Paraview and obtain the following grid. +/// We read the vtu output file in \a Paraview and obtain the following grid. /// \image html tutorial1.png /// \page tutorial1 From 5da5b2e325924e27c21375612088f508cdfde3c9 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 18 Apr 2012 14:33:44 +0200 Subject: [PATCH 15/51] Removed warnings from compiler. --- Doxyfile | 5 +++-- tutorials/tutorial3.cpp | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/Doxyfile b/Doxyfile index ed1b6bb6..b4ffcb96 100644 --- a/Doxyfile +++ b/Doxyfile @@ -51,7 +51,7 @@ PROJECT_LOGO = # If a relative path is entered, it will be relative to the location # where doxygen was started. If left blank the current directory will be used. -OUTPUT_DIRECTORY = +OUTPUT_DIRECTORY = ../Documentation # If the CREATE_SUBDIRS tag is set to YES, then doxygen will create # 4096 sub-directories (in 2 levels) under the output directory of each output @@ -610,7 +610,8 @@ WARN_LOGFILE = # directories like "/usr/src/myproject". Separate the files or directories # with spaces. -INPUT = opm/core tutorials examples +# INPUT = opm/core tutorials examples +INPUT = tutorials # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp index 85e8ac23..0ad1d606 100644 --- a/tutorials/tutorial3.cpp +++ b/tutorials/tutorial3.cpp @@ -220,7 +220,6 @@ int main () using namespace Opm; GridManager grid(nx, ny, nz, dx, dy, dz); int num_cells = grid.c_grid()->number_of_cells; - int num_faces = grid.c_grid()->number_of_faces; /// \endcode /// \page tutorial3 @@ -479,7 +478,7 @@ void TwophaseFluid::mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const template -void TwophaseFluid::pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const +void TwophaseFluid::pc(int /*c */, const Sat& /* s*/, Pcap& pcap, DPcap& dpcap) const { pcap = 0.; dpcap = 0.; @@ -514,7 +513,9 @@ double TwophaseFluid::s_max(int c) const /// /// + /// \page tutorial3 /// \details /// \section completecode3 Complete source code: /// \include tutorial3.cpp +/// \include generate_doc_figures.py From d3c0356b93592e9065baf56b4476e2a054518555 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 18 Apr 2012 14:38:04 +0200 Subject: [PATCH 16/51] Revert Doxyfile to correct input file setting. --- Doxyfile | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Doxyfile b/Doxyfile index b4ffcb96..ed1b6bb6 100644 --- a/Doxyfile +++ b/Doxyfile @@ -51,7 +51,7 @@ PROJECT_LOGO = # If a relative path is entered, it will be relative to the location # where doxygen was started. If left blank the current directory will be used. -OUTPUT_DIRECTORY = ../Documentation +OUTPUT_DIRECTORY = # If the CREATE_SUBDIRS tag is set to YES, then doxygen will create # 4096 sub-directories (in 2 levels) under the output directory of each output @@ -610,8 +610,7 @@ WARN_LOGFILE = # directories like "/usr/src/myproject". Separate the files or directories # with spaces. -# INPUT = opm/core tutorials examples -INPUT = tutorials +INPUT = opm/core tutorials examples # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is From 4aab2f18fbab370168410b86182de87621b03658 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 18 Apr 2012 14:46:01 +0200 Subject: [PATCH 17/51] Removed warning from compiler. --- opm/core/fluid/IncompPropertiesBasic.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/core/fluid/IncompPropertiesBasic.cpp b/opm/core/fluid/IncompPropertiesBasic.cpp index 8201dacd..1b93001b 100644 --- a/opm/core/fluid/IncompPropertiesBasic.cpp +++ b/opm/core/fluid/IncompPropertiesBasic.cpp @@ -50,12 +50,12 @@ namespace Opm const SaturationPropsBasic::RelPermFunc& relpermfunc, const std::vector& rho, const std::vector& mu, - const double porosity, - const double permeability, + const double por, //porosity + const double perm, const int dim, const int num_cells) { - rock_.init(dim, num_cells, porosity, permeability); + rock_.init(dim, num_cells, por, perm); pvt_.init(num_phases, rho, mu); satprops_.init(num_phases, relpermfunc); if (pvt_.numPhases() != satprops_.numPhases()) { From 6d04f0042c1f95d916fc1eda790a6215b7f899e1 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 18 Apr 2012 15:17:39 +0200 Subject: [PATCH 18/51] Added support for configurable input/output directories. --- generate_doc_figures.py | 36 ++++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/generate_doc_figures.py b/generate_doc_figures.py index 2829d816..20714318 100755 --- a/generate_doc_figures.py +++ b/generate_doc_figures.py @@ -1,13 +1,21 @@ from subprocess import call from paraview.simple import * from paraview import servermanager -from os import remove +from os import remove, mkdir, curdir +from os.path import join, isdir +figure_path = curdir +tutorial_data_path = curdir +tutorial_path = "tutorials" + +if not isdir(figure_path): + mkdir(figure_path) + connection = servermanager.Connect() # tutorial 1 -call("tutorials/tutorial1") -grid = servermanager.sources.XMLUnstructuredGridReader(FileName="tutorial1.vtu") +call(join(tutorial_path,"tutorial1")) +grid = servermanager.sources.XMLUnstructuredGridReader(FileName=join(tutorial_data_path,"tutorial1.vtu")) grid.UpdatePipeline() Show(grid) dp = GetDisplayProperties(grid) @@ -22,13 +30,13 @@ camera.SetViewUp(-0.19, 0.4, 0.9) camera.SetViewAngle(30) camera.SetFocalPoint(1.5, 1.5, 1) Render() -WriteImage("Figure/tutorial1.png") +WriteImage(join(figure_path,"tutorial1.png")) Hide(grid) -remove("tutorial1.vtu") +remove(join(tutorial_data_path,"tutorial1.vtu")) # tutorial 2 -call("tutorials/tutorial2") -grid = servermanager.sources.XMLUnstructuredGridReader(FileName="tutorial2.vtu") +call(join(tutorial_path,"tutorial2")) +grid = servermanager.sources.XMLUnstructuredGridReader(FileName=join(tutorial_data_path,"tutorial2.vtu")) grid.UpdatePipeline() Show(grid) dp = GetDisplayProperties(grid) @@ -45,15 +53,15 @@ camera.SetViewUp(0, 1, 0) camera.SetViewAngle(30) camera.SetFocalPoint(20, 20, 0.5) Render() -WriteImage("Figure/tutorial2.png") +WriteImage(join(figure_path,"tutorial2.png")) Hide(grid) -# remove("tutorial2.vtu") +remove(join(tutorial_data_path,"tutorial2.vtu")) -# # tutorial 3 -call("tutorials/tutorial3") +# tutorial 3 +call(join(tutorial_path,"tutorial3")) cases = ["000", "005", "010", "015", "019"] for case in cases: - grid = servermanager.sources.XMLUnstructuredGridReader(FileName="tutorial3-"+case+".vtu") + grid = servermanager.sources.XMLUnstructuredGridReader(FileName=join(tutorial_data_path,"tutorial3-"+case+".vtu")) grid.UpdatePipeline() Show(grid) dp = GetDisplayProperties(grid) @@ -69,7 +77,7 @@ for case in cases: camera.SetViewAngle(30) camera.SetFocalPoint(100, 100, 5) Render() - WriteImage("Figure/tutorial3-"+case+".png") + WriteImage(join(figure_path,"tutorial3-"+case+".png")) Hide(grid) -# remove("tutorial3.vtu") +# remove(join(tutorial_data_path,"tutorial3.vtu")) From 944f788b8804d0c3fcb85472ad06a3d26894e128 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 18 Apr 2012 15:42:44 +0200 Subject: [PATCH 19/51] Added collection and removal of temporary files in python script. --- Doxyfile | 7 ++++--- generate_doc_figures.py | 38 +++++++++++++++++++++++++------------- 2 files changed, 29 insertions(+), 16 deletions(-) diff --git a/Doxyfile b/Doxyfile index ed1b6bb6..4b1c43ee 100644 --- a/Doxyfile +++ b/Doxyfile @@ -51,7 +51,7 @@ PROJECT_LOGO = # If a relative path is entered, it will be relative to the location # where doxygen was started. If left blank the current directory will be used. -OUTPUT_DIRECTORY = +OUTPUT_DIRECTORY = ../Documentation # If the CREATE_SUBDIRS tag is set to YES, then doxygen will create # 4096 sub-directories (in 2 levels) under the output directory of each output @@ -610,7 +610,8 @@ WARN_LOGFILE = # directories like "/usr/src/myproject". Separate the files or directories # with spaces. -INPUT = opm/core tutorials examples +#INPUT = opm/core tutorials examples +INPUT = tutorials # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is @@ -688,7 +689,7 @@ EXAMPLE_RECURSIVE = NO # directories that contain image that are included in the documentation (see # the \image command). -IMAGE_PATH = Figure/ +IMAGE_PATH = ../Documentation/Figure/ # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/generate_doc_figures.py b/generate_doc_figures.py index 20714318..37edbf36 100755 --- a/generate_doc_figures.py +++ b/generate_doc_figures.py @@ -4,18 +4,22 @@ from paraview import servermanager from os import remove, mkdir, curdir from os.path import join, isdir -figure_path = curdir +figure_path = "../Documentation/Figure" tutorial_data_path = curdir tutorial_path = "tutorials" +collected_garbage_file = [] + if not isdir(figure_path): mkdir(figure_path) connection = servermanager.Connect() # tutorial 1 -call(join(tutorial_path,"tutorial1")) -grid = servermanager.sources.XMLUnstructuredGridReader(FileName=join(tutorial_data_path,"tutorial1.vtu")) +call(join(tutorial_path, "tutorial1")) +data_file_name = join(tutorial_data_path, "tutorial1.vtu") +grid = servermanager.sources.XMLUnstructuredGridReader(FileName = data_file_name) +collected_garbage_file.append(data_file_name) grid.UpdatePipeline() Show(grid) dp = GetDisplayProperties(grid) @@ -30,13 +34,14 @@ camera.SetViewUp(-0.19, 0.4, 0.9) camera.SetViewAngle(30) camera.SetFocalPoint(1.5, 1.5, 1) Render() -WriteImage(join(figure_path,"tutorial1.png")) +WriteImage(join(figure_path, "tutorial1.png")) Hide(grid) -remove(join(tutorial_data_path,"tutorial1.vtu")) # tutorial 2 -call(join(tutorial_path,"tutorial2")) -grid = servermanager.sources.XMLUnstructuredGridReader(FileName=join(tutorial_data_path,"tutorial2.vtu")) +call(join(tutorial_path, "tutorial2")) +data_file_name = join(tutorial_data_path, "tutorial2.vtu") +grid = servermanager.sources.XMLUnstructuredGridReader(FileName = data_file_name) +collected_garbage_file.append(data_file_name) grid.UpdatePipeline() Show(grid) dp = GetDisplayProperties(grid) @@ -53,15 +58,19 @@ camera.SetViewUp(0, 1, 0) camera.SetViewAngle(30) camera.SetFocalPoint(20, 20, 0.5) Render() -WriteImage(join(figure_path,"tutorial2.png")) +WriteImage(join(figure_path, "tutorial2.png")) Hide(grid) -remove(join(tutorial_data_path,"tutorial2.vtu")) # tutorial 3 -call(join(tutorial_path,"tutorial3")) +call(join(tutorial_path, "tutorial3")) +for case in range(0,20): + data_file_name = join(tutorial_data_path, "tutorial3-"+"%(case)03d"%{"case": case}+".vtu") + collected_garbage_file.append(data_file_name) + cases = ["000", "005", "010", "015", "019"] for case in cases: - grid = servermanager.sources.XMLUnstructuredGridReader(FileName=join(tutorial_data_path,"tutorial3-"+case+".vtu")) + data_file_name = join(tutorial_data_path, "tutorial3-"+case+".vtu") + grid = servermanager.sources.XMLUnstructuredGridReader(FileName = data_file_name) grid.UpdatePipeline() Show(grid) dp = GetDisplayProperties(grid) @@ -77,7 +86,10 @@ for case in cases: camera.SetViewAngle(30) camera.SetFocalPoint(100, 100, 5) Render() - WriteImage(join(figure_path,"tutorial3-"+case+".png")) + WriteImage(join(figure_path, "tutorial3-"+case+".png")) Hide(grid) -# remove(join(tutorial_data_path,"tutorial3.vtu")) + +# remove temporary files +for f in collected_garbage_file: + remove(f) From 0badf481b59b9173d696f7a5ec2741834d648fbc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 11:45:52 +0200 Subject: [PATCH 20/51] Moved UniformTableLinear and related func out of subnamespace utils. --- opm/core/utility/UniformTableLinear.hpp | 3 -- .../utility/buildUniformMonotoneTable.hpp | 34 +++++++++---------- 2 files changed, 16 insertions(+), 21 deletions(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index b7e8d751..058c5f68 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -29,8 +29,6 @@ #include namespace Opm { - namespace utils { - /// @brief This class uses linear interpolation to compute the value /// (and its derivative) of a function f sampled at uniform points. @@ -253,7 +251,6 @@ namespace Opm { return os; } - } // namespace utils } // namespace Opm #endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED diff --git a/opm/core/utility/buildUniformMonotoneTable.hpp b/opm/core/utility/buildUniformMonotoneTable.hpp index bf4ab69e..3ebad7c9 100644 --- a/opm/core/utility/buildUniformMonotoneTable.hpp +++ b/opm/core/utility/buildUniformMonotoneTable.hpp @@ -24,27 +24,25 @@ #include namespace Opm { - namespace utils { - template - void buildUniformMonotoneTable(const std::vector& xv, - const std::vector& yv, - const int samples, - UniformTableLinear& table) - { - MonotCubicInterpolator interp(xv, yv); - std::vector uniform_yv(samples); - double xmin = xv[0]; - double xmax = xv.back(); - for (int i = 0; i < samples; ++i) { - double w = double(i)/double(samples - 1); - double x = (1.0 - w)*xmin + w*xmax; - uniform_yv[i] = interp(x); - } - table = UniformTableLinear(xmin, xmax, uniform_yv); + template + void buildUniformMonotoneTable(const std::vector& xv, + const std::vector& yv, + const int samples, + UniformTableLinear& table) + { + MonotCubicInterpolator interp(xv, yv); + std::vector uniform_yv(samples); + double xmin = xv[0]; + double xmax = xv.back(); + for (int i = 0; i < samples; ++i) { + double w = double(i)/double(samples - 1); + double x = (1.0 - w)*xmin + w*xmax; + uniform_yv[i] = interp(x); } + table = UniformTableLinear(xmin, xmax, uniform_yv); + } - } // namespace utils } // namespace Opm From 74c8cc68e80ee662ca40edbf7c0d384b1f5d2181 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 11:46:41 +0200 Subject: [PATCH 21/51] Make column gravity Gauss-Seidel solver report average number of iterations. --- opm/core/transport/reorder/TransportModelTwophase.cpp | 10 ++++++---- opm/core/transport/reorder/TransportModelTwophase.hpp | 2 +- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index 33a91905..d52becc9 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -549,7 +549,7 @@ namespace Opm - void TransportModelTwophase::solveGravityColumn(const std::vector& cells) + int TransportModelTwophase::solveGravityColumn(const std::vector& cells) { // Set up column gravflux. const int nc = cells.size(); @@ -597,7 +597,7 @@ namespace Opm THROW("In solveGravityColumn(), we did not converge after " << num_iters << " iterations. Delta s = " << max_s_change); } - // Repeat if necessary. + return num_iters + 1; } @@ -629,11 +629,13 @@ namespace Opm saturation_ = &saturation[0]; // Solve on all columns. - + int num_iters = 0; for (std::vector >::size_type i = 0; i < columns.second.size(); i++) { // std::cout << "==== new column" << std::endl; - solveGravityColumn(columns.second[i]); + num_iters += solveGravityColumn(columns.second[i]); } + std::cout << "Gauss-Seidel column solver average iterations: " + << double(num_iters)/double(columns.second.size()) << std::endl; } } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 957a58f2..ba1062d6 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -52,7 +52,7 @@ namespace Opm void solveSingleCellGravity(const std::vector& cells, const int pos, const double* gravflux); - void solveGravityColumn(const std::vector& cells); + int solveGravityColumn(const std::vector& cells); void solveGravity(const std::pair, std::vector > >& columns, const double* porevolume, const double dt, From 0b09b8b5fdc78b52b928aa024aec8b52fcdb588e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 11:47:41 +0200 Subject: [PATCH 22/51] Make constructor take arguments by reference. Silence a warning. --- opm/core/fluid/PvtPropertiesBasic.cpp | 6 +++--- opm/core/fluid/PvtPropertiesBasic.hpp | 7 ++++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/opm/core/fluid/PvtPropertiesBasic.cpp b/opm/core/fluid/PvtPropertiesBasic.cpp index 53b6cbfd..3d7d7bb2 100644 --- a/opm/core/fluid/PvtPropertiesBasic.cpp +++ b/opm/core/fluid/PvtPropertiesBasic.cpp @@ -60,8 +60,8 @@ namespace Opm } void PvtPropertiesBasic::init(const int num_phases, - const std::vector rho, - const std::vector mu) + const std::vector& rho, + const std::vector& visc) { if (num_phases > 3 || num_phases < 1) { THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); @@ -70,7 +70,7 @@ namespace Opm formation_volume_factor_.clear(); formation_volume_factor_.resize(num_phases, 1.0); density_ = rho; - viscosity_ = mu; + viscosity_ = visc; } const double* PvtPropertiesBasic::surfaceDensities() const diff --git a/opm/core/fluid/PvtPropertiesBasic.hpp b/opm/core/fluid/PvtPropertiesBasic.hpp index 7bf234a7..d1a3e9f3 100644 --- a/opm/core/fluid/PvtPropertiesBasic.hpp +++ b/opm/core/fluid/PvtPropertiesBasic.hpp @@ -44,10 +44,11 @@ namespace Opm /// mu1 [mu2, mu3] (1.0) Viscosity in cP void init(const parameter::ParameterGroup& param); - /// Initialize from argument Basic multi phase fluid pvt properties. + /// Initialize from arguments. + /// Basic multi phase fluid pvt properties. void init(const int num_phases, - const std::vector rho, - const std::vector mu); + const std::vector& rho, + const std::vector& visc); /// Number of active phases. int numPhases() const; From 389a3d79dc021fc210a60c1097735939e72501fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 11:48:54 +0200 Subject: [PATCH 23/51] Follow previous removal of utils:: subnamespace. --- opm/core/fluid/blackoil/SinglePvtDead.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/fluid/blackoil/SinglePvtDead.hpp b/opm/core/fluid/blackoil/SinglePvtDead.hpp index 3a34a56c..137e2e30 100644 --- a/opm/core/fluid/blackoil/SinglePvtDead.hpp +++ b/opm/core/fluid/blackoil/SinglePvtDead.hpp @@ -74,8 +74,8 @@ namespace Opm double* output_dRdp) const; private: // PVT properties of dry gas or dead oil - utils::UniformTableLinear one_over_B_; - utils::UniformTableLinear viscosity_; + UniformTableLinear one_over_B_; + UniformTableLinear viscosity_; }; } From a41ab2a4a33acdbd03c866460afc8e61736ae741 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 11:49:59 +0200 Subject: [PATCH 24/51] Now saturation props read from deck may have multiple tables, and support SATNUM. --- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 8 +- opm/core/fluid/IncompPropertiesFromDeck.cpp | 12 +- opm/core/fluid/SaturationPropsFromDeck.cpp | 205 +++++++++++------- opm/core/fluid/SaturationPropsFromDeck.hpp | 42 ++-- 4 files changed, 168 insertions(+), 99 deletions(-) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 301a3974..dc90e0cf 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -199,11 +199,11 @@ namespace Opm /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void BlackoilPropertiesFromDeck::relperm(const int n, const double* s, - const int* /*cells*/, + const int* cells, double* kr, double* dkrds) const { - satprops_.relperm(n, s, kr, dkrds); + satprops_.relperm(n, s, cells, kr, dkrds); } @@ -218,11 +218,11 @@ namespace Opm /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void BlackoilPropertiesFromDeck::capPress(const int n, const double* s, - const int* /*cells*/, + const int* cells, double* pc, double* dpcds) const { - satprops_.capPress(n, s, pc, dpcds); + satprops_.capPress(n, s, cells, pc, dpcds); } diff --git a/opm/core/fluid/IncompPropertiesFromDeck.cpp b/opm/core/fluid/IncompPropertiesFromDeck.cpp index 519615ea..53f2370a 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.cpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.cpp @@ -101,11 +101,11 @@ namespace Opm /// 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*/, + const int* cells, double* kr, double* dkrds) const { - satprops_.relperm(n, s, kr, dkrds); + satprops_.relperm(n, s, cells, kr, dkrds); } @@ -120,11 +120,11 @@ namespace Opm /// 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*/, + const int* cells, double* pc, double* dpcds) const { - satprops_.capPress(n, s, pc, dpcds); + satprops_.capPress(n, s, cells, pc, dpcds); } @@ -136,11 +136,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*/, + const int* cells, double* smin, double* smax) const { - satprops_.satRange(n, smin, smax); + satprops_.satRange(n, cells, smin, smax); } } // namespace Opm diff --git a/opm/core/fluid/SaturationPropsFromDeck.cpp b/opm/core/fluid/SaturationPropsFromDeck.cpp index 8b614db1..f7c14abd 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.cpp +++ b/opm/core/fluid/SaturationPropsFromDeck.cpp @@ -41,49 +41,46 @@ namespace Opm if (!phase_usage_.phase_used[Liquid]) { THROW("SaturationPropsFromDeck::init() -- oil phase must be active."); } - const int samples = 200; - double swco = 0.0; - double swmax = 1.0; + + // Obtain SATNUM, if it exists. + // Otherwise, let the cell_to_func_ mapping be just empty. + int satfuncs_expected = 1; + if (deck.hasField("SATNUM")) { + cell_to_func_ = deck.getIntegerValue("SATNUM"); + satfuncs_expected = *std::max_element(cell_to_func_.begin(), cell_to_func_.end()); + // Must subtract 1 to get numbering from 0. + std::transform(cell_to_func_.begin(), cell_to_func_.end(), + cell_to_func_.begin(), std::bind2nd(std::plus(), -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_; - if (swof_table.size() != 1) { - THROW("We must have exactly one SWOF table."); + num_tables = swof_table.size(); + if (num_tables < satfuncs_expected) { + THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected); } - const std::vector& sw = swof_table[0][0]; - const std::vector& krw = swof_table[0][1]; - const std::vector& krow = swof_table[0][2]; - const std::vector& pcow = swof_table[0][3]; - buildUniformMonotoneTable(sw, krw, samples, krw_); - buildUniformMonotoneTable(sw, krow, samples, krow_); - buildUniformMonotoneTable(sw, pcow, samples, pcow_); - krocw_ = krow[0]; // At connate water -> ecl. SWOF - swco = sw[0]; - smin_[phase_usage_.phase_pos[Aqua]] = sw[0]; - swmax = sw.back(); - smax_[phase_usage_.phase_pos[Aqua]] = sw.back(); } if (phase_usage_.phase_used[Vapour]) { const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; - if (sgof_table.size() != 1) { - THROW("We must have exactly one SGOF table."); + 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."); } - const std::vector& sg = sgof_table[0][0]; - const std::vector& krg = sgof_table[0][1]; - const std::vector& krog = sgof_table[0][2]; - const std::vector& pcog = sgof_table[0][3]; - buildUniformMonotoneTable(sg, krg, samples, krg_); - buildUniformMonotoneTable(sg, krog, samples, krog_); - buildUniformMonotoneTable(sg, pcog, samples, pcog_); - smin_[phase_usage_.phase_pos[Vapour]] = sg[0]; - if (std::fabs(sg.back() + swco - 1.0) > 1e-2) { - THROW("Gas maximum saturation in SGOF table = " << sg.back() << - ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); - } - smax_[phase_usage_.phase_pos[Vapour]] = sg.back(); } - // These only consider water min/max sats. Consider gas sats? - smin_[phase_usage_.phase_pos[Liquid]] = 1.0 - swmax; - smax_[phase_usage_.phase_pos[Liquid]] = 1.0 - swco; + + // Initialize tables. + satfuncset_.resize(num_tables); + for (int table = 0; table < num_tables; ++table) { + satfuncset_[table].init(deck, table, phase_usage_); + } } @@ -101,6 +98,7 @@ namespace Opm /// 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. @@ -109,6 +107,7 @@ namespace Opm /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void SaturationPropsFromDeck::relperm(const int n, const double* s, + const int* cells, double* kr, double* dkrds) const { @@ -116,12 +115,12 @@ namespace Opm if (dkrds) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*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) { - evalKr(s + np*i, kr + np*i); + funcForCell(cells[i]).evalKr(s + np*i, kr + np*i); } } } @@ -132,6 +131,7 @@ namespace Opm /// 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. @@ -140,6 +140,7 @@ namespace Opm /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void SaturationPropsFromDeck::capPress(const int n, const double* s, + const int* cells, double* pc, double* dpcds) const { @@ -147,12 +148,12 @@ namespace Opm if (dpcds) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*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) { - evalPc(s + np*i, pc + np*i); + funcForCell(cells[i]).evalPc(s + np*i, pc + np*i); } } } @@ -162,30 +163,84 @@ namespace Opm /// 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. void SaturationPropsFromDeck::satRange(const int n, + const int* cells, double* smin, double* smax) const { 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] = smin_[p]; - smax[np*i + p] = smax_[p]; + smin[np*i + p] = funcForCell(cells[i]).smin_[p]; + smax[np*i + p] = funcForCell(cells[i]).smax_[p]; } } } - - - // Private methods below. - - - void SaturationPropsFromDeck::evalKr(const double* s, double* kr) const + // Map the cell number to the correct function set. + const SaturationPropsFromDeck::SatFuncSet& + SaturationPropsFromDeck::funcForCell(const int cell) const { - if (phase_usage_.num_phases == 3) { + return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]]; + } + + + + // ----------- Methods of SatFuncSet below ----------- + + + void SaturationPropsFromDeck::SatFuncSet::init(const EclipseGridParser& deck, + const int table_num, + const PhaseUsage phase_usg) + { + phase_usage = phase_usg; + const int samples = 200; + double swco = 0.0; + double swmax = 1.0; + if (phase_usage.phase_used[Aqua]) { + const SWOF::table_t& swof_table = deck.getSWOF().swof_; + const std::vector& sw = swof_table[table_num][0]; + const std::vector& krw = swof_table[table_num][1]; + const std::vector& krow = swof_table[table_num][2]; + const std::vector& pcow = swof_table[table_num][3]; + buildUniformMonotoneTable(sw, krw, samples, krw_); + buildUniformMonotoneTable(sw, krow, samples, krow_); + buildUniformMonotoneTable(sw, pcow, samples, pcow_); + krocw_ = krow[0]; // At connate water -> ecl. SWOF + swco = sw[0]; + smin_[phase_usage.phase_pos[Aqua]] = sw[0]; + swmax = sw.back(); + smax_[phase_usage.phase_pos[Aqua]] = sw.back(); + } + if (phase_usage.phase_used[Vapour]) { + const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; + const std::vector& sg = sgof_table[table_num][0]; + const std::vector& krg = sgof_table[table_num][1]; + const std::vector& krog = sgof_table[table_num][2]; + const std::vector& pcog = sgof_table[table_num][3]; + buildUniformMonotoneTable(sg, krg, samples, krg_); + buildUniformMonotoneTable(sg, krog, samples, krog_); + buildUniformMonotoneTable(sg, pcog, samples, pcog_); + smin_[phase_usage.phase_pos[Vapour]] = sg[0]; + if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { + THROW("Gas maximum saturation in SGOF table = " << sg.back() << + ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); + } + smax_[phase_usage.phase_pos[Vapour]] = sg.back(); + } + // These only consider water min/max sats. Consider gas sats? + smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; + } + + + void SaturationPropsFromDeck::SatFuncSet::evalKr(const double* s, double* kr) const + { + if (phase_usage.num_phases == 3) { // Stone-II relative permeability model. double sw = s[Aqua]; double sg = s[Vapour]; @@ -203,18 +258,18 @@ namespace Opm return; } // We have a two-phase situation. We know that oil is active. - if (phase_usage_.phase_used[Aqua]) { - int wpos = phase_usage_.phase_pos[Aqua]; - int opos = phase_usage_.phase_pos[Liquid]; + if (phase_usage.phase_used[Aqua]) { + int wpos = phase_usage.phase_pos[Aqua]; + int opos = phase_usage.phase_pos[Liquid]; double sw = s[wpos]; double krw = krw_(sw); double krow = krow_(sw); kr[wpos] = krw; kr[opos] = krow; } else { - ASSERT(phase_usage_.phase_used[Vapour]); - int gpos = phase_usage_.phase_pos[Vapour]; - int opos = phase_usage_.phase_pos[Liquid]; + ASSERT(phase_usage.phase_used[Vapour]); + int gpos = phase_usage.phase_pos[Vapour]; + int opos = phase_usage.phase_pos[Liquid]; double sg = s[gpos]; double krg = krg_(sg); double krog = krog_(sg); @@ -224,9 +279,9 @@ namespace Opm } - void SaturationPropsFromDeck::evalKrDeriv(const double* s, double* kr, double* dkrds) const + void SaturationPropsFromDeck::SatFuncSet::evalKrDeriv(const double* s, double* kr, double* dkrds) const { - const int np = phase_usage_.num_phases; + const int np = phase_usage.num_phases; std::fill(dkrds, dkrds + np*np, 0.0); if (np == 3) { @@ -256,9 +311,9 @@ namespace Opm return; } // We have a two-phase situation. We know that oil is active. - if (phase_usage_.phase_used[Aqua]) { - int wpos = phase_usage_.phase_pos[Aqua]; - int opos = phase_usage_.phase_pos[Liquid]; + if (phase_usage.phase_used[Aqua]) { + int wpos = phase_usage.phase_pos[Aqua]; + int opos = phase_usage.phase_pos[Liquid]; double sw = s[wpos]; double krw = krw_(sw); double dkrww = krw_.derivative(sw); @@ -269,9 +324,9 @@ namespace Opm dkrds[wpos + wpos*np] = dkrww; dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. } else { - ASSERT(phase_usage_.phase_used[Vapour]); - int gpos = phase_usage_.phase_pos[Vapour]; - int opos = phase_usage_.phase_pos[Liquid]; + ASSERT(phase_usage.phase_used[Vapour]); + int gpos = phase_usage.phase_pos[Vapour]; + int opos = phase_usage.phase_pos[Liquid]; double sg = s[gpos]; double krg = krg_(sg); double dkrgg = krg_.derivative(sg); @@ -286,20 +341,20 @@ namespace Opm } - void SaturationPropsFromDeck::evalPc(const double* s, double* pc) const + void SaturationPropsFromDeck::SatFuncSet::evalPc(const double* s, double* pc) const { - pc[phase_usage_.phase_pos[Liquid]] = 0.0; - if (phase_usage_.phase_used[Aqua]) { - int pos = phase_usage_.phase_pos[Aqua]; + pc[phase_usage.phase_pos[Liquid]] = 0.0; + if (phase_usage.phase_used[Aqua]) { + int pos = phase_usage.phase_pos[Aqua]; pc[pos] = pcow_(s[pos]); } - if (phase_usage_.phase_used[Vapour]) { - int pos = phase_usage_.phase_pos[Vapour]; + if (phase_usage.phase_used[Vapour]) { + int pos = phase_usage.phase_pos[Vapour]; pc[pos] = pcog_(s[pos]); } } - void SaturationPropsFromDeck::evalPcDeriv(const double* s, double* pc, double* dpcds) const + void SaturationPropsFromDeck::SatFuncSet::evalPcDeriv(const double* s, double* pc, double* dpcds) const { // The problem of determining three-phase capillary pressures // is very hard experimentally, usually one extends two-phase @@ -307,16 +362,16 @@ namespace Opm // In our approach the derivative matrix is quite sparse, only // the diagonal elements corresponding to non-oil phases are // (potentially) nonzero. - const int np = phase_usage_.num_phases; + const int np = phase_usage.num_phases; std::fill(dpcds, dpcds + np*np, 0.0); - pc[phase_usage_.phase_pos[Liquid]] = 0.0; - if (phase_usage_.phase_used[Aqua]) { - int pos = phase_usage_.phase_pos[Aqua]; + pc[phase_usage.phase_pos[Liquid]] = 0.0; + if (phase_usage.phase_used[Aqua]) { + int pos = phase_usage.phase_pos[Aqua]; pc[pos] = pcow_(s[pos]); dpcds[np*pos + pos] = pcow_.derivative(s[pos]); } - if (phase_usage_.phase_used[Vapour]) { - int pos = phase_usage_.phase_pos[Vapour]; + if (phase_usage.phase_used[Vapour]) { + int pos = phase_usage.phase_pos[Vapour]; pc[pos] = pcog_(s[pos]); dpcds[np*pos + pos] = pcog_.derivative(s[pos]); } diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index 0e5d3b2d..a43e14af 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -23,6 +23,7 @@ #include #include #include +#include namespace Opm { @@ -50,6 +51,7 @@ namespace Opm /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void relperm(const int n, const double* s, + const int* cells, double* kr, double* dkrds) const; @@ -64,6 +66,7 @@ namespace Opm /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void capPress(const int n, const double* s, + const int* cells, double* pc, double* dpcds) const; @@ -72,25 +75,36 @@ 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 satRange(const int n, + const int* cells, double* smin, double* smax) const; - private: - void evalKr(const double* s, double* kr) const; - void evalKrDeriv(const double* s, double* kr, double* dkrds) const; - void evalPc(const double* s, double* pc) const; - void evalPcDeriv(const double* s, double* pc, double* dpcds) const; private: PhaseUsage phase_usage_; - utils::UniformTableLinear krw_; - utils::UniformTableLinear krow_; - utils::UniformTableLinear pcow_; - utils::UniformTableLinear krg_; - utils::UniformTableLinear krog_; - utils::UniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; + class SatFuncSet + { + public: + void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); + void evalKr(const double* s, double* kr) const; + void evalKrDeriv(const double* s, double* kr, double* dkrds) const; + void evalPc(const double* s, double* pc) const; + void evalPcDeriv(const double* s, double* pc, double* dpcds) const; + double smin_[PhaseUsage::MaxNumPhases]; + double smax_[PhaseUsage::MaxNumPhases]; + private: + PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. + UniformTableLinear krw_; + UniformTableLinear krow_; + UniformTableLinear pcow_; + UniformTableLinear krg_; + UniformTableLinear krog_; + UniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; + std::vector satfuncset_; + std::vector cell_to_func_; // = SATNUM - 1 + + const SatFuncSet& funcForCell(const int cell) const; }; From 7b382a82e8dc8b646a7413ea7f01fe7c65813647 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 14:11:00 +0200 Subject: [PATCH 25/51] Reinject UniformTableLinear into utils subnamespace for backwards compatibility. --- opm/core/utility/UniformTableLinear.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 058c5f68..11c2d4ca 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -251,6 +251,12 @@ namespace Opm { return os; } + namespace utils + { + using Opm::UniformTableLinear; + } + + } // namespace Opm #endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED From b9b01326c639c64aa2078e74fc06374eb96fa4ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 15:00:56 +0200 Subject: [PATCH 26/51] Bugfix: in SATNUM treatment, allow for inactive cells. --- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 2 +- opm/core/fluid/IncompPropertiesFromDeck.cpp | 2 +- opm/core/fluid/SaturationPropsFromDeck.cpp | 17 ++++++++++------- opm/core/fluid/SaturationPropsFromDeck.hpp | 4 +++- 4 files changed, 15 insertions(+), 10 deletions(-) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index dc90e0cf..82063b69 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -27,7 +27,7 @@ namespace Opm { rock_.init(deck, global_cell); pvt_.init(deck); - satprops_.init(deck); + satprops_.init(deck, global_cell); if (pvt_.numPhases() != satprops_.numPhases()) { THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); diff --git a/opm/core/fluid/IncompPropertiesFromDeck.cpp b/opm/core/fluid/IncompPropertiesFromDeck.cpp index 53f2370a..9febc085 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.cpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.cpp @@ -31,7 +31,7 @@ namespace Opm { rock_.init(deck, global_cell); pvt_.init(deck); - satprops_.init(deck); + satprops_.init(deck, global_cell); if (pvt_.numPhases() != satprops_.numPhases()) { THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); diff --git a/opm/core/fluid/SaturationPropsFromDeck.cpp b/opm/core/fluid/SaturationPropsFromDeck.cpp index f7c14abd..fd75999e 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.cpp +++ b/opm/core/fluid/SaturationPropsFromDeck.cpp @@ -32,7 +32,8 @@ namespace Opm } /// Initialize from deck. - void SaturationPropsFromDeck::init(const EclipseGridParser& deck) + void SaturationPropsFromDeck::init(const EclipseGridParser& deck, + const std::vector& global_cell) { phase_usage_ = phaseUsageFromDeck(deck); @@ -42,15 +43,17 @@ namespace Opm THROW("SaturationPropsFromDeck::init() -- oil phase must be active."); } - // Obtain SATNUM, if it exists. + // 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")) { - cell_to_func_ = deck.getIntegerValue("SATNUM"); - satfuncs_expected = *std::max_element(cell_to_func_.begin(), cell_to_func_.end()); - // Must subtract 1 to get numbering from 0. - std::transform(cell_to_func_.begin(), cell_to_func_.end(), - cell_to_func_.begin(), std::bind2nd(std::plus(), -1)); + const std::vector& satnum = deck.getIntegerValue("SATNUM"); + satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); + int num_cells = global_cell.size(); + cell_to_func_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + cell_to_func_[cell] = satnum[global_cell[cell]] - 1; + } } // Find number of tables, check for consistency. diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index a43e14af..7c8be8d4 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -35,7 +35,9 @@ namespace Opm SaturationPropsFromDeck(); /// Initialize from deck. - void init(const EclipseGridParser& deck); + /// global_cell maps from grid cells to their original logical Cartesian indices. + void init(const EclipseGridParser& deck, + const std::vector& global_cell); /// \return P, the number of phases. int numPhases() const; From 0c3c195220218205e35168540dc0e19e8271a743 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 15:33:25 +0200 Subject: [PATCH 27/51] Pretty-formatting only. --- opm/core/pressure/tpfa/cfs_tpfa_residual.h | 46 +++++++++++----------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.h b/opm/core/pressure/tpfa/cfs_tpfa_residual.h index 3fdc38f9..8fdc8c5f 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.h +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.h @@ -54,15 +54,15 @@ struct cfs_tpfa_res_data { struct cfs_tpfa_res_data * -cfs_tpfa_res_construct(struct UnstructuredGrid *G , - struct WellCompletions *wconn , - int nphases); +cfs_tpfa_res_construct(struct UnstructuredGrid *G , + struct WellCompletions *wconn , + int nphases); void cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h); void -cfs_tpfa_res_assemble(struct UnstructuredGrid *G, +cfs_tpfa_res_assemble(struct UnstructuredGrid *G, double dt, struct cfs_tpfa_res_forces *forces, const double *zc, @@ -75,7 +75,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G, struct cfs_tpfa_res_data *h); void -cfs_tpfa_res_flux(struct UnstructuredGrid *G , +cfs_tpfa_res_flux(struct UnstructuredGrid *G , struct cfs_tpfa_res_forces *forces , int np , const double *trans , @@ -88,7 +88,7 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G , double *wflux ); void -cfs_tpfa_res_fpress(struct UnstructuredGrid *G, +cfs_tpfa_res_fpress(struct UnstructuredGrid *G, int np, const double *htrans, const double *pmobf, @@ -100,19 +100,19 @@ cfs_tpfa_res_fpress(struct UnstructuredGrid *G, #if 0 void -cfs_tpfa_retrieve_masstrans(struct UnstructuredGrid *G, - int np, - struct cfs_tpfa_data *h, - double *masstrans_f); +cfs_tpfa_retrieve_masstrans(struct UnstructuredGrid *G, + int np, + struct cfs_tpfa_data *h, + double *masstrans_f); void -cfs_tpfa_retrieve_gravtrans(struct UnstructuredGrid *G, - int np, - struct cfs_tpfa_data *h, - double *gravtrans_f); +cfs_tpfa_retrieve_gravtrans(struct UnstructuredGrid *G, + int np, + struct cfs_tpfa_data *h, + double *gravtrans_f); double -cfs_tpfa_impes_maxtime(struct UnstructuredGrid *G, +cfs_tpfa_impes_maxtime(struct UnstructuredGrid *G, struct compr_quantities *cq, const double *trans, const double *porevol, @@ -122,14 +122,14 @@ cfs_tpfa_impes_maxtime(struct UnstructuredGrid *G, const double *gravity); void -cfs_tpfa_expl_mass_transport(struct UnstructuredGrid *G, - well_t *W, - struct completion_data *wdata, - int np, - double dt, - const double *porevol, - struct cfs_tpfa_data *h, - double *surf_vol); +cfs_tpfa_expl_mass_transport(struct UnstructuredGrid *G, + well_t *W, + struct completion_data *wdata, + int np, + double dt, + const double *porevol, + struct cfs_tpfa_data *h, + double *surf_vol); #endif From db6b57f76a274a4977bde23c7611335515a7158e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 16:11:28 +0200 Subject: [PATCH 28/51] Silence warnings about unhandled cases. --- opm/core/WellsManager.cpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index f4838598..2e6eb3b6 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -539,9 +539,22 @@ namespace Opm well_data[i].target = guide_rate * parent_prod_spec.liquid_max_rate_; well_data[i].control = RATE; break; + case ProductionSpecification::NONE_CM: + case ProductionSpecification::ORAT: + case ProductionSpecification::WRAT: + case ProductionSpecification::REIN: + case ProductionSpecification::RESV: + case ProductionSpecification::VREP: + case ProductionSpecification::WGRA: + case ProductionSpecification::FLD: + case ProductionSpecification::GRUP: + THROW("Unhandled production specification control mode " << parent_prod_spec.control_mode_); + break; } - - } + } + case ProductionSpecification::RAT: + THROW("Unhandled production specification guide rate type " << ProductionSpecification::RAT); + break; } } From a085337dc05b2bdf9af1c9d8a65848094d3951f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 16:11:58 +0200 Subject: [PATCH 29/51] Initialise current well control in the appropriate place. --- opm/core/WellsManager.cpp | 1 + opm/core/utility/newwells.c | 3 --- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index 2e6eb3b6..dd386d69 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -598,6 +598,7 @@ namespace Opm // We only append a single control at this point. // TODO: Handle multiple controls. ok = well_controls_append(well_data[w].control, well_data[w].target, w_->ctrls[w]); + w_->ctrls[w]->current = 0; if (!ok) { THROW("Failed to add well controls."); } diff --git a/opm/core/utility/newwells.c b/opm/core/utility/newwells.c index bac7cc25..4e839505 100644 --- a/opm/core/utility/newwells.c +++ b/opm/core/utility/newwells.c @@ -407,9 +407,6 @@ well_controls_append(enum control_type type , ctrl->target[ctrl->num] = target; ctrl->num += 1; - - /* TODO: Review this: */ - ctrl->current = 0; } return ok; From 78006cff402ffd8483255d363f021cf9c84e696e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Apr 2012 09:50:36 +0200 Subject: [PATCH 30/51] Renamed some enums and functions relating to Wells, and started documenting it. --- opm/core/InjectionSpecification.hpp | 2 +- opm/core/WellsGroup.cpp | 4 +- opm/core/WellsManager.cpp | 22 +++--- opm/core/newwells.h | 108 +++++++++++++++++++--------- opm/core/utility/newwells.c | 16 ++--- 5 files changed, 97 insertions(+), 55 deletions(-) diff --git a/opm/core/InjectionSpecification.hpp b/opm/core/InjectionSpecification.hpp index 3cd164a8..772815b3 100644 --- a/opm/core/InjectionSpecification.hpp +++ b/opm/core/InjectionSpecification.hpp @@ -15,7 +15,7 @@ namespace Opm InjectionSpecification(); - surface_component injector_type_; + SurfaceComponent injector_type_; ControlMode control_mode_; double surface_flow_max_rate_; double reinjection_fraction_target_; diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 6bf82b39..18bcea2e 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -287,7 +287,7 @@ namespace Opm namespace { - surface_component toSurfaceComponent(std::string type) + SurfaceComponent toSurfaceComponent(std::string type) { if (type == "OIL") { return OIL; @@ -298,7 +298,7 @@ namespace Opm if (type == "GAS") { return GAS; } - THROW("Unknown type " << type << ", could not convert to surface_component"); + THROW("Unknown type " << type << ", could not convert to SurfaceComponent"); } InjectionSpecification::ControlMode toInjectionControlMode(std::string type) diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index dd386d69..f0933649 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -36,11 +36,11 @@ namespace struct WellData { - well_type type; - control_type control; + WellType type; + WellControlType control; double target; double reference_bhp_depth; - surface_component injected_phase; + SurfaceComponent injected_phase; }; @@ -573,7 +573,7 @@ namespace Opm // Set up the Wells struct. - w_ = wells_create(num_wells, num_perfs); + w_ = create_wells(num_wells, num_perfs); if (!w_) { THROW("Failed creating Wells struct."); } @@ -590,28 +590,26 @@ namespace Opm } const double* zfrac = (well_data[w].type == INJECTOR) ? fracs[well_data[w].injected_phase] : 0; - int ok = wells_add(well_data[w].type, well_data[w].reference_bhp_depth, nperf, + int ok = add_wells(well_data[w].type, well_data[w].reference_bhp_depth, nperf, zfrac, &cells[0], &wi[0], w_); if (!ok) { THROW("Failed to add a well."); } // We only append a single control at this point. // TODO: Handle multiple controls. - ok = well_controls_append(well_data[w].control, well_data[w].target, w_->ctrls[w]); + ok = append_well_controls(well_data[w].control, well_data[w].target, w_->ctrls[w]); w_->ctrls[w]->current = 0; if (!ok) { THROW("Failed to add well controls."); } } - - for(size_t i = 0; i < well_collection_.getLeafNodes().size(); i++) { + + // \TODO comment this. + for (size_t i = 0; i < well_collection_.getLeafNodes().size(); i++) { WellNode* node = static_cast(well_collection_.getLeafNodes()[i].get()); - // We know that getLeafNodes() is ordered the same way as they're indexed in w_ node->setWellsPointer(w_, i); } - - } @@ -619,7 +617,7 @@ namespace Opm /// Destructor. WellsManager::~WellsManager() { - wells_destroy(w_); + destroy_wells(w_); } diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 425eefd2..4be7ef4c 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -24,48 +24,92 @@ extern "C" { #endif -enum well_type { INJECTOR, PRODUCER }; -enum control_type { BHP , RATE }; -enum surface_component { WATER = 0, OIL = 1, GAS = 2 }; +/** Well type indicates desired/expected well behaviour. */ +enum WellType { INJECTOR, PRODUCER }; +/** Type of well control equation or inequality. + * BHP -> bottom hole pressure is specified. + * RATE -> flow rate is specified. + */ +enum WellControlType { BHP , RATE }; +/** Canonical component names and ordering. */ +enum SurfaceComponent { WATER = 0, OIL = 1, GAS = 2 }; -/* Control for single well */ -struct WellControls { - int num; - int cpty; - enum control_type *type; - double *target; - int current; +/** Controls for a single well. + * Each control specifies a well rate or bottom hole pressure. Only + * one control can be active at a time, indicated by current. The + * meaning of each control's target value depends on the control + * type, for BHP controls it is a pressure in Pascal, for RATE + * controls it is a volumetric rate in cubic(meter)/second. The + * active control should be interpreted as an equation, whereas the + * non-active controls should be interpreted as inequalities + * specifying constraints on the solution, where BHP controls yield + * minimum pressures, and RATE controls yield maximum rates. + */ +struct WellControls +{ + int num; /** Number of controls. */ + int cpty; /** Allocated capacity, for internal use only. */ + enum WellControlType *type; /** Array of control types. */ + double *target; /** Array of control targets. */ + int current; /** Index of current active control. */ }; -struct Wells { - int number_of_wells; - int well_cpty; - int perf_cpty; - enum well_type *type; - double *depth_ref; - double *zfrac; /* Surface volume fraction - * (3*number_of_wells) */ - int *well_connpos; - int *well_cells; - double *WI; /* Well index */ - - struct WellControls **ctrls; /* One struct for each well */ +/** Struct encapsulating static information about all wells in a scenario. */ +struct Wells +{ + int number_of_wells; /** Number of wells. */ + int well_cpty; /** Allocated well capacity, for internal use only. */ + int perf_cpty; /** Allocated perforation capacity, for internal use only. */ + enum WellType *type; /** Array of well types. */ + double *depth_ref; /** Array of well bhp reference depths. */ + double *zfrac; /** Component volume fractions for each well, size is (3*number_of_wells). + * This is intended to be used for injection wells. For production wells + * the component fractions will vary and cannot be specified a priori. + */ + int *well_connpos; /** Array of indices into well_cells (and WI). + * For a well w, well_connpos[w] and well_connpos[w+1] yield + * start and one-beyond-end indices into the well_cells array + * for accessing w's perforation cell indices. + */ + int *well_cells; /** Array of perforation cell indices. + * Size is number of perforations (== well_connpos[number_of_wells]). + */ + double *WI; /** Well productivity index, same size and structure as well_cells. */ + struct WellControls **ctrls; /** Well controls, one struct for each well. */ }; -struct CompletionData { - double *gpot; /* Gravity potential */ - double *A; /* RB^{-1} */ - double *phasemob; /* Phase mobility */ + +/** Struct encapsulating dynamic information about all wells in a scenario. + * All arrays in this struct contain data for each perforation, ordered + * the same as Wells::well_cells and Wells:WI. Letting NP be the number + * of perforations, the array sizes are: + * gpot 3*NP + * A 9*NP (matrix in Fortran order). + * phasemob 3*NP + * \TODO: Verify that the sizes are correct, check if we should refactor to handle two phases better. + */ +struct CompletionData +{ + double *gpot; /** Gravity potentials. */ + double *A; /** Volumes to surface-components matrix, A = RB^{-1}. */ + double *phasemob; /** Phase mobilities. */ }; +/** Contruction function initializing a Wells object. + * The arguments may be used to indicate expected capacity needed, + * they will be used internally for pre-allocation. + * \return NULL upon failure, otherwise a valid Wells object with 0 wells. + */ struct Wells * -wells_create(int nwells, int nperf); +create_wells(int nwells_reserve_cap, int nperf_reserve_cap); + +/** */ int -wells_add(enum well_type type , +add_wells(enum WellType type , double depth_ref, int nperf , const double *zfrac , /* Injection fraction or NULL */ @@ -74,15 +118,15 @@ wells_add(enum well_type type , struct Wells *W ); void -wells_destroy(struct Wells *W); +destroy_wells(struct Wells *W); int -well_controls_append(enum control_type type , +append_well_controls(enum WellControlType type , double target, struct WellControls *ctrl ); void -well_controls_clear(struct WellControls *ctrl); +clear_well_controls(struct WellControls *ctrl); #ifdef __cplusplus } diff --git a/opm/core/utility/newwells.c b/opm/core/utility/newwells.c index 4e839505..cf918abb 100644 --- a/opm/core/utility/newwells.c +++ b/opm/core/utility/newwells.c @@ -245,7 +245,7 @@ wells_reserve(int nwells, int nperf, struct Wells *W) /* ---------------------------------------------------------------------- */ struct Wells * -wells_create(int nwells, int nperf) +create_wells(int nwells_reserve_cap, int nperf_reserve_cap) /* ---------------------------------------------------------------------- */ { int ok; @@ -272,13 +272,13 @@ wells_create(int nwells, int nperf) if (ok) { W->well_connpos[0] = 0; - if ((nwells > 0) || (nperf > 0)) { - ok = wells_reserve(nwells, nperf, W); + if ((nwells_reserve_cap > 0) || (nperf_reserve_cap > 0)) { + ok = wells_reserve(nwells_reserve_cap, nperf_reserve_cap, W); } } if (! ok) { - wells_destroy(W); + destroy_wells(W); W = NULL; } } @@ -289,7 +289,7 @@ wells_create(int nwells, int nperf) /* ---------------------------------------------------------------------- */ void -wells_destroy(struct Wells *W) +destroy_wells(struct Wells *W) /* ---------------------------------------------------------------------- */ { int w; @@ -331,7 +331,7 @@ alloc_size(int n, int a, int cpty) /* ---------------------------------------------------------------------- */ int -wells_add(enum well_type type , +add_wells(enum WellType type , double depth_ref, int nperf , const double *zfrac , /* Injection fraction or NULL */ @@ -386,7 +386,7 @@ wells_add(enum well_type type , /* ---------------------------------------------------------------------- */ int -well_controls_append(enum control_type type , +append_well_controls(enum WellControlType type , double target, struct WellControls *ctrl ) /* ---------------------------------------------------------------------- */ @@ -415,7 +415,7 @@ well_controls_append(enum control_type type , /* ---------------------------------------------------------------------- */ void -well_controls_clear(struct WellControls *ctrl) +clear_well_controls(struct WellControls *ctrl) /* ---------------------------------------------------------------------- */ { if (ctrl != NULL) { From c09c660af89e388ebdccfb49227c46eeca82f6a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Apr 2012 11:34:58 +0200 Subject: [PATCH 31/51] Renamed add_wells() -> add_well(). Documented. --- opm/core/WellsManager.cpp | 4 +-- opm/core/newwells.h | 53 ++++++++++++++++++++++++++++++------- opm/core/utility/newwells.c | 14 +++++----- 3 files changed, 52 insertions(+), 19 deletions(-) diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index f0933649..8194ae32 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -590,8 +590,8 @@ namespace Opm } const double* zfrac = (well_data[w].type == INJECTOR) ? fracs[well_data[w].injected_phase] : 0; - int ok = add_wells(well_data[w].type, well_data[w].reference_bhp_depth, nperf, - zfrac, &cells[0], &wi[0], w_); + int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, nperf, + zfrac, &cells[0], &wi[0], w_); if (!ok) { THROW("Failed to add a well."); } diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 4be7ef4c..c28138a1 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -102,32 +102,65 @@ struct CompletionData * The arguments may be used to indicate expected capacity needed, * they will be used internally for pre-allocation. * \return NULL upon failure, otherwise a valid Wells object with 0 wells. + * Call add_well() to populate the Wells object. + * Call destroy_wells() to deallocate and clean up the Wells object. */ struct Wells * create_wells(int nwells_reserve_cap, int nperf_reserve_cap); -/** */ +/** Append a well to a Wells struct. + * If successful, W->number_of_wells is incremented by 1. + * The newly added well will have no controls associated with it, add + * controls using append_well_controls(). The current control index is set + * to -1 (invalid). + * \param[in] type Type of well. + * \param[in] depth_ref Reference depth for bhp. + * \param[in] nperf Number of perforations. + * \param[in] zfrac Injection fraction (three components) or NULL. + * \param[in] cells Perforation cell indices. + * \param[in] WI Well production index per perforation, or NULL. + * \param[inout] W The Wells object to be modified. + * \return 1 if successful, 0 if failed. + */ int -add_wells(enum WellType type , - double depth_ref, - int nperf , - const double *zfrac , /* Injection fraction or NULL */ - const int *cells , - const double *WI , /* Well index per perf (or NULL) */ - struct Wells *W ); +add_well(enum WellType type , + double depth_ref, + int nperf , + const double *zfrac , + const int *cells , + const double *WI , + struct Wells *W ); -void -destroy_wells(struct Wells *W); +/** Append a control to a well. + * If successful, ctrl->num is incremented by 1. + * Note that this function does not change ctrl->current. + * To append a control to a well with index w, pass its + * controls to this function via wellsptr->ctrls[w]. + * \param[in] type Control type. + * \param[in] target Target value for the control. + * \param[inout] ctrl The WellControls object to be modified. + * \return 1 if successful, 0 if failed. + */ int append_well_controls(enum WellControlType type , double target, struct WellControls *ctrl ); +/** Clear all controls from a well. */ void clear_well_controls(struct WellControls *ctrl); + +/** Destruction function for Wells objects. + * Assumes that create_wells() and add_wells() have been used to + * build the object. + */ +void +destroy_wells(struct Wells *W); + + #ifdef __cplusplus } #endif diff --git a/opm/core/utility/newwells.c b/opm/core/utility/newwells.c index cf918abb..55a82ffb 100644 --- a/opm/core/utility/newwells.c +++ b/opm/core/utility/newwells.c @@ -331,13 +331,13 @@ alloc_size(int n, int a, int cpty) /* ---------------------------------------------------------------------- */ int -add_wells(enum WellType type , - double depth_ref, - int nperf , - const double *zfrac , /* Injection fraction or NULL */ - const int *cells , - const double *WI , /* Well index per perf (or NULL) */ - struct Wells *W ) +add_well(enum WellType type , + double depth_ref, + int nperf , + const double *zfrac , /* Injection fraction or NULL */ + const int *cells , + const double *WI , /* Well index per perf (or NULL) */ + struct Wells *W ) /* ---------------------------------------------------------------------- */ { int ok, nw, nperf_tot, off; From 73c74e422839fea7cdaa98f858e8d003e55d30c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Apr 2012 11:38:10 +0200 Subject: [PATCH 32/51] Added file level docs. --- opm/core/newwells.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index c28138a1..64a24faa 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -20,6 +20,15 @@ #ifndef OPM_NEWWELLS_H_INCLUDED #define OPM_NEWWELLS_H_INCLUDED + +/** + * \file + * + * Main OPM-Core well data structure along with functions + * to create, populate and destroy it. + */ + + #ifdef __cplusplus extern "C" { #endif From b78553c10c245f93d7f4bda5399095914f0885fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Apr 2012 14:00:38 +0200 Subject: [PATCH 33/51] Accounting for (constant) formation volume factor in incompressible fluids. --- opm/core/fluid/IncompPropertiesFromDeck.cpp | 2 +- opm/core/fluid/PvtPropertiesIncompFromDeck.cpp | 18 +++++++++++++++--- opm/core/fluid/PvtPropertiesIncompFromDeck.hpp | 18 ++++++++++++++---- 3 files changed, 30 insertions(+), 8 deletions(-) diff --git a/opm/core/fluid/IncompPropertiesFromDeck.cpp b/opm/core/fluid/IncompPropertiesFromDeck.cpp index 9febc085..2789fe92 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.cpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.cpp @@ -87,7 +87,7 @@ namespace Opm /// \return Array of P density values. const double* IncompPropertiesFromDeck::density() const { - return pvt_.surfaceDensities(); + return pvt_.reservoirDensities(); } /// \param[in] n Number of data points. diff --git a/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp b/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp index dc5e8043..62cc1b23 100644 --- a/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp +++ b/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp @@ -51,18 +51,23 @@ namespace Opm if (deck.hasField("DENSITY")) { const std::vector& d = deck.getDENSITY().densities_[region_number]; enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; - density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water]; - density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil]; + 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& 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. @@ -76,6 +81,7 @@ namespace Opm 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"); @@ -84,7 +90,13 @@ namespace Opm const double* PvtPropertiesIncompFromDeck::surfaceDensities() const { - return density_.data(); + return surface_density_.data(); + } + + + const double* PvtPropertiesIncompFromDeck::reservoirDensities() const + { + return reservoir_density_.data(); } diff --git a/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp b/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp index 84f14374..acbabaa4 100644 --- a/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp +++ b/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp @@ -31,9 +31,6 @@ namespace Opm /// eclipse input (keywords DENSITY, PVTW, PVCDO). /// /// 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 PvtPropertiesIncompFromDeck { @@ -51,11 +48,24 @@ namespace Opm /// \return Array of size numPhases(). const double* surfaceDensities() const; + /// Densities of stock components at reservoir conditions. + /// Note: a reasonable question to ask is why there can be + /// different densities at surface and reservoir conditions, + /// when the phases are assumed incompressible. The answer is + /// that even if we approximate the phases as being + /// incompressible during simulation, the density difference + /// between surface and reservoir may be larger. For accurate + /// 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; + /// Viscosities. const double* viscosity() const; private: - std::tr1::array density_; + std::tr1::array surface_density_; + std::tr1::array reservoir_density_; std::tr1::array viscosity_; }; From 8cdabd9155c4ff628e7c96c865e2214c62c56eed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 20 Apr 2012 19:47:52 +0200 Subject: [PATCH 34/51] Remove an instance of EOL whitespace introduced in cset a09de38891c5. --- opm/core/newwells.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 64a24faa..48e4c559 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -35,7 +35,7 @@ extern "C" { /** Well type indicates desired/expected well behaviour. */ enum WellType { INJECTOR, PRODUCER }; -/** Type of well control equation or inequality. +/** Type of well control equation or inequality. * BHP -> bottom hole pressure is specified. * RATE -> flow rate is specified. */ From d0d45e0ccad4b5d03b5a9641f2ac969381671c61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 20 Apr 2012 19:53:09 +0200 Subject: [PATCH 35/51] Use more precise description of the WellControlType. --- opm/core/newwells.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 48e4c559..2884dfc2 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -35,7 +35,7 @@ extern "C" { /** Well type indicates desired/expected well behaviour. */ enum WellType { INJECTOR, PRODUCER }; -/** Type of well control equation or inequality. +/** Type of well control equation or inequality constraint. * BHP -> bottom hole pressure is specified. * RATE -> flow rate is specified. */ From 165c4308441bd5977009667bd418fa797b221dbc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 20 Apr 2012 20:54:58 +0200 Subject: [PATCH 36/51] Correct array sizes in struct CompletionData. --- opm/core/newwells.h | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 2884dfc2..ce04e6a3 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -92,13 +92,17 @@ struct Wells /** Struct encapsulating dynamic information about all wells in a scenario. - * All arrays in this struct contain data for each perforation, ordered - * the same as Wells::well_cells and Wells:WI. Letting NP be the number - * of perforations, the array sizes are: - * gpot 3*NP - * A 9*NP (matrix in Fortran order). - * phasemob 3*NP - * \TODO: Verify that the sizes are correct, check if we should refactor to handle two phases better. + * All arrays in this struct contain data for each perforation, + * ordered the same as Wells::well_cells and Wells:WI. The array + * sizes are, respectively, + * + * gpot n*NP + * A n²*NP (matrix in row-major (i.e., Fortran) order). + * phasemob n*NP + * + * in which "n" denotes the number of active fluid phases (and + * constituent components) and "NP" is the total number of + * perforations, well_connpos[ number_of_wells ]. */ struct CompletionData { From fd58342f40b66fae3b5be46a469ddacf984e3608 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 20 Apr 2012 23:53:04 +0200 Subject: [PATCH 37/51] Correct misunderstanding about row-major vs. column major ordering. --- opm/core/newwells.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index ce04e6a3..5b1234cc 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -97,7 +97,7 @@ struct Wells * sizes are, respectively, * * gpot n*NP - * A n²*NP (matrix in row-major (i.e., Fortran) order). + * A n²*NP (matrix in column-major (i.e., Fortran) order). * phasemob n*NP * * in which "n" denotes the number of active fluid phases (and From 868905ae24476b577e258323d2700598c5bd8a77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sat, 21 Apr 2012 10:10:11 +0200 Subject: [PATCH 38/51] Use complete sentences when describing individual WellControlTypes. --- opm/core/newwells.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 5b1234cc..3e29f058 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -36,8 +36,8 @@ extern "C" { /** Well type indicates desired/expected well behaviour. */ enum WellType { INJECTOR, PRODUCER }; /** Type of well control equation or inequality constraint. - * BHP -> bottom hole pressure is specified. - * RATE -> flow rate is specified. + * BHP -> Well constrained by bottom-hole pressure target. + * RATE -> Well constrained by total reservoir volume flow rate. */ enum WellControlType { BHP , RATE }; /** Canonical component names and ordering. */ From 75ee1494ec00d50cbebac1886dacd666e7dcc736 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sat, 21 Apr 2012 10:18:20 +0200 Subject: [PATCH 39/51] Formalise description of struct WellControls. --- opm/core/newwells.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 3e29f058..3e14fb9b 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -45,15 +45,16 @@ enum SurfaceComponent { WATER = 0, OIL = 1, GAS = 2 }; /** Controls for a single well. - * Each control specifies a well rate or bottom hole pressure. Only + * Each control specifies a well rate or bottom-hole pressure. Only * one control can be active at a time, indicated by current. The * meaning of each control's target value depends on the control * type, for BHP controls it is a pressure in Pascal, for RATE * controls it is a volumetric rate in cubic(meter)/second. The - * active control should be interpreted as an equation, whereas the - * non-active controls should be interpreted as inequalities - * specifying constraints on the solution, where BHP controls yield - * minimum pressures, and RATE controls yield maximum rates. + * active control as an equality constraint, whereas the + * non-active controls should be interpreted as inequality + * constraints (upper or lower bounds). For instance, a PRODUCER's BHP + * constraint defines a minimum acceptable bottom-hole pressure value + * for the well. */ struct WellControls { From 3807506aeb96dd46c9f656cdddddace4b3597bf0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sat, 21 Apr 2012 11:34:11 +0200 Subject: [PATCH 40/51] Hide memory management aspects of struct WellControls. --- opm/core/newwells.h | 3 ++- opm/core/utility/newwells.c | 54 ++++++++++++++++++++++++++++++++----- 2 files changed, 49 insertions(+), 8 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 3e14fb9b..381b31e5 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -59,10 +59,11 @@ enum SurfaceComponent { WATER = 0, OIL = 1, GAS = 2 }; struct WellControls { int num; /** Number of controls. */ - int cpty; /** Allocated capacity, for internal use only. */ enum WellControlType *type; /** Array of control types. */ double *target; /** Array of control targets. */ int current; /** Index of current active control. */ + + void *data; /** Internal management structure. */ }; diff --git a/opm/core/utility/newwells.c b/opm/core/utility/newwells.c index 55a82ffb..3960bb82 100644 --- a/opm/core/utility/newwells.c +++ b/opm/core/utility/newwells.c @@ -41,14 +41,42 @@ #include #include +struct WellControlMgmt { + int cpty; +}; + + +static void +destroy_ctrl_mgmt(struct WellControlMgmt *m) +{ + free(m); +} + + +static struct WellControlMgmt * +create_ctrl_mgmt(void) +{ + struct WellControlMgmt *m; + + m = malloc(1 * sizeof *m); + + if (m != NULL) { + m->cpty = 0; + } + + return m; +} + + /* ---------------------------------------------------------------------- */ static void well_controls_destroy(struct WellControls *ctrl) /* ---------------------------------------------------------------------- */ { if (ctrl != NULL) { - free(ctrl->target); - free(ctrl->type); + destroy_ctrl_mgmt(ctrl->data); + free (ctrl->target); + free (ctrl->type); } free(ctrl); @@ -67,10 +95,16 @@ well_controls_create(void) if (ctrl != NULL) { /* Initialise empty control set */ ctrl->num = 0; - ctrl->cpty = 0; ctrl->type = NULL; ctrl->target = NULL; ctrl->current = -1; + + ctrl->data = create_ctrl_mgmt(); + + if (ctrl->data == NULL) { + well_controls_destroy(ctrl); + ctrl = NULL; + } } return ctrl; @@ -85,6 +119,8 @@ well_controls_reserve(int nctrl, struct WellControls *ctrl) int c, ok; void *type, *target; + struct WellControlMgmt *m; + type = realloc(ctrl->type , nctrl * sizeof *ctrl->type ); target = realloc(ctrl->target, nctrl * sizeof *ctrl->target); @@ -93,12 +129,13 @@ well_controls_reserve(int nctrl, struct WellControls *ctrl) if (target != NULL) { ctrl->target = target; ok++; } if (ok == 2) { - for (c = ctrl->cpty; c < nctrl; c++) { + m = ctrl->data; + for (c = m->cpty; c < nctrl; c++) { ctrl->type [c] = BHP; ctrl->target[c] = -1.0; } - ctrl->cpty = nctrl; + m->cpty = nctrl; } return ok == 2; @@ -393,12 +430,15 @@ append_well_controls(enum WellControlType type , { int ok, alloc; + struct WellControlMgmt *m; + assert (ctrl != NULL); - ok = ctrl->num < ctrl->cpty; + m = ctrl->data; + ok = ctrl->num < m->cpty; if (! ok) { - alloc = alloc_size(ctrl->num, 1, ctrl->cpty); + alloc = alloc_size(ctrl->num, 1, m->cpty); ok = well_controls_reserve(alloc, ctrl); } From 1e9748cdb88b0e178b060fad9f8269f33f46c22e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sat, 21 Apr 2012 13:57:41 +0200 Subject: [PATCH 41/51] Hide memory management aspects of struct Wells. --- opm/core/newwells.h | 6 ++- opm/core/utility/newwells.c | 82 +++++++++++++++++++++++++++++-------- 2 files changed, 70 insertions(+), 18 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 381b31e5..e806782c 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -72,8 +72,7 @@ struct WellControls struct Wells { int number_of_wells; /** Number of wells. */ - int well_cpty; /** Allocated well capacity, for internal use only. */ - int perf_cpty; /** Allocated perforation capacity, for internal use only. */ + enum WellType *type; /** Array of well types. */ double *depth_ref; /** Array of well bhp reference depths. */ double *zfrac; /** Component volume fractions for each well, size is (3*number_of_wells). @@ -90,6 +89,9 @@ struct Wells */ double *WI; /** Well productivity index, same size and structure as well_cells. */ struct WellControls **ctrls; /** Well controls, one struct for each well. */ + + + void *data; /** Internal management structure. */ }; diff --git a/opm/core/utility/newwells.c b/opm/core/utility/newwells.c index 3960bb82..c241ea23 100644 --- a/opm/core/utility/newwells.c +++ b/opm/core/utility/newwells.c @@ -45,6 +45,11 @@ struct WellControlMgmt { int cpty; }; +struct WellMgmt { + int well_cpty; + int perf_cpty; +}; + static void destroy_ctrl_mgmt(struct WellControlMgmt *m) @@ -68,6 +73,29 @@ create_ctrl_mgmt(void) } +static void +destroy_well_mgmt(struct WellMgmt *m) +{ + free(m); +} + + +static struct WellMgmt * +create_well_mgmt(void) +{ + struct WellMgmt *m; + + m = malloc(1 * sizeof *m); + + if (m != NULL) { + m->well_cpty = 0; + m->perf_cpty = 0; + } + + return m; +} + + /* ---------------------------------------------------------------------- */ static void well_controls_destroy(struct WellControls *ctrl) @@ -197,7 +225,11 @@ initialise_new_wells(int nwells, struct Wells *W) { int ok, w; - for (w = W->well_cpty; w < nwells; w++) { + struct WellMgmt *m; + + m = W->data; + + for (w = m->well_cpty; w < nwells; w++) { W->type [w] = PRODUCER; W->depth_ref[w] = -1.0; @@ -208,7 +240,7 @@ initialise_new_wells(int nwells, struct Wells *W) W->well_connpos[w + 1] = W->well_connpos[w]; } - for (w = W->well_cpty, ok = 1; ok && (w < nwells); w++) { + for (w = m->well_cpty, ok = 1; ok && (w < nwells); w++) { W->ctrls[w] = well_controls_create(); ok = W->ctrls[w] != NULL; @@ -231,7 +263,11 @@ initialise_new_perfs(int nperf, struct Wells *W) { int k; - for (k = W->perf_cpty; k < nperf; k++) { + struct WellMgmt *m; + + m = W->data; + + for (k = m->perf_cpty; k < nperf; k++) { W->well_cells[k] = -1 ; W->WI [k] = 0.0; } @@ -245,12 +281,16 @@ wells_reserve(int nwells, int nperf, struct Wells *W) { int ok; - assert (nwells >= W->well_cpty); - assert (nperf >= W->perf_cpty); + struct WellMgmt *m; + + m = W->data; + + assert (nwells >= m->well_cpty); + assert (nperf >= m->perf_cpty); ok = 1; - if (nwells > W->well_cpty) { + if (nwells > m->well_cpty) { ok = wells_allocate(nwells, W); if (ok) { @@ -258,16 +298,16 @@ wells_reserve(int nwells, int nperf, struct Wells *W) } if (ok) { - W->well_cpty = nwells; + m->well_cpty = nwells; } } - if (ok && (nperf > W->perf_cpty)) { + if (ok && (nperf > m->perf_cpty)) { ok = perfs_allocate(nperf, W); if (ok) { initialise_new_perfs(nperf, W); - W->perf_cpty = nperf; + m->perf_cpty = nperf; } } @@ -292,8 +332,6 @@ create_wells(int nwells_reserve_cap, int nperf_reserve_cap) if (W != NULL) { W->number_of_wells = 0; - W->well_cpty = 0; - W->perf_cpty = 0; W->type = NULL; W->depth_ref = NULL; @@ -305,7 +343,9 @@ create_wells(int nwells_reserve_cap, int nperf_reserve_cap) W->ctrls = NULL; - ok = W->well_connpos != NULL; + W->data = create_well_mgmt(); + + ok = (W->well_connpos != NULL) && (W->data != NULL); if (ok) { W->well_connpos[0] = 0; @@ -331,11 +371,17 @@ destroy_wells(struct Wells *W) { int w; + struct WellMgmt *m; + if (W != NULL) { - for (w = 0; w < W->well_cpty; w++) { + m = W->data; + + for (w = 0; w < m->well_cpty; w++) { well_controls_destroy(W->ctrls[w]); } + destroy_well_mgmt(m); + free(W->ctrls); free(W->WI); free(W->well_cells); @@ -380,14 +426,18 @@ add_well(enum WellType type , int ok, nw, nperf_tot, off; int nwalloc, nperfalloc; + struct WellMgmt *m; + nw = W->number_of_wells; nperf_tot = W->well_connpos[nw]; - ok = (nw < W->well_cpty) && (nperf_tot + nperf <= W->perf_cpty); + m = W->data; + + ok = (nw < m->well_cpty) && (nperf_tot + nperf <= m->perf_cpty); if (! ok) { - nwalloc = alloc_size(nw , 1 , W->well_cpty); - nperfalloc = alloc_size(nperf_tot, nperf, W->perf_cpty); + nwalloc = alloc_size(nw , 1 , m->well_cpty); + nperfalloc = alloc_size(nperf_tot, nperf, m->perf_cpty); ok = wells_reserve(nwalloc, nperfalloc, W); } From 2a5d35a04b2ee3583ee2438ecd3af61e34f8db1d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sat, 21 Apr 2012 21:07:19 +0200 Subject: [PATCH 42/51] Prefer to speak of "data structures" rather than "struct"s in documentation. --- opm/core/newwells.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index e806782c..29ad6b73 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -68,7 +68,7 @@ struct WellControls -/** Struct encapsulating static information about all wells in a scenario. */ +/** Data structure aggregating static information about all wells in a scenario. */ struct Wells { int number_of_wells; /** Number of wells. */ @@ -88,15 +88,15 @@ struct Wells * Size is number of perforations (== well_connpos[number_of_wells]). */ double *WI; /** Well productivity index, same size and structure as well_cells. */ - struct WellControls **ctrls; /** Well controls, one struct for each well. */ + struct WellControls **ctrls; /** Well controls, one set of controls for each well. */ void *data; /** Internal management structure. */ }; -/** Struct encapsulating dynamic information about all wells in a scenario. - * All arrays in this struct contain data for each perforation, +/** Data structure aggregating dynamic information about all wells in a scenario. + * All arrays in this structure contain data for each perforation, * ordered the same as Wells::well_cells and Wells:WI. The array * sizes are, respectively, * @@ -126,7 +126,7 @@ struct Wells * create_wells(int nwells_reserve_cap, int nperf_reserve_cap); -/** Append a well to a Wells struct. +/** Append a new well to an existing Wells object. * If successful, W->number_of_wells is incremented by 1. * The newly added well will have no controls associated with it, add * controls using append_well_controls(). The current control index is set From ceeafefa52e9a264567083290a1a83502857a728 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 22 Apr 2012 11:36:49 +0200 Subject: [PATCH 43/51] Update add_well() comment to reflect intended pre- and post-condition. --- opm/core/newwells.h | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 29ad6b73..829494f8 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -126,19 +126,22 @@ struct Wells * create_wells(int nwells_reserve_cap, int nperf_reserve_cap); -/** Append a new well to an existing Wells object. - * If successful, W->number_of_wells is incremented by 1. - * The newly added well will have no controls associated with it, add - * controls using append_well_controls(). The current control index is set - * to -1 (invalid). - * \param[in] type Type of well. - * \param[in] depth_ref Reference depth for bhp. - * \param[in] nperf Number of perforations. - * \param[in] zfrac Injection fraction (three components) or NULL. - * \param[in] cells Perforation cell indices. - * \param[in] WI Well production index per perforation, or NULL. - * \param[inout] W The Wells object to be modified. - * \return 1 if successful, 0 if failed. +/** + * Append a new well to an existing Wells object. + * + * Increments W->number_of_wells by one if successful. The new well + * does not include operational constraints. Such information is + * specified using function append_well_controls(). The current + * control index is set to -1 (invalid). + * + * \param[in] type Type of well. + * \param[in] depth_ref Reference depth for bhp. + * \param[in] nperf Number of perforations. + * \param[in] zfrac Injection fraction (three components) or NULL. + * \param[in] cells Perforation cell indices. + * \param[in] WI Well production index per perforation, or NULL. + * \param[inout] W The Wells object to be modified. + * \return Non-zero (true) if successful and zero otherwise. */ int add_well(enum WellType type , From f2ab37c6c26bae661c2037526a507b99590c16a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 22 Apr 2012 11:50:41 +0200 Subject: [PATCH 44/51] Expand description of function create_wells(). Restore original parameter names in the process. --- opm/core/newwells.h | 29 ++++++++++++++++++++++------- opm/core/utility/newwells.c | 6 +++--- 2 files changed, 25 insertions(+), 10 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 829494f8..fdf8383e 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -115,15 +115,30 @@ struct CompletionData double *phasemob; /** Phase mobilities. */ }; -/** Contruction function initializing a Wells object. - * The arguments may be used to indicate expected capacity needed, - * they will be used internally for pre-allocation. - * \return NULL upon failure, otherwise a valid Wells object with 0 wells. - * Call add_well() to populate the Wells object. - * Call destroy_wells() to deallocate and clean up the Wells object. +/** + * Construct a Wells object initially capable of managing a given + * number of wells and total number of well connections + * (perforations). + * + * Function add_well() is used to populate the Wells object. No + * reallocation occurrs in function add_well() as long as the + * initially indicated capacites are sufficient. Call function + * destroy_wells() to dispose of the Wells object and its allocated + * memory resources. + * + * \param[in] nwells Expected number of wells in simulation scenario. + * Pass zero if the total number of wells is unknown. + * + * \param[in] nperf Expected total number of well connections + * (perforations) for all wells in simulation + * scenario. Pass zero if the total number of well + * connections is unknown. + * + * \return A valid Wells object with no wells if successful, and NULL + * otherwise. */ struct Wells * -create_wells(int nwells_reserve_cap, int nperf_reserve_cap); +create_wells(int nwells, int nperf); /** diff --git a/opm/core/utility/newwells.c b/opm/core/utility/newwells.c index c241ea23..411a5a58 100644 --- a/opm/core/utility/newwells.c +++ b/opm/core/utility/newwells.c @@ -322,7 +322,7 @@ wells_reserve(int nwells, int nperf, struct Wells *W) /* ---------------------------------------------------------------------- */ struct Wells * -create_wells(int nwells_reserve_cap, int nperf_reserve_cap) +create_wells(int nwells, int nperf) /* ---------------------------------------------------------------------- */ { int ok; @@ -349,8 +349,8 @@ create_wells(int nwells_reserve_cap, int nperf_reserve_cap) if (ok) { W->well_connpos[0] = 0; - if ((nwells_reserve_cap > 0) || (nperf_reserve_cap > 0)) { - ok = wells_reserve(nwells_reserve_cap, nperf_reserve_cap, W); + if ((nwells > 0) || (nperf > 0)) { + ok = wells_reserve(nwells, nperf, W); } } From 225acb5c267fcd217c559578a9a2aeb90788d9a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 22 Apr 2012 16:58:37 +0200 Subject: [PATCH 45/51] Document intended return code for append_well_controls(). --- opm/core/newwells.h | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index fdf8383e..fa9a737e 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -168,15 +168,17 @@ add_well(enum WellType type , struct Wells *W ); -/** Append a control to a well. - * If successful, ctrl->num is incremented by 1. - * Note that this function does not change ctrl->current. - * To append a control to a well with index w, pass its - * controls to this function via wellsptr->ctrls[w]. - * \param[in] type Control type. - * \param[in] target Target value for the control. - * \param[inout] ctrl The WellControls object to be modified. - * \return 1 if successful, 0 if failed. +/** + * Append operational constraint to an existing well. + * + * Increments ctrl->num by one if successful. Introducing a new + * operational constraint does not affect the well's notion of the + * currently active constraint represented by ctrl->current. + * + * \param[in] type Control type. + * \param[in] target Target value for the control. + * \param[inout] ctrl Existing set of well controls. + * \return Non-zero (true) if successful and zero (false) otherwise. */ int append_well_controls(enum WellControlType type , From c6c1f08bd75a1ca48ee4129e7384e84583b7a85d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 22 Apr 2012 17:03:13 +0200 Subject: [PATCH 46/51] Fix cross-reference misprint in destroy_wells() documentation. --- opm/core/newwells.h | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index fa9a737e..2df51a19 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -190,9 +190,13 @@ void clear_well_controls(struct WellControls *ctrl); -/** Destruction function for Wells objects. - * Assumes that create_wells() and add_wells() have been used to - * build the object. +/** + * Wells object destructor. + * + * Disposes of all resources managed by the Wells object. + * + * The Wells object must be built using function create_wells() and + * subsequently populated using function add_well(). */ void destroy_wells(struct Wells *W); From 375283b14183e762cfda5a6d14be5932dc3abaf0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 22 Apr 2012 21:15:28 +0200 Subject: [PATCH 47/51] Mention that clear_well_controls() does not affect capacity. --- opm/core/newwells.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 2df51a19..45dee397 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -185,7 +185,10 @@ append_well_controls(enum WellControlType type , double target, struct WellControls *ctrl ); -/** Clear all controls from a well. */ +/** + * Clear all controls from a well. + * + * Does not affect the control set capacity. */ void clear_well_controls(struct WellControls *ctrl); From 9b162de9791f95e5ffeccf2c59a20b403fb57550 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 22 Apr 2012 21:19:08 +0200 Subject: [PATCH 48/51] Fix spelling. --- opm/core/newwells.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 45dee397..26b213bc 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -121,7 +121,7 @@ struct CompletionData * (perforations). * * Function add_well() is used to populate the Wells object. No - * reallocation occurrs in function add_well() as long as the + * reallocation occurs in function add_well() as long as the * initially indicated capacites are sufficient. Call function * destroy_wells() to dispose of the Wells object and its allocated * memory resources. From 7bcf30f7f891dee86e07692c6dafc03cd5f383fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 22 Apr 2012 21:24:43 +0200 Subject: [PATCH 49/51] Further refinement of add_well() documentation. --- opm/core/newwells.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 26b213bc..3cce0338 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -150,12 +150,15 @@ create_wells(int nwells, int nperf); * control index is set to -1 (invalid). * * \param[in] type Type of well. - * \param[in] depth_ref Reference depth for bhp. + * \param[in] depth_ref Reference depth for well's BHP. * \param[in] nperf Number of perforations. * \param[in] zfrac Injection fraction (three components) or NULL. - * \param[in] cells Perforation cell indices. + * \param[in] cells Grid cells in which well is perforated. Should + * ideally be track ordered. * \param[in] WI Well production index per perforation, or NULL. - * \param[inout] W The Wells object to be modified. + * \param[in,out] W Existing set of wells to which new well will + * be added. + * * \return Non-zero (true) if successful and zero otherwise. */ int From c0a5879d53a6fac44a3cbafc67710276a34e31ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 23 Apr 2012 09:59:53 +0200 Subject: [PATCH 50/51] Spell-check comment. --- opm/core/newwells.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 3cce0338..4940c10c 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -122,7 +122,7 @@ struct CompletionData * * Function add_well() is used to populate the Wells object. No * reallocation occurs in function add_well() as long as the - * initially indicated capacites are sufficient. Call function + * initially indicated capacities are sufficient. Call function * destroy_wells() to dispose of the Wells object and its allocated * memory resources. * From 1f9fea92fa50dbc3e487515f9feb0265c416d3e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 23 Apr 2012 11:10:55 +0200 Subject: [PATCH 51/51] Catch exceptions from boost::create_directories() to give useful error message. --- examples/spu_2p.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index e5357758..4c877da6 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -266,7 +266,12 @@ main(int argc, char** argv) output_dir = param.getDefault("output_dir", std::string("output")); // Ensure that output dir exists boost::filesystem::path fpath(output_dir); - create_directories(fpath); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } output_interval = param.getDefault("output_interval", output_interval); } const int num_transport_substeps = param.getDefault("num_transport_substeps", 1);