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

Merge pull request #111 from atgeirr/dg-improvements

Browse Source 
Further improvements to limiters for DG time-of-flight computations
This commit is contained in:
Bård Skaflestad 2013-01-08 07:45:40 -08:00
commit a197dbdd0b
6 changed files with 378 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
examples/compute_tof.cpp
View File

@ -169,16 +169,19 @@ main(int argc, char** argv)
// Choice of tof solver.
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;
// 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);
use_cvi = param.getDefault("use_cvi", false);
use_limiter = param.getDefault("use_limiter", false);
dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(*grid->c_grid(), param));
} else {
use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
}
bool compute_tracer = param.getDefault("compute_tracer", false);
if (use_dg && compute_tracer) {
THROW("DG for tracer not yet implemented.");
}
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
@ -232,12 +235,16 @@ main(int argc, char** argv)
// Solve time-of-flight.
transport_timer.start();
std::vector<double> tof;
std::vector<double> tracer;
if (use_dg) {
Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid(), use_cvi, use_limiter);
tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
} else {
Opm::TransportModelTracerTof tofsolver(*grid->c_grid(), use_multidim_upwind);
tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof);
if (compute_tracer) {
tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tof, tracer);
} else {
tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof);
}
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
@ -249,7 +256,17 @@ main(int argc, char** argv)
if (output) {
std::string tof_filename = output_dir + "/tof.txt";
std::ofstream tof_stream(tof_filename.c_str());
tof_stream.precision(16);
std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
if (compute_tracer) {
std::string tracer_filename = output_dir + "/tracer.txt";
std::ofstream tracer_stream(tracer_filename.c_str());
tracer_stream.precision(16);
const int nt = tracer.size()/num_cells;
for (int i = 0; i < nt*num_cells; ++i) {
tracer_stream << tracer[i] << (((i + 1) % nt == 0) ? '\n' : ' ');
}
}
}
std::cout << "\n\n================ End of simulation ===============\n"

10
examples/compute_tof_from_files.cpp
View File

@ -123,13 +123,12 @@ main(int argc, char** argv)
// Choice of tof solver.
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;
// 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);
use_cvi = param.getDefault("use_cvi", false);
use_limiter = param.getDefault("use_limiter", false);
dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(grid, param));
} else {
use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
}
@ -164,8 +163,7 @@ main(int argc, char** argv)
std::vector<double> tof;
std::vector<double> tracer;
if (use_dg) {
Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi, use_limiter);
tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof);
dg_solver->solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof);
} else {
Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind);
if (compute_tracer) {

27
opm/core/transport/reorder/TransportModelInterface.cpp
View File

@ -29,15 +29,18 @@
void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux)
{
// Compute reordered sequence of single-cell problems
std::vector<int> sequence(grid.number_of_cells);
std::vector<int> components(grid.number_of_cells + 1);
sequence_.resize(grid.number_of_cells);
components_.resize(grid.number_of_cells + 1);
int ncomponents;
time::StopWatch clock;
clock.start();
compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
compute_sequence(&grid, darcyflux, &sequence_[0], &components_[0], &ncomponents);
clock.stop();
std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl;
// Make vector's size match actual used data.
components_.resize(ncomponents + 1);
// Invoke appropriate solve method for each interdependent component.
for (int comp = 0; comp < ncomponents; ++comp) {
#if 0
@ -50,11 +53,23 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g
}
#endif
#endif
const int comp_size = components[comp + 1] - components[comp];
const int comp_size = components_[comp + 1] - components_[comp];
if (comp_size == 1) {
solveSingleCell(sequence[components[comp]]);
solveSingleCell(sequence_[components_[comp]]);
} else {
solveMultiCell(comp_size, &sequence[components[comp]]);
solveMultiCell(comp_size, &sequence_[components_[comp]]);
}
}
}
const std::vector<int>& Opm::TransportModelInterface::sequence() const
{
return sequence_;
}
const std::vector<int>& Opm::TransportModelInterface::components() const
{
return components_;
}

10
opm/core/transport/reorder/TransportModelInterface.hpp
View File

