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

Merge pull request #109 from atgeirr/dg-improvements

Browse Source 
DG time-of-flight improvements
This commit is contained in:
Bård Skaflestad 2012-12-19 09:09:33 -08:00
commit 8f47facf82
8 changed files with 339 additions and 69 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
examples/compute_tof.cpp
View File

@ -170,10 +170,12 @@ main(int argc, char** argv)
bool use_dg = param.getDefault("use_dg", false);
int dg_degree = -1;
bool use_cvi = false;
bool use_limiter = false;
bool use_multidim_upwind = false;
if (use_dg) {
dg_degree = param.getDefault("dg_degree", 0);
use_cvi = param.getDefault("use_cvi", false);
use_limiter = param.getDefault("use_limiter", false);
} else {
use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
}
@ -231,7 +233,7 @@ main(int argc, char** argv)
transport_timer.start();
std::vector<double> tof;
if (use_dg) {
Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid(), use_cvi);
Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid(), use_cvi, use_limiter);
tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
} else {
Opm::TransportModelTracerTof tofsolver(*grid->c_grid(), use_multidim_upwind);

4
examples/compute_tof_from_files.cpp
View File

@ -124,10 +124,12 @@ main(int argc, char** argv)
bool use_dg = param.getDefault("use_dg", false);
int dg_degree = -1;
bool use_cvi = false;
bool use_limiter = false;
bool use_multidim_upwind = false;
if (use_dg) {
dg_degree = param.getDefault("dg_degree", 0);
use_cvi = param.getDefault("use_cvi", false);
use_limiter = param.getDefault("use_limiter", false);
} else {
use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
}
@ -157,7 +159,7 @@ main(int argc, char** argv)
transport_timer.start();
std::vector<double> tof;
if (use_dg) {
Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi);
Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi, use_limiter);
tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof);
} else {
Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind);

81
opm/core/grid/CellQuadrature.hpp
View File

