OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #116 from atgeirr/master

Browse Source 
More DG time-of-flight improvements.
This commit is contained in:
Bård Skaflestad 2013-01-16 06:54:19 -08:00
commit 8b482cadda
7 changed files with 633 additions and 130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile.am
View File

@ -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 \

4
examples/compute_tof.cpp
View File

@ -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<Opm::TransportModelTracerTofDiscGal> 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<double> tof;
std::vector<double> 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) {

4
examples/compute_tof_from_files.cpp
View File

@ -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<Opm::TransportModelTracerTofDiscGal> 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<double> tof;
std::vector<double> 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) {

310
opm/core/transport/reorder/DGBasis.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/DGBasis.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <numeric>
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

232
opm/core/transport/reorder/DGBasis.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#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

203
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
View File

@ -20,6 +20,7 @@
#include <opm/core/grid/CellQuadrature.hpp>
#include <opm/core/grid/FaceQuadrature.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <opm/core/transport/reorder/DGBasis.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/VelocityInterpolation.hpp>
@ -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 <grid.dimensions> 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<double>& 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<double> 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]);

8
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
View File

@ -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<double>& 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<DGBasisInterface> basis_func_;
double* tof_coeff_;
std::vector<double> rhs_; // single-cell right-hand-side
std::vector<double> jac_; // single-cell jacobian