@ -20,6 +20,8 @@
#ifndef OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED
#define OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED
#include <vector>
struct UnstructuredGrid;
namespace Opm
@ -31,7 +33,8 @@ namespace Opm
/// method that will have an interface geared to the model's
/// needs. (The solve() method is therefore not virtual in this
/// class.) The reorderAndTransport() method is provided as an aid
/// to implementing solve() in subclasses.
/// to implementing solve() in subclasses, together with the
/// sequence() and components() methods for accessing the ordering.
class TransportModelInterface
{
public:
@ -41,6 +44,11 @@ namespace Opm
virtual void solveMultiCell(const int num_cells, const int* cells) = 0;
protected:
void reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux);
const std::vector<int>& sequence() const;
const std::vector<int>& components() const;
private:
std::vector<int> sequence_;
std::vector<int> components_;
};

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

@ -23,6 +23,7 @@
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/VelocityInterpolation.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/linalg/blas_lapack.h>
#include <algorithm>
#include <cmath>
@ -125,27 +126,66 @@ namespace Opm
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
/// \param[in] param Parameters for the solver.
/// The following parameters are accepted (defaults):
/// 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.
/// limiter_method ("MinUpwindFace") Limiter method used. Accepted methods are:
/// MinUpwindFace Limit cell tof to >= inflow face tofs.
/// limiter_usage ("DuringComputations") Usage pattern for limiter. Accepted choices are:
/// DuringComputations Apply limiter to cells as they are computed,
/// so downstream cells' solutions may be affected
/// by limiting in upstream cells.
/// AsPostProcess Apply in dependency order, but only after
/// computing (unlimited) solution.
/// AsSimultaneousPostProcess Apply to each cell independently, using un-
/// limited solution in neighbouring cells.
TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi,
const bool use_limiter)
const parameter::ParameterGroup& param)
: grid_(grid),
use_cvi_(use_cvi),
use_limiter_(use_limiter),
use_cvi_(false),
use_limiter_(false),
limiter_relative_flux_threshold_(1e-3),
limiter_method_(MinUpwindAverage),
limiter_usage_(DuringComputations),
coord_(grid.dimensions),
velocity_(grid.dimensions)
{
use_cvi_ = param.getDefault("use_cvi", use_cvi_);
use_limiter_ = param.getDefault("use_limiter", use_limiter_);
if (use_limiter_) {
limiter_relative_flux_threshold_ = param.getDefault("limiter_relative_flux_threshold",
limiter_relative_flux_threshold_);
const std::string limiter_method_str = param.getDefault<std::string>("limiter_method", "MinUpwindAverage");
if (limiter_method_str == "MinUpwindFace") {
limiter_method_ = MinUpwindFace;
} else if (limiter_method_str == "MinUpwindAverage") {
limiter_method_ = MinUpwindAverage;
} else {
THROW("Unknown limiter method: " << limiter_method_str);
}
const std::string limiter_usage_str = param.getDefault<std::string>("limiter_usage", "DuringComputations");
if (limiter_usage_str == "DuringComputations") {
limiter_usage_ = DuringComputations;
} else if (limiter_usage_str == "AsPostProcess") {
limiter_usage_ = AsPostProcess;
} else if (limiter_usage_str == "AsSimultaneousPostProcess") {
limiter_usage_ = AsSimultaneousPostProcess;
} else {
THROW("Unknown limiter usage spec: " << limiter_usage_str);
}
}
// 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));
if (use_cvi_) {
velocity_interpolation_.reset(new VelocityInterpolationECVI(grid_));
} else {
velocity_interpolation_.reset(new VelocityInterpolationConstant(grid));
velocity_interpolation_.reset(new VelocityInterpolationConstant(grid_));
}
}
@ -194,6 +234,19 @@ namespace Opm
grad_basis_.resize(num_basis*grid_.dimensions);
velocity_interpolation_->setupFluxes(darcyflux);
reorderAndTransport(grid_, darcyflux);
switch (limiter_usage_) {
case AsPostProcess:
applyLimiterAsPostProcess();
break;
case AsSimultaneousPostProcess:
applyLimiterAsSimultaneousPostProcess();
break;
case DuringComputations:
// Do nothing.
break;
default:
THROW("Unknown limiter usage choice: " << limiter_usage_);
}
}
@ -370,8 +423,8 @@ namespace Opm
std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell);
// Apply limiter.
if (degree_ > 0 && use_limiter_) {
useLimiter(cell);
if (degree_ > 0 && use_limiter_ && limiter_usage_ == DuringComputations) {
applyLimiter(cell, tof_coeff_);
}
}
@ -389,7 +442,24 @@ namespace Opm
void TransportModelTracerTofDiscGal::useLimiter(const int cell)
void TransportModelTracerTofDiscGal::applyLimiter(const int cell, double* tof)
{
switch (limiter_method_) {
case MinUpwindFace:
applyMinUpwindFaceLimiter(cell, tof);
break;
case MinUpwindAverage:
applyMinUpwindAverageLimiter(cell, tof);
break;
default:
THROW("Limiter type not implemented: " << limiter_method_);
}
}
void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof)
{
if (degree_ != 1) {
THROW("This limiter only makes sense for our DG1 implementation.");
@ -399,15 +469,38 @@ namespace Opm
// 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.
// Upstream faces whose flux does not exceed the relative
// flux threshold are not considered for this minimum.
// 2. The TOF shall not be below zero in any point.
// Find total upstream/downstream fluxes.
double upstream_flux = 0.0;
double downstream_flux = 0.0;
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;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux < 0.0) {
upstream_flux += flux;
} else {
downstream_flux += flux;
}
}
// In the presence of sources, significant fluxes may be missing from the computed fluxes,
// setting the total flux to the (positive) maximum avoids this: since source is either
// inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
const double total_flux = std::max(-upstream_flux, downstream_flux);
// Find minimum tof on upstream faces and for this cell.
const int dim = grid_.dimensions;
const int num_basis = DGBasis::numBasisFunc(dim, degree_);
double max_slope_mult = 0.0;
double min_upstream_tof = 1e100;
double min_here_tof = 1e100;
int num_upstream_faces = 0;
// For inflow faces, ensure that cell tof does not dip below
// the minimum value from upstream (for all faces).
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;
@ -419,57 +512,206 @@ namespace Opm
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
if (flux >= 0.0) {
// This is a downstream face.
continue;
const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_);
if (upstream) {
++num_upstream_faces;
}
++num_upstream_faces;
bool interior = (upstream_cell >= 0);
// 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;
// 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]);
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);
min_here_tof = std::min(min_here_tof, 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 (interior) {
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_tof = std::min(min_upstream_tof, tof_upstream);
} else {
// Allow tof down to 0 on inflow boundaries.
min_upstream_tof = std::min(min_upstream_tof, 0.0);
}
}
}
// Compute maximum slope multiplier.
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;
max_slope_mult = 0.0;
tof_coeff_[num_basis*cell] = min_upstream;
break;
}
const double face_mult = (tof_c - min_upstream)/(tof_c - min_here);
max_slope_mult = std::max(max_slope_mult, face_mult);
}
ASSERT(max_slope_mult >= 0.0);
// Compute slope multiplier (limiter).
if (num_upstream_faces == 0) {
min_upstream_tof = 0.0;
min_here_tof = 0.0;
}
if (min_upstream_tof < 0.0) {
min_upstream_tof = 0.0;
}
const double tof_c = tof_coeff_[num_basis*cell];
double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof);
if (tof_c < min_upstream_tof) {
// Handle by setting a flat solution.
std::cout << "Trouble in cell " << cell << std::endl;
limiter = 0.0;
tof[num_basis*cell] = min_upstream_tof;
}
ASSERT(limiter >= 0.0);
// Actually do the limiting (if applicable).
const double limiter = max_slope_mult;
if (num_upstream_faces > 0 && limiter < 1.0) {
std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl;
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;
tof[i] *= limiter;
}
} else {
std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
// std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
}
}
void TransportModelTracerTofDiscGal::applyMinUpwindAverageLimiter(const int cell, double* tof)
{
if (degree_ != 1) {
THROW("This limiter only makes sense for our DG1 implementation.");
}
// Limiter principles:
// 1. Let M be the average TOF value of the upstream cells.
/// Then the value at all points in this cell shall be at least M.
// Upstream faces whose flux does not exceed the relative
// flux threshold are not considered for this minimum.
// 2. The TOF shall not be below zero in any point.
// Find total upstream/downstream fluxes.
double upstream_flux = 0.0;
double downstream_flux = 0.0;
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;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux < 0.0) {
upstream_flux += flux;
} else {
downstream_flux += flux;
}
}
// In the presence of sources, significant fluxes may be missing from the computed fluxes,
// setting the total flux to the (positive) maximum avoids this: since source is either
// inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
const double total_flux = std::max(-upstream_flux, downstream_flux);
// Find minimum tof on upstream faces and for this cell.
const int dim = grid_.dimensions;
const int num_basis = DGBasis::numBasisFunc(dim, degree_);
double min_upstream_tof = 1e100;
double min_here_tof = 1e100;
int num_upstream_faces = 0;
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];
}
const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_);
if (upstream) {
++num_upstream_faces;
}
bool interior = (upstream_cell >= 0);
// 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]);
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) {
min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]);
} else {
// Allow tof down to 0 on inflow boundaries.
min_upstream_tof = std::min(min_upstream_tof, 0.0);
}
}
}
}
// Compute slope multiplier (limiter).
if (num_upstream_faces == 0) {
min_upstream_tof = 0.0;
min_here_tof = 0.0;
}
if (min_upstream_tof < 0.0) {
min_upstream_tof = 0.0;
}
const double tof_c = tof_coeff_[num_basis*cell];
double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof);
if (tof_c < min_upstream_tof) {
// Handle by setting a flat solution.
std::cout << "Trouble in cell " << cell << std::endl;
limiter = 0.0;
tof[num_basis*cell] = min_upstream_tof;
}
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;
}
} else {
// std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
}
}
void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess()
{
// Apply the limiter sequentially to all cells.
// This means that a cell's limiting behaviour may be affected by
// any limiting applied to its upstream cells.
const std::vector<int>& seq = TransportModelInterface::sequence();
const int nc = seq.size();
ASSERT(nc == grid_.number_of_cells);
for (int i = 0; i < nc; ++i) {
const int cell = seq[i];
applyLimiter(cell, tof_coeff_);
}
}
void TransportModelTracerTofDiscGal::applyLimiterAsSimultaneousPostProcess()
{
// Apply the limiter simultaneously to all cells.
// 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_);
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]);
}
std::copy(tof_coeffs_new.begin(), tof_coeffs_new.end(), tof_coeff_);
}
} // namespace Opm

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