@ -54,6 +54,18 @@ namespace Opm
double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] };
return std::fabs(determinantOf(a, b, c) / 6.0);
}
/// Calculates the area of a triangle consisting of 3 vertices
/// with 2-dimensional coordinates
inline double triangleArea2d(const double* p0,
const double* p1,
const double* p2)
{
double a[2] = { p1[0] - p0[0], p1[1] - p0[1] };
double b[2] = { p2[0] - p0[0], p2[1] - p0[1] };
double a_cross_b = a[0]*b[1] - a[1]*b[0];
return 0.5*std::fabs(a_cross_b);
}
} // anonymous namespace
@ -61,10 +73,35 @@ namespace Opm
/// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by cell volume,
/// so weights always sum to cell volume.
///
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = cell volume, x_0 = cell centroid
/// Degree 2 method:
///
/// Degree 2 method for 2d (but see the note):
/// Based on subdivision of the cell into triangles,
/// with the centroid as a common vertex, and the triangle
/// edge midpoint rule.
/// Triangle i consists of the centroid C, nodes N_i and N_{i+1}.
/// Its area is A_i.
/// n = 2 * nn (nn = num nodes in face)
/// For i = 0..(nn-1):
/// w_i = 1/3 A_i.
/// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i
/// x_i = (N_i + N_{i+1})/2
/// x_{nn+i} = (C + N_i)/2
/// All N and A indices are interpreted cyclic, modulus nn.
/// Note: for simplicity of implementation, we currently use
/// n = 3 * nn
/// For i = 0..(nn-1):
/// w_{3*i + {0,1,2}} = 1/3 A_i
/// x_{3*i} = (N_i + N_{i+1})/2
/// x_{3*i + {1,2}} = (C + N_{i,i+1})/2
/// This is simpler, because we can implement it easily
/// based on iteration over faces without requiring any
/// particular (cyclic) ordering.
///
/// Degree 2 method for 3d:
/// Based on subdivision of each cell face into triangles
/// with the face centroid as a common vertex, and then
/// subdividing the cell into tetrahedra with the cell
@ -87,8 +124,8 @@ namespace Opm
const int degree)
: grid_(grid), cell_(cell), degree_(degree)
{
if (grid.dimensions != 3) {
THROW("CellQuadrature only implemented for 3D case.");
if (grid.dimensions > 3) {
THROW("CellQuadrature only implemented for up to 3 dimensions.");
}
if (degree > 2) {
THROW("CellQuadrature exact for polynomial degrees > 1 not implemented.");
@ -97,10 +134,14 @@ namespace Opm
int numQuadPts() const
{
if (degree_ < 2) {
if (degree_ < 2 || grid_.dimensions == 1) {
return 1;
}
// Degree 2 case.
if (grid_.dimensions == 2) {
return 3*(grid_.cell_facepos[cell_ + 1] - grid_.cell_facepos[cell_]);
}
ASSERT(grid_.dimensions == 3);
int sumnodes = 0;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
@ -118,6 +159,29 @@ namespace Opm
return;
}
// Degree 2 case.
if (dim == 2) {
if (index % 3 == 0) {
// Boundary midpoint. This is the face centroid.
const int hface = grid_.cell_facepos[cell_] + index/3;
const int face = grid_.cell_faces[hface];
const double* fc = grid_.face_centroids + dim*face;
std::copy(fc, fc + dim, coord);
} else {
// Interiour midpoint. This is the average of the
// cell centroid and a face node (they should
// always have two nodes in 2d).
const int hface = grid_.cell_facepos[cell_] + index/3;
const int face = grid_.cell_faces[hface];
const int nodeoff = (index % 3) - 1; // == 0 or 1
const int node = grid_.face_nodes[grid_.face_nodepos[face] + nodeoff];
const double* nc = grid_.node_coordinates + dim*node;
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dd] + cc[dd]);
}
}
return;
}
ASSERT(dim == 3);
int tetindex = index / 4;
const int subindex = index % 4;
const double* nc = grid_.node_coordinates;
@ -155,6 +219,15 @@ namespace Opm
// Degree 2 case.
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
if (dim == 2) {
const int hface = grid_.cell_facepos[cell_] + index/3;
const int face = grid_.cell_faces[hface];
const int* nptr = grid_.face_nodes + grid_.face_nodepos[face];
const double* nc0 = grid_.node_coordinates + dim*nptr[0];
const double* nc1 = grid_.node_coordinates + dim*nptr[1];
return triangleArea2d(nc0, nc1, cc)/3.0;
}
ASSERT(dim == 3);
int tetindex = index / 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {

31
opm/core/grid/FaceQuadrature.hpp
View File

@ -57,10 +57,15 @@ namespace Opm
/// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by face area,
/// so weights always sum to face area.
///
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = face area, x_0 = face centroid
/// Degree 2 method:
///
/// Degree 2 method for 2d:
/// Simpson's method (actually this is degree 3).
///
/// Degree 2 method for 3d:
/// Based on subdivision of the face into triangles,
/// with the centroid as a common vertex, and the triangle
/// edge midpoint rule.
@ -81,8 +86,8 @@ namespace Opm
const int degree)
: grid_(grid), face_(face), degree_(degree)
{
if (grid_.dimensions != 3) {
THROW("FaceQuadrature only implemented for 3D case.");
if (grid_.dimensions > 3) {
THROW("FaceQuadrature only implemented for up to 3 dimensions.");
}
if (degree_ > 2) {
THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented.");
@ -91,18 +96,22 @@ namespace Opm
int numQuadPts() const
{
if (degree_ < 2) {
if (degree_ < 2 || grid_.dimensions < 2) {
return 1;
}
// Degree 2 case.
if (grid_.dimensions == 2) {
return 3;
} else {
return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
}
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
if (degree_ < 2) {
if (degree_ < 2 || dim < 2) {
std::copy(fc, fc + dim, coord);
return;
}
@ -110,6 +119,13 @@ namespace Opm
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (dim == 2) {
ASSERT(nn == 2);
const double* pa[3] = { nc + dim*fnodes[0], fc, nc + dim*fnodes[1] };
std::copy(pa[index], pa[index] + dim, coord);
return;
}
ASSERT(dim == 3);
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
@ -134,6 +150,11 @@ namespace Opm
}
// Degree 2 case.
const int dim = grid_.dimensions;
if (dim == 2) {
const double simpsonw[3] = { 1.0/6.0, 4.0/6.0, 1.0/6.0 };
return grid_.face_areas[face_]*simpsonw[index];
}
ASSERT(dim == 3);
const double* fc = grid_.face_centroids + dim*face_;
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];

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

@ -128,11 +128,20 @@ namespace Opm
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi)
const bool use_cvi,
const bool use_limiter)
: grid_(grid),
use_cvi_(use_cvi),
use_limiter_(use_limiter),
coord_(grid.dimensions),
velocity_(grid.dimensions)
{
// A note about the use_cvi_ member variable:
// In principle, we should not need it, since the choice of velocity
// interpolation is made below, but we may need to use higher order
// quadrature to exploit CVI, so we store the choice.
// An alternative would be to add a virtual method isConstant() to
// the VelocityInterpolationInterface.
if (use_cvi) {
velocity_interpolation_.reset(new VelocityInterpolationECVI(grid));
} else {
@ -224,19 +233,23 @@ namespace Opm
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
if (upstream_cell < 0) {
// This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
continue;
}
if (flux >= 0.0) {
// This is an outflow boundary.
continue;
}
if (upstream_cell < 0) {
// This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
continue;
}
// Do quadrature over the face to compute
// \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds
// (where u_h^{ext} is the upstream unknown (tof)).
// Quadrature degree set to 2*D, since u_h^{ext} varies
// with degree D, and b_j too. We assume that the normal
// 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, degree_);
FaceQuadrature quad(grid_, face, 2*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]);
@ -253,7 +266,8 @@ namespace Opm
// Compute cell jacobian contribution. We use Fortran ordering
// for jac_, i.e. rows cycling fastest.
{
CellQuadrature quad(grid_, cell, 2*degree_ - 1);
const int deg_needed = use_cvi_ ? 2*degree_ : 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)
quad.quadPtCoord(quad_pt, &coord_[0]);
@ -351,8 +365,14 @@ namespace Opm
}
THROW("Lapack error: " << info << " encountered in cell " << cell);
}
// The solution ends up in rhs_, so we must copy it.
std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell);
// Apply limiter.
if (degree_ > 0 && use_limiter_) {
useLimiter(cell);
}
}
@ -369,4 +389,84 @@ namespace Opm
void TransportModelTracerTofDiscGal::useLimiter(const int cell)
{
if (degree_ != 1) {
THROW("This limiter only makes sense for our DG1 implementation.");
}
// Limiter principles:
// 1. Let M be the minimum TOF value on the upstream faces,
// evaluated in the upstream cells. Then the value at all
// points in this cell shall be at least M.
// 2. The TOF shall not be below zero in any point.
const int dim = grid_.dimensions;
const int num_basis = DGBasis::numBasisFunc(dim, degree_);
double limiter = 1e100;
// For inflow faces, ensure that cell tof does not dip below
// the minimum value from upstream (for that face).
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
int upstream_cell = -1;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
upstream_cell = grid_.face_cells[2*face+1];
} else {
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
// Evaluate the solution in all corners, and find the appropriate limiter.
bool upstream = (upstream_cell >= 0 && flux < 0.0);
double min_upstream = upstream ? 1e100 : 0.0;
double min_here = 1e100;
for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
DGBasis::eval(grid_, cell, degree_, nc, &basis_[0]);
const double tof_here = std::inner_product(basis_.begin(), basis_.end(),
tof_coeff_ + num_basis*cell, 0.0);
min_here = std::min(min_here, tof_here);
if (upstream) {
DGBasis::eval(grid_, upstream_cell, degree_, nc, &basis_nb_[0]);
const double tof_upstream
= std::inner_product(basis_nb_.begin(), basis_nb_.end(),
tof_coeff_ + num_basis*upstream_cell, 0.0);
min_upstream = std::min(min_upstream, tof_upstream);
}
}
if (min_here < min_upstream) {
// Must limit slope.
const double tof_c = tof_coeff_[num_basis*cell];
if (tof_c < min_upstream) {
// Handle by setting a flat solution.
std::cout << "Trouble in cell " << cell << std::endl;
limiter = 0.0;
tof_coeff_[num_basis*cell] = min_upstream;
break;
}
const double face_limit = (tof_c - min_upstream)/(tof_c - min_here);
limiter = std::min(limiter, face_limit);
}
}
if (limiter < 0.0) {
THROW("Error in limiter.");
}
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_coeff_[i] *= limiter;
}
} else {
std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
}
}
} // namespace Opm

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

