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

Moved DG basis to separate classes in its own file.

Browse Source 
Also introduced interface (base) class and changed the api, not used yet by the tof solver.
This commit is contained in:
Atgeirr Flø Rasmussen 2013-01-16 13:38:02 +01:00
parent 52c164040a
commit 33cc611ced
3 changed files with 537 additions and 0 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 \

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

@ -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 <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
{
// ---------------- 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

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

@ -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 <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() = 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