@ -33,6 +33,7 @@ namespace Opm
class IncompPropertiesInterface;
class VelocityInterpolationInterface;
namespace parameter { class ParameterGroup; }
/// Implements a discontinuous Galerkin solver for
/// (single-phase) time-of-flight using reordering.
@ -48,11 +49,23 @@ namespace Opm
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
/// \param[in] param Parameters for the solver.
/// The following parameters are accepted (defaults):
/// 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.
/// limiter_method ("MinUpwindFace") Limiter method used. Accepted methods are:
/// MinUpwindFace Limit cell tof to >= inflow face tofs.
/// limiter_usage ("DuringComputations") Usage pattern for limiter. Accepted choices are:
/// DuringComputations Apply limiter to cells as they are computed,
/// so downstream cells' solutions may be affected
/// by limiting in upstream cells.
/// AsPostProcess Apply in dependency order, but only after
/// computing (unlimited) solution.
/// AsSimultaneousPostProcess Apply to each cell independently, using un-
/// limited solution in neighbouring cells.
TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi,
const bool use_limiter = false);
const parameter::ParameterGroup& param);
/// Solve for time-of-flight.
@ -88,6 +101,11 @@ namespace Opm
boost::shared_ptr<VelocityInterpolationInterface> velocity_interpolation_;
bool use_cvi_;
bool use_limiter_;
double limiter_relative_flux_threshold_;
enum LimiterMethod { MinUpwindFace, MinUpwindAverage };
LimiterMethod limiter_method_;
enum LimiterUsage { DuringComputations, AsPostProcess, AsSimultaneousPostProcess };
LimiterUsage limiter_usage_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
@ -105,7 +123,15 @@ namespace Opm
std::vector<double> velocity_;
// Private methods
void useLimiter(const int cell);
// Apply some limiter, writing to array tof
// (will read data from tof_coeff_, it is ok to call
// with tof_coeff as tof argument.
void applyLimiter(const int cell, double* tof);
void applyMinUpwindFaceLimiter(const int cell, double* tof);
void applyMinUpwindAverageLimiter(const int cell, double* tof);
void applyLimiterAsPostProcess();
void applyLimiterAsSimultaneousPostProcess();
};
} // namespace Opm