From 13f945f081a21eb375b42e455dad473f93f5060c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 17 Jan 2013 14:35:50 +0100 Subject: [PATCH 01/14] Added test for DG basis classes. Only 2d so far. --- tests/test_dgbasis.cpp | 160 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 160 insertions(+) create mode 100644 tests/test_dgbasis.cpp diff --git a/tests/test_dgbasis.cpp b/tests/test_dgbasis.cpp new file mode 100644 index 000000000..d5eb8c8c7 --- /dev/null +++ b/tests/test_dgbasis.cpp @@ -0,0 +1,160 @@ +/* + Copyright 2013 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 // to suppress our messages when throwing + +#define BOOST_TEST_MODULE DGBasisTest +#include + +#include +#include +#include +#include + +using namespace Opm; + + +namespace +{ + + + bool aequal(double a, double b) + { + const double eps = 1e-15; + return std::fabs(a - b) < eps; + } + +} // anonymous namespace + + +namespace cart2d +{ + + static void test() + { + // Set up 2d 1-cell cartesian case. + GridManager g(1, 1); + const UnstructuredGrid& grid = *g.c_grid(); + + // Test DGBasisBoundedTotalDegree, degree 0. + { + DGBasisBoundedTotalDegree b(grid, 0); + BOOST_CHECK_EQUAL(b.numBasisFunc(), 1); + std::vector bx(b.numBasisFunc(), 0.0); + b.eval(0, grid.cell_centroids, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + double x[2] = { 0.123, 0.456 }; + b.eval(0, x, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + std::vector c(b.numBasisFunc(), 0.0); + b.addConstant(0.789, &c[0]); + BOOST_CHECK(aequal(c[0], 0.789)); + b.multiplyGradient(1.234, &c[0]); + BOOST_CHECK(aequal(c[0], 0.789)); + } + // Test DGBasisBoundedTotalDegree, degree 1. + { + DGBasisBoundedTotalDegree b(grid, 1); + BOOST_CHECK_EQUAL(b.numBasisFunc(), 3); + std::vector bx(b.numBasisFunc(), 0.0); + b.eval(0, grid.cell_centroids, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + BOOST_CHECK(aequal(bx[1], 0.0)); + BOOST_CHECK(aequal(bx[2], 0.0)); + double x[2] = { 0.123, 0.456 }; + b.eval(0, x, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + BOOST_CHECK(aequal(bx[1], 0.123 - 0.5)); + BOOST_CHECK(aequal(bx[2], 0.456 - 0.5)); + std::vector c(b.numBasisFunc(), 0.0); + c[0] = 1.0; c[1] = 2.0; c[2] = 3.0; + b.addConstant(0.789, &c[0]); + BOOST_CHECK(aequal(c[0], 1.789)); + BOOST_CHECK(aequal(c[1], 2.0)); + BOOST_CHECK(aequal(c[2], 3.0)); + const double fx = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2]; + b.multiplyGradient(1.234, &c[0]); + BOOST_CHECK(aequal(c[0], 1.789)); + BOOST_CHECK(aequal(c[1], 2.0*1.234)); + BOOST_CHECK(aequal(c[2], 3.0*1.234)); + const double fx2 = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2]; + BOOST_CHECK(aequal(fx2 - c[0], 1.234*(fx - c[0]))); + } + // Test DGBasisMultilin, degree 0. + { + DGBasisMultilin b(grid, 0); + BOOST_CHECK_EQUAL(b.numBasisFunc(), 1); + std::vector bx(b.numBasisFunc(), 0.0); + b.eval(0, grid.cell_centroids, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + double x[2] = { 0.123, 0.456 }; + b.eval(0, x, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + std::vector c(b.numBasisFunc(), 0.0); + b.addConstant(0.789, &c[0]); + BOOST_CHECK(aequal(c[0], 0.789)); + b.multiplyGradient(1.234, &c[0]); + BOOST_CHECK(aequal(c[0], 0.789)); + } + // Test DGBasisMultilin, degree 1. + { + DGBasisMultilin b(grid, 1); + BOOST_CHECK_EQUAL(b.numBasisFunc(), 4); + std::vector bx(b.numBasisFunc(), 0.0); + b.eval(0, grid.cell_centroids, &bx[0]); + BOOST_CHECK(aequal(bx[0], 0.25)); + BOOST_CHECK(aequal(bx[1], 0.25)); + BOOST_CHECK(aequal(bx[2], 0.25)); + BOOST_CHECK(aequal(bx[3], 0.25)); + double x[2] = { 0.123, 0.456 }; + b.eval(0, x, &bx[0]); + const double xm[2] = { 1.0 - x[0], x[0] }; + const double ym[2] = { 1.0 - x[1], x[1] }; + BOOST_CHECK(aequal(bx[0], xm[0]*ym[0])); + BOOST_CHECK(aequal(bx[1], xm[0]*ym[1])); + BOOST_CHECK(aequal(bx[2], xm[1]*ym[0])); + BOOST_CHECK(aequal(bx[3], xm[1]*ym[1])); + std::vector c(b.numBasisFunc(), 0.0); + c[0] = -1.567; c[1] = 1.42; c[2] = 0.59; c[3] = 3.225; + std::vector corig = c; + b.addConstant(0.789, &c[0]); + BOOST_CHECK(aequal(c[0], corig[0] + 0.25*0.789)); + BOOST_CHECK(aequal(c[1], corig[1] + 0.25*0.789)); + BOOST_CHECK(aequal(c[2], corig[2] + 0.25*0.789)); + BOOST_CHECK(aequal(c[3], corig[3] + 0.25*0.789)); + const double fx = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2] + c[3]*bx[3]; + const double fc = 0.25*(c[0] + c[1] + c[2] + c[3]); + b.multiplyGradient(1.234, &c[0]); + const double fx2 = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2] + c[3]*bx[3]; + BOOST_CHECK(aequal(fx2 - fc, 1.234*(fx - fc))); + } + + } +} // namespace cart2d + + +BOOST_AUTO_TEST_CASE(test_dgbasis) +{ + cart2d::test(); +} From d0fa32011bb1ddb6cd1e00ad1acd1147b7acd689 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 18 Jan 2013 13:23:37 +0100 Subject: [PATCH 02/14] Opm::WellState::init(): Handle shut wells This is a minor bugfix to account for the presence of shut wells (characterised by "ctrls[w]->current < 0"). The existing code would lead to indexing outside the "ctrls" array in the context of a shut well. --- opm/core/simulator/WellState.hpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 981ccd9ba..f0de92b07 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -42,12 +42,15 @@ namespace Opm // to pressure in first perforation cell. for (int w = 0; w < nw; ++w) { const WellControls* ctrl = wells->ctrls[w]; - if (ctrl->type[ctrl->current] == BHP) { - bhp_[w] = ctrl->target[ctrl->current]; - } else { + + if ((ctrl->current < 0) || // SHUT + (ctrl->type[ctrl->current] != BHP)) { const int cell = wells->well_cells[wells->well_connpos[w]]; bhp_[w] = state.pressure()[cell]; } + else { + bhp_[w] = ctrl->target[ctrl->current]; + } } perfrates_.resize(wells->well_connpos[nw], 0.0); perfpress_.resize(wells->well_connpos[nw], -1e100); From 763d14ce4d43b764c9f8a6d38bbdc0c620c1b30f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Jan 2013 10:46:06 +0100 Subject: [PATCH 03/14] Refactored limiter functions. The two limiter functions were very similar, and were unified in a single method, with a bool argument to choose method variety. --- .../TransportModelTracerTofDiscGal.cpp | 140 +++--------------- .../TransportModelTracerTofDiscGal.hpp | 3 +- 2 files changed, 24 insertions(+), 119 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 132041a4f..59acb1de5 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -411,10 +411,10 @@ namespace Opm { switch (limiter_method_) { case MinUpwindFace: - applyMinUpwindFaceLimiter(cell, tof); + applyMinUpwindLimiter(cell, true, tof); break; case MinUpwindAverage: - applyMinUpwindAverageLimiter(cell, tof); + applyMinUpwindLimiter(cell, false, tof); break; default: THROW("Limiter type not implemented: " << limiter_method_); @@ -424,18 +424,24 @@ namespace Opm - void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof) + void TransportModelTracerTofDiscGal::applyMinUpwindLimiter(const int cell, const bool face_min, double* tof) { if (basis_func_->degree() != 1) { THROW("This limiter only makes sense for our DG1 implementation."); } // Limiter principles: - // 1. Let M be the minimum TOF value on the upstream faces, - // evaluated in the upstream cells. Then the value at all - // points in this cell shall be at least M. - // Upstream faces whose flux does not exceed the relative - // flux threshold are not considered for this minimum. + // 1. Let M be either: + // - the minimum TOF value of all upstream faces, + // evaluated in the upstream cells + // (chosen if face_min is true). + // or: + // - the minimum average TOF value of all upstream cells + // (chosen if face_min is false). + // Then the value at all points in this cell shall be at + // least M. Upstream faces whose flux does not exceed the + // relative flux threshold are not considered for this + // minimum. // 2. The TOF shall not be below zero in any point. // Find total upstream/downstream fluxes. @@ -492,11 +498,15 @@ namespace Opm min_here_tof = std::min(min_here_tof, tof_here); if (upstream) { if (interior) { - basis_func_->eval(upstream_cell, nc, &basis_nb_[0]); - const double tof_upstream - = std::inner_product(basis_nb_.begin(), basis_nb_.end(), - tof_coeff_ + num_basis*upstream_cell, 0.0); - min_upstream_tof = std::min(min_upstream_tof, tof_upstream); + if (face_min) { + basis_func_->eval(upstream_cell, nc, &basis_nb_[0]); + const double tof_upstream + = std::inner_product(basis_nb_.begin(), basis_nb_.end(), + tof_coeff_ + num_basis*upstream_cell, 0.0); + min_upstream_tof = std::min(min_upstream_tof, tof_upstream); + } else { + min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]); + } } else { // Allow tof down to 0 on inflow boundaries. min_upstream_tof = std::min(min_upstream_tof, 0.0); @@ -535,110 +545,6 @@ namespace Opm - void TransportModelTracerTofDiscGal::applyMinUpwindAverageLimiter(const int cell, double* tof) - { - if (basis_func_->degree() != 1) { - THROW("This limiter only makes sense for our DG1 implementation."); - } - - // Limiter principles: - // 1. Let M be the average TOF value of the upstream cells. - /// Then the value at all points in this cell shall be at least M. - // Upstream faces whose flux does not exceed the relative - // flux threshold are not considered for this minimum. - // 2. The TOF shall not be below zero in any point. - - // Find total upstream/downstream fluxes. - double upstream_flux = 0.0; - double downstream_flux = 0.0; - 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) { - upstream_flux += flux; - } else { - downstream_flux += flux; - } - } - // In the presence of sources, significant fluxes may be missing from the computed fluxes, - // setting the total flux to the (positive) maximum avoids this: since source is either - // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. - const double total_flux = std::max(-upstream_flux, downstream_flux); - - // Find minimum tof on upstream faces and for this cell. - const int dim = grid_.dimensions; - const int num_basis = basis_func_->numBasisFunc(); - double min_upstream_tof = 1e100; - double min_here_tof = 1e100; - int num_upstream_faces = 0; - 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]; - } - const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_); - if (upstream) { - ++num_upstream_faces; - } - bool interior = (upstream_cell >= 0); - - // Evaluate the solution in all corners. - for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { - const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; - basis_func_->eval(cell, nc, &basis_[0]); - const double tof_here = std::inner_product(basis_.begin(), basis_.end(), - tof_coeff_ + num_basis*cell, 0.0); - min_here_tof = std::min(min_here_tof, tof_here); - if (upstream) { - if (interior) { - min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]); - } else { - // Allow tof down to 0 on inflow boundaries. - min_upstream_tof = std::min(min_upstream_tof, 0.0); - } - } - } - } - - // Compute slope multiplier (limiter). - if (num_upstream_faces == 0) { - min_upstream_tof = 0.0; - min_here_tof = 0.0; - } - if (min_upstream_tof < 0.0) { - min_upstream_tof = 0.0; - } - const double tof_c = tof_coeff_[num_basis*cell]; - double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof); - if (tof_c < min_upstream_tof) { - // Handle by setting a flat solution. - std::cout << "Trouble in cell " << cell << std::endl; - limiter = 0.0; - basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell); - } - ASSERT(limiter >= 0.0); - - // Actually do the limiting (if applicable). - if (limiter < 1.0) { - // std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl; - basis_func_->multiplyGradient(limiter, tof + num_basis*cell); - } else { - // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; - } - } - - void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess() diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index bcee06dfc..9ce8cec22 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -130,8 +130,7 @@ namespace Opm // (will read data from tof_coeff_, it is ok to call // with tof_coeff as tof argument. void applyLimiter(const int cell, double* tof); - void applyMinUpwindFaceLimiter(const int cell, double* tof); - void applyMinUpwindAverageLimiter(const int cell, double* tof); + void applyMinUpwindLimiter(const int cell, const bool face_min, double* tof); void applyLimiterAsPostProcess(); void applyLimiterAsSimultaneousPostProcess(); }; From 1093e09d650cab7771e297b5f5fb4cf5b1a1e4c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Jan 2013 11:09:07 +0100 Subject: [PATCH 04/14] Added functionAverage() method to DGBasis classes. --- opm/core/transport/reorder/DGBasis.cpp | 17 +++++++++++++++-- opm/core/transport/reorder/DGBasis.hpp | 12 ++++++++++++ 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/DGBasis.cpp b/opm/core/transport/reorder/DGBasis.cpp index 8527a1caa..93ca347bd 100644 --- a/opm/core/transport/reorder/DGBasis.cpp +++ b/opm/core/transport/reorder/DGBasis.cpp @@ -154,6 +154,12 @@ namespace Opm } } + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + double DGBasisBoundedTotalDegree::functionAverage(const double* coefficients) const + { + return coefficients[0]; + } @@ -300,11 +306,18 @@ namespace Opm double* coefficients) const { const int nb = numBasisFunc(); - const double average = std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb); + const double aver = functionAverage(coefficients); for (int ix = 0; ix < nb; ++ix) { - coefficients[ix] = factor*(coefficients[ix] - average) + average; + coefficients[ix] = factor*(coefficients[ix] - aver) + aver; } } + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + double DGBasisMultilin::functionAverage(const double* coefficients) const + { + const int nb = numBasisFunc(); + return std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb); + } } // namespace Opm diff --git a/opm/core/transport/reorder/DGBasis.hpp b/opm/core/transport/reorder/DGBasis.hpp index a74108148..596fefbd8 100644 --- a/opm/core/transport/reorder/DGBasis.hpp +++ b/opm/core/transport/reorder/DGBasis.hpp @@ -75,6 +75,10 @@ namespace Opm /// \param[out] coefficients Coefficients {c_i} for a single cell. virtual void multiplyGradient(const double factor, double* coefficients) const = 0; + + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + virtual double functionAverage(const double* coefficients) const = 0; }; @@ -144,6 +148,10 @@ namespace Opm virtual void multiplyGradient(const double factor, double* coefficients) const; + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + virtual double functionAverage(const double* coefficients) const; + private: const UnstructuredGrid& grid_; const int degree_; @@ -217,6 +225,10 @@ namespace Opm virtual void multiplyGradient(const double factor, double* coefficients) const; + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + virtual double functionAverage(const double* coefficients) const; + private: const UnstructuredGrid& grid_; const int degree_; From e02bf4f312b77c586d84e8d94f7acf15263b08cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Jan 2013 11:09:43 +0100 Subject: [PATCH 05/14] Use functionAverage() instead of direct coefficient access. This fixes the issue with limiters on multilinear basis. --- .../transport/reorder/TransportModelTracerTofDiscGal.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 59acb1de5..fdbec4409 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -498,14 +498,16 @@ namespace Opm min_here_tof = std::min(min_here_tof, tof_here); if (upstream) { if (interior) { + const double* upstream_coef = tof_coeff_ + num_basis*upstream_cell; if (face_min) { basis_func_->eval(upstream_cell, nc, &basis_nb_[0]); const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(), - tof_coeff_ + num_basis*upstream_cell, 0.0); + upstream_coef, 0.0); min_upstream_tof = std::min(min_upstream_tof, tof_upstream); } else { - min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]); + min_upstream_tof = std::min(min_upstream_tof, + basis_func_->functionAverage(upstream_coef)); } } else { // Allow tof down to 0 on inflow boundaries. @@ -523,7 +525,7 @@ namespace Opm if (min_upstream_tof < 0.0) { min_upstream_tof = 0.0; } - const double tof_c = tof_coeff_[num_basis*cell]; + const double tof_c = basis_func_->functionAverage(tof_coeff_ + num_basis*cell); double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof); if (tof_c < min_upstream_tof) { // Handle by setting a flat solution. From 604be4587133cd240d561621263b4033aefc5717 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Jan 2013 14:55:27 +0100 Subject: [PATCH 06/14] Add methods totalFlux() and minCornerVal(). Also started refactoring applyMinUpwindLimiter() using the added methods. --- .../TransportModelTracerTofDiscGal.cpp | 77 ++++++++++++------- .../TransportModelTracerTofDiscGal.hpp | 6 +- 2 files changed, 54 insertions(+), 29 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index fdbec4409..24846cf58 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -444,34 +444,13 @@ namespace Opm // minimum. // 2. The TOF shall not be below zero in any point. - // Find total upstream/downstream fluxes. - double upstream_flux = 0.0; - double downstream_flux = 0.0; - 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) { - upstream_flux += flux; - } else { - downstream_flux += flux; - } - } - // In the presence of sources, significant fluxes may be missing from the computed fluxes, - // setting the total flux to the (positive) maximum avoids this: since source is either - // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. - const double total_flux = std::max(-upstream_flux, downstream_flux); - - // Find minimum tof on upstream faces and for this cell. + // Find minimum tof on upstream faces/cells and for this cell. const int dim = grid_.dimensions; const int num_basis = basis_func_->numBasisFunc(); double min_upstream_tof = 1e100; double min_here_tof = 1e100; int num_upstream_faces = 0; + const double total_flux = totalFlux(cell); 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; @@ -490,12 +469,9 @@ namespace Opm bool interior = (upstream_cell >= 0); // Evaluate the solution in all corners. + min_here_tof = std::min(min_here_tof, minCornerVal(cell, face)); for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; - basis_func_->eval(cell, nc, &basis_[0]); - const double tof_here = std::inner_product(basis_.begin(), basis_.end(), - tof_coeff_ + num_basis*cell, 0.0); - min_here_tof = std::min(min_here_tof, tof_here); if (upstream) { if (interior) { const double* upstream_coef = tof_coeff_ + num_basis*upstream_cell; @@ -583,4 +559,51 @@ namespace Opm + double TransportModelTracerTofDiscGal::totalFlux(const int cell) const + { + // Find total upstream/downstream fluxes. + double upstream_flux = 0.0; + double downstream_flux = 0.0; + 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) { + upstream_flux += flux; + } else { + downstream_flux += flux; + } + } + // In the presence of sources, significant fluxes may be missing from the computed fluxes, + // setting the total flux to the (positive) maximum avoids this: since source is either + // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. + return std::max(-upstream_flux, downstream_flux); + } + + + + + double TransportModelTracerTofDiscGal::minCornerVal(const int cell, const int face) const + { + // Evaluate the solution in all corners. + const int dim = grid_.dimensions; + const int num_basis = basis_func_->numBasisFunc(); + double min_cornerval = 1e100; + for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { + const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; + basis_func_->eval(cell, nc, &basis_[0]); + const double tof_corner = std::inner_product(basis_.begin(), basis_.end(), + tof_coeff_ + num_basis*cell, 0.0); + min_cornerval = std::min(min_cornerval, tof_corner); + } + return min_cornerval; + } + + + + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 9ce8cec22..cf55ca8cb 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -119,8 +119,8 @@ namespace Opm std::vector orig_jac_; // single-cell jacobian (copy) // Below: storage for quantities needed by solveSingleCell(). std::vector coord_; - std::vector basis_; - std::vector basis_nb_; + mutable std::vector basis_; + mutable std::vector basis_nb_; std::vector grad_basis_; std::vector velocity_; @@ -133,6 +133,8 @@ namespace Opm void applyMinUpwindLimiter(const int cell, const bool face_min, double* tof); void applyLimiterAsPostProcess(); void applyLimiterAsSimultaneousPostProcess(); + double totalFlux(const int cell) const; + double minCornerVal(const int cell, const int face) const; }; } // namespace Opm From f99fd9254b7c69ea5cd0aecc6d428cadb590b90e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 23 Jan 2013 09:50:25 +0100 Subject: [PATCH 07/14] Refactored applyMinUpwindLimiter(). --- .../TransportModelTracerTofDiscGal.cpp | 34 ++++++------------- 1 file changed, 11 insertions(+), 23 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 24846cf58..54a6e5858 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -445,7 +445,6 @@ namespace Opm // 2. The TOF shall not be below zero in any point. // Find minimum tof on upstream faces/cells and for this cell. - const int dim = grid_.dimensions; const int num_basis = basis_func_->numBasisFunc(); double min_upstream_tof = 1e100; double min_here_tof = 1e100; @@ -463,33 +462,22 @@ namespace Opm upstream_cell = grid_.face_cells[2*face]; } const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_); + const bool interior = (upstream_cell >= 0); + + // Find minimum tof in this cell and upstream. + // The meaning of minimum upstream tof depends on method. + min_here_tof = std::min(min_here_tof, minCornerVal(cell, face)); if (upstream) { ++num_upstream_faces; - } - bool interior = (upstream_cell >= 0); - - // Evaluate the solution in all corners. - min_here_tof = std::min(min_here_tof, minCornerVal(cell, face)); - for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { - const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; - if (upstream) { - if (interior) { - const double* upstream_coef = tof_coeff_ + num_basis*upstream_cell; - if (face_min) { - basis_func_->eval(upstream_cell, nc, &basis_nb_[0]); - const double tof_upstream - = std::inner_product(basis_nb_.begin(), basis_nb_.end(), - upstream_coef, 0.0); - min_upstream_tof = std::min(min_upstream_tof, tof_upstream); - } else { - min_upstream_tof = std::min(min_upstream_tof, - basis_func_->functionAverage(upstream_coef)); - } + double upstream_tof = 0.0; + if (interior) { + if (face_min) { + upstream_tof = minCornerVal(upstream_cell, face); } else { - // Allow tof down to 0 on inflow boundaries. - min_upstream_tof = std::min(min_upstream_tof, 0.0); + upstream_tof = basis_func_->functionAverage(tof_coeff_ + num_basis*upstream_cell); } } + min_upstream_tof = std::min(min_upstream_tof, upstream_tof); } } From 977e8a19e2cf4150c20bdad6fc6076eb187814d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 23 Jan 2013 09:51:30 +0100 Subject: [PATCH 08/14] Added evalFunc() method. This public method is not virtual, and implemented in the base class using calls to the virtual methods. Not yet used by the DG solver. --- opm/core/transport/reorder/DGBasis.cpp | 18 ++++++++++++++++++ opm/core/transport/reorder/DGBasis.hpp | 15 +++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/opm/core/transport/reorder/DGBasis.cpp b/opm/core/transport/reorder/DGBasis.cpp index 93ca347bd..8ffb95ec4 100644 --- a/opm/core/transport/reorder/DGBasis.cpp +++ b/opm/core/transport/reorder/DGBasis.cpp @@ -25,12 +25,30 @@ namespace Opm { + // ---------------- Methods for class DGBasisInterface ---------------- + /// Virtual destructor. DGBasisInterface::~DGBasisInterface() { } + /// Evaluate function f = sum_i c_i b_i at the point x. + /// Note that this function is not virtual, but implemented in + /// terms of the virtual functions of the class. + /// \param[in] cell Cell index + /// \param[in] coefficients Coefficients {c_i} for a single cell. + /// \param[in] x Point at which to compute f(x). + double DGBasisInterface::evalFunc(const int cell, + const double* coefficients, + const double* x) const + { + bvals_.resize(numBasisFunc()); + eval(cell, x, &bvals_[0]); + return std::inner_product(bvals_.begin(), bvals_.end(), coefficients, 0.0); + } + + // ---------------- Methods for class DGBasisBoundedTotalDegree ---------------- diff --git a/opm/core/transport/reorder/DGBasis.hpp b/opm/core/transport/reorder/DGBasis.hpp index 596fefbd8..f1fc72341 100644 --- a/opm/core/transport/reorder/DGBasis.hpp +++ b/opm/core/transport/reorder/DGBasis.hpp @@ -20,6 +20,8 @@ #ifndef OPM_DGBASIS_HEADER_INCLUDED #define OPM_DGBASIS_HEADER_INCLUDED +#include + struct UnstructuredGrid; namespace Opm @@ -76,9 +78,22 @@ namespace Opm virtual void multiplyGradient(const double factor, double* coefficients) const = 0; + /// Evaluate function f = sum_i c_i b_i at the point x. + /// Note that this function is not virtual, but implemented in + /// terms of the virtual functions of the class. + /// \param[in] cell Cell index + /// \param[in] coefficients Coefficients {c_i} for a single cell. + /// \param[in] x Point at which to compute f(x). + double evalFunc(const int cell, + const double* coefficients, + const double* x) const; + /// Compute the average of the function f = sum_i c_i b_i. /// \param[in] coefficients Coefficients {c_i} for a single cell. virtual double functionAverage(const double* coefficients) const = 0; + + private: + mutable std::vector bvals_; // For evalFunc(). }; From 95fb0792e57dd9b878f09dd7c2b9ef839f8e0d7d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 23 Jan 2013 17:35:05 +0100 Subject: [PATCH 09/14] Bring ParameterGroup interface into scope. The constructor accepts a parameter::ParameterGroup reference and thus needs a valid interface in scope. Relying on header pollution is unwise. --- opm/core/fluid/BlackoilPropertiesBasic.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/core/fluid/BlackoilPropertiesBasic.hpp b/opm/core/fluid/BlackoilPropertiesBasic.hpp index 44a76013b..cc34a3cfe 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.hpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.hpp @@ -25,6 +25,7 @@ #include #include #include +#include namespace Opm { From 7dd4720057bc2a499284bf215315987658ba29bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 29 Jan 2013 13:42:24 +0100 Subject: [PATCH 10/14] Changed OpenRS->OPM in copyright notices and #include guards. --- opm/core/wells/WellCollection.cpp | 27 +++++++++++++-------------- opm/core/wells/WellCollection.hpp | 26 +++++++++++++------------- 2 files changed, 26 insertions(+), 27 deletions(-) diff --git a/opm/core/wells/WellCollection.cpp b/opm/core/wells/WellCollection.cpp index 727a7242b..f3920be84 100644 --- a/opm/core/wells/WellCollection.cpp +++ b/opm/core/wells/WellCollection.cpp @@ -1,22 +1,21 @@ /* -Copyright 2011 SINTEF ICT, Applied Mathematics. + Copyright 2011 SINTEF ICT, Applied Mathematics. + This file is part of the Open Porous Media project (OPM). -This file is part of The Open Reservoir Simulator Project (OpenRS). + 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. -OpenRS 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. -OpenRS 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 OpenRS. If not, see . - */ + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ #include diff --git a/opm/core/wells/WellCollection.hpp b/opm/core/wells/WellCollection.hpp index 1957576a6..4bda76d9b 100644 --- a/opm/core/wells/WellCollection.hpp +++ b/opm/core/wells/WellCollection.hpp @@ -1,22 +1,22 @@ /* -Copyright 2011 SINTEF ICT, Applied Mathematics. + Copyright 2011 SINTEF ICT, Applied Mathematics. + This file is part of the Open Porous Media project (OPM). -This file is part of The Open Reservoir Simulator Project (OpenRS). + 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. -OpenRS 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. -OpenRS 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 . +*/ -You should have received a copy of the GNU General Public License -along with OpenRS. If not, see . - */ #ifndef OPM_WELLCOLLECTION_HPP From 62a8a7f527caf837365dc5b72432378f8f980dd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 1 Feb 2013 16:06:39 +0100 Subject: [PATCH 11/14] Use portable method of zeroing vector of ints. The memset() technique is only applicable to platforms for which numerical zero is represented by all bits zero. --- opm/core/transport/reorder/tarjan.c | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/opm/core/transport/reorder/tarjan.c b/opm/core/transport/reorder/tarjan.c index 9150c2092..49d803847 100644 --- a/opm/core/transport/reorder/tarjan.c +++ b/opm/core/transport/reorder/tarjan.c @@ -20,8 +20,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ -#include #include +#include #ifdef MATLAB_MEX_FILE #include "tarjan.h" @@ -30,7 +30,13 @@ SOFTWARE. #endif +static void +clear_vector(size_t n, int *v) +{ + size_t i; + for (i = 0; i < n; i++) { v[i] = 0; } +} static int min(int a, int b){ return a < b? a : b;} @@ -81,11 +87,9 @@ tarjan (int nv, const int *ia, const int *ja, int *vert, int *comp, int *link = (int *) time + nv; int *status = (int *) link + nv; /* dual usage... */ - (void) cbottom; - - memset(work, 0, 3*nv * sizeof *work); - memset(vert, 0, nv * sizeof *vert ); - memset(comp, 0, (nv+1) * sizeof *comp ); + clear_vector(3 * ((size_t) nv), work); + clear_vector(1 * ((size_t) nv), vert); + clear_vector(1 + ((size_t) nv), comp); /* Init status all vertices */ for (i=0; i Date: Fri, 1 Feb 2013 16:15:35 +0100 Subject: [PATCH 12/14] Eliminate release-mode build warning. The 'cbottom' variable is only used within an assert(). Don't define the variable in release (i.e., "NDEBUG") mode. --- opm/core/transport/reorder/tarjan.c | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/core/transport/reorder/tarjan.c b/opm/core/transport/reorder/tarjan.c index 49d803847..8d02c4b54 100644 --- a/opm/core/transport/reorder/tarjan.c +++ b/opm/core/transport/reorder/tarjan.c @@ -79,7 +79,11 @@ tarjan (int nv, const int *ia, const int *ja, int *vert, int *comp, int *stack = comp + nv; int *bottom = stack; int *cstack = vert + nv-1; + +#if !defined(NDEBUG) int *cbottom = cstack; +#endif + int t = 0; int pos = 0; From 267ec64ae8610e792e0c961eeda9df59d2f20a3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 1 Feb 2013 16:25:46 +0100 Subject: [PATCH 13/14] Eliminate redundant explicit type conversion. The pointers in question are already type 'int *'. There is no need to explicitly convert them to that type too. --- opm/core/transport/reorder/tarjan.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/tarjan.c b/opm/core/transport/reorder/tarjan.c index 8d02c4b54..320244533 100644 --- a/opm/core/transport/reorder/tarjan.c +++ b/opm/core/transport/reorder/tarjan.c @@ -88,8 +88,8 @@ tarjan (int nv, const int *ia, const int *ja, int *vert, int *comp, int pos = 0; int *time = work; - int *link = (int *) time + nv; - int *status = (int *) link + nv; /* dual usage... */ + int *link = time + nv; + int *status = link + nv; /* dual usage... */ clear_vector(3 * ((size_t) nv), work); clear_vector(1 * ((size_t) nv), vert); From 967403996a881ac435312e17c2dbce31fe99755b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 1 Feb 2013 16:30:45 +0100 Subject: [PATCH 14/14] Replace an assignment with intended equality test This corrects a latent error that has been present since the inception of this module. --- opm/core/transport/reorder/tarjan.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/tarjan.c b/opm/core/transport/reorder/tarjan.c index 320244533..d3245610b 100644 --- a/opm/core/transport/reorder/tarjan.c +++ b/opm/core/transport/reorder/tarjan.c @@ -196,7 +196,7 @@ tarjan (int nv, const int *ia, const int *ja, int *vert, int *comp, } else { - assert(status[child] = DONE); + assert(status[child] == DONE); } } }