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

Refactored limiter functions.

Browse Source 
The two limiter functions were very similar, and were unified in a single
method, with a bool argument to choose method variety.
This commit is contained in:
Atgeirr Flø Rasmussen
2013-01-21 10:46:06 +01:00
parent 5521278813
commit 1113e65eec
2 changed files with 24 additions and 119 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -411,10 +411,10 @@ namespace Opm
{ {
switch (limiter_method_) { switch (limiter_method_) {
case MinUpwindFace: case MinUpwindFace:
applyMinUpwindFaceLimiter(cell, tof); applyMinUpwindLimiter(cell, true, tof);
break; break;
case MinUpwindAverage: case MinUpwindAverage:
applyMinUpwindAverageLimiter(cell, tof); applyMinUpwindLimiter(cell, false, tof);
break; break;
default: default:
THROW("Limiter type not implemented: " << limiter_method_); THROW("Limiter type not implemented: " << limiter_method_);
@@ -424,18 +424,24 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof) void TransportModelTracerTofDiscGal::applyMinUpwindLimiter(const int cell, const bool face_min, double* tof)
{ {
if (basis_func_->degree() != 1) { if (basis_func_->degree() != 1) {
THROW("This limiter only makes sense for our DG1 implementation."); THROW("This limiter only makes sense for our DG1 implementation.");
} }
// Limiter principles: // Limiter principles:
// 1. Let M be the minimum TOF value on the upstream faces, // 1. Let M be either:
// evaluated in the upstream cells. Then the value at all // - the minimum TOF value of all upstream faces,
// points in this cell shall be at least M. // evaluated in the upstream cells
// Upstream faces whose flux does not exceed the relative // (chosen if face_min is true).
// flux threshold are not considered for this minimum. // or:
// - the minimum average TOF value of all upstream cells
// (chosen if face_min is false).
// 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. // 2. The TOF shall not be below zero in any point.
// Find total upstream/downstream fluxes. // Find total upstream/downstream fluxes.
@@ -492,11 +498,15 @@ namespace Opm
min_here_tof = std::min(min_here_tof, tof_here); min_here_tof = std::min(min_here_tof, tof_here);
if (upstream) { if (upstream) {
if (interior) { if (interior) {
basis_func_->eval(upstream_cell, nc, &basis_nb_[0]); if (face_min) {
const double tof_upstream basis_func_->eval(upstream_cell, nc, &basis_nb_[0]);
= std::inner_product(basis_nb_.begin(), basis_nb_.end(), const double tof_upstream
tof_coeff_ + num_basis*upstream_cell, 0.0); = std::inner_product(basis_nb_.begin(), basis_nb_.end(),
min_upstream_tof = std::min(min_upstream_tof, tof_upstream); tof_coeff_ + num_basis*upstream_cell, 0.0);
min_upstream_tof = std::min(min_upstream_tof, tof_upstream);
} else {
min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]);
}
} else { } else {
// Allow tof down to 0 on inflow boundaries. // Allow tof down to 0 on inflow boundaries.
min_upstream_tof = std::min(min_upstream_tof, 0.0); min_upstream_tof = std::min(min_upstream_tof, 0.0);
@@ -535,110 +545,6 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyMinUpwindAverageLimiter(const int cell, double* tof)
{
if (basis_func_->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 = basis_func_->numBasisFunc();
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];
basis_func_->eval(cell, 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;
basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell);
}
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;
basis_func_->multiplyGradient(limiter, tof + num_basis*cell);
} else {
// std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
}
}
void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess() void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess()

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

@@ -130,8 +130,7 @@ namespace Opm
// (will read data from tof_coeff_, it is ok to call // (will read data from tof_coeff_, it is ok to call
// with tof_coeff as tof argument. // with tof_coeff as tof argument.
void applyLimiter(const int cell, double* tof); void applyLimiter(const int cell, double* tof);
void applyMinUpwindFaceLimiter(const int cell, double* tof); void applyMinUpwindLimiter(const int cell, const bool face_min, double* tof);
void applyMinUpwindAverageLimiter(const int cell, double* tof);
void applyLimiterAsPostProcess(); void applyLimiterAsPostProcess();
void applyLimiterAsSimultaneousPostProcess(); void applyLimiterAsSimultaneousPostProcess();
}; };