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 01/28] 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 02/28] 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 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 03/28] 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 04/28] 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 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 05/28] 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 06/28] 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 07/28] 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 08/28] 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 09/28] 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 10/28] 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 11/28] 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 12/28] 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 13/28] 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 14/28] 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 15/28] 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 16/28] 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 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 17/28] 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 18/28] 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 19/28] 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 20/28] 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 21/28] 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 22/28] 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 23/28] 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 24/28] 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 25/28] 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 26/28] 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 27/28] 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 28/28] 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;