From 3e723bc9650d18685a59e5dc9dd440bf22f75a96 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] 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_from_files.cpp | 4 +- .../TransportModelTracerTofDiscGal.cpp | 103 +++++++++++++++++- .../TransportModelTracerTofDiscGal.hpp | 8 +- 3 files changed, 108 insertions(+), 7 deletions(-) diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index 17ad7b85b..fd78e4045 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 20024d3aa..9738dc9cc 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 d088fff5a..3188b3716 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