From b25f81cf5800dbb32729e2b2efe545acaa72e653 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 13:10:09 +0100 Subject: [PATCH 01/13] Do not assume that the first basis function is the constant 1. Properly integrate b_i * porosity over the cell for all basis functions b_i. --- .../reorder/TransportModelTracerTofDiscGal.cpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index ae6ab28d..a6810e67 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -271,8 +271,18 @@ namespace Opm 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]; + { + CellQuadrature quad(grid_, cell, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // Integral of: b_i \phi + 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) { + rhs_[j] += w * basis_[j] * porevolume_[cell] / grid_.cell_volumes[cell]; + } + } + } // Compute upstream residual contribution. for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { From d670f3713143bcbf01960b36edc9be9fccf47677 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 13:17:43 +0100 Subject: [PATCH 02/13] Added (disabled) extra output for debugging purposes. Enable by defining EXTRA_VERBOSE. --- .../TransportModelTracerTofDiscGal.cpp | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index a6810e67..05755ab1 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -434,6 +434,29 @@ namespace Opm // Apply limiter. if (degree_ > 0 && use_limiter_ && limiter_usage_ == DuringComputations) { +#ifdef EXTRA_VERBOSE + std::cout << "Cell: " << cell << " "; + std::cout << "v = "; + for (int dd = 0; dd < dim; ++dd) { + std::cout << velocity_[dd] << ' '; + } + std::cout << " grad tau = "; + for (int dd = 0; dd < dim; ++dd) { + std::cout << tof_coeff_[num_basis*cell + dd + 1] << ' '; + } + const double prod = std::inner_product(velocity_.begin(), velocity_.end(), + tof_coeff_ + num_basis*cell + 1, 0.0); + const double vv = std::inner_product(velocity_.begin(), velocity_.end(), + velocity_.begin(), 0.0); + const double gg = std::inner_product(tof_coeff_ + num_basis*cell + 1, + tof_coeff_ + num_basis*cell + num_basis, + tof_coeff_ + num_basis*cell + 1, 0.0); + std::cout << " prod = " << std::inner_product(velocity_.begin(), velocity_.end(), + tof_coeff_ + num_basis*cell + 1, 0.0); + std::cout << " normalized = " << prod/std::sqrt(vv*gg); + std::cout << " angle = " << std::acos(prod/std::sqrt(vv*gg))*360.0/(2.0*M_PI); + std::cout << std::endl; +#endif applyLimiter(cell, tof_coeff_); } } From 6e716db0f91a27f4faec62237ce3f5b5f8fee7c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 13:33:55 +0100 Subject: [PATCH 03/13] Do not require unnecessary high quadrature precision. --- opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 05755ab1..6f56f0cb 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -329,7 +329,9 @@ namespace Opm // Compute cell jacobian contribution. We use Fortran ordering // for jac_, i.e. rows cycling fastest. { - const int deg_needed = use_cvi_ ? 2*degree_ : 2*degree_ - 1; + // Even with ECVI velocity interpolation, degree of precision 1 + // is sufficient for optimal convergence order for DG1. + const int deg_needed = 2*degree_ - 1; CellQuadrature quad(grid_, cell, deg_needed); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // b_i (v \cdot \grad b_j) From a2f2fcdfbb42ba3cf9f830ce9cf434397ebafd73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 13:37:39 +0100 Subject: [PATCH 04/13] Renamed class DGBasis -> DGBasisBoundedTotalDegree. --- .../transport/reorder/TransportModelTracerTofDiscGal.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 6f56f0cb..2bb690b1 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -37,8 +37,9 @@ namespace Opm - /// A class providing discontinuous Galerkin basis functions. - struct DGBasis + /// A class providing discontinuous Galerkin basis functions + /// of bounded total degree. + struct DGBasisBoundedTotalDegree { static int numBasisFunc(const int dimensions, const int degree) @@ -117,6 +118,8 @@ namespace Opm + typedef DGBasisBoundedTotalDegree DGBasis; + // typedef DGBasisMultilin DGBasis; From 4000eeaef695b38ac545b942d8296c50bf6375d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 13:39:35 +0100 Subject: [PATCH 05/13] Reduce source sum failure from error to warning. Should extend to sum over sources and boundaries before reinstating. --- 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 2bb690b1..0f1474eb 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -221,7 +221,8 @@ namespace Opm // 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); + // THROW("Sources do not sum to zero: " << cum_src); + MESSAGE("Warning: sources do not sum to zero: " << cum_src); } #endif degree_ = degree; From b5a2556df772bcee140f47224a538e1ccd3fb658 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 13:42:05 +0100 Subject: [PATCH 06/13] Added class DGBasisMultilin. Multilinear DG1 basis functions added. Will not work with current limiter. --- .../TransportModelTracerTofDiscGal.cpp | 98 +++++++++++++++++++ 1 file changed, 98 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 0f1474eb..19a64cb3 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -113,6 +113,101 @@ namespace Opm + /// A class providing discontinuous Galerkin basis functions of + /// multi-degree 1 (bilinear or trilinear functions). + struct DGBasisMultilin + { + static int numBasisFunc(const int dimensions, + const int degree) + { + switch (dimensions) { + case 1: + return degree + 1; + case 2: + return (degree + 1)*(degree + 1); + case 3: + return (degree + 1)*(degree + 1)*(degree + 1); + 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. + /// (for 2 dims:) + /// (Bi)degree 1: (x-)(y-), (x-)(y+), (x+)(y-), (x+)(y+) + /// where (x-) = (1/2 - x + xc), (x+) = (1/2 + x - xc) + /// and similar for (y-), (y+). + static void eval(const UnstructuredGrid& grid, + const int cell, + const int degree, + const double* x, + double* f_x) + { + const int dim = grid.dimensions; + const int num_basis = numBasisFunc(dim, degree); + const double* cc = grid.cell_centroids + dim*cell; + switch (degree) { + case 0: + f_x[0] = 1; + break; + case 1: + std::fill(f_x, f_x + num_basis, 1.0); + for (int dd = 0; dd < dim; ++dd) { + const double f[2] = { 0.5 - x[dd] + cc[dd], 0.5 + x[dd] - cc[dd] }; + const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example. + for (int ix = 0; ix < num_basis; ++ix) { + f_x[ix] *= f[(ix/divi) % 2]; + } + } + 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); + const double* cc = grid.cell_centroids + dim*cell; + switch (degree) { + case 0: + std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0); + break; + case 1: + std::fill(grad_f_x, grad_f_x + num_basis*dim, 1.0); + for (int dd = 0; dd < dim; ++dd) { + const double f[2] = { 0.5 - x[dd] + cc[dd], 0.5 + x[dd] - cc[dd] }; + const double fder[2] = { -1.0, 1.0 }; + const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example. + for (int ix = 0; ix < num_basis; ++ix) { + const int ind = (ix/divi) % 2; + for (int dder = 0; dder < dim; ++dder) { + grad_f_x[ix*dim + dder] *= (dder == dd ? fder[ind] : f[ind]); + } + } + } + break; + default: + THROW("Maximum degree is 1 for now."); + } + } + }; + @@ -123,6 +218,9 @@ namespace Opm + + + // --------------- Methods of TransportModelTracerTofDiscGal --------------- From ae2016e0630dec6b31317ed47ad24492ddf4dd0e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 13:44:18 +0100 Subject: [PATCH 07/13] Added addConstant() and multiplyGradient() methods. So far only to class DGBasisBoundedTotalDegree. --- .../TransportModelTracerTofDiscGal.cpp | 36 +++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 19a64cb3..44585673 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -108,6 +108,42 @@ namespace Opm } } } + + /// Modify basis coefficients to add to the function value. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to (f + increment) by modifying the c_i. This is done without + /// modifying its gradient. + /// \param[in] increment Add this value to the function. + /// \param[in] dim Number of space dimensions. + /// \param[in] degree Polynomial degree of basis. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + static void addConstant(const double increment, + const int /*dim*/, + const int /*degree*/, + double* coefficients) + { + coefficients[0] += increment; + } + + /// Modify basis coefficients to change the function's slope. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to a function g with the property that grad g = factor * grad f + /// by modifying the c_i. This is done without modifying the average, + /// i.e. the integrals of g and f over the cell are the same. + /// \param[in] factor Multiply gradient by this factor. + /// \param[in] dim Number of space dimensions. + /// \param[in] degree Polynomial degree of basis. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + static void multiplyGradient(const double factor, + const int dim, + const int degree, + double* coefficients) + { + const int nb = numBasisFunc(dim, degree); + for (int ix = 1; ix < nb; ++ix) { + coefficients[ix] *= factor; + } + } }; From d822b92e1cd412c7413f17a27125ddc593a04e2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 13:52:44 +0100 Subject: [PATCH 08/13] Implement limiters with addConstant() and multiplyGradient(). This is instead of directly manipulating the coefficients, requiring assumptions on the basis used. --- .../reorder/TransportModelTracerTofDiscGal.cpp | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 44585673..3025308b 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -727,16 +727,14 @@ namespace Opm // Handle by setting a flat solution. std::cout << "Trouble in cell " << cell << std::endl; limiter = 0.0; - tof[num_basis*cell] = min_upstream_tof; + DGBasis::addConstant(min_upstream_tof - tof_c, dim, degree_, 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; - for (int i = num_basis*cell + 1; i < num_basis*(cell+1); ++i) { - tof[i] *= limiter; - } + DGBasis::multiplyGradient(limiter, dim, degree_, tof + num_basis*cell); } else { // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; } @@ -835,16 +833,14 @@ namespace Opm // Handle by setting a flat solution. std::cout << "Trouble in cell " << cell << std::endl; limiter = 0.0; - tof[num_basis*cell] = min_upstream_tof; + DGBasis::addConstant(min_upstream_tof - tof_c, dim, degree_, 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; - for (int i = num_basis*cell + 1; i < num_basis*(cell+1); ++i) { - tof[i] *= limiter; - } + DGBasis::multiplyGradient(limiter, dim, degree_, tof + num_basis*cell); } else { // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; } From 52c164040ab3de632f5618f3662a4d016bd1bce2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2013 20:44:28 +0100 Subject: [PATCH 09/13] Add methods addConstant() and multiplyGradient() to multilinear basis. Not tested yet, though. --- .../TransportModelTracerTofDiscGal.cpp | 42 +++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 3025308b..6118849f 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -242,6 +242,48 @@ namespace Opm THROW("Maximum degree is 1 for now."); } } + + /// Modify basis coefficients to add to the function value. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to (f + increment) by modifying the c_i. This is done without + /// modifying its gradient. + /// \param[in] increment Add this value to the function. + /// \param[in] dim Number of space dimensions. + /// \param[in] degree Polynomial degree of basis. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + static void addConstant(const double increment, + const int dim, + const int degree, + double* coefficients) + { + const int nb = numBasisFunc(dim, degree); + const double term = increment/double(nb); + for (int ix = 0; ix < nb; ++ix) { + coefficients[ix] += term; + } + } + + /// Modify basis coefficients to change the function's slope. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to a function g with the property that grad g = factor * grad f + /// by modifying the c_i. This is done without modifying the average, + /// i.e. the integrals of g and f over the cell are the same. + /// \param[in] factor Multiply gradient by this factor. + /// \param[in] dim Number of space dimensions. + /// \param[in] degree Polynomial degree of basis. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + static void multiplyGradient(const double factor, + const int dim, + const int degree, + double* coefficients) + { + const int nb = numBasisFunc(dim, degree); + const double average = std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb); + for (int ix = 0; ix < nb; ++ix) { + coefficients[ix] = factor*(coefficients[ix] - average) + average; + } + } + }; From 33cc611ced33d937fe8f42d9b34c20e38b7646f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 Jan 2013 13:38:02 +0100 Subject: [PATCH 10/13] Moved DG basis to separate classes in its own file. Also introduced interface (base) class and changed the api, not used yet by the tof solver. --- Makefile.am | 2 + opm/core/transport/reorder/DGBasis.cpp | 304 +++++++++++++++++++++++++ opm/core/transport/reorder/DGBasis.hpp | 231 +++++++++++++++++++ 3 files changed, 537 insertions(+) create mode 100644 opm/core/transport/reorder/DGBasis.cpp create mode 100644 opm/core/transport/reorder/DGBasis.hpp diff --git a/Makefile.am b/Makefile.am index 556d2bb9..f6d0b5ea 100644 --- a/Makefile.am +++ b/Makefile.am @@ -108,6 +108,7 @@ opm/core/simulator/SimulatorCompressibleTwophase.cpp \ opm/core/simulator/SimulatorIncompTwophase.cpp \ opm/core/simulator/SimulatorReport.cpp \ opm/core/simulator/SimulatorTimer.cpp \ +opm/core/transport/reorder/DGBasis.cpp \ opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \ opm/core/transport/reorder/TransportModelInterface.cpp \ opm/core/transport/reorder/TransportModelTracerTof.cpp \ @@ -237,6 +238,7 @@ opm/core/transport/JacobianSystem.hpp \ opm/core/transport/NormSupport.hpp \ opm/core/transport/SimpleFluid2pWrapper.hpp \ opm/core/transport/SinglePointUpwindTwoPhase.hpp \ +opm/core/transport/reorder/DGBasis.hpp \ opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp \ opm/core/transport/reorder/TransportModelInterface.hpp \ opm/core/transport/reorder/TransportModelTracerTof.hpp \ diff --git a/opm/core/transport/reorder/DGBasis.cpp b/opm/core/transport/reorder/DGBasis.cpp new file mode 100644 index 00000000..2462f635 --- /dev/null +++ b/opm/core/transport/reorder/DGBasis.cpp @@ -0,0 +1,304 @@ +/* + 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 +#include +#include +#include + +namespace Opm +{ + + + // ---------------- Methods for class DGBasisBoundedTotalDegree ---------------- + + + /// Constructor. + /// \param[in] grid grid on which basis is used (cell-wise) + /// \param[in] degree polynomial degree of basis + DGBasisBoundedTotalDegree::DGBasisBoundedTotalDegree(const UnstructuredGrid& grid, + const int degree_arg) + : grid_(grid), + degree_(degree_arg) + { + if (grid_.dimensions > 3) { + THROW("Grid dimension must be 1, 2 or 3."); + } + if (degree_ > 1 || degree_ < 0) { + THROW("Degree must be 0 or 1."); + } + } + + /// Destructor. + DGBasisBoundedTotalDegree::~DGBasisBoundedTotalDegree() + { + } + + /// The number of basis functions per cell. + int DGBasisBoundedTotalDegree::numBasisFunc() const + { + 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."); + } + } + + /// The number of space dimensions. + int DGBasisBoundedTotalDegree::dimensions() const + { + return grid_.dimensions; + } + + /// The polynomial degree of the basis functions. + int DGBasisBoundedTotalDegree::degree() const + { + return degree_; + } + + /// Evaluate all basis functions associated with cell at x, + /// writing to f_x. The array f_x must have size equal to + /// numBasisFunc(). + void DGBasisBoundedTotalDegree::eval(const int cell, + const double* x, + double* f_x) const + { + const int dim = 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 basis functions associated with + /// cell at x, writing to grad_f_x. The array grad_f_x must + /// have size numBasisFunc() * dimensions(). The dimensions() + /// components of the first basis function gradient come + /// before the components of the second etc. + void DGBasisBoundedTotalDegree::evalGrad(const int /*cell*/, + const double* /*x*/, + double* grad_f_x) const + { + const int dim = dimensions(); + const int num_basis = numBasisFunc(); + std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0); + if (degree_ == 1) { + for (int ix = 0; ix < dim; ++ix) { + grad_f_x[dim*(ix + 1) + ix] = 1.0; + } + } + } + + /// Modify basis coefficients to add to the function value. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to (f + increment) by modifying the c_i. This is done without + /// modifying its gradient. + /// \param[in] increment Add this value to the function. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + void DGBasisBoundedTotalDegree::addConstant(const double increment, + double* coefficients) const + { + coefficients[0] += increment; + } + + /// Modify basis coefficients to change the function's slope. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to a function g with the property that grad g = factor * grad f + /// by modifying the c_i. This is done without modifying the average, + /// i.e. the integrals of g and f over the cell are the same. + /// \param[in] factor Multiply gradient by this factor. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + void DGBasisBoundedTotalDegree::multiplyGradient(const double factor, + double* coefficients) const + { + const int nb = numBasisFunc(); + for (int ix = 1; ix < nb; ++ix) { + coefficients[ix] *= factor; + } + } + + + + + // ---------------- Methods for class DGBasisMultilin ---------------- + + + /// Constructor. + /// \param[in] grid grid on which basis is used (cell-wise) + /// \param[in] degree polynomial degree of basis + DGBasisMultilin::DGBasisMultilin(const UnstructuredGrid& grid, + const int degree_arg) + : grid_(grid), + degree_(degree_arg) + { + if (grid_.dimensions > 3) { + THROW("Grid dimension must be 1, 2 or 3."); + } + if (degree_ > 1 || degree_ < 0) { + THROW("Degree must be 0 or 1."); + } + } + + /// Destructor. + DGBasisMultilin::~DGBasisMultilin() + { + } + + /// The number of basis functions per cell. + int DGBasisMultilin::numBasisFunc() const + { + switch (dimensions()) { + case 1: + return degree_ + 1; + case 2: + return (degree_ + 1)*(degree_ + 1); + case 3: + return (degree_ + 1)*(degree_ + 1)*(degree_ + 1); + default: + THROW("Dimensions must be 1, 2 or 3."); + } + } + + /// The number of space dimensions. + int DGBasisMultilin::dimensions() const + { + return grid_.dimensions; + } + + /// The polynomial degree of the basis functions. + int DGBasisMultilin::degree() const + { + return degree_; + } + + /// Evaluate all basis functions associated with cell at x, + /// writing to f_x. The array f_x must have size equal to + /// numBasisFunc(). + void DGBasisMultilin::eval(const int cell, + const double* x, + double* f_x) const + { + const int dim = dimensions(); + const int num_basis = numBasisFunc(); + const double* cc = grid_.cell_centroids + dim*cell; + switch (degree_) { + case 0: + f_x[0] = 1; + break; + case 1: + std::fill(f_x, f_x + num_basis, 1.0); + for (int dd = 0; dd < dim; ++dd) { + const double f[2] = { 0.5 - x[dd] + cc[dd], 0.5 + x[dd] - cc[dd] }; + const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example. + for (int ix = 0; ix < num_basis; ++ix) { + f_x[ix] *= f[(ix/divi) % 2]; + } + } + break; + default: + THROW("Maximum degree is 1 for now."); + } + } + + + /// Evaluate gradients of all basis functions associated with + /// cell at x, writing to grad_f_x. The array grad_f_x must + /// have size numBasisFunc() * dimensions(). The dimensions() + /// components of the first basis function gradient come + /// before the components of the second etc. + void DGBasisMultilin::evalGrad(const int cell, + const double* x, + double* grad_f_x) const + { + const int dim = dimensions(); + const int num_basis = numBasisFunc(); + const double* cc = grid_.cell_centroids + dim*cell; + switch (degree_) { + case 0: + std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0); + break; + case 1: + std::fill(grad_f_x, grad_f_x + num_basis*dim, 1.0); + for (int dd = 0; dd < dim; ++dd) { + const double f[2] = { 0.5 - x[dd] + cc[dd], 0.5 + x[dd] - cc[dd] }; + const double fder[2] = { -1.0, 1.0 }; + const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example. + for (int ix = 0; ix < num_basis; ++ix) { + const int ind = (ix/divi) % 2; + for (int dder = 0; dder < dim; ++dder) { + grad_f_x[ix*dim + dder] *= (dder == dd ? fder[ind] : f[ind]); + } + } + } + break; + default: + THROW("Maximum degree is 1 for now."); + } + } + + /// Modify basis coefficients to add to the function value. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to (f + increment) by modifying the c_i. This is done without + /// modifying its gradient. + /// \param[in] increment Add this value to the function. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + void DGBasisMultilin::addConstant(const double increment, + double* coefficients) const + { + const int nb = numBasisFunc(); + const double term = increment/double(nb); + for (int ix = 0; ix < nb; ++ix) { + coefficients[ix] += term; + } + } + + /// Modify basis coefficients to change the function's slope. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to a function g with the property that grad g = factor * grad f + /// by modifying the c_i. This is done without modifying the average, + /// i.e. the integrals of g and f over the cell are the same. + /// \param[in] factor Multiply gradient by this factor. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + void DGBasisMultilin::multiplyGradient(const double factor, + double* coefficients) const + { + const int nb = numBasisFunc(); + const double average = std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb); + for (int ix = 0; ix < nb; ++ix) { + coefficients[ix] = factor*(coefficients[ix] - average) + average; + } + } + + +} // namespace Opm diff --git a/opm/core/transport/reorder/DGBasis.hpp b/opm/core/transport/reorder/DGBasis.hpp new file mode 100644 index 00000000..4219b6c9 --- /dev/null +++ b/opm/core/transport/reorder/DGBasis.hpp @@ -0,0 +1,231 @@ +/* + 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 . +*/ + +#ifndef OPM_DGBASIS_HEADER_INCLUDED +#define OPM_DGBASIS_HEADER_INCLUDED + +struct UnstructuredGrid; + +namespace Opm +{ + + /// Base class for Discontinuous Galerkin bases, intended for time-of-flight computations. + class DGBasisInterface + { + public: + /// Virtual destructor. + virtual ~DGBasisInterface() = 0; + + /// The number of basis functions per cell. + virtual int numBasisFunc() const = 0; + + /// The number of space dimensions. + virtual int dimensions() const = 0; + + /// The polynomial degree of the basis functions. + virtual int degree() const = 0; + + /// Evaluate all basis functions associated with cell at x, + /// writing to f_x. The array f_x must have size equal to + /// numBasisFunc(). + virtual void eval(const int cell, + const double* x, + double* f_x) const = 0; + + /// Evaluate gradients of all basis functions associated with + /// cell at x, writing to grad_f_x. The array grad_f_x must + /// have size numBasisFunc() * dimensions(). The dimensions() + /// components of the first basis function gradient come + /// before the components of the second etc. + virtual void evalGrad(const int cell, + const double* x, + double* grad_f_x) const = 0; + + /// Modify basis coefficients to add to the function value. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to (f + increment) by modifying the c_i. This is done without + /// modifying its gradient. + /// \param[in] increment Add this value to the function. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + virtual void addConstant(const double increment, + double* coefficients) const = 0; + + /// Modify basis coefficients to change the function's slope. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to a function g with the property that grad g = factor * grad f + /// by modifying the c_i. This is done without modifying the average, + /// i.e. the integrals of g and f over the cell are the same. + /// \param[in] factor Multiply gradient by this factor. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + virtual void multiplyGradient(const double factor, + double* coefficients) const = 0; + }; + + + + + + /// A class providing discontinuous Galerkin basis functions + /// of bounded total degree. + /// + /// The basis functions are the following for each cell (example for 3d): + /// Degree 0: 1. + /// Degree 1: 1, x - xc, y - yc, z - zc + /// where (xc, yc, zc) are the coordinates of the cell centroid. + /// Further degrees await development. + class DGBasisBoundedTotalDegree : public DGBasisInterface + { + public: + /// Constructor. + /// \param[in] grid grid on which basis is used (cell-wise) + /// \param[in] degree polynomial degree of basis + DGBasisBoundedTotalDegree(const UnstructuredGrid& grid, const int degree); + + /// Destructor. + virtual ~DGBasisBoundedTotalDegree(); + + /// The number of basis functions per cell. + virtual int numBasisFunc() const; + + /// The number of space dimensions. + virtual int dimensions() const; + + /// The polynomial degree of the basis functions. + virtual int degree() const; + + /// Evaluate all basis functions associated with cell at x, + /// writing to f_x. The array f_x must have size equal to + /// numBasisFunc(). + virtual void eval(const int cell, + const double* x, + double* f_x) const; + + /// Evaluate gradients of all basis functions associated with + /// cell at x, writing to grad_f_x. The array grad_f_x must + /// have size numBasisFunc() * dimensions(). The dimensions() + /// components of the first basis function gradient come + /// before the components of the second etc. + virtual void evalGrad(const int cell, + const double* x, + double* grad_f_x) const; + + /// Modify basis coefficients to add to the function value. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to (f + increment) by modifying the c_i. This is done without + /// modifying its gradient. + /// \param[in] increment Add this value to the function. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + virtual void addConstant(const double increment, + double* coefficients) const; + + /// Modify basis coefficients to change the function's slope. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to a function g with the property that grad g = factor * grad f + /// by modifying the c_i. This is done without modifying the average, + /// i.e. the integrals of g and f over the cell are the same. + /// \param[in] factor Multiply gradient by this factor. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + virtual void multiplyGradient(const double factor, + double* coefficients) const; + + private: + const UnstructuredGrid& grid_; + const int degree_; + }; + + + + + /// A class providing discontinuous Galerkin basis functions of + /// multi-degree 1 (bilinear or trilinear functions). + /// + /// The basis functions for a cell are the following + /// Degree 0: 1. + /// (for 2 dims:) + /// (Bi)degree 1: (x-)(y-), (x-)(y+), (x+)(y-), (x+)(y+) + /// where (x-) = (1/2 - x + xc), (x+) = (1/2 + x - xc) + /// and xc is the x-coordinate of the cell centroid. + /// Similar for (y-), (y+). + class DGBasisMultilin : public DGBasisInterface + { + /// Constructor. + /// \param[in] grid grid on which basis is used (cell-wise) + /// \param[in] degree polynomial degree of basis (in each coordinate) + DGBasisMultilin(const UnstructuredGrid& grid, const int degree); + + /// Destructor. + virtual ~DGBasisMultilin(); + + /// The number of basis functions per cell. + virtual int numBasisFunc() const; + + /// The number of space dimensions. + virtual int dimensions() const; + + /// The polynomial degree of the basis functions. + virtual int degree() const; + + /// Evaluate all basis functions associated with cell at x, + /// writing to f_x. The array f_x must have size equal to + /// numBasisFunc(). + virtual void eval(const int cell, + const double* x, + double* f_x) const; + + /// Evaluate gradients of all basis functions associated with + /// cell at x, writing to grad_f_x. The array grad_f_x must + /// have size numBasisFunc() * dimensions(). The dimensions() + /// components of the first basis function gradient come + /// before the components of the second etc. + virtual void evalGrad(const int cell, + const double* x, + double* grad_f_x) const; + + /// Modify basis coefficients to add to the function value. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to (f + increment) by modifying the c_i. This is done without + /// modifying its gradient. + /// \param[in] increment Add this value to the function. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + virtual void addConstant(const double increment, + double* coefficients) const; + + /// Modify basis coefficients to change the function's slope. + /// A function f = sum_i c_i b_i is assumed, and we change + /// it to a function g with the property that grad g = factor * grad f + /// by modifying the c_i. This is done without modifying the average, + /// i.e. the integrals of g and f over the cell are the same. + /// \param[in] factor Multiply gradient by this factor. + /// \param[out] coefficients Coefficients {c_i} for a single cell. + virtual void multiplyGradient(const double factor, + double* coefficients) const; + + private: + const UnstructuredGrid& grid_; + const int degree_; + + }; + + + + +} // namespace Opm + + +#endif // OPM_DGBASIS_HEADER_INCLUDED From 564db5df87c6229066b778344e349d684688891b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 Jan 2013 13:45:15 +0100 Subject: [PATCH 11/13] Increased quadrature order in a term. This is for the benefit of bi/tri-linear basis functions, as stated in the comments. --- .../reorder/TransportModelTracerTofDiscGal.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 6118849f..c1853eae 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -510,8 +510,15 @@ namespace Opm // for jac_, i.e. rows cycling fastest. { // Even with ECVI velocity interpolation, degree of precision 1 - // is sufficient for optimal convergence order for DG1. - const int deg_needed = 2*degree_ - 1; + // is sufficient for optimal convergence order for DG1 when we + // use linear (total degree 1) basis functions. + // With bi(tri)-linear basis functions, it still seems sufficient + // for convergence order 2, but the solution looks much better and + // has significantly lower error with degree of precision 2. + // For now, we err on the side of caution, and use 2*degree, even + // though this is wasteful for the pure linear basis functions. + // const int deg_needed = 2*degree_ - 1; + const int deg_needed = 2*degree_; CellQuadrature quad(grid_, cell, deg_needed); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // b_i (v \cdot \grad b_j) From fb060da933adc5b79e698612e85e92c4c00f0367 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 Jan 2013 15:11:46 +0100 Subject: [PATCH 12/13] Implement (empty) virtual destructor, make methods public. --- opm/core/transport/reorder/DGBasis.cpp | 6 ++++++ opm/core/transport/reorder/DGBasis.hpp | 3 ++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/DGBasis.cpp b/opm/core/transport/reorder/DGBasis.cpp index 2462f635..8527a1ca 100644 --- a/opm/core/transport/reorder/DGBasis.cpp +++ b/opm/core/transport/reorder/DGBasis.cpp @@ -26,6 +26,12 @@ namespace Opm { + /// Virtual destructor. + DGBasisInterface::~DGBasisInterface() + { + } + + // ---------------- Methods for class DGBasisBoundedTotalDegree ---------------- diff --git a/opm/core/transport/reorder/DGBasis.hpp b/opm/core/transport/reorder/DGBasis.hpp index 4219b6c9..a7410814 100644 --- a/opm/core/transport/reorder/DGBasis.hpp +++ b/opm/core/transport/reorder/DGBasis.hpp @@ -30,7 +30,7 @@ namespace Opm { public: /// Virtual destructor. - virtual ~DGBasisInterface() = 0; + virtual ~DGBasisInterface(); /// The number of basis functions per cell. virtual int numBasisFunc() const = 0; @@ -164,6 +164,7 @@ namespace Opm /// Similar for (y-), (y+). class DGBasisMultilin : public DGBasisInterface { + public: /// Constructor. /// \param[in] grid grid on which basis is used (cell-wise) /// \param[in] degree polynomial degree of basis (in each coordinate) From c3f9e64c9c3deb472b2b3d013e90fcf03912703c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 Jan 2013 15:13:45 +0100 Subject: [PATCH 13/13] Complete separation of basis func classes. Also: - Add use_tensorial_basis parameter allowing run-time choice of basis. - Remove degree argument from solveTof() method, degree is instead obtained from parameters in constructors. Modified compute_tof* programs to match. --- examples/compute_tof.cpp | 4 +- examples/compute_tof_from_files.cpp | 4 +- .../TransportModelTracerTofDiscGal.cpp | 335 ++---------------- .../TransportModelTracerTofDiscGal.hpp | 8 +- 4 files changed, 46 insertions(+), 305 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 1dfc6ec5..9189befb 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -168,12 +168,10 @@ main(int argc, char** argv) // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); - int dg_degree = -1; bool use_multidim_upwind = false; // Need to initialize dg solver here, since it uses parameters now. boost::scoped_ptr dg_solver; if (use_dg) { - dg_degree = param.getDefault("dg_degree", 0); dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(*grid->c_grid(), param)); } else { use_multidim_upwind = param.getDefault("use_multidim_upwind", false); @@ -237,7 +235,7 @@ main(int argc, char** argv) std::vector tof; std::vector tracer; if (use_dg) { - dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof); + dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); } else { Opm::TransportModelTracerTof tofsolver(*grid->c_grid(), use_multidim_upwind); if (compute_tracer) { diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index bf5c077d..b8421396 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -122,12 +122,10 @@ main(int argc, char** argv) // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); - int dg_degree = -1; bool use_multidim_upwind = false; // Need to initialize dg solver here, since it uses parameters now. boost::scoped_ptr dg_solver; if (use_dg) { - dg_degree = param.getDefault("dg_degree", 0); dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(grid, param)); } else { use_multidim_upwind = param.getDefault("use_multidim_upwind", false); @@ -163,7 +161,7 @@ main(int argc, char** argv) std::vector tof; std::vector tracer; if (use_dg) { - dg_solver->solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof); + dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof); } else { Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind); if (compute_tracer) { diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index c1853eae..132041a4 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 @@ -33,272 +34,6 @@ namespace Opm { - // --------------- Helpers for TransportModelTracerTofDiscGal --------------- - - - - /// A class providing discontinuous Galerkin basis functions - /// of bounded total degree. - struct DGBasisBoundedTotalDegree - { - 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; - } - } - } - - /// Modify basis coefficients to add to the function value. - /// A function f = sum_i c_i b_i is assumed, and we change - /// it to (f + increment) by modifying the c_i. This is done without - /// modifying its gradient. - /// \param[in] increment Add this value to the function. - /// \param[in] dim Number of space dimensions. - /// \param[in] degree Polynomial degree of basis. - /// \param[out] coefficients Coefficients {c_i} for a single cell. - static void addConstant(const double increment, - const int /*dim*/, - const int /*degree*/, - double* coefficients) - { - coefficients[0] += increment; - } - - /// Modify basis coefficients to change the function's slope. - /// A function f = sum_i c_i b_i is assumed, and we change - /// it to a function g with the property that grad g = factor * grad f - /// by modifying the c_i. This is done without modifying the average, - /// i.e. the integrals of g and f over the cell are the same. - /// \param[in] factor Multiply gradient by this factor. - /// \param[in] dim Number of space dimensions. - /// \param[in] degree Polynomial degree of basis. - /// \param[out] coefficients Coefficients {c_i} for a single cell. - static void multiplyGradient(const double factor, - const int dim, - const int degree, - double* coefficients) - { - const int nb = numBasisFunc(dim, degree); - for (int ix = 1; ix < nb; ++ix) { - coefficients[ix] *= factor; - } - } - }; - - - - - /// A class providing discontinuous Galerkin basis functions of - /// multi-degree 1 (bilinear or trilinear functions). - struct DGBasisMultilin - { - static int numBasisFunc(const int dimensions, - const int degree) - { - switch (dimensions) { - case 1: - return degree + 1; - case 2: - return (degree + 1)*(degree + 1); - case 3: - return (degree + 1)*(degree + 1)*(degree + 1); - 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. - /// (for 2 dims:) - /// (Bi)degree 1: (x-)(y-), (x-)(y+), (x+)(y-), (x+)(y+) - /// where (x-) = (1/2 - x + xc), (x+) = (1/2 + x - xc) - /// and similar for (y-), (y+). - static void eval(const UnstructuredGrid& grid, - const int cell, - const int degree, - const double* x, - double* f_x) - { - const int dim = grid.dimensions; - const int num_basis = numBasisFunc(dim, degree); - const double* cc = grid.cell_centroids + dim*cell; - switch (degree) { - case 0: - f_x[0] = 1; - break; - case 1: - std::fill(f_x, f_x + num_basis, 1.0); - for (int dd = 0; dd < dim; ++dd) { - const double f[2] = { 0.5 - x[dd] + cc[dd], 0.5 + x[dd] - cc[dd] }; - const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example. - for (int ix = 0; ix < num_basis; ++ix) { - f_x[ix] *= f[(ix/divi) % 2]; - } - } - 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); - const double* cc = grid.cell_centroids + dim*cell; - switch (degree) { - case 0: - std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0); - break; - case 1: - std::fill(grad_f_x, grad_f_x + num_basis*dim, 1.0); - for (int dd = 0; dd < dim; ++dd) { - const double f[2] = { 0.5 - x[dd] + cc[dd], 0.5 + x[dd] - cc[dd] }; - const double fder[2] = { -1.0, 1.0 }; - const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example. - for (int ix = 0; ix < num_basis; ++ix) { - const int ind = (ix/divi) % 2; - for (int dder = 0; dder < dim; ++dder) { - grad_f_x[ix*dim + dder] *= (dder == dd ? fder[ind] : f[ind]); - } - } - } - break; - default: - THROW("Maximum degree is 1 for now."); - } - } - - /// Modify basis coefficients to add to the function value. - /// A function f = sum_i c_i b_i is assumed, and we change - /// it to (f + increment) by modifying the c_i. This is done without - /// modifying its gradient. - /// \param[in] increment Add this value to the function. - /// \param[in] dim Number of space dimensions. - /// \param[in] degree Polynomial degree of basis. - /// \param[out] coefficients Coefficients {c_i} for a single cell. - static void addConstant(const double increment, - const int dim, - const int degree, - double* coefficients) - { - const int nb = numBasisFunc(dim, degree); - const double term = increment/double(nb); - for (int ix = 0; ix < nb; ++ix) { - coefficients[ix] += term; - } - } - - /// Modify basis coefficients to change the function's slope. - /// A function f = sum_i c_i b_i is assumed, and we change - /// it to a function g with the property that grad g = factor * grad f - /// by modifying the c_i. This is done without modifying the average, - /// i.e. the integrals of g and f over the cell are the same. - /// \param[in] factor Multiply gradient by this factor. - /// \param[in] dim Number of space dimensions. - /// \param[in] degree Polynomial degree of basis. - /// \param[out] coefficients Coefficients {c_i} for a single cell. - static void multiplyGradient(const double factor, - const int dim, - const int degree, - double* coefficients) - { - const int nb = numBasisFunc(dim, degree); - const double average = std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb); - for (int ix = 0; ix < nb; ++ix) { - coefficients[ix] = factor*(coefficients[ix] - average) + average; - } - } - - }; - - - - - - - typedef DGBasisBoundedTotalDegree DGBasis; - // typedef DGBasisMultilin DGBasis; - - - - - - // --------------- Methods of TransportModelTracerTofDiscGal --------------- @@ -331,6 +66,14 @@ namespace Opm coord_(grid.dimensions), velocity_(grid.dimensions) { + const int dg_degree = param.getDefault("dg_degree", 0); + const bool use_tensorial_basis = param.getDefault("use_tensorial_basis", false); + if (use_tensorial_basis) { + basis_func_.reset(new DGBasisMultilin(grid_, dg_degree)); + } else { + basis_func_.reset(new DGBasisBoundedTotalDegree(grid_, dg_degree)); + } + use_cvi_ = param.getDefault("use_cvi", use_cvi_); use_limiter_ = param.getDefault("use_limiter", use_limiter_); if (use_limiter_) { @@ -387,7 +130,6 @@ namespace Opm void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux, const double* porevolume, const double* source, - const int degree, std::vector& tof_coeff) { darcyflux_ = darcyflux; @@ -401,8 +143,7 @@ namespace Opm MESSAGE("Warning: sources do not sum to zero: " << cum_src); } #endif - degree_ = degree; - const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_); + const int num_basis = basis_func_->numBasisFunc(); tof_coeff.resize(num_basis*grid_.number_of_cells); std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0); tof_coeff_ = &tof_coeff[0]; @@ -445,18 +186,19 @@ namespace Opm // equal to Res = Jac*c - rhs, and we compute rhs directly. const int dim = grid_.dimensions; - const int num_basis = DGBasis::numBasisFunc(dim, degree_); + const int num_basis = basis_func_->numBasisFunc(); std::fill(rhs_.begin(), rhs_.end(), 0.0); std::fill(jac_.begin(), jac_.end(), 0.0); // Compute cell residual contribution. { - CellQuadrature quad(grid_, cell, degree_); + const int deg_needed = basis_func_->degree(); + CellQuadrature quad(grid_, cell, deg_needed); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // Integral of: b_i \phi quad.quadPtCoord(quad_pt, &coord_[0]); - DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + basis_func_->eval(cell, &coord_[0], &basis_[0]); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { rhs_[j] += w * basis_[j] * porevolume_[cell] / grid_.cell_volumes[cell]; @@ -492,11 +234,12 @@ namespace Opm // velocity is constant (this assumption may have to go // for higher order than DG1). const double normal_velocity = flux / grid_.face_areas[face]; - FaceQuadrature quad(grid_, face, 2*degree_); + const int deg_needed = 2*basis_func_->degree(); + FaceQuadrature quad(grid_, face, deg_needed); 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]); - DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]); + basis_func_->eval(cell, &coord_[0], &basis_[0]); + basis_func_->eval(upstream_cell, &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); @@ -517,14 +260,14 @@ namespace Opm // has significantly lower error with degree of precision 2. // For now, we err on the side of caution, and use 2*degree, even // though this is wasteful for the pure linear basis functions. - // const int deg_needed = 2*degree_ - 1; - const int deg_needed = 2*degree_; + // const int deg_needed = 2*basis_func_->degree() - 1; + const int deg_needed = 2*basis_func_->degree(); CellQuadrature quad(grid_, cell, deg_needed); 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]); + basis_func_->eval(cell, &coord_[0], &basis_[0]); + basis_func_->evalGrad(cell, &coord_[0], &grad_basis_[0]); velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { @@ -553,11 +296,11 @@ namespace Opm // Do quadrature over the face to compute // \int_{\partial K} b_i (v(x) \cdot n) b_j ds const double normal_velocity = flux / grid_.face_areas[face]; - FaceQuadrature quad(grid_, face, 2*degree_); + FaceQuadrature quad(grid_, face, 2*basis_func_->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]); + basis_func_->eval(cell, &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) { @@ -578,10 +321,10 @@ namespace Opm 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, 2*degree_); + CellQuadrature quad(grid_, cell, 2*basis_func_->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]); + basis_func_->eval(cell, &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) { @@ -622,7 +365,7 @@ namespace Opm std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); // Apply limiter. - if (degree_ > 0 && use_limiter_ && limiter_usage_ == DuringComputations) { + if (basis_func_->degree() > 0 && use_limiter_ && limiter_usage_ == DuringComputations) { #ifdef EXTRA_VERBOSE std::cout << "Cell: " << cell << " "; std::cout << "v = "; @@ -683,7 +426,7 @@ namespace Opm void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof) { - if (degree_ != 1) { + if (basis_func_->degree() != 1) { THROW("This limiter only makes sense for our DG1 implementation."); } @@ -719,7 +462,7 @@ namespace Opm // Find minimum tof on upstream faces and for this cell. const int dim = grid_.dimensions; - const int num_basis = DGBasis::numBasisFunc(dim, degree_); + const int num_basis = basis_func_->numBasisFunc(); double min_upstream_tof = 1e100; double min_here_tof = 1e100; int num_upstream_faces = 0; @@ -743,13 +486,13 @@ namespace Opm // 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]; - DGBasis::eval(grid_, cell, degree_, nc, &basis_[0]); + 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) { - DGBasis::eval(grid_, upstream_cell, degree_, nc, &basis_nb_[0]); + 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); @@ -776,14 +519,14 @@ namespace Opm // Handle by setting a flat solution. std::cout << "Trouble in cell " << cell << std::endl; limiter = 0.0; - DGBasis::addConstant(min_upstream_tof - tof_c, dim, degree_, tof + num_basis*cell); + 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; - DGBasis::multiplyGradient(limiter, dim, degree_, tof + num_basis*cell); + basis_func_->multiplyGradient(limiter, tof + num_basis*cell); } else { // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; } @@ -794,7 +537,7 @@ namespace Opm void TransportModelTracerTofDiscGal::applyMinUpwindAverageLimiter(const int cell, double* tof) { - if (degree_ != 1) { + if (basis_func_->degree() != 1) { THROW("This limiter only makes sense for our DG1 implementation."); } @@ -829,7 +572,7 @@ namespace Opm // Find minimum tof on upstream faces and for this cell. const int dim = grid_.dimensions; - const int num_basis = DGBasis::numBasisFunc(dim, degree_); + const int num_basis = basis_func_->numBasisFunc(); double min_upstream_tof = 1e100; double min_here_tof = 1e100; int num_upstream_faces = 0; @@ -853,7 +596,7 @@ namespace Opm // 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]; - DGBasis::eval(grid_, cell, degree_, nc, &basis_[0]); + 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); @@ -882,14 +625,14 @@ namespace Opm // Handle by setting a flat solution. std::cout << "Trouble in cell " << cell << std::endl; limiter = 0.0; - DGBasis::addConstant(min_upstream_tof - tof_c, dim, degree_, tof + num_basis*cell); + 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; - DGBasis::multiplyGradient(limiter, dim, degree_, tof + num_basis*cell); + basis_func_->multiplyGradient(limiter, tof + num_basis*cell); } else { // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; } @@ -921,7 +664,7 @@ namespace Opm // This means that each cell is limited independently from all other cells, // we write the resulting dofs to a new array instead of writing to tof_coeff_. // Afterwards we copy the results back to tof_coeff_. - const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_); + const int num_basis = basis_func_->numBasisFunc(); std::vector tof_coeffs_new(tof_coeff_, tof_coeff_ + num_basis*grid_.number_of_cells); for (int c = 0; c < grid_.number_of_cells; ++c) { applyLimiter(c, &tof_coeffs_new[0]); diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 529f75b1..bcee06df 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -33,6 +33,7 @@ namespace Opm class IncompPropertiesInterface; class VelocityInterpolationInterface; + class DGBasisInterface; namespace parameter { class ParameterGroup; } /// Implements a discontinuous Galerkin solver for @@ -51,6 +52,9 @@ namespace Opm /// \param[in] grid A 2d or 3d grid. /// \param[in] param Parameters for the solver. /// The following parameters are accepted (defaults): + /// dg_degree (0) Polynomial degree of basis functions. + /// use_tensorial_basis (false) Use tensor-product basis, interpreting dg_degree as + /// bi/tri-degree not total degree. /// use_cvi (false) Use ECVI velocity interpolation. /// use_limiter (false) Use a slope limiter. If true, the next three parameters are used. /// limiter_relative_flux_threshold (1e-3) Ignore upstream fluxes below this threshold, relative to total cell flux. @@ -74,7 +78,6 @@ namespace Opm /// \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 /// the K coefficients corresponding to the first @@ -84,7 +87,6 @@ namespace Opm void solveTof(const double* darcyflux, const double* porevolume, const double* source, - const int degree, std::vector& tof_coeff); private: @@ -109,7 +111,7 @@ namespace Opm 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_; + boost::shared_ptr basis_func_; double* tof_coeff_; std::vector rhs_; // single-cell right-hand-side std::vector jac_; // single-cell jacobian