From 91f9b3c2debaf47db8adbc7f58c8ef2470b96f9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 6 Dec 2012 10:18:57 +0100 Subject: [PATCH 1/6] Modify precision for quadrature depending on use of CVI method. --- .../transport/reorder/TransportModelTracerTofDiscGal.cpp | 9 ++++++++- .../transport/reorder/TransportModelTracerTofDiscGal.hpp | 1 + 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 073018bf..f65fe91d 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -133,6 +133,12 @@ namespace Opm 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 { @@ -253,7 +259,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]); diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 006c37fe..d088fff5 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -84,6 +84,7 @@ namespace Opm const UnstructuredGrid& grid_; boost::shared_ptr velocity_interpolation_; + bool use_cvi_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell const double* source_; // one volumetric source term per cell From ff1bfc810af75d9326c18d6abcea05cae890fa6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Dec 2012 12:22:56 +0100 Subject: [PATCH 2/6] Implemented 2d quadrature rules for degree 1 and 2 (order 2 and 3). Also added cartesian 2d cell testcase. --- opm/core/grid/CellQuadrature.hpp | 81 +++++++++++++++- opm/core/grid/FaceQuadrature.hpp | 33 +++++-- tests/test_quadratures.cpp | 159 +++++++++++++++++++++---------- 3 files changed, 215 insertions(+), 58 deletions(-) diff --git a/opm/core/grid/CellQuadrature.hpp b/opm/core/grid/CellQuadrature.hpp index 4aeceb50..8fa45429 100644 --- a/opm/core/grid/CellQuadrature.hpp +++ b/opm/core/grid/CellQuadrature.hpp @@ -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) { diff --git a/opm/core/grid/FaceQuadrature.hpp b/opm/core/grid/FaceQuadrature.hpp index 23a29d1a..26a7016a 100644 --- a/opm/core/grid/FaceQuadrature.hpp +++ b/opm/core/grid/FaceQuadrature.hpp @@ -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. - return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]); + 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 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_]; diff --git a/tests/test_quadratures.cpp b/tests/test_quadratures.cpp index 36afe3d6..b659a21b 100644 --- a/tests/test_quadratures.cpp +++ b/tests/test_quadratures.cpp @@ -66,22 +66,6 @@ namespace mutable std::vector 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 void testSingleCase(const UnstructuredGrid& grid, @@ -98,45 +82,124 @@ namespace } // anonymous namespace - -static void test3dCart() +namespace cart2d { - // Set up 3d 1-cell cartesian case. - GridManager g(1, 1, 1); - const UnstructuredGrid& grid = *g.c_grid(); + struct ConstantFunc + { + double operator()(const double*) const + { + return 1.234; + } + }; - // CellQuadrature tests. - testSingleCase(grid, 0, 1, 5.0); - testSingleCase(grid, 0, 2, 5.0); - testSingleCase(grid, 0, 2, 6.25); + struct LinearFunc + { + double operator()(const double* x) const + { + return 1.0*x[0] + 2.0*x[1] + 3.0; + } + }; - // FaceQuadrature tests, degree 1 precision. - testSingleCase(grid, 0, 1, 4.5); - testSingleCase(grid, 1, 1, 5.5); - testSingleCase(grid, 2, 1, 4.0); - testSingleCase(grid, 3, 1, 6.0); - testSingleCase(grid, 4, 1, 4.5); - testSingleCase(grid, 5, 1, 5.5); + 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; + } + }; - // FaceQuadrature tests, degree 2 precision. - testSingleCase(grid, 0, 2, 4.5); - testSingleCase(grid, 1, 2, 5.5); - testSingleCase(grid, 2, 2, 4.0); - testSingleCase(grid, 3, 2, 6.0); - testSingleCase(grid, 4, 2, 4.5); - testSingleCase(grid, 5, 2, 5.5); + static void test() + { + // Set up 2d 1-cell cartesian case. + GridManager g(1, 1); + const UnstructuredGrid& grid = *g.c_grid(); - // FaceQuadrature tests, quadratic function, degree 2 precision. - testSingleCase(grid, 0, 2, 6.0); - testSingleCase(grid, 1, 2, 6.5); - testSingleCase(grid, 2, 2, 5.0); - testSingleCase(grid, 3, 2, 7.5); - testSingleCase(grid, 4, 2, 4.25); - testSingleCase(grid, 5, 2, 8.25); -} + // CellQuadrature tests. + testSingleCase(grid, 0, 1, 1.234); + testSingleCase(grid, 0, 1, 4.5); + testSingleCase(grid, 0, 2, 1.234); + testSingleCase(grid, 0, 2, 4.5); + testSingleCase(grid, 0, 2, 5.25); + // FaceQuadrature tests, degree 1 precision. + testSingleCase(grid, 0, 1, 4); + testSingleCase(grid, 1, 1, 5); + testSingleCase(grid, 2, 1, 3.5); + testSingleCase(grid, 3, 1, 5.5); + + // FaceQuadrature tests, degree 2 precision. + testSingleCase(grid, 0, 2, 4); + testSingleCase(grid, 1, 2, 5); + testSingleCase(grid, 2, 2, 3.5); + testSingleCase(grid, 3, 2, 5.5); + + // FaceQuadrature tests, quadratic function, degree 2 precision. + testSingleCase(grid, 0, 2, 4.0); + testSingleCase(grid, 1, 2, 7.5); + testSingleCase(grid, 2, 2, 4.0); + testSingleCase(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); + const UnstructuredGrid& grid = *g.c_grid(); + + // CellQuadrature tests. + testSingleCase(grid, 0, 1, 5.0); + testSingleCase(grid, 0, 2, 5.0); + testSingleCase(grid, 0, 2, 6.25); + + // FaceQuadrature tests, degree 1 precision. + testSingleCase(grid, 0, 1, 4.5); + testSingleCase(grid, 1, 1, 5.5); + testSingleCase(grid, 2, 1, 4.0); + testSingleCase(grid, 3, 1, 6.0); + testSingleCase(grid, 4, 1, 4.5); + testSingleCase(grid, 5, 1, 5.5); + + // FaceQuadrature tests, degree 2 precision. + testSingleCase(grid, 0, 2, 4.5); + testSingleCase(grid, 1, 2, 5.5); + testSingleCase(grid, 2, 2, 4.0); + testSingleCase(grid, 3, 2, 6.0); + testSingleCase(grid, 4, 2, 4.5); + testSingleCase(grid, 5, 2, 5.5); + + // FaceQuadrature tests, quadratic function, degree 2 precision. + testSingleCase(grid, 0, 2, 6.0); + testSingleCase(grid, 1, 2, 6.5); + testSingleCase(grid, 2, 2, 5.0); + testSingleCase(grid, 3, 2, 7.5); + testSingleCase(grid, 4, 2, 4.25); + testSingleCase(grid, 5, 2, 8.25); + } + +} // namespace cart3d BOOST_AUTO_TEST_CASE(test_quadratures) { - test3dCart(); + cart2d::test(); + cart3d::test(); } From 389a943f341ca31af6fbc0dd75dc56ef51d967c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Dec 2012 13:17:36 +0100 Subject: [PATCH 3/6] Bugfix: multiply by edge length in 2d face quadrature weights. --- opm/core/grid/FaceQuadrature.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/grid/FaceQuadrature.hpp b/opm/core/grid/FaceQuadrature.hpp index 26a7016a..9ed06353 100644 --- a/opm/core/grid/FaceQuadrature.hpp +++ b/opm/core/grid/FaceQuadrature.hpp @@ -152,7 +152,7 @@ namespace Opm 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 simpsonw[index]; + return grid_.face_areas[face_]*simpsonw[index]; } ASSERT(dim == 3); const double* fc = grid_.face_centroids + dim*face_; From b66920203d3bf2c772368f49b5ded976e003b4db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 17 Dec 2012 11:07:40 +0100 Subject: [PATCH 4/6] Ensure correct output in case of Lapack error. Make sure that we output the original values, since Lapack may overwrite those used in the call. --- opm/core/utility/VelocityInterpolation.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/opm/core/utility/VelocityInterpolation.cpp b/opm/core/utility/VelocityInterpolation.cpp index f895d88f..2714b37c 100644 --- a/opm/core/utility/VelocityInterpolation.cpp +++ b/opm/core/utility/VelocityInterpolation.cpp @@ -101,6 +101,7 @@ namespace Opm std::vector N(dim*dim); // Normals matrix. Fortran ordering! std::vector orig_N(dim*dim); // Normals matrix. Fortran ordering! std::vector f(dim); // Flux vector. + std::vector orig_f(dim); // Flux vector. std::vector piv(dim); // For LAPACK solve const SparseTable& all_ci = bcmethod_.cornerInfo(); const std::vector& 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); } From d8de94a5dab929a24010e1df246e4cc0929d7666 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 18 Dec 2012 14:12:56 +0100 Subject: [PATCH 5/6] Bugfix: use sufficient quadrature order. --- opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index f65fe91d..20024d3a 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -242,7 +242,7 @@ namespace Opm // \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds // (where u_h^{ext} is the upstream unknown (tof)). 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]); From 668aba99423dc5b7a8e874f1490c43e816eacf12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 18 Dec 2012 14:15:31 +0100 Subject: [PATCH 6/6] Added limiter for DG1, parameter 'use_limiter'. The limiter is experimental and unfinished, untested work in progress. Limiter is therefore inactive by default. Also fixed a minor bug: use_cvi_ was not initialized. --- examples/compute_tof.cpp | 4 +- examples/compute_tof_from_files.cpp | 4 +- .../TransportModelTracerTofDiscGal.cpp | 103 +++++++++++++++++- .../TransportModelTracerTofDiscGal.hpp | 8 +- 4 files changed, 111 insertions(+), 8 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 259bd253..e308782b 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -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 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); diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index 17ad7b85..fd78e404 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -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 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); diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 20024d3a..9738dc9c 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -128,8 +128,11 @@ 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) { @@ -230,17 +233,21 @@ 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, 2*degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { @@ -358,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); + } } @@ -376,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 diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index d088fff5..3188b371 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -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,9 +83,11 @@ namespace Opm TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&); TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&); + // Data members const UnstructuredGrid& grid_; boost::shared_ptr 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 @@ -100,6 +103,9 @@ namespace Opm std::vector basis_nb_; std::vector grad_basis_; std::vector velocity_; + + // Private methods + void useLimiter(const int cell); }; } // namespace Opm