@ -51,7 +51,8 @@ namespace Opm
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi);
const bool use_cvi,
const bool use_limiter = false);
/// Solve for time-of-flight.
@ -82,8 +83,11 @@ namespace Opm
TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&);
TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&);
// Data members
const UnstructuredGrid& grid_;
boost::shared_ptr<VelocityInterpolationInterface> velocity_interpolation_;
bool use_cvi_;
bool use_limiter_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
@ -99,6 +103,9 @@ namespace Opm
std::vector<double> basis_nb_;
std::vector<double> grad_basis_;
std::vector<double> velocity_;
// Private methods
void useLimiter(const int cell);
};
} // namespace Opm

4
opm/core/utility/VelocityInterpolation.cpp
View File

@ -101,6 +101,7 @@ namespace Opm
std::vector<double> N(dim*dim); // Normals matrix. Fortran ordering!
std::vector<double> orig_N(dim*dim); // Normals matrix. Fortran ordering!
std::vector<double> f(dim); // Flux vector.
std::vector<double> orig_f(dim); // Flux vector.
std::vector<MAT_SIZE_T> piv(dim); // For LAPACK solve
const SparseTable<WachspressCoord::CornerInfo>& all_ci = bcmethod_.cornerInfo();
const std::vector<int>& adj_faces = bcmethod_.adjacentFaces();
@ -129,6 +130,7 @@ namespace Opm
MAT_SIZE_T ldb = n;
MAT_SIZE_T info = 0;
orig_N = N;
orig_f = f;
dgesv_(&n, &nrhs, &N[0], &lda, &piv[0], &f[0], &ldb, &info);
if (info != 0) {
// Print the local matrix and rhs.
@ -142,7 +144,7 @@ namespace Opm
}
std::cerr << "and f = \n";
for (int row = 0; row < n; ++row) {
std::cerr << " " << f[row] << '\n';
std::cerr << " " << orig_f[row] << '\n';
}
THROW("Lapack error: " << info << " encountered in cell " << cell);
}

