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/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/DGBasis.cpp b/opm/core/transport/reorder/DGBasis.cpp new file mode 100644 index 00000000..8527a1ca --- /dev/null +++ b/opm/core/transport/reorder/DGBasis.cpp @@ -0,0 +1,310 @@ +/* + 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 +{ + + + /// Virtual destructor. + DGBasisInterface::~DGBasisInterface() + { + } + + + // ---------------- 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..a7410814 --- /dev/null +++ b/opm/core/transport/reorder/DGBasis.hpp @@ -0,0 +1,232 @@ +/* + 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(); + + /// 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 + { + public: + /// 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 diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index ae6ab28d..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,93 +34,6 @@ namespace Opm { - // --------------- Helpers for TransportModelTracerTofDiscGal --------------- - - - - /// A class providing discontinuous Galerkin basis functions. - struct DGBasis - { - static int numBasisFunc(const int dimensions, - const int degree) - { - switch (dimensions) { - case 1: - return degree + 1; - case 2: - return (degree + 2)*(degree + 1)/2; - case 3: - return (degree + 3)*(degree + 2)*(degree + 1)/6; - default: - THROW("Dimensions must be 1, 2 or 3."); - } - } - - /// Evaluate all nonzero basis functions at x, - /// writing to f_x. The array f_x must have - /// size numBasisFunc(grid.dimensions, degree). - /// - /// The basis functions are the following - /// Degree 0: 1. - /// Degree 1: x - xc, y - yc, z - zc etc. - /// Further degrees await development. - static void eval(const UnstructuredGrid& grid, - const int cell, - const int degree, - const double* x, - double* f_x) - { - const int dim = grid.dimensions; - const double* cc = grid.cell_centroids + dim*cell; - // Note intentional fallthrough in this switch statement! - switch (degree) { - case 1: - for (int ix = 0; ix < dim; ++ix) { - f_x[1 + ix] = x[ix] - cc[ix]; - } - case 0: - f_x[0] = 1; - break; - default: - THROW("Maximum degree is 1 for now."); - } - } - - /// Evaluate gradients of all nonzero basis functions at x, - /// writing to grad_f_x. The array grad_f_x must have size - /// numBasisFunc(grid.dimensions, degree) * grid.dimensions. - /// The components of the first basis function - /// gradient come before the components of the second etc. - static void evalGrad(const UnstructuredGrid& grid, - const int /*cell*/, - const int degree, - const double* /*x*/, - double* grad_f_x) - { - const int dim = grid.dimensions; - const int num_basis = numBasisFunc(dim, degree); - std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0); - if (degree > 1) { - THROW("Maximum degree is 1 for now."); - } else if (degree == 1) { - for (int ix = 0; ix < dim; ++ix) { - grad_f_x[dim*(ix + 1) + ix] = 1.0; - } - } - } - }; - - - - - - - - - - - - // --------------- Methods of TransportModelTracerTofDiscGal --------------- @@ -152,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_) { @@ -208,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; @@ -218,11 +139,11 @@ 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; - 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]; @@ -265,14 +186,25 @@ 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. - // Note: Assumes that \int_K b_j = 0 for all j > 0 - rhs_[0] += porevolume_[cell]; + { + 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]); + 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]; + } + } + } // Compute upstream residual contribution. for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { @@ -302,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); @@ -319,13 +252,22 @@ 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 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*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) { @@ -354,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) { @@ -379,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) { @@ -423,7 +365,30 @@ 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 = "; + 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_); } } @@ -461,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."); } @@ -497,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; @@ -521,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); @@ -554,16 +519,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; + 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; - for (int i = num_basis*cell + 1; i < num_basis*(cell+1); ++i) { - tof[i] *= limiter; - } + basis_func_->multiplyGradient(limiter, tof + num_basis*cell); } else { // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; } @@ -574,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."); } @@ -609,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; @@ -633,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); @@ -662,16 +625,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; + 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; - for (int i = num_basis*cell + 1; i < num_basis*(cell+1); ++i) { - tof[i] *= limiter; - } + basis_func_->multiplyGradient(limiter, tof + num_basis*cell); } else { // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; } @@ -703,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