From 368e1813a17fee85f84e8b0bdab17b7df2c16e51 Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Thu, 28 Jun 2012 13:20:18 +0200 Subject: [PATCH 01/59] Using Datamapper.h in writeVtkData --- opm/core/utility/writeVtkData.hpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/opm/core/utility/writeVtkData.hpp b/opm/core/utility/writeVtkData.hpp index 200b2b261..96cdfac22 100644 --- a/opm/core/utility/writeVtkData.hpp +++ b/opm/core/utility/writeVtkData.hpp @@ -26,25 +26,23 @@ #include #include #include +#include struct UnstructuredGrid; namespace Opm { - /// Intended to map strings (giving the output field names) to data. - typedef std::map*> DataMap; - /// Vtk output for cartesian grids. void writeVtkData(const std::tr1::array& dims, - const std::tr1::array& cell_size, - const DataMap& data, - std::ostream& os); + const std::tr1::array& cell_size, + const DataMap& data, + std::ostream& os); /// Vtk output for general grids. void writeVtkData(const UnstructuredGrid& grid, - const DataMap& data, - std::ostream& os); + const DataMap& data, + std::ostream& os); } // namespace Opm #endif // OPM_WRITEVTKDATA_HEADER_INCLUDED From 3e9cc7492231568263d214161754934e7099c692 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Aug 2012 16:08:26 +0200 Subject: [PATCH 02/59] Initial commit of tof computation by reordering. --- .../reorder/TransportModelTracerTof.cpp | 100 ++++++++++++++++++ .../reorder/TransportModelTracerTof.hpp | 63 +++++++++++ 2 files changed, 163 insertions(+) create mode 100644 opm/core/transport/reorder/TransportModelTracerTof.cpp create mode 100644 opm/core/transport/reorder/TransportModelTracerTof.hpp diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp new file mode 100644 index 000000000..683b732c3 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -0,0 +1,100 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include +#include + +namespace Opm +{ + + + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid) + : grid_(grid) + { + } + + + + + /// Solve for time-of-flight at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in, out] tof Array of time-of-flight values. + void TransportModelTracerTof::solveTof(const double* darcyflux, + const double* porevolume, + std::vector& tof) + { + darcyflux_ = darcyflux; + porevolume_ = porevolume; + tof.resize(grid_.number_of_cells); + std::fill(tof.begin(), tof.end(), 0.0); + tof_ = &tof[0]; + reorderAndTransport(grid_, darcyflux); + } + + + + + void TransportModelTracerTof::solveSingleCell(const int cell) + { + // Compute flux terms. + double upwind_term = 0.0; + double downwind_flux = 0.0; + for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { + int f = grid_.cell_faces[i]; + double flux; + int other; + // Compute cell flux + if (cell == grid_.face_cells[2*f]) { + flux = darcyflux_[f]; + other = grid_.face_cells[2*f+1]; + } else { + flux =-darcyflux_[f]; + other = grid_.face_cells[2*f]; + } + // Add flux to upwind_term or downwind_flux, if interior. + if (other != -1) { + if (flux < 0.0) { + upwind_term += flux*tof_[other]; + } else { + downwind_flux += flux; + } + } + } + + // Compute tof. + tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux; + } + + + + + void TransportModelTracerTof::solveMultiCell(const int /*num_cells*/, const int* /*cells*/) + { + THROW("Not yet implemented."); + } + + + + +} // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp new file mode 100644 index 000000000..1ab7ef3b4 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -0,0 +1,63 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED +#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED + +#include +#include +#include +#include +struct UnstructuredGrid; + +namespace Opm +{ + + class IncompPropertiesInterface; + + /// Implements a reordering transport solver for incompressible two-phase flow. + class TransportModelTracerTof : public TransportModelInterface + { + public: + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTof(const UnstructuredGrid& grid); + + /// Solve for time-of-flight at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in, out] tof Array of time-of-flight values. + void solveTof(const double* darcyflux, + const double* porevolume, + std::vector& tof); + + private: + virtual void solveSingleCell(const int cell); + virtual void solveMultiCell(const int num_cells, const int* cells); + + private: + const UnstructuredGrid& grid_; + const double* darcyflux_; // one flux per grid face + const double* porevolume_; // one volume per cell + double* tof_; + }; + +} // namespace Opm + +#endif // OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED From 931dcc4a3de845393d129783f5dbe0c0cf5a7cc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 11:26:51 +0200 Subject: [PATCH 03/59] Implemented rudimentary solveMultiCell(). Simply calls solveSingleCell() once for each cell in block. --- opm/core/transport/reorder/TransportModelTracerTof.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 683b732c3..0825706e3 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -89,9 +89,12 @@ namespace Opm - void TransportModelTracerTof::solveMultiCell(const int /*num_cells*/, const int* /*cells*/) + void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells) { - THROW("Not yet implemented."); + std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; + for (int i = 0; i < num_cells; ++i) { + solveSingleCell(cells[i]); + } } From d9cff689b49f382a20ed61c28ef64180a2cea5f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 16:48:21 +0200 Subject: [PATCH 04/59] Added class SinglePvtDead, add parameter 'props_use_spline' to simulators. Recall that the class that used to be called SinglePvtDead has been renamed to SinglePvtDeadSpline. If 'props_use_spline' is true, that class is used (this is the default), which makes a monotone spline that is uniformly, densely sampled. The new class simply uses linear interpolation in the input tables. --- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 5 +++-- opm/core/fluid/BlackoilPropertiesFromDeck.hpp | 3 ++- opm/core/fluid/blackoil/BlackoilPvtProperties.hpp | 2 +- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 004e6f88f..e7244d87c 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -23,10 +23,11 @@ namespace Opm { BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid) + const UnstructuredGrid& grid, + const bool use_spline) { rock_.init(deck, grid); - pvt_.init(deck); + pvt_.init(deck, use_spline); satprops_.init(deck, grid); if (pvt_.numPhases() != satprops_.numPhases()) { THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index a2a6c3de9..c434b7f4b 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -43,7 +43,8 @@ namespace Opm /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + const UnstructuredGrid& grid, + const bool use_spline); /// Destructor. virtual ~BlackoilPropertiesFromDeck(); diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp index 627f266b7..ca193a4d4 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp @@ -47,7 +47,7 @@ namespace Opm BlackoilPvtProperties(); /// Initialize from deck. - void init(const EclipseGridParser& deck); + void init(const EclipseGridParser& deck, const bool use_spline); /// Number of active phases. int numPhases() const; From 7a79bd1872ff5e3207c96d7dfcb98dae0194521e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 17:56:01 +0200 Subject: [PATCH 05/59] Enable choice of spline-smoothed saturation props (or not). --- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 22 +- opm/core/fluid/BlackoilPropertiesFromDeck.hpp | 3 +- opm/core/fluid/IncompPropertiesFromDeck.hpp | 2 +- opm/core/fluid/SaturationPropsFromDeck.hpp | 45 ++- .../fluid/SaturationPropsFromDeck_impl.hpp | 271 ++++++++++++++++++ opm/core/fluid/SaturationPropsInterface.hpp | 85 ++++++ 6 files changed, 397 insertions(+), 31 deletions(-) create mode 100644 opm/core/fluid/SaturationPropsFromDeck_impl.hpp create mode 100644 opm/core/fluid/SaturationPropsInterface.hpp diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index e7244d87c..fdc9519cf 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -28,10 +28,20 @@ namespace Opm { rock_.init(deck, grid); pvt_.init(deck, use_spline); - satprops_.init(deck, grid); - if (pvt_.numPhases() != satprops_.numPhases()) { + // Unfortunate lack of pointer smartness here... + if (use_spline) { + SaturationPropsFromDeck<>* ptr = new SaturationPropsFromDeck<>(); + ptr->init(deck, grid); + satprops_.reset(ptr); + } else { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + ptr->init(deck, grid); + satprops_.reset(ptr); + } + if (pvt_.numPhases() != satprops_->numPhases()) { THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); } } @@ -236,7 +246,7 @@ namespace Opm double* kr, double* dkrds) const { - satprops_.relperm(n, s, cells, kr, dkrds); + satprops_->relperm(n, s, cells, kr, dkrds); } @@ -255,7 +265,7 @@ namespace Opm double* pc, double* dpcds) const { - satprops_.capPress(n, s, cells, pc, dpcds); + satprops_->capPress(n, s, cells, pc, dpcds); } @@ -271,7 +281,7 @@ namespace Opm double* smin, double* smax) const { - satprops_.satRange(n, cells, smin, smax); + satprops_->satRange(n, cells, smin, smax); } diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index c434b7f4b..60dc1e754 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -26,6 +26,7 @@ #include #include #include +#include struct UnstructuredGrid; @@ -166,7 +167,7 @@ namespace Opm private: RockFromDeck rock_; BlackoilPvtProperties pvt_; - SaturationPropsFromDeck satprops_; + boost::scoped_ptr satprops_; mutable std::vector B_; mutable std::vector dB_; mutable std::vector R_; diff --git a/opm/core/fluid/IncompPropertiesFromDeck.hpp b/opm/core/fluid/IncompPropertiesFromDeck.hpp index 68623e5cc..5820bb668 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.hpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.hpp @@ -135,7 +135,7 @@ namespace Opm private: RockFromDeck rock_; PvtPropertiesIncompFromDeck pvt_; - SaturationPropsFromDeck satprops_; + SaturationPropsFromDeck<> satprops_; }; diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index b83a12711..a978bbac0 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -20,8 +20,8 @@ #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED +#include #include -#include #include #include @@ -30,7 +30,23 @@ struct UnstructuredGrid; namespace Opm { - class SaturationPropsFromDeck : public BlackoilPhases + + /// Class storing saturation functions in a uniform table, + /// densely sampled from a monotone spline, + /// using linear interpolation. + class SatFuncSetUniform; + + /// Class storing saturation functions in a nonuniform table + /// using linear interpolation. + class SatFuncSetNonuniform; + + + /// Interface to saturation functions from deck. + /// Possible values for template argument: + /// SatFuncSetNonuniform, + /// SatFuncSetUniform (default). + template + class SaturationPropsFromDeck : public SaturationPropsInterface { public: /// Default constructor. @@ -88,30 +104,12 @@ namespace Opm private: PhaseUsage phase_usage_; - 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; + typedef SatFuncSet Funcs; + + const Funcs& funcForCell(const int cell) const; }; @@ -119,6 +117,7 @@ namespace Opm } // namespace Opm +#include #endif // OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED diff --git a/opm/core/fluid/SaturationPropsFromDeck_impl.hpp b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp new file mode 100644 index 000000000..f64f12bcf --- /dev/null +++ b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp @@ -0,0 +1,271 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED +#define OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED + + +#include +#include +#include +#include + +namespace Opm +{ + + /// Class storing saturation functions in a uniform table, + /// densely sampled from a monotone spline, + /// using linear interpolation. + class SatFuncSetUniform : public BlackoilPhases + { + 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) + }; + + + + + /// Class storing saturation functions in a nonuniform table + /// using linear interpolation. + class SatFuncSetNonuniform : public BlackoilPhases + { + 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_. + NonuniformTableLinear krw_; + NonuniformTableLinear krow_; + NonuniformTableLinear pcow_; + NonuniformTableLinear krg_; + NonuniformTableLinear krog_; + NonuniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; + + + + // ----------- Methods of SaturationPropsFromDeck --------- + + + /// Default constructor. + template + SaturationPropsFromDeck::SaturationPropsFromDeck() + { + } + + /// Initialize from deck. + template + void SaturationPropsFromDeck::init(const EclipseGridParser& deck, + const UnstructuredGrid& grid) + { + phase_usage_ = phaseUsageFromDeck(deck); + + // Extract input data. + // Oil phase should be active. + if (!phase_usage_.phase_used[Liquid]) { + THROW("SaturationPropsFromDeck::init() -- oil phase must be active."); + } + + // Obtain SATNUM, if it exists, and create cell_to_func_. + // Otherwise, let the cell_to_func_ mapping be just empty. + int satfuncs_expected = 1; + if (deck.hasField("SATNUM")) { + const std::vector& satnum = deck.getIntegerValue("SATNUM"); + satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); + const int num_cells = grid.number_of_cells; + cell_to_func_.resize(num_cells); + const int* gc = grid.global_cell; + for (int cell = 0; cell < num_cells; ++cell) { + const int deck_pos = (gc == NULL) ? cell : gc[cell]; + cell_to_func_[cell] = satnum[deck_pos] - 1; + } + } + + // Find number of tables, check for consistency. + enum { Uninitialized = -1 }; + int num_tables = Uninitialized; + if (phase_usage_.phase_used[Aqua]) { + const SWOF::table_t& swof_table = deck.getSWOF().swof_; + num_tables = swof_table.size(); + if (num_tables < satfuncs_expected) { + THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected); + } + } + if (phase_usage_.phase_used[Vapour]) { + const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; + int num_sgof_tables = sgof_table.size(); + if (num_sgof_tables < satfuncs_expected) { + THROW("Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected); + } + if (num_tables == Uninitialized) { + num_tables = num_sgof_tables; + } else if (num_tables != num_sgof_tables) { + THROW("Inconsistent number of tables in SWOF and SGOF."); + } + } + + // Initialize tables. + satfuncset_.resize(num_tables); + for (int table = 0; table < num_tables; ++table) { + satfuncset_[table].init(deck, table, phase_usage_); + } + } + + + + + /// \return P, the number of phases. + template + int SaturationPropsFromDeck::numPhases() const + { + return phase_usage_.num_phases; + } + + + + + /// Relative permeability. + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + /// \param[in] cells Array of n cell indices to be associated with the s values. + /// \param[out] kr Array of nP relperm values, array must be valid before calling. + /// \param[out] dkrds If non-null: array of nP^2 relperm derivative values, + /// array must be valid before calling. + /// The P^2 derivative matrix is + /// m_{ij} = \frac{dkr_i}{ds^j}, + /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) + template + void SaturationPropsFromDeck::relperm(const int n, + const double* s, + const int* cells, + double* kr, + double* dkrds) const + { + ASSERT (cells != 0); + + const int np = phase_usage_.num_phases; + if (dkrds) { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i); + } + } else { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + funcForCell(cells[i]).evalKr(s + np*i, kr + np*i); + } + } + } + + + + + /// Capillary pressure. + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + /// \param[in] cells Array of n cell indices to be associated with the s values. + /// \param[out] pc Array of nP capillary pressure values, array must be valid before calling. + /// \param[out] dpcds If non-null: array of nP^2 derivative values, + /// array must be valid before calling. + /// The P^2 derivative matrix is + /// m_{ij} = \frac{dpc_i}{ds^j}, + /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) + template + void SaturationPropsFromDeck::capPress(const int n, + const double* s, + const int* cells, + double* pc, + double* dpcds) const + { + ASSERT (cells != 0); + + const int np = phase_usage_.num_phases; + if (dpcds) { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i); + } + } else { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + funcForCell(cells[i]).evalPc(s + np*i, pc + np*i); + } + } + } + + + + + /// Obtain the range of allowable saturation values. + /// \param[in] n Number of data points. + /// \param[in] cells Array of n cell indices. + /// \param[out] smin Array of nP minimum s values, array must be valid before calling. + /// \param[out] smax Array of nP maximum s values, array must be valid before calling. + template + void SaturationPropsFromDeck::satRange(const int n, + const int* cells, + double* smin, + double* smax) const + { + ASSERT (cells != 0); + + const int np = phase_usage_.num_phases; + for (int i = 0; i < n; ++i) { + for (int p = 0; p < np; ++p) { + smin[np*i + p] = funcForCell(cells[i]).smin_[p]; + smax[np*i + p] = funcForCell(cells[i]).smax_[p]; + } + } + } + + + // Map the cell number to the correct function set. + template + const typename SaturationPropsFromDeck::Funcs& + SaturationPropsFromDeck::funcForCell(const int cell) const + { + return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]]; + } + + + +} // namespace Opm + +#endif // OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED diff --git a/opm/core/fluid/SaturationPropsInterface.hpp b/opm/core/fluid/SaturationPropsInterface.hpp new file mode 100644 index 000000000..e55995677 --- /dev/null +++ b/opm/core/fluid/SaturationPropsInterface.hpp @@ -0,0 +1,85 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED +#define OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED + +#include + + +namespace Opm +{ + + class SaturationPropsInterface : public BlackoilPhases + { + public: + /// Virtual destructor. + virtual ~SaturationPropsInterface() {}; + + /// \return P, the number of phases. + virtual int numPhases() const = 0; + + /// Relative permeability. + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + /// \param[out] kr Array of nP relperm values, array must be valid before calling. + /// \param[out] dkrds If non-null: array of nP^2 relperm derivative values, + /// array must be valid before calling. + /// The P^2 derivative matrix is + /// m_{ij} = \frac{dkr_i}{ds^j}, + /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) + virtual void relperm(const int n, + const double* s, + const int* cells, + double* kr, + double* dkrds) const = 0; + + /// Capillary pressure. + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + /// \param[out] pc Array of nP capillary pressure values, array must be valid before calling. + /// \param[out] dpcds If non-null: array of nP^2 derivative values, + /// array must be valid before calling. + /// The P^2 derivative matrix is + /// m_{ij} = \frac{dpc_i}{ds^j}, + /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) + virtual void capPress(const int n, + const double* s, + const int* cells, + double* pc, + double* dpcds) const = 0; + + /// Obtain the range of allowable saturation values. + /// \param[in] n Number of data points. + /// \param[out] smin Array of nP minimum s values, array must be valid before calling. + /// \param[out] smax Array of nP maximum s values, array must be valid before calling. + virtual void satRange(const int n, + const int* cells, + double* smin, + double* smax) const = 0; + + }; + + + +} // namespace Opm + + + +#endif // OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED From 0101e1f575b0ce50222a5b1d61bcb84d875d4196 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 4 Sep 2012 11:49:05 +0200 Subject: [PATCH 06/59] Updated doc to match new parameter names. New parameters are (default): pvt_tab_size (200) sat_tab_size (200) threephase_model ("simple") [also accepts "stone2"]. --- opm/core/fluid/BlackoilPropertiesFromDeck.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index 97a29d0af..80a8d8192 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -53,8 +53,9 @@ namespace Opm /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. /// \param[in] param Parameters. Accepted parameters include: - /// dead_tab_size (1025) number of uniform sample points for dead-oil pvt tables. - /// tab_size_kr (200) number of uniform sample points for saturation tables. + /// pvt_tab_size (200) number of uniform sample points for dead-oil pvt tables. + /// sat_tab_size (200) number of uniform sample points for saturation tables. + /// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2"). /// For both parameters, a 0 or negative value indicates that no spline fitting is to /// be done, and the input fluid data used directly for linear interpolation. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, From 876d23942c619854c81239454409957e3ec38741 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 4 Sep 2012 14:21:51 +0200 Subject: [PATCH 07/59] Use porevolume of last step in computation of gravity residual. --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 01882053b..7749218f2 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -349,7 +349,7 @@ namespace Opm gf[1] = gravflux[pos]; } s0 = tm.saturation_[cell]; - dtpv = tm.dt_/tm.porevolume0_[cell]; + dtpv = tm.dt_/tm.porevolume_[cell]; } double operator()(double s) const From ff78e358f09960c179ad83a7c16fb6c6cc2df0eb Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 4 Sep 2012 14:22:56 +0200 Subject: [PATCH 08/59] Removed unappropriate conversion from water saturation to complete saturation. --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 7749218f2..9bfd8b7bb 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -522,9 +522,7 @@ namespace Opm cells[c] = c; } mob_.resize(2*nc); - std::vector boths; - Opm::toBothSat(saturation, boths); - props_.relperm(cells.size(), &boths[0], &cells[0], &mob_[0], 0); + props_.relperm(cells.size(), &saturation[0], &cells[0], &mob_[0], 0); props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); for (int c = 0; c < nc; ++c) { From 55793cc9093d1352c464e9ca25f694ec48ca60a6 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 4 Sep 2012 15:01:14 +0200 Subject: [PATCH 09/59] Use same search interval in computation of flux and gravity residuals. --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 9bfd8b7bb..0364c4d3a 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -444,9 +444,8 @@ namespace Opm GravityResidual res(*this, cells, pos, gravflux); if (std::fabs(res(saturation_[cell])) > tol_) { int iters_used; - saturation_[cell] = RootFinder::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used); + saturation_[cell] = RootFinder::solve(res, 0.0, 1.0, maxit_, tol_, iters_used); } - saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]); mobility(saturation_[cell], cell, &mob_[2*cell]); } From 8c68fd03730f2f10ea98a532cccc27f208bb88f7 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 4 Sep 2012 15:13:55 +0200 Subject: [PATCH 10/59] Add initial guess in root solver for gravity residual. --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 0364c4d3a..9ea6271de 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -444,7 +444,7 @@ namespace Opm GravityResidual res(*this, cells, pos, gravflux); if (std::fabs(res(saturation_[cell])) > tol_) { int iters_used; - saturation_[cell] = RootFinder::solve(res, 0.0, 1.0, maxit_, tol_, iters_used); + saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); } mobility(saturation_[cell], cell, &mob_[2*cell]); } From 3d01a6009938861bd565b8426a9a88409a6f08fa Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 5 Sep 2012 10:10:02 +0200 Subject: [PATCH 11/59] Changed solvegravity interface. Fixed bug. --- .../TransportModelCompressibleTwophase.cpp | 30 +++++++++++-------- .../TransportModelCompressibleTwophase.hpp | 2 -- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 9ea6271de..d43f444a2 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -380,8 +380,8 @@ namespace Opm { double sat[2] = { s, 1.0 - s }; props_.relperm(1, sat, &cell, mob, 0); - mob[0] /= visc_[0]; - mob[1] /= visc_[1]; + mob[0] /= visc_[2*cell + 0]; + mob[1] /= visc_[2*cell + 1]; } @@ -407,7 +407,7 @@ namespace Opm { // Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f) // + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ] - // But b_w,i * rho_w,S = rho_w,i, which we conmpute with a call to props_.density(). + // But b_w,i * rho_w,S = rho_w,i, which we compute with a call to props_.density(). // We assume that we already have stored T_ij in trans_. // We also assume that the A_ matrices are updated from an earlier call to solve(). const int nc = grid_.number_of_cells; @@ -437,8 +437,8 @@ namespace Opm void TransportModelCompressibleTwophase::solveSingleCellGravity(const std::vector& cells, - const int pos, - const double* gravflux) + const int pos, + const double* gravflux) { const int cell = cells[pos]; GravityResidual res(*this, cells, pos, gravflux); @@ -505,8 +505,6 @@ namespace Opm void TransportModelCompressibleTwophase::solveGravity(const std::vector >& columns, - const double* pressure, - const double* porevolume0, const double dt, std::vector& saturation, std::vector& surfacevol) @@ -521,16 +519,22 @@ namespace Opm cells[c] = c; } mob_.resize(2*nc); - props_.relperm(cells.size(), &saturation[0], &cells[0], &mob_[0], 0); - props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); - for (int c = 0; c < nc; ++c) { - mob_[2*c + 0] /= visc_[2*c + 0]; - mob_[2*c + 1] /= visc_[2*c + 1]; + + // props_.relperm(cells.size(), &saturation[0], &cells[0], &mob_[0], 0); + + // props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); + // for (int c = 0; c < nc; ++c) { + // mob_[2*c + 0] /= visc_[2*c + 0]; + // mob_[2*c + 1] /= visc_[2*c + 1]; + // } + + const int np = props_.numPhases(); + for (int cell = 0; cell < nc; ++cell) { + mobility(saturation_[cell], cell, &mob_[np*cell]); } // Set up other variables. - porevolume0_ = porevolume0; dt_ = dt; toWaterSat(saturation, saturation_); diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 4ee93e0af..c0be99a9c 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -79,8 +79,6 @@ namespace Opm /// \param[in, out] saturation Phase saturations. /// \param[in, out] surfacevol Surface volume densities for each phase. void solveGravity(const std::vector >& columns, - const double* pressure, - const double* porevolume0, const double dt, std::vector& saturation, std::vector& surfacevol); From 0f91bc6a349a83fdcd3966e8a03eb62950f54085 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 5 Sep 2012 11:28:54 +0200 Subject: [PATCH 12/59] Massive whitespace cleanup: entire fluid subdir. --- opm/core/fluid/BlackoilPropertiesBasic.cpp | 72 +++---- opm/core/fluid/BlackoilPropertiesBasic.hpp | 20 +- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 10 +- opm/core/fluid/BlackoilPropertiesFromDeck.hpp | 12 +- .../fluid/BlackoilPropertiesInterface.hpp | 6 +- opm/core/fluid/IncompPropertiesBasic.cpp | 70 +++--- opm/core/fluid/IncompPropertiesBasic.hpp | 34 +-- opm/core/fluid/IncompPropertiesFromDeck.cpp | 42 ++-- opm/core/fluid/IncompPropertiesFromDeck.hpp | 14 +- opm/core/fluid/IncompPropertiesInterface.hpp | 6 +- opm/core/fluid/PvtPropertiesBasic.cpp | 112 +++++----- opm/core/fluid/PvtPropertiesBasic.hpp | 18 +- .../fluid/PvtPropertiesIncompFromDeck.cpp | 70 +++--- .../fluid/PvtPropertiesIncompFromDeck.hpp | 12 +- opm/core/fluid/RockBasic.hpp | 8 +- opm/core/fluid/RockCompressibility.cpp | 6 +- opm/core/fluid/RockFromDeck.cpp | 2 +- opm/core/fluid/RockFromDeck.hpp | 2 +- opm/core/fluid/SaturationPropsBasic.cpp | 204 +++++++++--------- opm/core/fluid/SaturationPropsBasic.hpp | 22 +- opm/core/fluid/SaturationPropsFromDeck.hpp | 10 +- .../fluid/SaturationPropsFromDeck_impl.hpp | 20 +- opm/core/fluid/SaturationPropsInterface.hpp | 4 +- .../fluid/blackoil/BlackoilPvtProperties.hpp | 6 +- 24 files changed, 391 insertions(+), 391 deletions(-) diff --git a/opm/core/fluid/BlackoilPropertiesBasic.cpp b/opm/core/fluid/BlackoilPropertiesBasic.cpp index be8211c15..2496d98ef 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.cpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.cpp @@ -26,20 +26,20 @@ namespace Opm { BlackoilPropertiesBasic::BlackoilPropertiesBasic(const parameter::ParameterGroup& param, - const int dim, - const int num_cells) + const int dim, + const int num_cells) { - double poro = param.getDefault("porosity", 1.0); - using namespace Opm::unit; - using namespace Opm::prefix; - double perm = param.getDefault("permeability", 100.0)*milli*darcy; + double poro = param.getDefault("porosity", 1.0); + using namespace Opm::unit; + using namespace Opm::prefix; + double perm = param.getDefault("permeability", 100.0)*milli*darcy; rock_.init(dim, num_cells, poro, perm); - pvt_.init(param); + pvt_.init(param); satprops_.init(param); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } } BlackoilPropertiesBasic::~BlackoilPropertiesBasic() @@ -90,11 +90,11 @@ namespace Opm /// \param[out] dmudp If non-null: array of nP viscosity derivative values, /// array must be valid before calling. void BlackoilPropertiesBasic::viscosity(const int n, - const double* p, - const double* z, - const int* /*cells*/, - double* mu, - double* dmudp) const + const double* p, + const double* z, + const int* /*cells*/, + double* mu, + double* dmudp) const { if (dmudp) { THROW("BlackoilPropertiesBasic::viscosity() -- derivatives of viscosity not yet implemented."); @@ -114,16 +114,16 @@ namespace Opm /// array must be valid before calling. The matrices are output /// in Fortran order. void BlackoilPropertiesBasic::matrix(const int n, - const double* /*p*/, - const double* /*z*/, - const int* /*cells*/, - double* A, - double* dAdp) const + const double* /*p*/, + const double* /*z*/, + const int* /*cells*/, + double* A, + double* dAdp) const { - const int np = numPhases(); - ASSERT(np <= 2); - double B[2]; // Must be enough since component classes do not handle more than 2. - pvt_.B(1, 0, 0, B); + const int np = numPhases(); + ASSERT(np <= 2); + double B[2]; // Must be enough since component classes do not handle more than 2. + pvt_.B(1, 0, 0, B); // Compute A matrix // #pragma omp parallel for for (int i = 0; i < n; ++i) { @@ -152,8 +152,8 @@ namespace Opm /// of a call to the method matrix(). /// \param[out] rho Array of nP density values, array must be valid before calling. void BlackoilPropertiesBasic::density(const int n, - const double* A, - double* rho) const + const double* A, + double* rho) const { const int np = numPhases(); const double* sdens = pvt_.surfaceDensities(); @@ -186,10 +186,10 @@ namespace Opm /// m_{ij} = \frac{dkr_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void BlackoilPropertiesBasic::relperm(const int n, - const double* s, - const int* /*cells*/, - double* kr, - double* dkrds) const + const double* s, + const int* /*cells*/, + double* kr, + double* dkrds) const { satprops_.relperm(n, s, kr, dkrds); } @@ -205,10 +205,10 @@ namespace Opm /// m_{ij} = \frac{dpc_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void BlackoilPropertiesBasic::capPress(const int n, - const double* s, - const int* /*cells*/, - double* pc, - double* dpcds) const + const double* s, + const int* /*cells*/, + double* pc, + double* dpcds) const { satprops_.capPress(n, s, pc, dpcds); } @@ -226,7 +226,7 @@ namespace Opm double* smin, double* smax) const { - satprops_.satRange(n, smin, smax); + satprops_.satRange(n, smin, smax); } diff --git a/opm/core/fluid/BlackoilPropertiesBasic.hpp b/opm/core/fluid/BlackoilPropertiesBasic.hpp index d879036ef..44a76013b 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.hpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.hpp @@ -35,16 +35,16 @@ namespace Opm { public: /// Construct from parameters. - /// The following parameters are accepted (defaults): - /// num_phases (2) Must be 1 or 2. - /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". - /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 - /// mu1 [mu2, mu3] (1.0) Viscosity in cP - /// porosity (1.0) Porosity - /// permeability (100.0) Permeability in mD + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1 or 2. + /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". + /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 + /// mu1 [mu2, mu3] (1.0) Viscosity in cP + /// porosity (1.0) Porosity + /// permeability (100.0) Permeability in mD BlackoilPropertiesBasic(const parameter::ParameterGroup& param, - const int dim, - const int num_cells); + const int dim, + const int num_cells); /// Destructor. virtual ~BlackoilPropertiesBasic(); @@ -151,7 +151,7 @@ namespace Opm double* dpcds) const; - /// Obtain the range of allowable saturation values. +/// Obtain the range of allowable saturation values. /// In cell cells[i], saturation of phase p is allowed to be /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \param[in] n Number of data points. diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 7cdae4ae2..e2457f39b 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -33,10 +33,10 @@ namespace Opm satprops_.reset(ptr); ptr->init(deck, grid, 200); - if (pvt_.numPhases() != satprops_->numPhases()) { - THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } + if (pvt_.numPhases() != satprops_->numPhases()) { + THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); + } } BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, @@ -319,7 +319,7 @@ namespace Opm double* smin, double* smax) const { - satprops_->satRange(n, cells, smin, smax); + satprops_->satRange(n, cells, smin, smax); } diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index aceb534b8..eb8b947b4 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -41,15 +41,15 @@ namespace Opm public: /// Initialize from deck and grid. /// \param[in] deck Deck input parser - /// \param[in] grid Grid to which property object applies, needed for the + /// \param[in] grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + const UnstructuredGrid& grid); /// Initialize from deck, grid and parameters. /// \param[in] deck Deck input parser - /// \param[in] grid Grid to which property object applies, needed for the + /// \param[in] grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. /// \param[in] param Parameters. Accepted parameters include: @@ -167,9 +167,9 @@ namespace Opm double* dpcds) const; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \param[in] n Number of data points. /// \param[in] cells Array of n cell indices. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. diff --git a/opm/core/fluid/BlackoilPropertiesInterface.hpp b/opm/core/fluid/BlackoilPropertiesInterface.hpp index dd8a857fb..501c8a406 100644 --- a/opm/core/fluid/BlackoilPropertiesInterface.hpp +++ b/opm/core/fluid/BlackoilPropertiesInterface.hpp @@ -138,9 +138,9 @@ namespace Opm double* dpcds) const = 0; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \param[in] n Number of data points. /// \param[in] cells Array of n cell indices. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. diff --git a/opm/core/fluid/IncompPropertiesBasic.cpp b/opm/core/fluid/IncompPropertiesBasic.cpp index 0ce3c88aa..ab68017c1 100644 --- a/opm/core/fluid/IncompPropertiesBasic.cpp +++ b/opm/core/fluid/IncompPropertiesBasic.cpp @@ -28,22 +28,22 @@ namespace Opm { IncompPropertiesBasic::IncompPropertiesBasic(const parameter::ParameterGroup& param, - const int dim, - const int num_cells) + const int dim, + const int num_cells) { - double poro = param.getDefault("porosity", 1.0); - using namespace Opm::unit; - using namespace Opm::prefix; - double perm = param.getDefault("permeability", 100.0)*milli*darcy; + double poro = param.getDefault("porosity", 1.0); + using namespace Opm::unit; + using namespace Opm::prefix; + double perm = param.getDefault("permeability", 100.0)*milli*darcy; rock_.init(dim, num_cells, poro, perm); - pvt_.init(param); + pvt_.init(param); satprops_.init(param); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } - viscosity_.resize(pvt_.numPhases()); - pvt_.mu(1, 0, 0, &viscosity_[0]); + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } + viscosity_.resize(pvt_.numPhases()); + pvt_.mu(1, 0, 0, &viscosity_[0]); } IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases, @@ -56,14 +56,14 @@ namespace Opm const int num_cells) { rock_.init(dim, num_cells, por, perm); - pvt_.init(num_phases, rho, mu); + pvt_.init(num_phases, rho, mu); satprops_.init(num_phases, relpermfunc); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } - viscosity_.resize(pvt_.numPhases()); - pvt_.mu(1, 0, 0, &viscosity_[0]); + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } + viscosity_.resize(pvt_.numPhases()); + pvt_.mu(1, 0, 0, &viscosity_[0]); } IncompPropertiesBasic::~IncompPropertiesBasic() @@ -109,7 +109,7 @@ namespace Opm /// \return Array of P viscosity values. const double* IncompPropertiesBasic::viscosity() const { - return &viscosity_[0]; + return &viscosity_[0]; } /// \return Array of P density values. @@ -117,7 +117,7 @@ namespace Opm { // No difference between reservoir and surface densities // modelled by this class. - return pvt_.surfaceDensities(); + return pvt_.surfaceDensities(); } /// \return Array of P density values. @@ -125,7 +125,7 @@ namespace Opm { // No difference between reservoir and surface densities // modelled by this class. - return pvt_.surfaceDensities(); + return pvt_.surfaceDensities(); } /// \param[in] n Number of data points. @@ -138,10 +138,10 @@ namespace Opm /// m_{ij} = \frac{dkr_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) void IncompPropertiesBasic::relperm(const int n, - const double* s, - const int* /*cells*/, - double* kr, - double* dkrds) const + const double* s, + const int* /*cells*/, + double* kr, + double* dkrds) const { satprops_.relperm(n, s, kr, dkrds); } @@ -157,10 +157,10 @@ namespace Opm /// m_{ij} = \frac{dpc_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) void IncompPropertiesBasic::capPress(const int n, - const double* s, - const int* /*cells*/, - double* pc, - double* dpcds) const + const double* s, + const int* /*cells*/, + double* pc, + double* dpcds) const { satprops_.capPress(n, s, pc, dpcds); } @@ -174,11 +174,11 @@ namespace Opm /// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling. void IncompPropertiesBasic::satRange(const int n, - const int* /*cells*/, - double* smin, - double* smax) const + const int* /*cells*/, + double* smin, + double* smax) const { - satprops_.satRange(n, smin, smax); + satprops_.satRange(n, smin, smax); } } // namespace Opm diff --git a/opm/core/fluid/IncompPropertiesBasic.hpp b/opm/core/fluid/IncompPropertiesBasic.hpp index b90b0876f..72c39b1b0 100644 --- a/opm/core/fluid/IncompPropertiesBasic.hpp +++ b/opm/core/fluid/IncompPropertiesBasic.hpp @@ -42,29 +42,29 @@ namespace Opm { public: /// Construct from parameters. - /// The following parameters are accepted (defaults): - /// num_phases (2) Must be 1 or 2. - /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". - /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 - /// mu1 [mu2, mu3] (1.0) Viscosity in cP - /// porosity (1.0) Porosity - /// permeability (100.0) Permeability in mD + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1 or 2. + /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". + /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 + /// mu1 [mu2, mu3] (1.0) Viscosity in cP + /// porosity (1.0) Porosity + /// permeability (100.0) Permeability in mD IncompPropertiesBasic(const parameter::ParameterGroup& param, - const int dim, - const int num_cells); + const int dim, + const int num_cells); /// Construct from arguments a basic two phase fluid. IncompPropertiesBasic(const int num_phases, const SaturationPropsBasic::RelPermFunc& relpermfunc, const std::vector& rho, - const std::vector& mu, + const std::vector& mu, const double porosity, const double permeability, const int dim, - const int num_cells); + const int num_cells); - /// Destructor. + /// Destructor. virtual ~IncompPropertiesBasic(); // ---- Rock interface ---- @@ -132,9 +132,9 @@ namespace Opm double* dpcds) const; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \param[in] n Number of data points. /// \param[in] cells Array of n cell indices. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. @@ -145,9 +145,9 @@ namespace Opm double* smax) const; private: RockBasic rock_; - PvtPropertiesBasic pvt_; + PvtPropertiesBasic pvt_; SaturationPropsBasic satprops_; - std::vector viscosity_; + std::vector viscosity_; }; diff --git a/opm/core/fluid/IncompPropertiesFromDeck.cpp b/opm/core/fluid/IncompPropertiesFromDeck.cpp index 45e49791f..9aa690980 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.cpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.cpp @@ -27,15 +27,15 @@ namespace Opm { IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid) + const UnstructuredGrid& grid) { rock_.init(deck, grid); - pvt_.init(deck); + pvt_.init(deck); satprops_.init(deck, grid, 200); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } } IncompPropertiesFromDeck::~IncompPropertiesFromDeck() @@ -81,19 +81,19 @@ namespace Opm /// \return Array of P viscosity values. const double* IncompPropertiesFromDeck::viscosity() const { - return pvt_.viscosity(); + return pvt_.viscosity(); } /// \return Array of P density values. const double* IncompPropertiesFromDeck::density() const { - return pvt_.reservoirDensities(); + return pvt_.reservoirDensities(); } /// \return Array of P density values. const double* IncompPropertiesFromDeck::surfaceDensity() const { - return pvt_.surfaceDensities(); + return pvt_.surfaceDensities(); } /// \param[in] n Number of data points. @@ -106,10 +106,10 @@ namespace Opm /// m_{ij} = \frac{dkr_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) void IncompPropertiesFromDeck::relperm(const int n, - const double* s, - const int* cells, - double* kr, - double* dkrds) const + const double* s, + const int* cells, + double* kr, + double* dkrds) const { satprops_.relperm(n, s, cells, kr, dkrds); } @@ -125,10 +125,10 @@ namespace Opm /// m_{ij} = \frac{dpc_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) void IncompPropertiesFromDeck::capPress(const int n, - const double* s, - const int* cells, - double* pc, - double* dpcds) const + const double* s, + const int* cells, + double* pc, + double* dpcds) const { satprops_.capPress(n, s, cells, pc, dpcds); } @@ -142,11 +142,11 @@ namespace Opm /// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling. void IncompPropertiesFromDeck::satRange(const int n, - const int* cells, - double* smin, - double* smax) const + const int* cells, + double* smin, + double* smax) const { - satprops_.satRange(n, cells, smin, smax); + satprops_.satRange(n, cells, smin, smax); } } // namespace Opm diff --git a/opm/core/fluid/IncompPropertiesFromDeck.hpp b/opm/core/fluid/IncompPropertiesFromDeck.hpp index 8f917d3b9..6680eaa49 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.hpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.hpp @@ -47,13 +47,13 @@ namespace Opm public: /// Initialize from deck and grid. /// \param deck Deck input parser - /// \param grid Grid to which property object applies, needed for the + /// \param grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. IncompPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + const UnstructuredGrid& grid); - /// Destructor. + /// Destructor. virtual ~IncompPropertiesFromDeck(); // ---- Rock interface ---- @@ -121,9 +121,9 @@ namespace Opm double* dpcds) const; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \param[in] n Number of data points. /// \param[in] cells Array of n cell indices. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. @@ -134,7 +134,7 @@ namespace Opm double* smax) const; private: RockFromDeck rock_; - PvtPropertiesIncompFromDeck pvt_; + PvtPropertiesIncompFromDeck pvt_; SaturationPropsFromDeck satprops_; }; diff --git a/opm/core/fluid/IncompPropertiesInterface.hpp b/opm/core/fluid/IncompPropertiesInterface.hpp index a5025e8a1..4f8bfb516 100644 --- a/opm/core/fluid/IncompPropertiesInterface.hpp +++ b/opm/core/fluid/IncompPropertiesInterface.hpp @@ -109,9 +109,9 @@ namespace Opm double* pc, double* dpcds) const = 0; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \param[in] n Number of data points. /// \param[in] cells Array of n cell indices. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. diff --git a/opm/core/fluid/PvtPropertiesBasic.cpp b/opm/core/fluid/PvtPropertiesBasic.cpp index 3d7d7bb2d..5220502ad 100644 --- a/opm/core/fluid/PvtPropertiesBasic.cpp +++ b/opm/core/fluid/PvtPropertiesBasic.cpp @@ -34,41 +34,41 @@ namespace Opm void PvtPropertiesBasic::init(const parameter::ParameterGroup& param) { - int num_phases = param.getDefault("num_phases", 2); - if (num_phases > 3 || num_phases < 1) { - THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); - } - density_.resize(num_phases); - viscosity_.resize(num_phases); - // We currently do not allow the user to set B. - formation_volume_factor_.clear(); - formation_volume_factor_.resize(num_phases, 1.0); + int num_phases = param.getDefault("num_phases", 2); + if (num_phases > 3 || num_phases < 1) { + THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); + } + density_.resize(num_phases); + viscosity_.resize(num_phases); + // We currently do not allow the user to set B. + formation_volume_factor_.clear(); + formation_volume_factor_.resize(num_phases, 1.0); - // Setting mu and rho from parameters - using namespace Opm::prefix; - using namespace Opm::unit; - const double kgpm3 = kilogram/cubic(meter); - const double cP = centi*Poise; - std::string rname[3] = { "rho1", "rho2", "rho3" }; - double rdefault[3] = { 1.0e3, 1.0e3, 1.0e3 }; - std::string vname[3] = { "mu1", "mu2", "mu3" }; - double vdefault[3] = { 1.0, 1.0, 1.0 }; - for (int phase = 0; phase < num_phases; ++phase) { - density_[phase] = kgpm3*param.getDefault(rname[phase], rdefault[phase]); - viscosity_[phase] = cP*param.getDefault(vname[phase], vdefault[phase]); - } + // Setting mu and rho from parameters + using namespace Opm::prefix; + using namespace Opm::unit; + const double kgpm3 = kilogram/cubic(meter); + const double cP = centi*Poise; + std::string rname[3] = { "rho1", "rho2", "rho3" }; + double rdefault[3] = { 1.0e3, 1.0e3, 1.0e3 }; + std::string vname[3] = { "mu1", "mu2", "mu3" }; + double vdefault[3] = { 1.0, 1.0, 1.0 }; + for (int phase = 0; phase < num_phases; ++phase) { + density_[phase] = kgpm3*param.getDefault(rname[phase], rdefault[phase]); + viscosity_[phase] = cP*param.getDefault(vname[phase], vdefault[phase]); + } } void PvtPropertiesBasic::init(const int num_phases, const std::vector& rho, const std::vector& visc) { - if (num_phases > 3 || num_phases < 1) { - THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); - } - // We currently do not allow the user to set B. - formation_volume_factor_.clear(); - formation_volume_factor_.resize(num_phases, 1.0); + if (num_phases > 3 || num_phases < 1) { + THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); + } + // We currently do not allow the user to set B. + formation_volume_factor_.clear(); + formation_volume_factor_.resize(num_phases, 1.0); density_ = rho; viscosity_ = visc; } @@ -87,69 +87,69 @@ namespace Opm void PvtPropertiesBasic::mu(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_mu) const + const double* /*p*/, + const double* /*z*/, + double* output_mu) const { - const int np = numPhases(); + const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_mu[np*i + phase] = viscosity_[phase]; } - } + } } void PvtPropertiesBasic::B(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_B) const + const double* /*p*/, + const double* /*z*/, + double* output_B) const { - const int np = numPhases(); + const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[np*i + phase] = formation_volume_factor_[phase]; } - } + } } void PvtPropertiesBasic::dBdp(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_B, - double* output_dBdp) const + const double* /*p*/, + const double* /*z*/, + double* output_B, + double* output_dBdp) const { - const int np = numPhases(); + const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[np*i + phase] = formation_volume_factor_[phase]; output_dBdp[np*i + phase] = 0.0; } - } + } } void PvtPropertiesBasic::R(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_R) const + const double* /*p*/, + const double* /*z*/, + double* output_R) const { - const int np = numPhases(); - std::fill(output_R, output_R + n*np, 0.0); + const int np = numPhases(); + std::fill(output_R, output_R + n*np, 0.0); } void PvtPropertiesBasic::dRdp(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_R, - double* output_dRdp) const + const double* /*p*/, + const double* /*z*/, + double* output_R, + double* output_dRdp) const { - const int np = numPhases(); - std::fill(output_R, output_R + n*np, 0.0); - std::fill(output_dRdp, output_dRdp + n*np, 0.0); + const int np = numPhases(); + std::fill(output_R, output_R + n*np, 0.0); + std::fill(output_dRdp, output_dRdp + n*np, 0.0); } } // namespace Opm diff --git a/opm/core/fluid/PvtPropertiesBasic.hpp b/opm/core/fluid/PvtPropertiesBasic.hpp index d1a3e9f3d..3e30e807c 100644 --- a/opm/core/fluid/PvtPropertiesBasic.hpp +++ b/opm/core/fluid/PvtPropertiesBasic.hpp @@ -38,11 +38,11 @@ namespace Opm PvtPropertiesBasic(); /// Initialize from parameters. - /// The following parameters are accepted (defaults): - /// num_phases (2) Must be 1, 2 or 3. - /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 - /// mu1 [mu2, mu3] (1.0) Viscosity in cP - void init(const parameter::ParameterGroup& param); + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1, 2 or 3. + /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 + /// mu1 [mu2, mu3] (1.0) Viscosity in cP + void init(const parameter::ParameterGroup& param); /// Initialize from arguments. /// Basic multi phase fluid pvt properties. @@ -55,7 +55,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of size numPhases(). - const double* surfaceDensities() const; + const double* surfaceDensities() const; /// Viscosity as a function of p and z. void mu(const int n, @@ -90,9 +90,9 @@ namespace Opm double* output_dRdp) const; private: - std::vector density_; - std::vector viscosity_; - std::vector formation_volume_factor_; + std::vector density_; + std::vector viscosity_; + std::vector formation_volume_factor_; }; } diff --git a/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp b/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp index 62cc1b231..34a4ba002 100644 --- a/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp +++ b/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp @@ -38,54 +38,54 @@ namespace Opm { typedef std::vector > > table_t; // If we need multiple regions, this class and the SinglePvt* classes must change. - int region_number = 0; + int region_number = 0; PhaseUsage phase_usage = phaseUsageFromDeck(deck); - if (phase_usage.phase_used[PhaseUsage::Vapour] || - !phase_usage.phase_used[PhaseUsage::Aqua] || - !phase_usage.phase_used[PhaseUsage::Liquid]) { - THROW("PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n"); - } + if (phase_usage.phase_used[PhaseUsage::Vapour] || + !phase_usage.phase_used[PhaseUsage::Aqua] || + !phase_usage.phase_used[PhaseUsage::Liquid]) { + THROW("PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n"); + } - // Surface densities. Accounting for different orders in eclipse and our code. - if (deck.hasField("DENSITY")) { - const std::vector& d = deck.getDENSITY().densities_[region_number]; - enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; - surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water]; - surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil]; - } else { - THROW("Input is missing DENSITY\n"); - } + // Surface densities. Accounting for different orders in eclipse and our code. + if (deck.hasField("DENSITY")) { + const std::vector& d = deck.getDENSITY().densities_[region_number]; + enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; + surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water]; + surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil]; + } else { + THROW("Input is missing DENSITY\n"); + } // Make reservoir densities the same as surface densities initially. // We will modify them with formation volume factors if found. reservoir_density_ = surface_density_; // Water viscosity. - if (deck.hasField("PVTW")) { - const std::vector& pvtw = deck.getPVTW().pvtw_[region_number]; - if (pvtw[2] != 0.0 || pvtw[4] != 0.0) { - MESSAGE("Compressibility effects in PVTW are ignored."); - } + if (deck.hasField("PVTW")) { + const std::vector& pvtw = deck.getPVTW().pvtw_[region_number]; + if (pvtw[2] != 0.0 || pvtw[4] != 0.0) { + MESSAGE("Compressibility effects in PVTW are ignored."); + } reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1]; - viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3]; - } else { - // Eclipse 100 default. - // viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise; - THROW("Input is missing PVTW\n"); - } + viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3]; + } else { + // Eclipse 100 default. + // viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise; + THROW("Input is missing PVTW\n"); + } // Oil viscosity. - if (deck.hasField("PVCDO")) { - const std::vector& pvcdo = deck.getPVCDO().pvcdo_[region_number]; - if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) { - MESSAGE("Compressibility effects in PVCDO are ignored."); - } + if (deck.hasField("PVCDO")) { + const std::vector& pvcdo = deck.getPVCDO().pvcdo_[region_number]; + if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) { + MESSAGE("Compressibility effects in PVCDO are ignored."); + } reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1]; - viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3]; - } else { - THROW("Input is missing PVCDO\n"); - } + viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3]; + } else { + THROW("Input is missing PVCDO\n"); + } } const double* PvtPropertiesIncompFromDeck::surfaceDensities() const diff --git a/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp b/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp index acbabaa49..01bbeed1d 100644 --- a/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp +++ b/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp @@ -39,14 +39,14 @@ namespace Opm PvtPropertiesIncompFromDeck(); /// Initialize from deck. - void init(const EclipseGridParser& deck); + void init(const EclipseGridParser& deck); /// Number of active phases. int numPhases() const; /// Densities of stock components at surface conditions. /// \return Array of size numPhases(). - const double* surfaceDensities() const; + const double* surfaceDensities() const; /// Densities of stock components at reservoir conditions. /// Note: a reasonable question to ask is why there can be @@ -58,15 +58,15 @@ namespace Opm /// reporting and using data given in terms of surface values, /// we need to handle this difference. /// \return Array of size numPhases(). - const double* reservoirDensities() const; + const double* reservoirDensities() const; /// Viscosities. const double* viscosity() const; private: - std::tr1::array surface_density_; - std::tr1::array reservoir_density_; - std::tr1::array viscosity_; + std::tr1::array surface_density_; + std::tr1::array reservoir_density_; + std::tr1::array viscosity_; }; } diff --git a/opm/core/fluid/RockBasic.hpp b/opm/core/fluid/RockBasic.hpp index 965637807..2adc3ffd6 100644 --- a/opm/core/fluid/RockBasic.hpp +++ b/opm/core/fluid/RockBasic.hpp @@ -35,9 +35,9 @@ namespace Opm /// Initialize with homogenous porosity and permeability. void init(const int dimensions, - const int num_cells, - const double poro, - const double perm); + const int num_cells, + const double poro, + const double perm); /// \return D, the number of spatial dimensions. int numDimensions() const @@ -66,7 +66,7 @@ namespace Opm } private: - int dimensions_; + int dimensions_; std::vector porosity_; std::vector permeability_; }; diff --git a/opm/core/fluid/RockCompressibility.cpp b/opm/core/fluid/RockCompressibility.cpp index 104861013..af65d852b 100644 --- a/opm/core/fluid/RockCompressibility.cpp +++ b/opm/core/fluid/RockCompressibility.cpp @@ -69,8 +69,8 @@ namespace Opm const double cpnorm = rock_comp_*(pressure - pref_); return (1.0 + cpnorm + 0.5*cpnorm*cpnorm); } else { - // return Opm::linearInterpolation(p_, poromult_, pressure); - return Opm::linearInterpolationExtrap(p_, poromult_, pressure); + // return Opm::linearInterpolation(p_, poromult_, pressure); + return Opm::linearInterpolationExtrap(p_, poromult_, pressure); } } @@ -81,7 +81,7 @@ namespace Opm } else { //const double poromult = Opm::linearInterpolation(p_, poromult_, pressure); //const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure); - const double poromult = Opm::linearInterpolationExtrap(p_, poromult_, pressure); + const double poromult = Opm::linearInterpolationExtrap(p_, poromult_, pressure); const double dporomultdp = Opm::linearInterpolationDerivativeExtrap(p_, poromult_, pressure); return dporomultdp/poromult; diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 7ba8903fd..ec7db9a6a 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -51,7 +51,7 @@ namespace Opm /// Initialize from deck and cell mapping. /// \param deck Deck input parser - /// \param grid grid to which property object applies, needed for the + /// \param grid grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. void RockFromDeck::init(const EclipseGridParser& deck, diff --git a/opm/core/fluid/RockFromDeck.hpp b/opm/core/fluid/RockFromDeck.hpp index 65027e425..63e71a6e4 100644 --- a/opm/core/fluid/RockFromDeck.hpp +++ b/opm/core/fluid/RockFromDeck.hpp @@ -37,7 +37,7 @@ namespace Opm /// Initialize from deck and grid. /// \param deck Deck input parser - /// \param grid Grid to which property object applies, needed for the + /// \param grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. void init(const EclipseGridParser& deck, diff --git a/opm/core/fluid/SaturationPropsBasic.cpp b/opm/core/fluid/SaturationPropsBasic.cpp index 760f934d4..10ddf23dd 100644 --- a/opm/core/fluid/SaturationPropsBasic.cpp +++ b/opm/core/fluid/SaturationPropsBasic.cpp @@ -29,64 +29,64 @@ namespace Opm namespace { - struct KrFunConstant - { - double kr(double) - { - return 1.0; - } - double dkrds(double) - { - return 0.0; - } - }; + struct KrFunConstant + { + double kr(double) + { + return 1.0; + } + double dkrds(double) + { + return 0.0; + } + }; - struct KrFunLinear - { - double kr(double s) - { - return s; - } - double dkrds(double) - { - return 1.0; - } - }; + struct KrFunLinear + { + double kr(double s) + { + return s; + } + double dkrds(double) + { + return 1.0; + } + }; - struct KrFunQuadratic - { - double kr(double s) - { - return s*s; - } - double dkrds(double s) - { - return 2.0*s; - } - }; + struct KrFunQuadratic + { + double kr(double s) + { + return s*s; + } + double dkrds(double s) + { + return 2.0*s; + } + }; - template - static inline void evalAllKrDeriv(const int n, const int np, - const double* s, double* kr, double* dkrds, Fun fun) - { - if (dkrds == 0) { + template + static inline void evalAllKrDeriv(const int n, const int np, + const double* s, double* kr, double* dkrds, Fun fun) + { + if (dkrds == 0) { // #pragma omp parallel for - for (int i = 0; i < n*np; ++i) { - kr[i] = fun.kr(s[i]); - } - return; - } + for (int i = 0; i < n*np; ++i) { + kr[i] = fun.kr(s[i]); + } + return; + } // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0); - for (int phase = 0; phase < np; ++phase) { - kr[i*np + phase] = fun.kr(s[i*np + phase]); - // Only diagonal elements in derivative. - dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]); - } - } - } + for (int i = 0; i < n; ++i) { + std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0); + for (int phase = 0; phase < np; ++phase) { + kr[i*np + phase] = fun.kr(s[i*np + phase]); + // Only diagonal elements in derivative. + dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]); + } + } + } } // anon namespace @@ -109,25 +109,25 @@ namespace Opm /// Initialize from parameters. void SaturationPropsBasic::init(const parameter::ParameterGroup& param) { - int num_phases = param.getDefault("num_phases", 2); - if (num_phases > 2 || num_phases < 1) { - THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases); - } + int num_phases = param.getDefault("num_phases", 2); + if (num_phases > 2 || num_phases < 1) { + THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases); + } num_phases_ = num_phases; - //std::string rpf = param.getDefault("relperm_func", std::string("Unset")); - std::string rpf = param.getDefault("relperm_func", std::string("Linear")); - if (rpf == "Constant") { - relperm_func_ = Constant; - if(num_phases!=1){ - THROW("Constant relperm with more than one phase???"); - } - } else if (rpf == "Linear") { - relperm_func_ = Linear; - } else if (rpf == "Quadratic") { - relperm_func_ = Quadratic; - } else { - THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf); - } + //std::string rpf = param.getDefault("relperm_func", std::string("Unset")); + std::string rpf = param.getDefault("relperm_func", std::string("Linear")); + if (rpf == "Constant") { + relperm_func_ = Constant; + if(num_phases!=1){ + THROW("Constant relperm with more than one phase???"); + } + } else if (rpf == "Linear") { + relperm_func_ = Linear; + } else if (rpf == "Quadratic") { + relperm_func_ = Quadratic; + } else { + THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf); + } } @@ -136,7 +136,7 @@ namespace Opm /// \return P, the number of phases. int SaturationPropsBasic::numPhases() const { - return num_phases_; + return num_phases_; } @@ -152,29 +152,29 @@ namespace Opm /// m_{ij} = \frac{dkr_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void SaturationPropsBasic::relperm(const int n, - const double* s, - double* kr, - double* dkrds) const + const double* s, + double* kr, + double* dkrds) const { - switch (relperm_func_) { - case Constant: - { - evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant()); - break; - } - case Linear: - { - evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear()); - break; - } - case Quadratic: - { - evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic()); - break; - } - default: - THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_); - } + switch (relperm_func_) { + case Constant: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant()); + break; + } + case Linear: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear()); + break; + } + case Quadratic: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic()); + break; + } + default: + THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_); + } } @@ -190,13 +190,13 @@ namespace Opm /// m_{ij} = \frac{dpc_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void SaturationPropsBasic::capPress(const int n, - const double* /*s*/, - double* pc, - double* dpcds) const + const double* /*s*/, + double* pc, + double* dpcds) const { - std::fill(pc, pc + num_phases_*n, 0.0); + std::fill(pc, pc + num_phases_*n, 0.0); if (dpcds) { - std::fill(dpcds, dpcds + num_phases_*num_phases_*n, 0.0); + std::fill(dpcds, dpcds + num_phases_*num_phases_*n, 0.0); } } @@ -207,11 +207,11 @@ namespace Opm /// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling. void SaturationPropsBasic::satRange(const int n, - double* smin, - double* smax) const + double* smin, + double* smax) const { - std::fill(smin, smin + num_phases_*n, 0.0); - std::fill(smax, smax + num_phases_*n, 1.0); + std::fill(smin, smin + num_phases_*n, 0.0); + std::fill(smax, smax + num_phases_*n, 1.0); } diff --git a/opm/core/fluid/SaturationPropsBasic.hpp b/opm/core/fluid/SaturationPropsBasic.hpp index 789e0a351..2f90672d8 100644 --- a/opm/core/fluid/SaturationPropsBasic.hpp +++ b/opm/core/fluid/SaturationPropsBasic.hpp @@ -40,16 +40,16 @@ namespace Opm SaturationPropsBasic(); /// Initialize from parameters. - /// The following parameters are accepted (defaults): - /// num_phases (2) Must be 1 or 2. - /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1 or 2. + /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". void init(const parameter::ParameterGroup& param); - enum RelPermFunc { Constant, Linear, Quadratic }; + enum RelPermFunc { Constant, Linear, Quadratic }; /// Initialize from arguments a basic Saturation property. void init(const int num_phases, - const RelPermFunc& relperm_func) + const RelPermFunc& relperm_func) { num_phases_ = num_phases; relperm_func_ = relperm_func; @@ -86,18 +86,18 @@ namespace Opm double* pc, double* dpcds) const; - /// Obtain the range of allowable saturation values. + /// Obtain the range of allowable saturation values. /// \param[in] n Number of data points. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling. - void satRange(const int n, - double* smin, - double* smax) const; + void satRange(const int n, + double* smin, + double* smax) const; private: - int num_phases_; - RelPermFunc relperm_func_; + int num_phases_; + RelPermFunc relperm_func_; }; diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index 42461985c..0a0599ec8 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -50,7 +50,7 @@ namespace Opm /// Initialize from deck and grid. /// \param[in] deck Deck input parser - /// \param[in] grid Grid to which property object applies, needed for the + /// \param[in] grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. /// \param[in] samples Number of uniform sample points for saturation tables. @@ -92,14 +92,14 @@ namespace Opm double* pc, double* dpcds) const; - /// Obtain the range of allowable saturation values. + /// Obtain the range of allowable saturation values. /// \param[in] n Number of data points. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling. - void satRange(const int n, + void satRange(const int n, const int* cells, - double* smin, - double* smax) const; + double* smin, + double* smax) const; private: PhaseUsage phase_usage_; diff --git a/opm/core/fluid/SaturationPropsFromDeck_impl.hpp b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp index a4df7f425..f7a80b453 100644 --- a/opm/core/fluid/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp @@ -105,7 +105,7 @@ namespace Opm template int SaturationPropsFromDeck::numPhases() const { - return phase_usage_.num_phases; + return phase_usage_.num_phases; } @@ -191,18 +191,18 @@ namespace Opm template void SaturationPropsFromDeck::satRange(const int n, const int* cells, - double* smin, - double* smax) const + double* smin, + double* smax) const { ASSERT (cells != 0); - const int np = phase_usage_.num_phases; - for (int i = 0; i < n; ++i) { - for (int p = 0; p < np; ++p) { - smin[np*i + p] = funcForCell(cells[i]).smin_[p]; - smax[np*i + p] = funcForCell(cells[i]).smax_[p]; - } - } + const int np = phase_usage_.num_phases; + for (int i = 0; i < n; ++i) { + for (int p = 0; p < np; ++p) { + smin[np*i + p] = funcForCell(cells[i]).smin_[p]; + smax[np*i + p] = funcForCell(cells[i]).smax_[p]; + } + } } diff --git a/opm/core/fluid/SaturationPropsInterface.hpp b/opm/core/fluid/SaturationPropsInterface.hpp index e55995677..30f277cc1 100644 --- a/opm/core/fluid/SaturationPropsInterface.hpp +++ b/opm/core/fluid/SaturationPropsInterface.hpp @@ -65,11 +65,11 @@ namespace Opm double* pc, double* dpcds) const = 0; - /// Obtain the range of allowable saturation values. + /// Obtain the range of allowable saturation values. /// \param[in] n Number of data points. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling. - virtual void satRange(const int n, + virtual void satRange(const int n, const int* cells, double* smin, double* smax) const = 0; diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp index d3a42e321..8633fd9ee 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp @@ -70,7 +70,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of size numPhases(). - const double* surfaceDensities() const; + const double* surfaceDensities() const; /// Viscosity as a function of p and z. void mu(const int n, @@ -111,11 +111,11 @@ namespace Opm PhaseUsage phase_usage_; - int region_number_; + int region_number_; std::vector > props_; - double densities_[MaxNumPhases]; + double densities_[MaxNumPhases]; mutable std::vector data1_; mutable std::vector data2_; }; From 1d98e043a5aaa90ada79a6867be4478ef8a3dffc Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 5 Sep 2012 13:34:25 +0200 Subject: [PATCH 13/59] Fixed source term (measured at reservoir conditions). --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index d43f444a2..7b4a72e68 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -152,8 +152,8 @@ namespace Opm B_cell = 1.0/tm.A_[np*np*cell + 0]; double src_flux = -tm.source_[cell]; bool src_is_inflow = src_flux < 0.0; - influx = src_is_inflow ? B_cell*src_flux : 0.0; - outflux = !src_is_inflow ? B_cell*src_flux : 0.0; + influx = src_is_inflow ? src_flux : 0.0; + outflux = !src_is_inflow ? src_flux : 0.0; comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell]; dtpv = tm.dt_/tm.porevolume0_[cell]; for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { From 680276debfafb346fe5620deab0c0cf71dbed41c Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 5 Sep 2012 14:07:51 +0200 Subject: [PATCH 14/59] Fixed documentation. --- .../transport/reorder/TransportModelCompressibleTwophase.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index c0be99a9c..147d7da10 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -74,7 +74,6 @@ namespace Opm /// vertical stack, that do not interact with other columns (for /// gravity segregation. /// \param[in] columns Vector of cell-columns. - /// \param[in] porevolume0 Array of pore volumes at start of timestep. /// \param[in] dt Time step. /// \param[in, out] saturation Phase saturations. /// \param[in, out] surfacevol Surface volume densities for each phase. From bdcf0291e07be19c0abe1dd908b741ad2667b929 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Sep 2012 09:52:13 +0200 Subject: [PATCH 15/59] Fix error message. --- opm/core/wells/WellsManager.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 0215b053e..c31345fa4 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -544,7 +544,7 @@ namespace Opm cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; } else if (wci_line.injector_type_[0] == 'G') { if (!pu.phase_used[BlackoilPhases::Vapour]) { - THROW("Water phase not used, yet found water-injecting well."); + THROW("Gas phase not used, yet found gas-injecting well."); } cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; } From fa6b7729725a4d89d79699c7d0e32224b67367aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Sep 2012 09:53:11 +0200 Subject: [PATCH 16/59] Changed well initialization and property calculation. Bhp is now initialized to bhp target for bhp-controlled wells. Mobilities and pvt properties are now calculated from well perforation pressure and injection specifications for injectors, producers still use cell properties as before. --- opm/core/pressure/CompressibleTpfa.cpp | 39 +++++++++++++++++++++----- opm/core/simulator/WellState.hpp | 13 +++++++-- 2 files changed, 42 insertions(+), 10 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 954e07fa4..859585991 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -36,6 +36,7 @@ #include #include #include +#include namespace Opm { @@ -446,7 +447,7 @@ namespace Opm /// Compute per-iteration dynamic properties for wells. void CompressibleTpfa::computeWellDynamicData(const double /*dt*/, const BlackoilState& /*state*/, - const WellState& /*well_state*/) + const WellState& well_state) { // These are the variables that get computed by this function: // @@ -458,18 +459,42 @@ namespace Opm wellperf_A_.resize(nperf*np*np); wellperf_phasemob_.resize(nperf*np); // The A matrix is set equal to the perforation grid cells' - // matrix, for both injectors and producers. + // matrix for producers, computed from bhp and injection + // component fractions from // The mobilities are set equal to the perforation grid cells' - // mobilities, for both injectors and producers. + // mobilities for producers. + std::vector mu(np); for (int w = 0; w < nw; ++w) { + bool producer = (wells_->type[w] == PRODUCER); + const double* comp_frac = &wells_->comp_frac[np*w]; for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) { const int c = wells_->well_cells[j]; - const double* cA = &cell_A_[np*np*c]; double* wpA = &wellperf_A_[np*np*j]; - std::copy(cA, cA + np*np, wpA); - const double* cM = &cell_phasemob_[np*c]; double* wpM = &wellperf_phasemob_[np*j]; - std::copy(cM, cM + np, wpM); + if (producer) { + const double* cA = &cell_A_[np*np*c]; + std::copy(cA, cA + np*np, wpA); + const double* cM = &cell_phasemob_[np*c]; + std::copy(cM, cM + np, wpM); + } else { + const double bhp = well_state.bhp()[w]; + double perf_p = bhp; + for (int phase = 0; phase < np; ++phase) { + perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase]; + } + // Hack warning: comp_frac is used as a component + // surface-volume variable in calls to matrix() and + // viscosity(), but as a saturation in the call to + // relperm(). This is probably ok as long as injectors + // only inject pure fluids. + props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL); + props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL); + ASSERT(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6); + props_.relperm (1, comp_frac, &c, wpM , NULL); + for (int phase = 0; phase < np; ++phase) { + wpM[phase] /= mu[phase]; + } + } } } } diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 59d9f0602..7ddc886e2 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -37,10 +37,17 @@ namespace Opm if (wells) { const int nw = wells->number_of_wells; bhp_.resize(nw); - // Initialize bhp to be pressure in first perforation cell. + // Initialize bhp to be target pressure + // if bhp-controlled well, otherwise set + // to pressure in first perforation cell. for (int w = 0; w < nw; ++w) { - const int cell = wells->well_cells[wells->well_connpos[w]]; - bhp_[w] = state.pressure()[cell]; + const WellControls* ctrl = wells->ctrls[w]; + if (ctrl->type[ctrl->current] == BHP) { + bhp_[w] = ctrl->target[ctrl->current]; + } else { + const int cell = wells->well_cells[wells->well_connpos[w]]; + bhp_[w] = state.pressure()[cell]; + } } perfrates_.resize(wells->well_connpos[nw]); } From 67b5f007fde0ae8b5351e5e4cb4e9dacf4cda32d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Sep 2012 10:40:36 +0200 Subject: [PATCH 17/59] Made initialization from SWAT/SGAS etc. more robust and general. --- opm/core/utility/initState_impl.hpp | 66 ++++++++++++++++++++++------- 1 file changed, 51 insertions(+), 15 deletions(-) diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index 7d01ade0e..c139bf0e0 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include namespace Opm @@ -512,15 +513,23 @@ namespace Opm State& state) { const int num_phases = props.numPhases(); + const PhaseUsage pu = phaseUsageFromDeck(deck); + if (num_phases != pu.num_phases) { + THROW("initStateFromDeck(): Uuser specified property object with " << num_phases << " phases, " + "found " << pu.num_phases << " phases in deck."); + } state.init(grid, num_phases); if (deck.hasField("EQUIL")) { if (num_phases != 2) { THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios."); } + if (pu.phase_used[BlackoilPhases::Vapour]) { + THROW("initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas)."); + } // Set saturations depending on oil-water contact. const EQUIL& equil= deck.getEQUIL(); if (equil.equil.size() != 1) { - THROW("No region support yet."); + THROW("initStateFromDeck(): No region support yet."); } const double woc = equil.equil[0].water_oil_contact_depth_; initWaterOilContact(grid, props, woc, WaterBelow, state); @@ -528,37 +537,64 @@ namespace Opm const double datum_z = equil.equil[0].datum_depth_; const double datum_p = equil.equil[0].datum_depth_pressure_; initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state); - } else if (deck.hasField("SWAT") && deck.hasField("PRESSURE")) { - // Set saturations from SWAT, pressure from PRESSURE. + } else if (deck.hasField("PRESSURE")) { + // Set saturations from SWAT/SGAS, pressure from PRESSURE. std::vector& s = state.saturation(); std::vector& p = state.pressure(); - const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& p_deck = deck.getFloatingPointValue("PRESSURE"); const int num_cells = grid.number_of_cells; if (num_phases == 2) { - for (int c = 0; c < num_cells; ++c) { - int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; - s[2*c] = sw_deck[c_deck]; - s[2*c + 1] = 1.0 - s[2*c]; - p[c] = p_deck[c_deck]; + if (!pu.phase_used[BlackoilPhases::Aqua]) { + // oil-gas: we require SGAS + if (!deck.hasField("SGAS")) { + THROW("initStateFromDeck(): missing SGAS keyword in 2-phase init"); + } + const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c + gpos] = sg_deck[c_deck]; + s[2*c + opos] = 1.0 - sg_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } else { + // water-oil or water-gas: we require SWAT + if (!deck.hasField("SWAT")) { + THROW("initStateFromDeck(): missing SWAT keyword in 2-phase init"); + } + const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int nwpos = (wpos + 1) % 2; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c + wpos] = sw_deck[c_deck]; + s[2*c + nwpos] = 1.0 - sw_deck[c_deck]; + p[c] = p_deck[c_deck]; + } } } else if (num_phases == 3) { - if (!deck.hasField("SGAS")) { - THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found)."); + const bool has_swat_sgas = deck.hasField("SWAT") && deck.hasField("SGAS"); + if (!has_swat_sgas) { + THROW("initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init."); } + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); for (int c = 0; c < num_cells; ++c) { int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; - s[3*c] = sw_deck[c_deck]; - s[3*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); - s[3*c + 2] = sg_deck[c_deck]; + s[3*c + wpos] = sw_deck[c_deck]; + s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); + s[3*c + gpos] = sg_deck[c_deck]; p[c] = p_deck[c_deck]; } } else { THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases."); } } else { - THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE."); + THROW("initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS."); } // Finally, init face pressures. From 65447604ae680821eb960dbeb068d5dc972406de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Sep 2012 20:56:08 +0200 Subject: [PATCH 18/59] Typo fix. --- opm/core/utility/initState_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index c139bf0e0..575e18a41 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -515,7 +515,7 @@ namespace Opm const int num_phases = props.numPhases(); const PhaseUsage pu = phaseUsageFromDeck(deck); if (num_phases != pu.num_phases) { - THROW("initStateFromDeck(): Uuser specified property object with " << num_phases << " phases, " + THROW("initStateFromDeck(): user specified property object with " << num_phases << " phases, " "found " << pu.num_phases << " phases in deck."); } state.init(grid, num_phases); From 2017481a58d650f1b2b182a543848879eda60f32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 20 Sep 2012 14:33:57 +0200 Subject: [PATCH 19/59] Improve diagnostic output if crossflow is detected. --- opm/core/utility/miscUtilities.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index aeafbfca7..53e997731 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -374,7 +374,10 @@ namespace Opm if (perf_rate > 0.0) { // perf_rate is a total inflow rate, we want a water rate. if (wells->type[w] != INJECTOR) { - std::cout << "**** Warning: crossflow in well with index " << w << " ignored." << std::endl; + std::cout << "**** Warning: crossflow in well " + << w << " perf " << perf - wells->well_connpos[w] + << " ignored. Rate was " + << perf_rate/Opm::unit::day << " m^3/day." << std::endl; perf_rate = 0.0; } else { ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6); From 50a23c0f5d7e0d2a254585b84b55f72ad2ac400f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 20 Sep 2012 14:34:33 +0200 Subject: [PATCH 20/59] Support shut wells in compressible tpfa solver. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 60 ++++++++++++++++------ 1 file changed, 44 insertions(+), 16 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index eb8dc8f34..8eecbfed6 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -741,6 +741,29 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt, } +/* ---------------------------------------------------------------------- */ +static void +welleq_coeff_shut(int np, struct cfs_tpfa_res_data *h, + double *res, double *w2c, double *w2w) +/* ---------------------------------------------------------------------- */ +{ + int p; + double fwi; + const double *dpflux_w; + + /* Sum reservoir phase flux derivatives set by + * compute_darcyflux_and_deriv(). */ + dpflux_w = h->pimpl->flux_work + (1 * np); + for (p = 0, fwi = 0.0; p < np; p++) { + fwi += dpflux_w[ p ]; + } + + *res = 0.0; + *w2c = 0.0; + *w2w = fwi; +} + + /* ---------------------------------------------------------------------- */ static void welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h, @@ -815,23 +838,25 @@ assemble_completion_to_well(int w, int c, int nc, int np, ctrl = W->ctrls[ w ]; if (ctrl->current < 0) { - assert (0); /* Shut wells currently not supported */ + /* Interpreting a negative current control index to mean a shut well */ + welleq_coeff_shut(np, h, &res, &w2c, &w2w); } + else { + switch (ctrl->type[ ctrl->current ]) { + case BHP : + welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ], + h, &res, &w2c, &w2w); + break; - switch (ctrl->type[ ctrl->current ]) { - case BHP : - welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ], - h, &res, &w2c, &w2w); - break; + case RESERVOIR_RATE: + assert (W->number_of_phases == np); + welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w); + break; - case RESERVOIR_RATE: - assert (W->number_of_phases == np); - welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w); - break; - - case SURFACE_RATE: - assert (0); /* Surface rate currently not supported */ - break; + case SURFACE_RATE: + assert (0); /* Surface rate currently not supported */ + break; + } } /* Assemble completion contributions */ @@ -854,7 +879,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells , struct cfs_tpfa_res_data *h ) { int w, i, c, np, np2, nc; - int is_neumann; + int is_neumann, is_open; double pw, dp; double *WI, *gpot, *pmobp; const double *Ac, *dAc; @@ -876,6 +901,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells , for (w = i = 0; w < W->number_of_wells; w++) { pw = wpress[ w ]; + is_open = W->ctrls[w]->current >= 0; for (; i < W->well_connpos[w + 1]; i++, gpot += np, pmobp += np) { @@ -888,7 +914,9 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells , init_completion_contrib(i, np, Ac, dAc, h->pimpl); - assemble_completion_to_cell(c, nc + w, np, dt, h); + if (is_open) { + assemble_completion_to_cell(c, nc + w, np, dt, h); + } /* Prepare for RESV controls */ compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot, From e3388575d620f05946405e4b5a4152cff7d9c7bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 20 Sep 2012 14:35:03 +0200 Subject: [PATCH 21/59] Fix treatment of WELOPEN keyword. Now you can actually shut and open wells with WELOPEN. The following caveats apply: - this may interact improperly with group controls, - dynamic usage of WCONINJE/WCONPROD should not be mixed with WELOPEN. --- opm/core/wells/WellsManager.cpp | 51 ++++++++++++++++----------------- 1 file changed, 25 insertions(+), 26 deletions(-) diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index c31345fa4..95cb0369a 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -716,32 +716,31 @@ namespace Opm #endif if (deck.hasField("WELOPEN")) { - const WELOPEN& welopen = deck.getWELOPEN(); - - for (size_t i = 0; i < welopen.welopen.size(); ++i) { - WelopenLine line = welopen.welopen[i]; - std::string wellname = line.well_; - std::map::const_iterator it = well_names_to_index.find(wellname); - if (it == well_names_to_index.end()) { - THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS."); - } - int index = it->second; - if (line.openshutflag_ == "SHUT") { - // We currently don't care if the well is open or not. - /// \TODO Should this perhaps be allowed? I.e. should it be if(well_shut) { shutwell(); } else { /* do nothing*/ }? - //ASSERT(w_->ctrls[index]->current < 0); - } else if (line.openshutflag_ == "OPEN") { - //ASSERT(w_->ctrls[index]->current >= 0); - } else { - THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT."); - } - - // We revert back to it's original control. - // Note that this is OK as ~~ = id. - w_->ctrls[index]->current = ~w_->ctrls[index]->current; - - - } + const WELOPEN& welopen = deck.getWELOPEN(); + for (size_t i = 0; i < welopen.welopen.size(); ++i) { + WelopenLine line = welopen.welopen[i]; + std::string wellname = line.well_; + std::map::const_iterator it = well_names_to_index.find(wellname); + if (it == well_names_to_index.end()) { + THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS."); + } + const int index = it->second; + if (line.openshutflag_ == "SHUT") { + int& cur_ctrl = w_->ctrls[index]->current; + if (cur_ctrl >= 0) { + cur_ctrl = ~cur_ctrl; + } + ASSERT(w_->ctrls[index]->current < 0); + } else if (line.openshutflag_ == "OPEN") { + int& cur_ctrl = w_->ctrls[index]->current; + if (cur_ctrl < 0) { + cur_ctrl = ~cur_ctrl; + } + ASSERT(w_->ctrls[index]->current >= 0); + } else { + THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT."); + } + } } // Build the well_collection_ well group hierarchy. From bdcd5236bd9d592410d3caa480236de168226bbf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 20 Sep 2012 15:03:38 +0200 Subject: [PATCH 22/59] Don't crash on models without wells. The user will legitimately want to run models that do not specify wells (e.g., using boundary conditions). While we do not yet fully support that configuration (no wells), we absolutely must not crash by dereferencing null pointers or generating pointers into ::empty() std::vector<>s. This commit installs the required guards needed to avoid said failure mode. --- opm/core/pressure/CompressibleTpfa.cpp | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 859585991..42b3fba75 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -123,7 +123,7 @@ namespace Opm WellState& well_state) { const int nc = grid_.number_of_cells; - const int nw = wells_->number_of_wells; + const int nw = (wells_ != 0) ? wells_->number_of_wells : 0; // Set up dynamic data. computePerSolveDynamicData(dt, state, well_state); @@ -454,8 +454,8 @@ namespace Opm // std::vector wellperf_A_; // std::vector wellperf_phasemob_; const int np = props_.numPhases(); - const int nw = wells_->number_of_wells; - const int nperf = wells_->well_connpos[nw]; + const int nw = (wells_ != 0) ? wells_->number_of_wells : 0; + const int nperf = (wells_ != 0) ? wells_->well_connpos[nw] : 0; wellperf_A_.resize(nperf*np*np); wellperf_phasemob_.resize(nperf*np); // The A matrix is set equal to the perforation grid cells' @@ -512,9 +512,9 @@ namespace Opm const double* z = &state.surfacevol()[0]; UnstructuredGrid* gg = const_cast(&grid_); CompletionData completion_data; - completion_data.gpot = &wellperf_gpot_[0]; - completion_data.A = &wellperf_A_[0]; - completion_data.phasemob = &wellperf_phasemob_[0]; + completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0; + completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0; + completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0; cfs_tpfa_res_wells wells_tmp; wells_tmp.W = const_cast(wells_); wells_tmp.data = &completion_data; @@ -599,9 +599,9 @@ namespace Opm { UnstructuredGrid* gg = const_cast(&grid_); CompletionData completion_data; - completion_data.gpot = const_cast(&wellperf_gpot_[0]); - completion_data.A = const_cast(&wellperf_A_[0]); - completion_data.phasemob = const_cast(&wellperf_phasemob_[0]); + completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast(&wellperf_gpot_[0]) : 0; + completion_data.A = ! wellperf_A_.empty() ? const_cast(&wellperf_A_[0]) : 0; + completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast(&wellperf_phasemob_[0]) : 0; cfs_tpfa_res_wells wells_tmp; wells_tmp.W = const_cast(wells_); wells_tmp.data = &completion_data; @@ -609,6 +609,9 @@ namespace Opm forces.wells = &wells_tmp; forces.src = NULL; + double* wpress = ! well_state.bhp ().empty() ? & well_state.bhp ()[0] : 0; + double* wflux = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0; + cfs_tpfa_res_flux(gg, &forces, props_.numPhases(), @@ -617,9 +620,9 @@ namespace Opm &face_phasemob_[0], &face_gravcap_[0], &state.pressure()[0], - &well_state.bhp()[0], + wpress, &state.faceflux()[0], - &well_state.perfRates()[0]); + wflux); cfs_tpfa_res_fpress(gg, props_.numPhases(), &htrans_[0], From b137eb0b65dfaecf01db6bef9f2ba5cbc84b28b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 20 Sep 2012 15:48:48 +0200 Subject: [PATCH 23/59] Interpret ``wells != 0 && wells->W == 0'' as ``no wells''. The CompressibleTpfa class always passes a non-null `forces->wells' object to the constructor, assembly, and reconstruction routines but uses ``forces->wells->W == 0'' to signify a simulation model without wells. This is, arguably, an error in the CompressibleTpfa class but one that does not require a lot of work to support in the cfs_tpfa_residual module. Insert the extra tests in an effort to honour the ``liberal in what you accept, strict in what you produce'' principle. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 24 +++++++++++++--------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index eb8dc8f34..023e3c2d1 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -142,7 +142,7 @@ impl_allocate(struct UnstructuredGrid *G , nnu = G->number_of_cells; nwperf = 0; - if (wells != NULL) { assert (wells->W != NULL); + if ((wells != NULL) && (wells->W != NULL)) { nnu += wells->W->number_of_wells; nwperf = wells->W->well_connpos[ wells->W->number_of_wells ]; } @@ -185,13 +185,15 @@ construct_matrix(struct UnstructuredGrid *G , /* ---------------------------------------------------------------------- */ { int f, c1, c2, w, i, nc, nnu; + int wells_present; size_t nnz; struct CSRMatrix *A; nc = nnu = G->number_of_cells; - if (wells != NULL) { - assert (wells->W != NULL); + + wells_present = (wells != NULL) && (wells->W != NULL); + if (wells_present) { nnu += wells->W->number_of_wells; } @@ -214,7 +216,7 @@ construct_matrix(struct UnstructuredGrid *G , } } - if (wells != NULL) { + if (wells_present) { /* Well <-> cell connections */ struct Wells *W = wells->W; @@ -252,7 +254,7 @@ construct_matrix(struct UnstructuredGrid *G , } } - if (wells != NULL) { + if (wells_present) { /* Fill well <-> cell connections */ struct Wells *W = wells->W; @@ -1127,8 +1129,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G , nf = G->number_of_faces; nwperf = 0; - if (wells != NULL) { - assert (wells->W != NULL); + if ((wells != NULL) && (wells->W != NULL)) { nwperf = wells->W->well_connpos[ wells->W->number_of_wells ]; } @@ -1194,7 +1195,9 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , assemble_cell_contrib(G, c, h); } - if ((forces != NULL) && (forces->wells != NULL)) { + if ((forces != NULL) && + (forces->wells != NULL) && + (forces->wells->W != NULL)) { compute_well_compflux_and_deriv(forces->wells, cq->nphases, cpress, wpress, h->pimpl); @@ -1297,8 +1300,9 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G , { compute_flux(G, np, trans, pmobf, gravcap_f, cpress, fflux); - if ((forces != NULL) && (forces->wells != NULL) && - (wpress != NULL) && (wflux != NULL)) { + if ((forces != NULL) && + (forces->wells != NULL) && + (forces->wells->W != NULL) && (wpress != NULL) && (wflux != NULL)) { compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux); } From ab21d44c9ab40a25db5041b892cf0c61de7197e0 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Fri, 21 Sep 2012 15:06:47 +0200 Subject: [PATCH 24/59] Disable warning for using DUNE's FieldVector::size In DUNE 2.2 FieldVector::size changed from being a member to being a method. A compatibility warning is issued if you include the relevant headers. This warning can be silenced for DUNE modules by using passing the option --enable-fieldvector-size-is-method to ./configure. This patch effectively does the same, but through a macro definition. --- opm/core/linalg/LinearSolverIstl.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index 3c4435e8e..a3f9583ee 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -26,6 +26,10 @@ #include +// Silence compatibility warning from DUNE headers since we don't use +// the deprecated member anyway (in this compilation unit) +#define DUNE_COMMON_FIELDVECTOR_SIZE_IS_METHOD 1 + // TODO: clean up includes. #include #include From 56e81968e3965ee8a7828b0faea799f664df22a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 24 Sep 2012 16:43:00 +0200 Subject: [PATCH 25/59] Add support for new three-phase relperm option to BlackoilPropertiesFromDeck. New parameter option added: 'threephase_model' can now be 'gwseg'. --- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 23 +++++++++++++------ opm/core/fluid/SaturationPropsFromDeck.hpp | 1 + 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index e2457f39b..78e3c5db2 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -28,8 +28,8 @@ namespace Opm { rock_.init(deck, grid); pvt_.init(deck, 200); - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, 200); @@ -50,30 +50,39 @@ namespace Opm // Unfortunate lack of pointer smartness here... const int sat_samples = param.getDefault("sat_tab_size", 200); std::string threephase_model = param.getDefault("threephase_model", "simple"); - bool use_stone2 = (threephase_model == "stone2"); if (sat_samples > 1) { - if (use_stone2) { + if (threephase_model == "stone2") { SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); - } else { + } else if (threephase_model == "simple") { SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); } } else { - if (use_stone2) { + if (threephase_model == "stone2") { SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); - } else { + } else if (threephase_model == "simple") { SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); } } diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index 0a0599ec8..ea2acb342 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -26,6 +26,7 @@ #include #include #include +#include #include struct UnstructuredGrid; From 9f69e9fa51ba7a2c8c411d700a0ca1d48a931eb6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 24 Sep 2012 17:09:50 +0200 Subject: [PATCH 26/59] Guard against input error. If no valid threephase_model is input, throw instead of crashing. --- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 78e3c5db2..271f93951 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -66,6 +66,8 @@ namespace Opm = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); + } else { + THROW("Unknown threephase_model: " << threephase_model); } } else { if (threephase_model == "stone2") { @@ -83,6 +85,8 @@ namespace Opm = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); + } else { + THROW("Unknown threephase_model: " << threephase_model); } } From b9a2c141136b1c87e5acf84a7fcc59258f070de7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Sep 2012 14:00:17 +0200 Subject: [PATCH 27/59] Add proper support for source terms. This fixes the problem with infinite tofs at sinks. --- .../reorder/TransportModelTracerTof.cpp | 25 ++++++++++++++++--- .../reorder/TransportModelTracerTof.hpp | 9 ++++--- 2 files changed, 27 insertions(+), 7 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 0825706e3..70c2fc400 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -21,6 +21,8 @@ #include #include #include +#include +#include namespace Opm { @@ -37,15 +39,25 @@ namespace Opm /// Solve for time-of-flight at next timestep. - /// \param[in] darcyflux Array of signed face fluxes. - /// \param[in] porevolume Array of pore volumes. - /// \param[in, out] tof Array of time-of-flight values. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[out] tof Array of time-of-flight values. void TransportModelTracerTof::solveTof(const double* darcyflux, const double* porevolume, + const double* source, std::vector& tof) { darcyflux_ = darcyflux; porevolume_ = porevolume; + source_ = source; +#ifndef NDEBUG + // Sanity check for sources. + const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0); + if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) { + THROW("Sources do not sum to zero: " << cum_src); + } +#endif tof.resize(grid_.number_of_cells); std::fill(tof.begin(), tof.end(), 0.0); tof_ = &tof[0]; @@ -58,8 +70,13 @@ namespace Opm void TransportModelTracerTof::solveSingleCell(const int cell) { // Compute flux terms. + // Sources have zero tof, and therefore do not contribute + // to upwind_term. Sinks on the other hand, must be added + // to the downwind_flux (note sign change resulting from + // different sign conventions: pos. source is injection, + // pos. flux is outflow). double upwind_term = 0.0; - double downwind_flux = 0.0; + double downwind_flux = std::max(-source_[cell], 0.0); for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { int f = grid_.cell_faces[i]; double flux; diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp index 1ab7ef3b4..21375d895 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -40,11 +40,13 @@ namespace Opm TransportModelTracerTof(const UnstructuredGrid& grid); /// Solve for time-of-flight at next timestep. - /// \param[in] darcyflux Array of signed face fluxes. - /// \param[in] porevolume Array of pore volumes. - /// \param[in, out] tof Array of time-of-flight values. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[out] tof Array of time-of-flight values. void solveTof(const double* darcyflux, const double* porevolume, + const double* source, std::vector& tof); private: @@ -55,6 +57,7 @@ namespace Opm const UnstructuredGrid& grid_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell + const double* source_; // one volumetric source term per cell double* tof_; }; From 1a227dbf86eb6395233ce41848ab55fb49e0c32d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 08:58:03 +0200 Subject: [PATCH 28/59] Added skeleton of general order DG tof solver. --- .../TransportModelTracerTofDiscGal.cpp | 188 ++++++++++++++++++ .../TransportModelTracerTofDiscGal.hpp | 74 +++++++ 2 files changed, 262 insertions(+) create mode 100644 opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp create mode 100644 opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp new file mode 100644 index 000000000..05ba59b08 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -0,0 +1,188 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + + /// A class providing discontinuous Galerkin basis functions. + struct DGBasis + { + static int numBasisFunc(const int dimensions, + const int degree) + { + switch (dimensions) { + case 1: + return degree + 1; + case 2: + return (degree + 2)*(degree + 1)/2; + case 3: + return (degree + 3)*(degree + 2)*(degree + 1)/6; + default: + THROW("Dimensions must be 1, 2 or 3."); + } + } + + /// Evaluate all nonzero basis functions at x, + /// writing to f_x. The array f_x must have + /// size numBasisFunc(grid.dimensions, degree). + /// + /// The basis functions are the following + /// Degree 0: 1. + /// Degree 1: x - xc, y - yc, z - zc etc. + /// Further degrees await development. + static void eval(const UnstructuredGrid& grid, + const int cell, + const int degree, + const double* x, + double* f_x) + { + const int dim = grid.dimensions; + const double* cc = grid.cell_centroids + dim*cell; + // Note intentional fallthrough in this switch statement! + switch (degree) { + case 1: + for (int ix = 0; ix < dim; ++ix) { + f_x[1 + ix] = x[ix] - cc[ix]; + } + case 0: + f_x[0] = 1; + break; + default: + THROW("Maximum degree is 1 for now."); + } + } + + /// Evaluate gradients of all nonzero basis functions at x, + /// writing to grad_f_x. The array grad_f_x must have size + /// numBasisFunc(grid.dimensions, degree) * grid.dimensions. + /// The components of the first basis function + /// gradient come before the components of the second etc. + static void evalGrad(const UnstructuredGrid& grid, + const int /*cell*/, + const int degree, + const double* /*x*/, + double* grad_f_x) + { + const int dim = grid.dimensions; + const int num_basis = numBasisFunc(dim, degree); + std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0); + if (degree > 1) { + THROW("Maximum degree is 1 for now."); + } else if (degree == 1) { + for (int ix = 0; ix < dim; ++ix) { + grad_f_x[dim*(ix + 1) + ix] = 1.0; + } + } + } + }; + + + + + /// A class providing numerical quadrature functions. + struct Quadrature + { + static void x(){} + }; + + + + + + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid) + : grid_(grid) + { + } + + + + + /// Solve for time-of-flight at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[in] degree Polynomial degree of DG basis functions used. + /// \param[out] tof_coeff Array of time-of-flight solution coefficients. + /// The values are ordered by cell, meaning that + /// the K coefficients corresponding to the first + /// cell comes before the K coefficients corresponding + /// to the second cell etc. + /// K depends on degree and grid dimension. + void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux, + const double* porevolume, + const double* source, + const int degree, + std::vector& tof_coeff) + { + darcyflux_ = darcyflux; + porevolume_ = porevolume; + source_ = source; +#ifndef NDEBUG + // Sanity check for sources. + const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0); + if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) { + THROW("Sources do not sum to zero: " << cum_src); + } +#endif + degree_ = degree; + tof_coeff.resize(grid_.number_of_cells); + std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0); + tof_coeff_ = &tof_coeff[0]; + reorderAndTransport(grid_, darcyflux); + } + + + + + void TransportModelTracerTofDiscGal::solveSingleCell(const int cell) + { + // Residual: + // For each cell K, basis function b_j (spanning V_h), + // writing the solution u_h|K = \sum_i c_i b_i + // Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx + // + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds + // - \int_K \phi b_j + THROW("Not implemented yet!"); + } + + + + + void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells) + { + std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; + for (int i = 0; i < num_cells; ++i) { + solveSingleCell(cells[i]); + } + } + + + + +} // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp new file mode 100644 index 000000000..94d5edf17 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -0,0 +1,74 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED +#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED + +#include +#include +#include +#include +struct UnstructuredGrid; + +namespace Opm +{ + + class IncompPropertiesInterface; + + /// Implements a reordering transport solver for incompressible two-phase flow. + class TransportModelTracerTofDiscGal : public TransportModelInterface + { + public: + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTofDiscGal(const UnstructuredGrid& grid); + + /// Solve for time-of-flight at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[in] degree Polynomial degree of DG basis functions used. + /// \param[out] tof_coeff Array of time-of-flight solution coefficients. + /// The values are ordered by cell, meaning that + /// the K coefficients corresponding to the first + /// cell comes before the K coefficients corresponding + /// to the second cell etc. + /// K depends on degree and grid dimension. + void solveTof(const double* darcyflux, + const double* porevolume, + const double* source, + const int degree, + std::vector& tof_coeff); + + private: + virtual void solveSingleCell(const int cell); + virtual void solveMultiCell(const int num_cells, const int* cells); + + private: + const UnstructuredGrid& grid_; + const double* darcyflux_; // one flux per grid face + const double* porevolume_; // one volume per cell + const double* source_; // one volumetric source term per cell + int degree_; + double* tof_coeff_; + }; + +} // namespace Opm + +#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED From fb8e79857cf572eed6f51a498e5f857addba7fac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 26 Sep 2012 10:14:45 +0200 Subject: [PATCH 29/59] Add tentative implementation of surface-rate targets. Not tested at present. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 33 ++++++++++++++++++++-- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index 0e9b6f7d2..bcbacec28 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -820,7 +820,34 @@ welleq_coeff_resv(int np, struct cfs_tpfa_res_data *h, /* ---------------------------------------------------------------------- */ static void -assemble_completion_to_well(int w, int c, int nc, int np, +welleq_coeff_surfrate(int i, int np, struct cfs_tpfa_res_data *h, + struct WellControls *ctrl, + double *res, double *w2c, double *w2w) +/* ---------------------------------------------------------------------- */ +{ + int p; + double f; + const double *pflux, *dpflux_w, *dpflux_c, *distr; + + pflux = h->pimpl->compflux_p + (i * (1 * np)); + dpflux_w = h->pimpl->compflux_deriv_p + (i * (2 * np)); + dpflux_c = dpflux_w + (1 * (1 * np)); + distr = ctrl->distr + (ctrl->current * (1 * np)); + + *res = *w2c = *w2w = 0.0; + for (p = 0; p < np; p++) { + f = distr[ p ]; + + *res += f * pflux [ p ]; + *w2w += f * dpflux_w[ p ]; + *w2c += f * dpflux_c[ p ]; + } +} + + +/* ---------------------------------------------------------------------- */ +static void +assemble_completion_to_well(int i, int w, int c, int nc, int np, double pw, double dt, struct cfs_tpfa_res_wells *wells, struct cfs_tpfa_res_data *h ) @@ -856,7 +883,7 @@ assemble_completion_to_well(int w, int c, int nc, int np, break; case SURFACE_RATE: - assert (0); /* Surface rate currently not supported */ + welleq_coeff_surfrate(i, np, h, ctrl, &res, &w2c, &w2w); break; } } @@ -925,7 +952,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells , h->pimpl->flux_work, h->pimpl->flux_work + np); - assemble_completion_to_well(w, c, nc, np, pw, dt, wells, h); + assemble_completion_to_well(i, w, c, nc, np, pw, dt, wells, h); } ctrl = W->ctrls[ w ]; From 99a12a7edc2c42d6a6d2579b56788177ff7ad66d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 13:30:54 +0200 Subject: [PATCH 30/59] Initial version of DG(1) for tof implemented. Basis functions, quadratures and velocity interpolation are basic versions, not handling any higher than DG(1) for now. These are currently in helper classes and functions. The code in the main solver class is written with the aim of supporting DG(n) generally. --- .../TransportModelTracerTofDiscGal.cpp | 228 +++++++++++++++++- .../TransportModelTracerTofDiscGal.hpp | 8 + 2 files changed, 230 insertions(+), 6 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 05ba59b08..e626a2b3c 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -103,20 +104,115 @@ namespace Opm - /// A class providing numerical quadrature functions. - struct Quadrature + /// A class providing numerical quadrature for cells. + class CellQuadrature { - static void x(){} + public: + CellQuadrature(const UnstructuredGrid& grid, + const int cell, + const int degree) + : grid_(grid), cell_(cell), degree_(degree) + { + if (degree > 1) { + THROW("Only quadrature degree up to 1 for now."); + } + } + + int numQuadPts() const + { + return 1; + } + + void quadPtCoord(int /*index*/, double* coord) const + { + const double* cc = grid_.cell_centroids + grid_.dimensions*cell_; + std::copy(cc, cc + grid_.dimensions, coord); + } + + double quadPtWeight(int /*index*/) const + { + return 1.0; + } + + private: + const UnstructuredGrid& grid_; + const int cell_; + const int degree_; }; + /// A class providing numerical quadrature for faces. + class FaceQuadrature + { + public: + FaceQuadrature(const UnstructuredGrid& grid, + const int face, + const int degree) + : grid_(grid), face_(face), degree_(degree) + { + if (degree > 1) { + THROW("Only quadrature degree up to 1 for now."); + } + } + + int numQuadPts() const + { + return 1; + } + + void quadPtCoord(int /*index*/, double* coord) const + { + const double* fc = grid_.face_centroids + grid_.dimensions*face_; + std::copy(fc, fc + grid_.dimensions, coord); + } + + double quadPtWeight(int /*index*/) const + { + return 1.0; + } + + private: + const UnstructuredGrid& grid_; + const int face_; + const int degree_; + }; + + + // Initial version: only a constant interpolation. + static void interpolateVelocity(const UnstructuredGrid& grid, + const int cell, + const double* darcyflux, + const double* /*x*/, + double* v) + { + const int dim = grid.dimensions; + std::fill(v, v + dim, 0.0); + const double* cc = grid.cell_centroids + cell*dim; + for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell+1]; ++hface) { + const int face = grid.cell_faces[hface]; + const double* fc = grid.face_centroids + face*dim; + double flux = 0.0; + if (cell == grid.face_cells[2*face]) { + flux = darcyflux[face]; + } else { + flux = -darcyflux[face]; + } + for (int dd = 0; dd < dim; ++dd) { + v[dd] += flux * (fc[dd] - cc[dd]) / grid.cell_volumes[cell]; + } + } + } + + /// Construct solver. /// \param[in] grid A 2d or 3d grid. TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid) - : grid_(grid) + : grid_(grid), + coord_(grid.dimensions), + velocity_(grid.dimensions) { } @@ -151,9 +247,15 @@ namespace Opm } #endif degree_ = degree; - tof_coeff.resize(grid_.number_of_cells); + const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_); + tof_coeff.resize(num_basis*grid_.number_of_cells); std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0); tof_coeff_ = &tof_coeff[0]; + rhs_.resize(num_basis); + jac_.resize(num_basis*num_basis); + basis_.resize(num_basis); + basis_nb_.resize(num_basis); + grad_basis_.resize(num_basis*grid_.dimensions); reorderAndTransport(grid_, darcyflux); } @@ -168,7 +270,121 @@ namespace Opm // Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx // + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds // - \int_K \phi b_j - THROW("Not implemented yet!"); + // This is linear in c_i, so we do not need any nonlinear iterations. + // We assemble the jacobian and the right-hand side. The residual is + // equal to Res = Jac*c - rhs, and we compute rhs directly. + + const int dim = grid_.dimensions; + const int num_basis = DGBasis::numBasisFunc(degree_, dim); + + std::fill(rhs_.begin(), rhs_.end(), 0.0); + std::fill(jac_.begin(), jac_.end(), 0.0); + + // Compute cell residual contribution. + // Note: Assumes that \int_K b_j = 0 for all j > 0 + rhs_[0] += porevolume_[cell]; + + // Compute upstream residual contribution. + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + int upstream_cell = -1; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + upstream_cell = grid_.face_cells[2*face+1]; + } else { + flux = -darcyflux_[face]; + upstream_cell = grid_.face_cells[2*face]; + } + if (upstream_cell < 0) { + // This is an outer boundary. Assumed tof = 0 on inflow, so no contribution. + continue; + } + if (flux >= 0.0) { + // This is an outflow boundary. + continue; + } + // Do quadrature over the face to compute + // \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds + // (where u_h^{ext} is the upstream unknown (tof)). + FaceQuadrature quad(grid_, face, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // u^ext flux B (B = {b_j}) + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]); + const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(), + tof_coeff_ + num_basis*upstream_cell, 0.0); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + rhs_[j] -= w * tof_upstream * flux * basis_[j]; + } + } + } + + // Compute cell jacobian contribution. We use Fortran ordering + // for jac_, i.e. rows cycling fastest. + { + CellQuadrature quad(grid_, cell, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // b_i (v \cdot \grad b_j) + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]); + interpolateVelocity(grid_, cell, darcyflux_, &coord_[0], &velocity_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + for (int dd = 0; dd < dim; ++dd) { + jac_[j*num_basis + i] += w * basis_[i] * grad_basis_[dim*j + dd] * velocity_[dd]; + } + } + } + } + } + + // Compute downstream jacobian contribution. + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + } else { + flux = -darcyflux_[face]; + } + if (flux <= 0.0) { + // This is an inflow boundary. + continue; + } + // Do quadrature over the face to compute + // \int_{\partial K} b_i (v(x) \cdot n) b_j ds + FaceQuadrature quad(grid_, face, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // u^ext flux B (B = {b_j}) + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + jac_[j*num_basis + i] += w * basis_[i] * flux * basis_[j]; + } + } + } + } + + // Solve linear equation. + MAT_SIZE_T n = num_basis; + MAT_SIZE_T nrhs = 1; + MAT_SIZE_T lda = num_basis; + std::vector piv(num_basis); + MAT_SIZE_T ldb = num_basis; + MAT_SIZE_T info = 0; + dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info); + if (info != 0) { + THROW("Lapack error: " << info); + } + // The solution ends up in rhs_, so we must copy it. + std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); } diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 94d5edf17..b9163f8ed 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -67,6 +67,14 @@ namespace Opm const double* source_; // one volumetric source term per cell int degree_; double* tof_coeff_; + std::vector rhs_; // single-cell right-hand-side + std::vector jac_; // single-cell jacobian + // Below: storage for quantities needed by solveSingleCell(). + std::vector coord_; + std::vector basis_; + std::vector basis_nb_; + std::vector grad_basis_; + std::vector velocity_; }; } // namespace Opm From 93094ebeec7dca44eff7bbf331cf31e749ad9c3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 14:21:00 +0200 Subject: [PATCH 31/59] Fix argument order in call. --- opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index e626a2b3c..c042614a2 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -275,7 +275,7 @@ namespace Opm // equal to Res = Jac*c - rhs, and we compute rhs directly. const int dim = grid_.dimensions; - const int num_basis = DGBasis::numBasisFunc(degree_, dim); + const int num_basis = DGBasis::numBasisFunc(dim, degree_); std::fill(rhs_.begin(), rhs_.end(), 0.0); std::fill(jac_.begin(), jac_.end(), 0.0); From 04a76988d9874132808242897528ddfe3edd5f36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 15:15:04 +0200 Subject: [PATCH 32/59] Add sink term contribution. --- .../TransportModelTracerTofDiscGal.cpp | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index c042614a2..22e60a36c 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -343,7 +343,7 @@ namespace Opm } } - // Compute downstream jacobian contribution. + // Compute downstream jacobian contribution from faces. for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { const int face = grid_.cell_faces[hface]; double flux = 0.0; @@ -372,6 +372,25 @@ namespace Opm } } + // Compute downstream jacobian contribution from sink terms. + if (source_[cell] < 0.0) { + // A sink. + const double flux = -source_[cell]; // Sign convention for flux: outflux > 0. + // Do quadrature over the cell to compute + // \int_{K} b_i flux b_j dx + CellQuadrature quad(grid_, cell, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + jac_[j*num_basis + i] += w * basis_[i] * flux * basis_[j]; + } + } + } + } + // Solve linear equation. MAT_SIZE_T n = num_basis; MAT_SIZE_T nrhs = 1; From d83ab5856da0cc99de260ed36d60ecb955082501 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Sep 2012 09:49:36 +0200 Subject: [PATCH 33/59] Fix: forgotten multiply by cell volume in a quadrature. --- opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 22e60a36c..cdf909b3f 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -325,6 +325,7 @@ namespace Opm // Compute cell jacobian contribution. We use Fortran ordering // for jac_, i.e. rows cycling fastest. { + const double cvol = grid_.cell_volumes[cell]; CellQuadrature quad(grid_, cell, degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // b_i (v \cdot \grad b_j) @@ -336,7 +337,7 @@ namespace Opm for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { for (int dd = 0; dd < dim; ++dd) { - jac_[j*num_basis + i] += w * basis_[i] * grad_basis_[dim*j + dd] * velocity_[dd]; + jac_[j*num_basis + i] += w * basis_[i] * grad_basis_[dim*j + dd] * velocity_[dd] * cvol; } } } From f6b2306dab3cc5f520ba672c270ca9a402bb5898 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 28 Sep 2012 14:44:04 +0200 Subject: [PATCH 34/59] Work in progress: degree 2 quadratures. Also, changed quadrature degrees used to get exact quadratures for all terms. --- .../TransportModelTracerTofDiscGal.cpp | 282 ++++++++++++++++-- 1 file changed, 257 insertions(+), 25 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index cdf909b3f..e9e812ef3 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -29,6 +29,10 @@ namespace Opm { + // --------------- Helpers for TransportModelTracerTofDiscGal --------------- + + + /// A class providing discontinuous Galerkin basis functions. struct DGBasis { @@ -103,8 +107,83 @@ namespace Opm + static void cross(const double* a, const double* b, double* res) + { + res[0] = a[1]*b[2] - a[2]*b[1]; + res[1] = a[2]*b[0] - a[0]*b[2]; + res[2] = a[0]*b[1] - a[1]*b[0]; + } + + + + + static double triangleArea3d(const double* p0, + const double* p1, + const double* p2) + { + double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; + double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; + double cr[3]; + cross(a, b, cr); + return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]); + } + + + + + /// Calculates the determinant of a 3 x 3 matrix, represented as + /// three three-dimensional arrays. + static double determinantOf(const double* a0, + const double* a1, + const double* a2) + { + return + a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - + a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + + a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); + } + + + + + /// Computes the volume of a tetrahedron consisting of 4 vertices + /// with 3-dimensional coordinates + static double tetVolume(const double* p0, + const double* p1, + const double* p2, + const double* p3) + { + double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; + double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; + double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] }; + return std::fabs(determinantOf(a, b, c) / 6.0); + } + + + /// A class providing numerical quadrature for cells. + /// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i). + /// Note that this class does multiply weights by cell volume, + /// so weights always sum to cell volume. + /// Degree 1 method: + /// Midpoint (centroid) method. + /// n = 1, w_0 = cell volume, x_0 = cell centroid + /// Degree 2 method: + /// Based on subdivision of each cell face into triangles + /// with the face centroid as a common vertex, and then + /// subdividing the cell into tetrahedra with the cell + /// centroid as a common vertex. Then apply the tetrahedron + /// rule with the following 4 nodes (uniform weights): + /// a = 0.138196601125010515179541316563436 + /// x_i has all barycentric coordinates = a, except for + /// the i'th coordinate which is = 1 - 3a. + /// This rule is from http://nines.cs.kuleuven.be/ecf, + /// it is the second degree 2 4-point rule for tets, + /// referenced to Stroud(1971). + /// The tetrahedra are numbered T_{i,j}, and are given by the + /// cell centroid C, the face centroid FC_i, and two nodes + /// of face i: FN_{i,j}, FN_{i,j+1}. class CellQuadrature { public: @@ -113,25 +192,93 @@ namespace Opm const int degree) : grid_(grid), cell_(cell), degree_(degree) { - if (degree > 1) { - THROW("Only quadrature degree up to 1 for now."); + if (degree > 2) { + THROW("CellQuadrature exact for polynomial degrees > 1 not implemented."); + } + if (degree == 2) { + // Prepare subdivision. } } int numQuadPts() const { - return 1; + if (degree_ < 2) { + return 1; + } + // Degree 2 case. + int sumnodes = 0; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + } + return 4*sumnodes; } - void quadPtCoord(int /*index*/, double* coord) const + void quadPtCoord(const int index, double* coord) const { - const double* cc = grid_.cell_centroids + grid_.dimensions*cell_; - std::copy(cc, cc + grid_.dimensions, coord); + const int dim = grid_.dimensions; + const double* cc = grid_.cell_centroids + dim*cell_; + if (degree_ < 2) { + std::copy(cc, cc + dim, coord); + return; + } + // Degree 2 case. + int tetindex = index / 4; + const int subindex = index % 4; + const double* nc = grid_.node_coordinates; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + if (nfn <= tetindex) { + // Our tet is not associated with this face. + tetindex -= nfn; + continue; + } + const double* fc = grid_.face_centroids + dim*face; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face]; + const int node0 = fnodes[tetindex]; + const int node1 = fnodes[(tetindex + 1) % nfn]; + const double* n0c = nc + dim*node0; + const double* n1c = nc + dim*node1; + const double a = 0.138196601125010515179541316563436; + // Barycentric coordinates of our point in the tet. + double baryc[4] = { a, a, a, a }; + baryc[subindex] = 1.0 - 3.0*a; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd]; + } + return; + } + THROW("Should never reach this point."); } - double quadPtWeight(int /*index*/) const + double quadPtWeight(const int index) const { - return 1.0; + if (degree_ < 2) { + return grid_.cell_volumes[cell_]; + } + // Degree 2 case. + const int dim = grid_.dimensions; + const double* cc = grid_.cell_centroids + dim*cell_; + int tetindex = index / 4; + const double* nc = grid_.node_coordinates; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + if (nfn <= tetindex) { + // Our tet is not associated with this face. + tetindex -= nfn; + continue; + } + const double* fc = grid_.face_centroids + dim*face; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face]; + const int node0 = fnodes[tetindex]; + const int node1 = fnodes[(tetindex + 1) % nfn]; + const double* n0c = nc + dim*node0; + const double* n1c = nc + dim*node1; + return 0.25*tetVolume(cc, fc, n0c, n1c); + } + THROW("Should never reach this point."); } private: @@ -142,7 +289,28 @@ namespace Opm + + /// A class providing numerical quadrature for faces. + /// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i). + /// Note that this class does multiply weights by face area, + /// so weights always sum to face area. + /// Degree 1 method: + /// Midpoint (centroid) method. + /// n = 1, w_0 = face area, x_0 = face centroid + /// Degree 2 method: + /// Based on subdivision of the face into triangles, + /// with the centroid as a common vertex, and the triangle + /// edge midpoint rule. + /// Triangle i consists of the centroid C, nodes N_i and N_{i+1}. + /// Its area is A_i. + /// n = 2 * nn (nn = num nodes in face) + /// For i = 0..(nn-1): + /// w_i = 1/3 A_i. + /// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i + /// x_i = (N_i + N_{i+1})/2 + /// x_{nn+i} = (C + N_i)/2 + /// All N and A indices are interpreted cyclic, modulus nn. class FaceQuadrature { public: @@ -151,25 +319,79 @@ namespace Opm const int degree) : grid_(grid), face_(face), degree_(degree) { - if (degree > 1) { - THROW("Only quadrature degree up to 1 for now."); + if (grid_.dimensions != 3) { + THROW("FaceQuadrature only implemented for 3D case."); + } + if (degree_ > 2) { + THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented."); } } int numQuadPts() const { - return 1; + if (degree_ < 2) { + return 1; + } + // Degree 2 case. + return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]); } - void quadPtCoord(int /*index*/, double* coord) const + void quadPtCoord(const int index, double* coord) const { - const double* fc = grid_.face_centroids + grid_.dimensions*face_; - std::copy(fc, fc + grid_.dimensions, coord); + const int dim = grid_.dimensions; + const double* fc = grid_.face_centroids + dim*face_; + if (degree_ < 2) { + std::copy(fc, fc + dim, coord); + return; + } + // Degree 2 case. + const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_]; + const double* nc = grid_.node_coordinates; + if (index < nn) { + // Boundary edge midpoint. + const int node0 = fnodes[index]; + const int node1 = fnodes[(index + 1)%nn]; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]); + } + } else { + // Interiour edge midpoint. + // Recall that index is now in [nn, 2*nn). + const int node = fnodes[index - nn]; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]); + } + } } - double quadPtWeight(int /*index*/) const + double quadPtWeight(const int index) const { - return 1.0; + if (degree_ < 2) { + return grid_.face_areas[face_]; + } + // Degree 2 case. + const int dim = grid_.dimensions; + const double* fc = grid_.face_centroids + dim*face_; + const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_]; + const double* nc = grid_.node_coordinates; + if (index < nn) { + // Boundary edge midpoint. + const int node0 = fnodes[index]; + const int node1 = fnodes[(index + 1)%nn]; + const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc); + return area / 3.0; + } else { + // Interiour edge midpoint. + // Recall that index is now in [nn, 2*nn). + const int node0 = fnodes[(index - 1) % nn]; + const int node1 = fnodes[index - nn]; + const int node2 = fnodes[(index + 1) % nn]; + const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc); + const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc); + return (area0 + area1) / 3.0; + } } private: @@ -179,6 +401,8 @@ namespace Opm }; + + // Initial version: only a constant interpolation. static void interpolateVelocity(const UnstructuredGrid& grid, const int cell, @@ -206,6 +430,9 @@ namespace Opm + // --------------- Methods of TransportModelTracerTofDiscGal --------------- + + /// Construct solver. /// \param[in] grid A 2d or 3d grid. @@ -307,9 +534,9 @@ namespace Opm // Do quadrature over the face to compute // \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds // (where u_h^{ext} is the upstream unknown (tof)). + const double normal_velocity = flux / grid_.face_areas[face]; FaceQuadrature quad(grid_, face, degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { - // u^ext flux B (B = {b_j}) quad.quadPtCoord(quad_pt, &coord_[0]); DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]); @@ -317,7 +544,7 @@ namespace Opm tof_coeff_ + num_basis*upstream_cell, 0.0); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { - rhs_[j] -= w * tof_upstream * flux * basis_[j]; + rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j]; } } } @@ -325,8 +552,7 @@ namespace Opm // Compute cell jacobian contribution. We use Fortran ordering // for jac_, i.e. rows cycling fastest. { - const double cvol = grid_.cell_volumes[cell]; - CellQuadrature quad(grid_, cell, degree_); + CellQuadrature quad(grid_, cell, 2*degree_ - 1); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // b_i (v \cdot \grad b_j) quad.quadPtCoord(quad_pt, &coord_[0]); @@ -337,7 +563,7 @@ namespace Opm for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { for (int dd = 0; dd < dim; ++dd) { - jac_[j*num_basis + i] += w * basis_[i] * grad_basis_[dim*j + dd] * velocity_[dd] * cvol; + jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd]; } } } @@ -359,7 +585,8 @@ namespace Opm } // Do quadrature over the face to compute // \int_{\partial K} b_i (v(x) \cdot n) b_j ds - FaceQuadrature quad(grid_, face, degree_); + const double normal_velocity = flux / grid_.face_areas[face]; + FaceQuadrature quad(grid_, face, 2*degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // u^ext flux B (B = {b_j}) quad.quadPtCoord(quad_pt, &coord_[0]); @@ -367,26 +594,31 @@ namespace Opm const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { - jac_[j*num_basis + i] += w * basis_[i] * flux * basis_[j]; + jac_[j*num_basis + i] += w * basis_[i] * normal_velocity * basis_[j]; } } } } // Compute downstream jacobian contribution from sink terms. + // Contribution from inflow sources would be + // similar to the contribution from upstream faces, but + // it is zero since we let all external inflow be associated + // with a zero tof. if (source_[cell] < 0.0) { // A sink. const double flux = -source_[cell]; // Sign convention for flux: outflux > 0. + const double flux_density = flux / grid_.cell_volumes[cell]; // Do quadrature over the cell to compute // \int_{K} b_i flux b_j dx - CellQuadrature quad(grid_, cell, degree_); + CellQuadrature quad(grid_, cell, 2*degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { quad.quadPtCoord(quad_pt, &coord_[0]); DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { - jac_[j*num_basis + i] += w * basis_[i] * flux * basis_[j]; + jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j]; } } } From 0eb54ca90a6a47b81b44fd01e8ce26ae036f77f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 1 Oct 2012 16:39:35 +0200 Subject: [PATCH 35/59] Added perfPress() to WellState. --- opm/core/simulator/WellState.hpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 7ddc886e2..981ccd9ba 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -49,7 +49,8 @@ namespace Opm bhp_[w] = state.pressure()[cell]; } } - perfrates_.resize(wells->well_connpos[nw]); + perfrates_.resize(wells->well_connpos[nw], 0.0); + perfpress_.resize(wells->well_connpos[nw], -1e100); } } @@ -61,9 +62,14 @@ namespace Opm std::vector& perfRates() { return perfrates_; } const std::vector& perfRates() const { return perfrates_; } + /// One pressure per well connection. + std::vector& perfPress() { return perfpress_; } + const std::vector& perfPress() const { return perfpress_; } + private: std::vector bhp_; std::vector perfrates_; + std::vector perfpress_; }; } // namespace Opm From 6fb248d403ef80a1bbd15a81f5a83c3dcf95030c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 1 Oct 2012 16:40:10 +0200 Subject: [PATCH 36/59] Update WellState::perfPress() after pressure solve. --- opm/core/pressure/CompressibleTpfa.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 42b3fba75..8edb41a0d 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -632,6 +632,23 @@ namespace Opm &state.pressure()[0], &state.faceflux()[0], &state.facepressure()[0]); + + // Compute well perforation pressures (not done by the C code). + if (wells_ != 0) { + const int nw = wells_->number_of_wells; + const int np = props_.numPhases(); + for (int w = 0; w < nw; ++w) { + const double* comp_frac = &wells_->comp_frac[np*w]; + for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) { + const double bhp = well_state.bhp()[w]; + double perf_p = bhp; + for (int phase = 0; phase < np; ++phase) { + perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase]; + } + well_state.perfPress()[j] = perf_p; + } + } + } } } // namespace Opm From 17c1be65418acd6959ecd9f0ca3cfd3c3d876a3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 1 Oct 2012 16:40:47 +0200 Subject: [PATCH 37/59] Modified functions dealing with transport source. In preparation for switching to new convention for inflow sources in the compressible case: source being surface volumes, not reservoir volumes. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 78 +++++++++++++++++++++- opm/core/utility/miscUtilitiesBlackoil.hpp | 19 +++++- 2 files changed, 93 insertions(+), 4 deletions(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 6f4f75565..b4877e7e2 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -21,7 +21,10 @@ #include #include #include +#include #include +#include +#include #include #include #include @@ -31,16 +34,20 @@ namespace Opm { - /// @brief Computes injected and produced volumes of all phases. + /// @brief Computes injected and produced surface volumes of all phases. /// Note 1: assumes that only the first phase is injected. /// Note 2: assumes that transport has been done with an /// implicit method, i.e. that the current state /// gives the mobilities used for the preceding timestep. + /// Note 3: Gives surface volume values, not reservoir volumes + /// (as the incompressible version of the function does). + /// Also, assumes that src is given in surface volumes + /// for injector terms! /// @param[in] props fluid and rock properties. /// @param[in] p pressure (one value per cell) /// @param[in] z surface-volume values (for all P phases) /// @param[in] s saturation values (for all P phases) - /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. + /// @param[in] src if < 0: total outflow, if > 0: first phase inflow /// @param[in] dt timestep used /// @param[out] injected must point to a valid array with P elements, /// where P = s.size()/src.size(). @@ -63,6 +70,9 @@ namespace Opm std::fill(produced, produced + np, 0.0); std::vector visc(np); std::vector mob(np); + std::vector A(np*np); + std::vector prod_resv_phase(np); + std::vector prod_surfvol(np); for (int c = 0; c < num_cells; ++c) { if (src[c] > 0.0) { injected[0] += src[c]*dt; @@ -71,13 +81,21 @@ namespace Opm const double* sat = &s[np*c]; props.relperm(1, sat, &c, &mob[0], 0); props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0); + props.matrix(1, &press[c], &z[np*c], &c, &A[0], 0); double totmob = 0.0; for (int p = 0; p < np; ++p) { mob[p] /= visc[p]; totmob += mob[p]; } + std::fill(prod_surfvol.begin(), prod_surfvol.end(), 0.0); for (int p = 0; p < np; ++p) { - produced[p] += (mob[p]/totmob)*flux; + prod_resv_phase[p] = (mob[p]/totmob)*flux; + for (int q = 0; q < np; ++q) { + prod_surfvol[q] += prod_resv_phase[p]*A[q + np*p]; + } + } + for (int p = 0; p < np; ++p) { + produced[p] += prod_surfvol[p]; } } } @@ -251,4 +269,58 @@ namespace Opm } } + + /// Compute two-phase transport source terms from well terms. + /// Note: Unlike the incompressible version of this function, + /// this version computes surface volume injection rates, + /// production rates are still total reservoir volumes. + /// \param[in] props Fluid and rock properties. + /// \param[in] wells Wells data structure. + /// \param[in] well_state Well pressures and fluxes. + /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign: + /// (+) positive inflow of first (water) phase (surface volume), + /// (-) negative total outflow of both phases (reservoir volume). + void computeTransportSource(const BlackoilPropertiesInterface& props, + const Wells* wells, + const WellState& well_state, + std::vector& transport_src) + { + int nc = props.numCells(); + transport_src.clear(); + transport_src.resize(nc, 0.0); + // Well contributions. + if (wells) { + const int nw = wells->number_of_wells; + const int np = wells->number_of_phases; + if (np != 2) { + THROW("computeTransportSource() requires a 2 phase case."); + } + std::vector A(np, np); + for (int w = 0; w < nw; ++w) { + const double* comp_frac = wells->comp_frac + np*w; + for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) { + const int perf_cell = wells->well_cells[perf]; + double perf_rate = well_state.perfRates()[perf]; + if (perf_rate > 0.0) { + // perf_rate is a total inflow reservoir rate, we want a surface water rate. + if (wells->type[w] != INJECTOR) { + std::cout << "**** Warning: crossflow in well " + << w << " perf " << perf - wells->well_connpos[w] + << " ignored. Reservoir rate was " + << perf_rate/Opm::unit::day << " m^3/day." << std::endl; + perf_rate = 0.0; + } else { + ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6); + perf_rate *= comp_frac[0]; // Water reservoir volume rate. + props.matrix(1, &well_state.perfPress()[perf], comp_frac, &perf_cell, &A[0], 0); + perf_rate *= A[0]; // Water surface volume rate. + } + } + transport_src[perf_cell] += perf_rate; + } + } + } + } + + } // namespace Opm diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 565acc6ef..0225787ac 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -22,12 +22,13 @@ #include -struct UnstructuredGrid; +struct Wells; namespace Opm { class BlackoilPropertiesInterface; + class WellState; /// @brief Computes injected and produced volumes of all phases. /// Note 1: assumes that only the first phase is injected. @@ -131,6 +132,22 @@ namespace Opm const double* saturation, double* surfacevol); + + /// Compute two-phase transport source terms from well terms. + /// Note: Unlike the incompressible version of this function, + /// this version computes surface volume injection rates, + /// production rates are still total reservoir volumes. + /// \param[in] props Fluid and rock properties. + /// \param[in] wells Wells data structure. + /// \param[in] well_state Well pressures and fluxes. + /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign: + /// (+) positive inflow of first (water) phase (surface volume), + /// (-) negative total outflow of both phases (reservoir volume). + void computeTransportSource(const BlackoilPropertiesInterface& props, + const Wells* wells, + const WellState& well_state, + std::vector& transport_src); + } // namespace Opm #endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED From 017663cc5b0f84b91edcd34ed6aca55010b08317 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 2 Oct 2012 11:12:23 +0200 Subject: [PATCH 38/59] Change interface for (blackoil) computeInjectedProduced(). Also use new computeInjectedProduced() and computeTransportSource() functions in SimulatorCompressibleTwophase. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 46 +++++++++++----------- opm/core/utility/miscUtilitiesBlackoil.hpp | 32 ++++++++------- 2 files changed, 42 insertions(+), 36 deletions(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index b4877e7e2..4f49cd6aa 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -41,31 +41,33 @@ namespace Opm /// gives the mobilities used for the preceding timestep. /// Note 3: Gives surface volume values, not reservoir volumes /// (as the incompressible version of the function does). - /// Also, assumes that src is given in surface volumes + /// Also, assumes that transport_src is given in surface volumes /// for injector terms! - /// @param[in] props fluid and rock properties. - /// @param[in] p pressure (one value per cell) - /// @param[in] z surface-volume values (for all P phases) - /// @param[in] s saturation values (for all P phases) - /// @param[in] src if < 0: total outflow, if > 0: first phase inflow - /// @param[in] dt timestep used - /// @param[out] injected must point to a valid array with P elements, - /// where P = s.size()/src.size(). - /// @param[out] produced must also point to a valid array with P elements. + /// @param[in] props fluid and rock properties. + /// @param[in] state state variables (pressure, sat, surfvol) + /// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow + /// @param[in] dt timestep used + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/src.size(). + /// @param[out] produced must also point to a valid array with P elements. void computeInjectedProduced(const BlackoilPropertiesInterface& props, - const std::vector& press, - const std::vector& z, - const std::vector& s, - const std::vector& src, + const BlackoilState& state, + const std::vector& transport_src, const double dt, double* injected, double* produced) { - const int num_cells = src.size(); - const int np = s.size()/src.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and src vectors do not match."); + const int num_cells = transport_src.size(); + if (props.numCells() != num_cells) { + THROW("Size of transport_src vector does not match number of cells in props."); } + const int np = props.numPhases(); + if (int(state.saturation().size()) != num_cells*np) { + THROW("Sizes of state vectors do not match number of cells."); + } + const std::vector& press = state.pressure(); + const std::vector& s = state.saturation(); + const std::vector& z = state.surfacevol(); std::fill(injected, injected + np, 0.0); std::fill(produced, produced + np, 0.0); std::vector visc(np); @@ -74,10 +76,10 @@ namespace Opm std::vector prod_resv_phase(np); std::vector prod_surfvol(np); for (int c = 0; c < num_cells; ++c) { - if (src[c] > 0.0) { - injected[0] += src[c]*dt; - } else if (src[c] < 0.0) { - const double flux = -src[c]*dt; + if (transport_src[c] > 0.0) { + injected[0] += transport_src[c]*dt; + } else if (transport_src[c] < 0.0) { + const double flux = -transport_src[c]*dt; const double* sat = &s[np*c]; props.relperm(1, sat, &c, &mob[0], 0); props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0); diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 0225787ac..f4949ce2d 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -28,31 +28,34 @@ namespace Opm { class BlackoilPropertiesInterface; + class BlackoilState; class WellState; - /// @brief Computes injected and produced volumes of all phases. + + /// @brief Computes injected and produced surface volumes of all phases. /// Note 1: assumes that only the first phase is injected. /// Note 2: assumes that transport has been done with an /// implicit method, i.e. that the current state /// gives the mobilities used for the preceding timestep. - /// @param[in] props fluid and rock properties. - /// @param[in] p pressure (one value per cell) - /// @param[in] z surface-volume values (for all P phases) - /// @param[in] s saturation values (for all P phases) - /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. - /// @param[in] dt timestep used - /// @param[out] injected must point to a valid array with P elements, - /// where P = s.size()/src.size(). - /// @param[out] produced must also point to a valid array with P elements. + /// Note 3: Gives surface volume values, not reservoir volumes + /// (as the incompressible version of the function does). + /// Also, assumes that transport_src is given in surface volumes + /// for injector terms! + /// @param[in] props fluid and rock properties. + /// @param[in] state state variables (pressure, sat, surfvol) + /// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow + /// @param[in] dt timestep used + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/src.size(). + /// @param[out] produced must also point to a valid array with P elements. void computeInjectedProduced(const BlackoilPropertiesInterface& props, - const std::vector& p, - const std::vector& z, - const std::vector& s, - const std::vector& src, + const BlackoilState& state, + const std::vector& transport_src, const double dt, double* injected, double* produced); + /// @brief Computes total mobility for a set of saturation values. /// @param[in] props rock and fluid properties /// @param[in] cells cells with which the saturation values are associated @@ -67,6 +70,7 @@ namespace Opm const std::vector& s, std::vector& totmob); + /// @brief Computes total mobility and omega for a set of saturation values. /// @param[in] props rock and fluid properties /// @param[in] cells cells with which the saturation values are associated From 4f276a887056ddf22fd1264f0717a9c9b96061a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 2 Oct 2012 14:35:28 +0200 Subject: [PATCH 39/59] Bugfix: size of vector for A should be np*np. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 4f49cd6aa..62f179932 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -297,7 +297,7 @@ namespace Opm if (np != 2) { THROW("computeTransportSource() requires a 2 phase case."); } - std::vector A(np, np); + std::vector A(np*np); for (int w = 0; w < nw; ++w) { const double* comp_frac = wells->comp_frac + np*w; for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) { From 5acbf00b14b4dfb781b509ceac5d75bfa7a84dcd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 2 Oct 2012 14:37:18 +0200 Subject: [PATCH 40/59] Update compressible transport solver for new src convention. Namely, that inflowing transport sources are water *surface* volumes, not water *reservoir* volumes. --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 7b4a72e68..bad712630 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -152,7 +152,7 @@ namespace Opm B_cell = 1.0/tm.A_[np*np*cell + 0]; double src_flux = -tm.source_[cell]; bool src_is_inflow = src_flux < 0.0; - influx = src_is_inflow ? src_flux : 0.0; + influx = src_is_inflow ? B_cell* src_flux : 0.0; outflux = !src_is_inflow ? src_flux : 0.0; comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell]; dtpv = tm.dt_/tm.porevolume0_[cell]; From a9783eefc7e22cc77eb88cd13d11e8346c47e819 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 2 Oct 2012 15:46:33 +0200 Subject: [PATCH 41/59] Add explanatory comment. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 62f179932..ae806e9a4 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -77,8 +77,12 @@ namespace Opm std::vector prod_surfvol(np); for (int c = 0; c < num_cells; ++c) { if (transport_src[c] > 0.0) { + // Inflowing transport source is a surface volume flux + // for the first phase. injected[0] += transport_src[c]*dt; } else if (transport_src[c] < 0.0) { + // Outflowing transport source is a total reservoir + // volume flux. const double flux = -transport_src[c]*dt; const double* sat = &s[np*c]; props.relperm(1, sat, &c, &mob[0], 0); From a92549290f50b30cb646be10a59d54f1b1005153 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 4 Oct 2012 21:09:47 +0200 Subject: [PATCH 42/59] New function clone_wells() Used to create a deep copy (clone) of an existing Wells object. While here, add test case for common Wells object operations. --- opm/core/newwells.h | 13 +++ tests/test_wells.cpp | 205 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 218 insertions(+) create mode 100644 tests/test_wells.cpp diff --git a/opm/core/newwells.h b/opm/core/newwells.h index f16675b6c..20cfce61d 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -235,6 +235,19 @@ void destroy_wells(struct Wells *W); +/** + * Create a deep-copy (i.e., clone) of an existing Wells object, including its + * controls. + * + * @param[in] W Existing Wells object. + * @return Complete clone of the input object. Dispose of resources using + * function destroy_wells() when no longer needed. Returns @c NULL in case of + * allocation failure. + */ +struct Wells * +clone_wells(const struct Wells *W); + + #ifdef __cplusplus } #endif diff --git a/tests/test_wells.cpp b/tests/test_wells.cpp new file mode 100644 index 000000000..4f8ad7bb7 --- /dev/null +++ b/tests/test_wells.cpp @@ -0,0 +1,205 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE WellsModuleTest +#include + +#include + +#include +#include +#include + +BOOST_AUTO_TEST_CASE(Construction) +{ + const int nphases = 2; + const int nwells = 2; + const int nperfs = 2; + + boost::shared_ptr W(create_wells(nphases, nwells, nperfs), + destroy_wells); + + if (W) { + int cells[] = { 0, 9 }; + double WI = 1.0; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0], + &WI, "INJECTOR", W.get()); + + const double pfrac[] = { 0.0, 0.0 }; + const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1], + &WI, "PRODUCER", W.get()); + + if (ok0 && ok1) { + BOOST_CHECK_EQUAL(W->number_of_phases, nphases); + BOOST_CHECK_EQUAL(W->number_of_wells , nwells ); + + BOOST_CHECK_EQUAL(W->well_connpos[0], 0); + BOOST_CHECK_EQUAL(W->well_connpos[1], 1); + BOOST_CHECK_EQUAL(W->well_connpos[W->number_of_wells], nperfs); + + BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[0]], cells[0]); + BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[1]], cells[1]); + + BOOST_CHECK_EQUAL(W->WI[W->well_connpos[0]], WI); + BOOST_CHECK_EQUAL(W->WI[W->well_connpos[1]], WI); + + using std::string; + BOOST_CHECK_EQUAL(string(W->name[0]), string("INJECTOR")); + BOOST_CHECK_EQUAL(string(W->name[1]), string("PRODUCER")); + } + } +} + + +BOOST_AUTO_TEST_CASE(Controls) +{ + const int nphases = 2; + const int nwells = 1; + const int nperfs = 2; + + boost::shared_ptr W(create_wells(nphases, nwells, nperfs), + destroy_wells); + + if (W) { + int cells[] = { 0 , 9 }; + double WI [] = { 1.0, 1.0 }; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok = add_well(INJECTOR, 0.0, nperfs, &ifrac[0], &cells[0], + &WI[0], "INJECTOR", W.get()); + + if (ok) { + const double distr[] = { 1.0, 0.0 }; + const bool ok1 = append_well_controls(BHP, 1, &distr[0], + 0, W.get()); + const bool ok2 = append_well_controls(SURFACE_RATE, 1, + &distr[0], 0, W.get()); + + if (ok1 && ok2) { + WellControls* ctrls = W->ctrls[0]; + + BOOST_CHECK_EQUAL(ctrls->num , 2); + BOOST_CHECK_EQUAL(ctrls->current, -1); + + set_current_control(0, 0, W.get()); + BOOST_CHECK_EQUAL(ctrls->current, 0); + + set_current_control(0, 1, W.get()); + BOOST_CHECK_EQUAL(ctrls->current, 1); + + BOOST_CHECK_EQUAL(ctrls->type[0], BHP); + BOOST_CHECK_EQUAL(ctrls->type[1], SURFACE_RATE); + + BOOST_CHECK_EQUAL(ctrls->target[0], 1.0); + BOOST_CHECK_EQUAL(ctrls->target[1], 1.0); + } + } + } +} + + +BOOST_AUTO_TEST_CASE(Copy) +{ + const int nphases = 2; + const int nwells = 2; + const int nperfs = 2; + + boost::shared_ptr W1(create_wells(nphases, nwells, nperfs), + destroy_wells); + boost::shared_ptr W2; + + if (W1) { + int cells[] = { 0, 9 }; + const double WI = 1.0; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0], + &WI, "INJECTOR", W1.get()); + + const double pfrac[] = { 0.0, 0.0 }; + const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1], + &WI, "PRODUCER", W1.get()); + + bool ok = ok0 && ok1; + for (int w = 0; ok && (w < W1->number_of_wells); ++w) { + const double distr[] = { 1.0, 0.0 }; + const bool okc1 = append_well_controls(BHP, 1, &distr[0], + w, W1.get()); + const bool okc2 = append_well_controls(SURFACE_RATE, 1, + &distr[0], w, + W1.get()); + + ok = okc1 && okc2; + } + + if (ok) { + W2.reset(clone_wells(W1.get()), destroy_wells); + } + } + + if (W2) { + BOOST_CHECK_EQUAL(W2->number_of_phases, W1->number_of_phases); + BOOST_CHECK_EQUAL(W2->number_of_wells , W1->number_of_wells ); + BOOST_CHECK_EQUAL(W2->well_connpos[0] , W1->well_connpos[0] ); + + for (int w = 0; w < W1->number_of_wells; ++w) { + using std::string; + BOOST_CHECK_EQUAL(string(W2->name[w]), string(W1->name[w])); + BOOST_CHECK_EQUAL( W2->type[w] , W1->type[w] ); + + BOOST_CHECK_EQUAL(W2->well_connpos[w + 1], + W1->well_connpos[w + 1]); + + for (int j = W1->well_connpos[w]; + j < W1->well_connpos[w + 1]; ++j) { + BOOST_CHECK_EQUAL(W2->well_cells[j], W1->well_cells[j]); + BOOST_CHECK_EQUAL(W2->WI [j], W1->WI [j]); + } + + BOOST_CHECK(W1->ctrls[w] != 0); + BOOST_CHECK(W2->ctrls[w] != 0); + + WellControls* c1 = W1->ctrls[w]; + WellControls* c2 = W2->ctrls[w]; + + BOOST_CHECK_EQUAL(c2->num , c1->num ); + BOOST_CHECK_EQUAL(c2->current, c1->current); + + for (int c = 0; c < c1->num; ++c) { + BOOST_CHECK_EQUAL(c2->type [c], c1->type [c]); + BOOST_CHECK_EQUAL(c2->target[c], c1->target[c]); + + for (int p = 0; p < W1->number_of_phases; ++p) { + BOOST_CHECK_EQUAL(c2->distr[c*W1->number_of_phases + p], + c1->distr[c*W1->number_of_phases + p]); + } + } + } + } +} From 4d488c98a79cbc38e4832d982961b56acbe4accc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 24 Aug 2012 16:35:34 +0200 Subject: [PATCH 43/59] Unequivocally exclude MATLAB timing printing. It is not actually needed and prevents building when symbol MATLAB_MEX_FILE is defined. --- opm/core/transport/reorder/TransportModelInterface.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp index c52b528d6..f8e138a34 100644 --- a/opm/core/transport/reorder/TransportModelInterface.cpp +++ b/opm/core/transport/reorder/TransportModelInterface.cpp @@ -35,6 +35,7 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g // Invoke appropriate solve method for each interdependent component. for (int comp = 0; comp < ncomponents; ++comp) { +#if 0 #ifdef MATLAB_MEX_FILE // \TODO replace this with general signal handling code, check if it costs performance. if (interrupt_signal) { @@ -42,6 +43,7 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g "cells finished.\n", i, grid.number_of_cells); break; } +#endif #endif const int comp_size = components[comp + 1] - components[comp]; if (comp_size == 1) { From 24504875ea18adca065faacdca1a24e7203c2f7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 24 Aug 2012 19:16:31 +0200 Subject: [PATCH 44/59] Reference from canonical location. The header was removed from this directory upon import from the preexisting "opmtransport" repository. --- opm/core/transport/reorder/reordersequence.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/reordersequence.cpp b/opm/core/transport/reorder/reordersequence.cpp index 29f3a0aa2..0e7b9f946 100644 --- a/opm/core/transport/reorder/reordersequence.cpp +++ b/opm/core/transport/reorder/reordersequence.cpp @@ -1,11 +1,11 @@ /* Copyright 2011 (c) Jostein R. Natvig */ +#include + #ifdef MATLAB_MEX_FILE -#include "grid.h" #include "reordersequence.h" #include "tarjan.h" #else -#include #include #include #endif From e5c5a64a4adf20395a6401b692fb6865393cb121 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 4 Oct 2012 21:09:47 +0200 Subject: [PATCH 45/59] New function clone_wells() Used to create a deep copy (clone) of an existing Wells object. While here, add test case for common Wells object operations. --- opm/core/newwells.h | 13 +++ tests/test_wells.cpp | 205 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 218 insertions(+) create mode 100644 tests/test_wells.cpp diff --git a/opm/core/newwells.h b/opm/core/newwells.h index f16675b6c..20cfce61d 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -235,6 +235,19 @@ void destroy_wells(struct Wells *W); +/** + * Create a deep-copy (i.e., clone) of an existing Wells object, including its + * controls. + * + * @param[in] W Existing Wells object. + * @return Complete clone of the input object. Dispose of resources using + * function destroy_wells() when no longer needed. Returns @c NULL in case of + * allocation failure. + */ +struct Wells * +clone_wells(const struct Wells *W); + + #ifdef __cplusplus } #endif diff --git a/tests/test_wells.cpp b/tests/test_wells.cpp new file mode 100644 index 000000000..4f8ad7bb7 --- /dev/null +++ b/tests/test_wells.cpp @@ -0,0 +1,205 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE WellsModuleTest +#include + +#include + +#include +#include +#include + +BOOST_AUTO_TEST_CASE(Construction) +{ + const int nphases = 2; + const int nwells = 2; + const int nperfs = 2; + + boost::shared_ptr W(create_wells(nphases, nwells, nperfs), + destroy_wells); + + if (W) { + int cells[] = { 0, 9 }; + double WI = 1.0; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0], + &WI, "INJECTOR", W.get()); + + const double pfrac[] = { 0.0, 0.0 }; + const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1], + &WI, "PRODUCER", W.get()); + + if (ok0 && ok1) { + BOOST_CHECK_EQUAL(W->number_of_phases, nphases); + BOOST_CHECK_EQUAL(W->number_of_wells , nwells ); + + BOOST_CHECK_EQUAL(W->well_connpos[0], 0); + BOOST_CHECK_EQUAL(W->well_connpos[1], 1); + BOOST_CHECK_EQUAL(W->well_connpos[W->number_of_wells], nperfs); + + BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[0]], cells[0]); + BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[1]], cells[1]); + + BOOST_CHECK_EQUAL(W->WI[W->well_connpos[0]], WI); + BOOST_CHECK_EQUAL(W->WI[W->well_connpos[1]], WI); + + using std::string; + BOOST_CHECK_EQUAL(string(W->name[0]), string("INJECTOR")); + BOOST_CHECK_EQUAL(string(W->name[1]), string("PRODUCER")); + } + } +} + + +BOOST_AUTO_TEST_CASE(Controls) +{ + const int nphases = 2; + const int nwells = 1; + const int nperfs = 2; + + boost::shared_ptr W(create_wells(nphases, nwells, nperfs), + destroy_wells); + + if (W) { + int cells[] = { 0 , 9 }; + double WI [] = { 1.0, 1.0 }; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok = add_well(INJECTOR, 0.0, nperfs, &ifrac[0], &cells[0], + &WI[0], "INJECTOR", W.get()); + + if (ok) { + const double distr[] = { 1.0, 0.0 }; + const bool ok1 = append_well_controls(BHP, 1, &distr[0], + 0, W.get()); + const bool ok2 = append_well_controls(SURFACE_RATE, 1, + &distr[0], 0, W.get()); + + if (ok1 && ok2) { + WellControls* ctrls = W->ctrls[0]; + + BOOST_CHECK_EQUAL(ctrls->num , 2); + BOOST_CHECK_EQUAL(ctrls->current, -1); + + set_current_control(0, 0, W.get()); + BOOST_CHECK_EQUAL(ctrls->current, 0); + + set_current_control(0, 1, W.get()); + BOOST_CHECK_EQUAL(ctrls->current, 1); + + BOOST_CHECK_EQUAL(ctrls->type[0], BHP); + BOOST_CHECK_EQUAL(ctrls->type[1], SURFACE_RATE); + + BOOST_CHECK_EQUAL(ctrls->target[0], 1.0); + BOOST_CHECK_EQUAL(ctrls->target[1], 1.0); + } + } + } +} + + +BOOST_AUTO_TEST_CASE(Copy) +{ + const int nphases = 2; + const int nwells = 2; + const int nperfs = 2; + + boost::shared_ptr W1(create_wells(nphases, nwells, nperfs), + destroy_wells); + boost::shared_ptr W2; + + if (W1) { + int cells[] = { 0, 9 }; + const double WI = 1.0; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0], + &WI, "INJECTOR", W1.get()); + + const double pfrac[] = { 0.0, 0.0 }; + const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1], + &WI, "PRODUCER", W1.get()); + + bool ok = ok0 && ok1; + for (int w = 0; ok && (w < W1->number_of_wells); ++w) { + const double distr[] = { 1.0, 0.0 }; + const bool okc1 = append_well_controls(BHP, 1, &distr[0], + w, W1.get()); + const bool okc2 = append_well_controls(SURFACE_RATE, 1, + &distr[0], w, + W1.get()); + + ok = okc1 && okc2; + } + + if (ok) { + W2.reset(clone_wells(W1.get()), destroy_wells); + } + } + + if (W2) { + BOOST_CHECK_EQUAL(W2->number_of_phases, W1->number_of_phases); + BOOST_CHECK_EQUAL(W2->number_of_wells , W1->number_of_wells ); + BOOST_CHECK_EQUAL(W2->well_connpos[0] , W1->well_connpos[0] ); + + for (int w = 0; w < W1->number_of_wells; ++w) { + using std::string; + BOOST_CHECK_EQUAL(string(W2->name[w]), string(W1->name[w])); + BOOST_CHECK_EQUAL( W2->type[w] , W1->type[w] ); + + BOOST_CHECK_EQUAL(W2->well_connpos[w + 1], + W1->well_connpos[w + 1]); + + for (int j = W1->well_connpos[w]; + j < W1->well_connpos[w + 1]; ++j) { + BOOST_CHECK_EQUAL(W2->well_cells[j], W1->well_cells[j]); + BOOST_CHECK_EQUAL(W2->WI [j], W1->WI [j]); + } + + BOOST_CHECK(W1->ctrls[w] != 0); + BOOST_CHECK(W2->ctrls[w] != 0); + + WellControls* c1 = W1->ctrls[w]; + WellControls* c2 = W2->ctrls[w]; + + BOOST_CHECK_EQUAL(c2->num , c1->num ); + BOOST_CHECK_EQUAL(c2->current, c1->current); + + for (int c = 0; c < c1->num; ++c) { + BOOST_CHECK_EQUAL(c2->type [c], c1->type [c]); + BOOST_CHECK_EQUAL(c2->target[c], c1->target[c]); + + for (int p = 0; p < W1->number_of_phases; ++p) { + BOOST_CHECK_EQUAL(c2->distr[c*W1->number_of_phases + p], + c1->distr[c*W1->number_of_phases + p]); + } + } + } + } +} From 044c13e9fb98a731c2138e23e2c24fecd9759516 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Wed, 5 Sep 2012 11:47:38 +0200 Subject: [PATCH 46/59] Added constuctor to WellManager which used Wells struct. --- opm/core/wells/WellsManager.cpp | 4 ++++ opm/core/wells/WellsManager.hpp | 5 ++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 95cb0369a..fa142309e 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -223,6 +223,10 @@ namespace Opm : w_(0) { } + WellsManager::WellsManager(struct Wells* W) + : w_(W) + { + } diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index 6df65ec41..650eaf371 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -41,7 +41,10 @@ namespace Opm { public: /// Default constructor -- no wells. - WellsManager(); + WellsManager(); + /// Construct from mrst type output. + /// Wellmanger is not properly initialized + WellsManager(struct Wells* W); /// Construct from input deck and grid. /// The permeability argument may be zero if the input contain From d6154d89618724a367055fc6e1bf2c651e1be2ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 4 Oct 2012 23:46:20 +0200 Subject: [PATCH 47/59] Clone wells object when constructing from existing. This installs a measure of safety on the part of the interface in that the caller is free to dispose of the wells object upon returning from the WellsManager constructor. --- opm/core/wells/WellsManager.cpp | 6 +++++- opm/core/wells/WellsManager.hpp | 8 ++++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index fa142309e..8aaca9e2e 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -223,8 +223,12 @@ namespace Opm : w_(0) { } + + + + /// Construct from existing wells object. WellsManager::WellsManager(struct Wells* W) - : w_(W) + : w_(clone_wells(W)) { } diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index 650eaf371..19289c5d9 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -42,8 +42,12 @@ namespace Opm public: /// Default constructor -- no wells. WellsManager(); - /// Construct from mrst type output. - /// Wellmanger is not properly initialized + + /// Construct from existing wells object. + /// WellsManager is not properly initialised in the sense that the logic to + /// manage control switching does not exist. + /// + /// @param[in] W Existing wells object. WellsManager(struct Wells* W); /// Construct from input deck and grid. From 9e90dcebe5c727fb4fb9a8f32833a2917986164e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 8 Oct 2012 14:27:56 +0200 Subject: [PATCH 48/59] Fix sign of production rate controls. In the Wells struct, production rate control targets must be negative (and injection rate control targets are always positive). In the WellsGroup classes, there are separate variables for injection and production, and all rates are positive. Therefore, upon adding or modification of a production rate control, the negated value must be used. --- opm/core/wells/WellsGroup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index afa17f12d..c340c1aec 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -863,7 +863,7 @@ namespace Opm return; } // We're a producer, so we need to negate the input - double ntarget = target; + double ntarget = -target; double distr[3] = { 0.0, 0.0, 0.0 }; const int* phase_pos = phaseUsage().phase_pos; From 36721602b293fa501ff40edba27f72f76bb4d67e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 9 Oct 2012 09:54:26 +0200 Subject: [PATCH 49/59] Add timing of topological sort. --- opm/core/transport/reorder/TransportModelInterface.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp index f8e138a34..8a90ae34f 100644 --- a/opm/core/transport/reorder/TransportModelInterface.cpp +++ b/opm/core/transport/reorder/TransportModelInterface.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -31,7 +32,11 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g std::vector sequence(grid.number_of_cells); std::vector components(grid.number_of_cells + 1); int ncomponents; + time::StopWatch clock; + clock.start(); compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents); + clock.stop(); + std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl; // Invoke appropriate solve method for each interdependent component. for (int comp = 0; comp < ncomponents; ++comp) { From 76259dcd0fef1dcde9d760a1e9371ab561df71e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 9 Oct 2012 09:54:54 +0200 Subject: [PATCH 50/59] Improve docs and give more info on error. --- .../reorder/TransportModelTracerTof.cpp | 6 ++++-- .../reorder/TransportModelTracerTof.hpp | 14 ++++++++++--- .../TransportModelTracerTofDiscGal.cpp | 21 ++++++++++++++++--- .../TransportModelTracerTofDiscGal.hpp | 17 ++++++++++++--- 4 files changed, 47 insertions(+), 11 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 70c2fc400..8d93f2480 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -38,10 +38,12 @@ namespace Opm - /// Solve for time-of-flight at next timestep. + /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. /// \param[out] tof Array of time-of-flight values. void TransportModelTracerTof::solveTof(const double* darcyflux, const double* porevolume, diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp index 21375d895..270ac328a 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -31,7 +31,13 @@ namespace Opm class IncompPropertiesInterface; - /// Implements a reordering transport solver for incompressible two-phase flow. + /// Implements a first-order finite volume solver for + /// (single-phase) time-of-flight using reordering. + /// The equation solved is: + /// v \cdot \grad\tau = \phi + /// where v is the fluid velocity, \tau is time-of-flight and + /// \phi is the porosity. This is a boundary value problem, where + /// \tau is specified to be zero on all inflow boundaries. class TransportModelTracerTof : public TransportModelInterface { public: @@ -39,10 +45,12 @@ namespace Opm /// \param[in] grid A 2d or 3d grid. TransportModelTracerTof(const UnstructuredGrid& grid); - /// Solve for time-of-flight at next timestep. + /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. /// \param[out] tof Array of time-of-flight values. void solveTof(const double* darcyflux, const double* porevolume, diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index e9e812ef3..4bc2f1efb 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -446,10 +446,12 @@ namespace Opm - /// Solve for time-of-flight at next timestep. + /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. /// \param[in] degree Polynomial degree of DG basis functions used. /// \param[out] tof_coeff Array of time-of-flight solution coefficients. /// The values are ordered by cell, meaning that @@ -633,7 +635,20 @@ namespace Opm MAT_SIZE_T info = 0; dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info); if (info != 0) { - THROW("Lapack error: " << info); + // Print the local matrix and rhs. + std::cerr << "Failed solving single-cell system Ax = b in cell " << cell + << " with A = \n"; + for (int row = 0; row < n; ++row) { + for (int col = 0; col < n; ++col) { + std::cerr << " " << jac_[row + n*col]; + } + std::cerr << '\n'; + } + std::cerr << "and b = \n"; + for (int row = 0; row < n; ++row) { + std::cerr << " " << rhs_[row] << '\n'; + } + THROW("Lapack error: " << info << " encountered in cell " << cell); } // The solution ends up in rhs_, so we must copy it. std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index b9163f8ed..0285e1f02 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -31,7 +31,15 @@ namespace Opm class IncompPropertiesInterface; - /// Implements a reordering transport solver for incompressible two-phase flow. + /// Implements a discontinuous Galerkin solver for + /// (single-phase) time-of-flight using reordering. + /// The equation solved is: + /// v \cdot \grad\tau = \phi + /// where v is the fluid velocity, \tau is time-of-flight and + /// \phi is the porosity. This is a boundary value problem, where + /// \tau is specified to be zero on all inflow boundaries. + /// The user may specify the polynomial degree of the basis function space + /// used, but only degrees 0 and 1 are supported so far. class TransportModelTracerTofDiscGal : public TransportModelInterface { public: @@ -39,10 +47,13 @@ namespace Opm /// \param[in] grid A 2d or 3d grid. TransportModelTracerTofDiscGal(const UnstructuredGrid& grid); - /// Solve for time-of-flight at next timestep. + + /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. /// \param[in] degree Polynomial degree of DG basis functions used. /// \param[out] tof_coeff Array of time-of-flight solution coefficients. /// The values are ordered by cell, meaning that From 31d17a0dcdda0dcbd37e0558247487f093732c0b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 9 Oct 2012 12:21:17 +0200 Subject: [PATCH 51/59] Whitespace cleanup. --- .../reorder/TransportModelTracerTofDiscGal.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 4bc2f1efb..2d58ec918 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -109,9 +109,9 @@ namespace Opm static void cross(const double* a, const double* b, double* res) { - res[0] = a[1]*b[2] - a[2]*b[1]; - res[1] = a[2]*b[0] - a[0]*b[2]; - res[2] = a[0]*b[1] - a[1]*b[0]; + res[0] = a[1]*b[2] - a[2]*b[1]; + res[1] = a[2]*b[0] - a[0]*b[2]; + res[2] = a[0]*b[1] - a[1]*b[0]; } @@ -135,12 +135,12 @@ namespace Opm /// three three-dimensional arrays. static double determinantOf(const double* a0, const double* a1, - const double* a2) + const double* a2) { - return - a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - - a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + - a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); + return + a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - + a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + + a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); } @@ -156,7 +156,7 @@ namespace Opm double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] }; - return std::fabs(determinantOf(a, b, c) / 6.0); + return std::fabs(determinantOf(a, b, c) / 6.0); } From a740bcfe5f03ae9edebcb50b885395d1a0319ea8 Mon Sep 17 00:00:00 2001 From: kristinf Date: Tue, 9 Oct 2012 15:17:46 +0200 Subject: [PATCH 52/59] sign error in total_produced --- opm/core/wells/WellsGroup.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index c340c1aec..d407c0335 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -572,6 +572,7 @@ namespace Opm break; } const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase); + const double total_reinjected = - total_produced; // Production negative, injection positive const double my_guide_rate = injectionGuideRate(true); for (size_t i = 0; i < children_.size(); ++i) { // Apply for all children. @@ -580,11 +581,11 @@ namespace Opm const double children_guide_rate = children_[i]->injectionGuideRate(true); #ifdef DIRTY_WELLCTRL_HACK children_[i]->applyInjGroupControl(InjectionSpecification::RESV, - (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, + (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_, false); #else children_[i]->applyInjGroupControl(InjectionSpecification::RATE, - (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, + (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_, false); #endif } @@ -600,15 +601,15 @@ namespace Opm if (phaseUsage().phase_used[BlackoilPhases::Vapour]) { total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour); } - - const double my_guide_rate = injectionGuideRate(true); + const double total_reinjected = - total_produced; // Production negative, injection positive + const double my_guide_rate = injectionGuideRate(true); for (size_t i = 0; i < children_.size(); ++i) { // Apply for all children. // Note, we do _not_ want to call the applyProdGroupControl in this object, // as that would check if we're under group control, something we're not. const double children_guide_rate = children_[i]->injectionGuideRate(true); children_[i]->applyInjGroupControl(InjectionSpecification::RESV, - (children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_, + (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().voidage_replacment_fraction_, false); } From 89eee7e22039836c0c6a540228fdab9cfe87ef41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 10 Oct 2012 14:09:09 +0200 Subject: [PATCH 53/59] Bugfix: order of function arguments. Order of arguments for computePhaseFlowRatesPerWell() was wrong. This fix was done previously for SimulatorCompressibleTwophase, but the incompressible sim was ignored. Also added an ASSERT that may help catch some misuse. --- opm/core/utility/miscUtilities.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 53e997731..73e10ac90 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -553,6 +553,7 @@ namespace Opm { const int np = wells.number_of_phases; const int nw = wells.number_of_wells; + ASSERT(int(flow_rates_per_well_cell.size()) == wells.well_connpos[nw]); phase_flow_per_well.resize(nw * np); for (int wix = 0; wix < nw; ++wix) { for (int phase = 0; phase < np; ++phase) { From ca15ce6eec2d5215ac4f9a0e442a8c7a1add0a02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 10 Oct 2012 14:12:38 +0200 Subject: [PATCH 54/59] Fix comment. --- opm/core/fluid/IncompPropertiesInterface.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/fluid/IncompPropertiesInterface.hpp b/opm/core/fluid/IncompPropertiesInterface.hpp index 4f8bfb516..d3aaaa7b0 100644 --- a/opm/core/fluid/IncompPropertiesInterface.hpp +++ b/opm/core/fluid/IncompPropertiesInterface.hpp @@ -62,7 +62,7 @@ namespace Opm /// \return Array of P viscosity values. virtual const double* viscosity() const = 0; - /// Densities of fluid phases at surface conditions. + /// Densities of fluid phases at reservoir conditions. /// \return Array of P density values. virtual const double* density() const = 0; From 49af5dcd7b6d705ea66d8a709d0339c5707c46a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Oct 2012 11:11:33 +0200 Subject: [PATCH 55/59] Use new velocity interpolation interface. The class TransportModelTracerTofDiscGal now uses VelocityInterpolationInterface, and acts as a factory internally, choosing an interpolation method depending on the parameter 'use_cvi'. --- .../TransportModelTracerTofDiscGal.cpp | 42 ++++++------------- .../TransportModelTracerTofDiscGal.hpp | 13 +++++- 2 files changed, 25 insertions(+), 30 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 2d58ec918..718f2cf6b 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -20,10 +20,11 @@ #include #include #include +#include #include #include -#include #include +#include namespace Opm { @@ -403,32 +404,6 @@ namespace Opm - // Initial version: only a constant interpolation. - static void interpolateVelocity(const UnstructuredGrid& grid, - const int cell, - const double* darcyflux, - const double* /*x*/, - double* v) - { - const int dim = grid.dimensions; - std::fill(v, v + dim, 0.0); - const double* cc = grid.cell_centroids + cell*dim; - for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell+1]; ++hface) { - const int face = grid.cell_faces[hface]; - const double* fc = grid.face_centroids + face*dim; - double flux = 0.0; - if (cell == grid.face_cells[2*face]) { - flux = darcyflux[face]; - } else { - flux = -darcyflux[face]; - } - for (int dd = 0; dd < dim; ++dd) { - v[dd] += flux * (fc[dd] - cc[dd]) / grid.cell_volumes[cell]; - } - } - } - - // --------------- Methods of TransportModelTracerTofDiscGal --------------- @@ -436,11 +411,19 @@ namespace Opm /// Construct solver. /// \param[in] grid A 2d or 3d grid. - TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid) + /// \param[in] use_cvi If true, use corner point velocity interpolation. + /// Otherwise, use the basic constant interpolation. + TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, + const bool use_cvi) : grid_(grid), coord_(grid.dimensions), velocity_(grid.dimensions) { + if (use_cvi) { + velocity_interpolation_.reset(new VelocityInterpolationECVI(grid)); + } else { + velocity_interpolation_.reset(new VelocityInterpolationConstant(grid)); + } } @@ -485,6 +468,7 @@ namespace Opm basis_.resize(num_basis); basis_nb_.resize(num_basis); grad_basis_.resize(num_basis*grid_.dimensions); + velocity_interpolation_->setupFluxes(darcyflux); reorderAndTransport(grid_, darcyflux); } @@ -560,7 +544,7 @@ namespace Opm quad.quadPtCoord(quad_pt, &coord_[0]); DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]); - interpolateVelocity(grid_, cell, darcyflux_, &coord_[0], &velocity_[0]); + velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 0285e1f02..68307e1c6 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -21,15 +21,18 @@ #define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED #include +#include #include #include #include + struct UnstructuredGrid; namespace Opm { class IncompPropertiesInterface; + class VelocityInterpolationInterface; /// Implements a discontinuous Galerkin solver for /// (single-phase) time-of-flight using reordering. @@ -45,7 +48,10 @@ namespace Opm public: /// Construct solver. /// \param[in] grid A 2d or 3d grid. - TransportModelTracerTofDiscGal(const UnstructuredGrid& grid); + /// \param[in] use_cvi If true, use corner point velocity interpolation. + /// Otherwise, use the basic constant interpolation. + TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, + const bool use_cvi); /// Solve for time-of-flight. @@ -72,7 +78,12 @@ namespace Opm virtual void solveMultiCell(const int num_cells, const int* cells); private: + // Disable copying and assignment. + TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&); + TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&); + const UnstructuredGrid& grid_; + boost::shared_ptr velocity_interpolation_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell const double* source_; // one volumetric source term per cell From 88af4a4ce3e8db10228ad40b1951adfc3e37de02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 17 Oct 2012 12:40:43 +0200 Subject: [PATCH 56/59] Fix output in case of LAPACK error. Make copy of matrix before calling dgesv, since it will overwrite it. --- opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 4 +++- opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 718f2cf6b..404452e0e 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -465,6 +465,7 @@ namespace Opm tof_coeff_ = &tof_coeff[0]; rhs_.resize(num_basis); jac_.resize(num_basis*num_basis); + orig_jac_.resize(num_basis*num_basis); basis_.resize(num_basis); basis_nb_.resize(num_basis); grad_basis_.resize(num_basis*grid_.dimensions); @@ -617,6 +618,7 @@ namespace Opm std::vector piv(num_basis); MAT_SIZE_T ldb = num_basis; MAT_SIZE_T info = 0; + orig_jac_ = jac_; dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info); if (info != 0) { // Print the local matrix and rhs. @@ -624,7 +626,7 @@ namespace Opm << " with A = \n"; for (int row = 0; row < n; ++row) { for (int col = 0; col < n; ++col) { - std::cerr << " " << jac_[row + n*col]; + std::cerr << " " << orig_jac_[row + n*col]; } std::cerr << '\n'; } diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 68307e1c6..9a2e3b9c2 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -91,6 +91,7 @@ namespace Opm double* tof_coeff_; std::vector rhs_; // single-cell right-hand-side std::vector jac_; // single-cell jacobian + std::vector orig_jac_; // single-cell jacobian (copy) // Below: storage for quantities needed by solveSingleCell(). std::vector coord_; std::vector basis_; From d8dd982408002ebf6b82461c1dfac92e075df6e0 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Wed, 24 Oct 2012 09:30:20 +0200 Subject: [PATCH 57/59] Make GCC 4.6.3 happy in C++0x mode It complains about not finding a match for the pair<> template class, because the first parameter (this) is allegedly const. However, this isn't a const method, so I suspect it is a compiler bug. In order to move on, I slap on a harmless cast which will make this particular compiler happy, and which should have no effects elsewhere, but put it in a #if..#else..#endif macro to avoid warnings on others; hopefully this also makes it easier to spot and remove in the future. --- opm/core/wells/WellsGroup.cpp | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index d407c0335..e1bed0cb0 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -761,16 +761,28 @@ namespace Opm wells_->ctrls[self_index_]->current = ~ wells_->ctrls[self_index_]->current; } } - + +// macro to insert const_cast to get a round bug in GCC 4.6.3 where it +// suddenly believes that "this" is a const pointer, although we are not +// in a const method. +#if ( __GNUC__ == 4 ) && ( __GNUC_MINOR__ == 6 ) && ( __GNUC_PATCHLEVEL__ == 3 ) +#define CONST_CAST(T,v) const_cast(v) +#else +#define CONST_CAST(T,v) v +#endif + std::pair WellNode::getWorstOffending(const std::vector& well_reservoirrates_phase, const std::vector& well_surfacerates_phase, ProductionSpecification::ControlMode mode) { const int np = phaseUsage().num_phases; const int index = self_index_*np; - return std::make_pair(this, rateByMode(&well_reservoirrates_phase[index], - &well_surfacerates_phase[index], - mode)); + // note: CONST_CAST is just to work around a bug in GCC 4.6.3; it + // is not really needed, and should be a harmless no-op. + return std::make_pair(CONST_CAST(WellNode*,this), + rateByMode(&well_reservoirrates_phase[index], + &well_surfacerates_phase[index], + mode)); } void WellNode::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode, From 0c599b88685262be6a8b2ee16e8db681cc0d5f6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 24 Oct 2012 21:12:29 +0200 Subject: [PATCH 58/59] Fix well classification that was only correct by accident Specifically, the tests if (!wells->type[self_index] == INJECTOR) if (!wells->type[self_index] == PRODUCER) produced the expected results *only* because INJECTOR==0 and PRODUCER==1 in the WellType enumeration, thus (!INJECTOR == PRODUCER) and (!PRODUCER == INJECTOR). Installing the (much) more appropriate if (wells->type[self_index] != INJECTOR) if (wells->type[self_index] != PRODUCER) is not only more readable, it is also future-proof and scales better if we ever introduce new WellTypes (e.g., a MONITOR). --- opm/core/wells/WellsGroup.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index e1bed0cb0..d34d36982 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -794,7 +794,7 @@ namespace Opm && (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) { return; } - if (!wells_->type[self_index_] == INJECTOR) { + if (wells_->type[self_index_] != INJECTOR) { ASSERT(target == 0.0); return; } @@ -871,7 +871,7 @@ namespace Opm std::cout << "Returning" << std::endl; return; } - if (!wells_->type[self_index_] == PRODUCER) { + if (wells_->type[self_index_] != PRODUCER) { ASSERT(target == 0.0); return; } From 28661b73427e55568e421735b2cdc02b1abfa5b3 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Wed, 24 Oct 2012 21:14:24 +0200 Subject: [PATCH 59/59] Remove superfluous construction by std::make_pair Since we know the type of the components, we may just as well create the pair directly! (Using make_pair invokes compiler bugs in GCC). --- opm/core/wells/WellsGroup.cpp | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index e1bed0cb0..8c5414313 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -762,27 +762,16 @@ namespace Opm } } -// macro to insert const_cast to get a round bug in GCC 4.6.3 where it -// suddenly believes that "this" is a const pointer, although we are not -// in a const method. -#if ( __GNUC__ == 4 ) && ( __GNUC_MINOR__ == 6 ) && ( __GNUC_PATCHLEVEL__ == 3 ) -#define CONST_CAST(T,v) const_cast(v) -#else -#define CONST_CAST(T,v) v -#endif - std::pair WellNode::getWorstOffending(const std::vector& well_reservoirrates_phase, const std::vector& well_surfacerates_phase, ProductionSpecification::ControlMode mode) { const int np = phaseUsage().num_phases; const int index = self_index_*np; - // note: CONST_CAST is just to work around a bug in GCC 4.6.3; it - // is not really needed, and should be a harmless no-op. - return std::make_pair(CONST_CAST(WellNode*,this), - rateByMode(&well_reservoirrates_phase[index], - &well_surfacerates_phase[index], - mode)); + return std::pair(this, + rateByMode(&well_reservoirrates_phase[index], + &well_surfacerates_phase[index], + mode)); } void WellNode::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,