99
tests/test_quadratures.cpp
View File

@ -66,22 +66,6 @@ namespace
mutable std::vector<double> pt;
};
struct LinearFunc
{
double operator()(const double* x) const
{
return 1.0*x[0] + 2.0*x[1] + 1.0*x[2] + 3.0;
}
};
struct QuadraticFunc
{
double operator()(const double* x) const
{
return 1.0*x[0]*x[1] + 2.0*x[1] + 4.0*x[2] + 3.0;
}
};
template <class Quadrature, class Func>
void testSingleCase(const UnstructuredGrid& grid,
@ -98,8 +82,85 @@ namespace
} // anonymous namespace
namespace cart2d
{
struct ConstantFunc
{
double operator()(const double*) const
{
return 1.234;
}
};
static void test3dCart()
struct LinearFunc
{
double operator()(const double* x) const
{
return 1.0*x[0] + 2.0*x[1] + 3.0;
}
};
struct QuadraticFunc
{
double operator()(const double* x) const
{
return 3.0*x[0]*x[0] + 1.0*x[0]*x[1] + 2.0*x[1] + 3.0;
}
};
static void test()
{
// Set up 2d 1-cell cartesian case.
GridManager g(1, 1);
const UnstructuredGrid& grid = *g.c_grid();
// CellQuadrature tests.
testSingleCase<CellQuadrature, ConstantFunc>(grid, 0, 1, 1.234);
testSingleCase<CellQuadrature, LinearFunc>(grid, 0, 1, 4.5);
testSingleCase<CellQuadrature, ConstantFunc>(grid, 0, 2, 1.234);
testSingleCase<CellQuadrature, LinearFunc>(grid, 0, 2, 4.5);
testSingleCase<CellQuadrature, QuadraticFunc>(grid, 0, 2, 5.25);
// FaceQuadrature tests, degree 1 precision.
testSingleCase<FaceQuadrature, LinearFunc>(grid, 0, 1, 4);
testSingleCase<FaceQuadrature, LinearFunc>(grid, 1, 1, 5);
testSingleCase<FaceQuadrature, LinearFunc>(grid, 2, 1, 3.5);
testSingleCase<FaceQuadrature, LinearFunc>(grid, 3, 1, 5.5);
// FaceQuadrature tests, degree 2 precision.
testSingleCase<FaceQuadrature, LinearFunc>(grid, 0, 2, 4);
testSingleCase<FaceQuadrature, LinearFunc>(grid, 1, 2, 5);
testSingleCase<FaceQuadrature, LinearFunc>(grid, 2, 2, 3.5);
testSingleCase<FaceQuadrature, LinearFunc>(grid, 3, 2, 5.5);
// FaceQuadrature tests, quadratic function, degree 2 precision.
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 0, 2, 4.0);
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 1, 2, 7.5);
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 2, 2, 4.0);
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 3, 2, 6.5);
}
} // namespace cart2d
namespace cart3d
{
struct LinearFunc
{
double operator()(const double* x) const
{
return 1.0*x[0] + 2.0*x[1] + 1.0*x[2] + 3.0;
}
};
struct QuadraticFunc
{
double operator()(const double* x) const
{
return 1.0*x[0]*x[1] + 2.0*x[1] + 4.0*x[2] + 3.0;
}
};
static void test()
{
// Set up 3d 1-cell cartesian case.
GridManager g(1, 1, 1);
@ -135,8 +196,10 @@ static void test3dCart()
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 5, 2, 8.25);
}
} // namespace cart3d
BOOST_AUTO_TEST_CASE(test_quadratures)
{
test3dCart();
cart2d::test();
cart3d::test();
}