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

Merge pull request #605 from atgeirr/fix_mdu

Browse Source 
Fixes for tof and tracer solvers
This commit is contained in:
Atgeirr Flø Rasmussen 2014-07-03 16:03:40 +02:00
commit 81c47c07db
4 changed files with 315 additions and 227 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

252
opm/core/tof/TofDiscGalReorder.cpp
View File

@ -59,6 +59,8 @@ namespace Opm
basis_func_.reset(new DGBasisBoundedTotalDegree(grid_, dg_degree));
}
tracers_ensure_unity_ = param.getDefault("tracers_ensure_unity", true);
use_cvi_ = param.getDefault("use_cvi", use_cvi_);
use_limiter_ = param.getDefault("use_limiter", use_limiter_);
if (use_limiter_) {
@ -273,13 +275,66 @@ namespace Opm
// right-hand-sides, first for tof and then (optionally) for
// all tracers.
const int dim = grid_.dimensions;
const int num_basis = basis_func_->numBasisFunc();
++num_singlesolves_;
std::fill(rhs_.begin(), rhs_.end(), 0.0);
std::fill(jac_.begin(), jac_.end(), 0.0);
// Add cell contributions to res_ and jac_.
cellContribs(cell);
// Add face contributions to res_ and jac_.
faceContribs(cell);
// Solve linear equation.
solveLinearSystem(cell);
// The solution ends up in rhs_, so we must copy it.
std::copy(rhs_.begin(), rhs_.begin() + num_basis, tof_coeff_ + num_basis*cell);
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
std::copy(rhs_.begin() + num_basis, rhs_.end(), tracer_coeff_ + num_tracers_*num_basis*cell);
}
// Apply limiter.
if (basis_func_->degree() > 0 && use_limiter_ && limiter_usage_ == DuringComputations) {
applyLimiter(cell, tof_coeff_);
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
for (int tr = 0; tr < num_tracers_; ++tr) {
applyTracerLimiter(cell, tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis);
}
}
}
// Ensure that tracer averages sum to 1.
if (num_tracers_ && tracers_ensure_unity_ && tracerhead_by_cell_[cell] == NoTracerHead) {
std::vector<double> tr_aver(num_tracers_);
double tr_sum = 0.0;
for (int tr = 0; tr < num_tracers_; ++tr) {
const double* local_basis = tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis;
tr_aver[tr] = basis_func_->functionAverage(local_basis);
tr_sum += tr_aver[tr];
}
if (tr_sum == 0.0) {
std::cout << "Tracer sum is zero in cell " << cell << std::endl;
} else {
for (int tr = 0; tr < num_tracers_; ++tr) {
const double increment = tr_aver[tr]/tr_sum - tr_aver[tr];
double* local_basis = tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis;
basis_func_->addConstant(increment, local_basis);
}
}
}
}
void TofDiscGalReorder::cellContribs(const int cell)
{
const int num_basis = basis_func_->numBasisFunc();
const int dim = grid_.dimensions;
// Compute cell residual contribution.
{
const int deg_needed = basis_func_->degree();
@ -296,7 +351,70 @@ namespace Opm
}
}
// Compute upstream residual contribution.
// Compute cell jacobian contribution. We use Fortran ordering
// for jac_, i.e. rows cycling fastest.
{
// Even with ECVI velocity interpolation, degree of precision 1
// is sufficient for optimal convergence order for DG1 when we
// use linear (total degree 1) basis functions.
// With bi(tri)-linear basis functions, it still seems sufficient
// for convergence order 2, but the solution looks much better and
// has significantly lower error with degree of precision 2.
// For now, we err on the side of caution, and use 2*degree, even
// though this is wasteful for the pure linear basis functions.
// const int deg_needed = 2*basis_func_->degree() - 1;
const int deg_needed = 2*basis_func_->degree();
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]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
basis_func_->evalGrad(cell, &coord_[0], &grad_basis_[0]);
velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
for (int dd = 0; dd < dim; ++dd) {
jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd];
}
}
}
}
}
// Compute downstream jacobian contribution from sink terms.
// Contribution from inflow sources would be
// similar to the contribution from upstream faces, but
// it is zero since we let all external inflow be associated
// with a zero tof.
if (source_[cell] < 0.0) {
// A sink.
const double flux = -source_[cell]; // Sign convention for flux: outflux > 0.
const double flux_density = flux / grid_.cell_volumes[cell];
// Do quadrature over the cell to compute
// \int_{K} b_i flux b_j dx
CellQuadrature quad(grid_, cell, 2*basis_func_->degree());
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j];
}
}
}
}
}
void TofDiscGalReorder::faceContribs(const int cell)
{
const int num_basis = basis_func_->numBasisFunc();
// Compute upstream residual contribution from 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;
@ -352,37 +470,6 @@ namespace Opm
}
}
// Compute cell jacobian contribution. We use Fortran ordering
// for jac_, i.e. rows cycling fastest.
{
// Even with ECVI velocity interpolation, degree of precision 1
// is sufficient for optimal convergence order for DG1 when we
// use linear (total degree 1) basis functions.
// With bi(tri)-linear basis functions, it still seems sufficient
// for convergence order 2, but the solution looks much better and
// has significantly lower error with degree of precision 2.
// For now, we err on the side of caution, and use 2*degree, even
// though this is wasteful for the pure linear basis functions.
// const int deg_needed = 2*basis_func_->degree() - 1;
const int deg_needed = 2*basis_func_->degree();
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]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
basis_func_->evalGrad(cell, &coord_[0], &grad_basis_[0]);
velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
for (int dd = 0; dd < dim; ++dd) {
jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd];
}
}
}
}
}
// Compute downstream jacobian contribution from faces.
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
@ -412,33 +499,17 @@ namespace Opm
}
}
}
// Compute downstream jacobian contribution from sink terms.
// Contribution from inflow sources would be
// similar to the contribution from upstream faces, but
// it is zero since we let all external inflow be associated
// with a zero tof.
if (source_[cell] < 0.0) {
// A sink.
const double flux = -source_[cell]; // Sign convention for flux: outflux > 0.
const double flux_density = flux / grid_.cell_volumes[cell];
// Do quadrature over the cell to compute
// \int_{K} b_i flux b_j dx
CellQuadrature quad(grid_, cell, 2*basis_func_->degree());
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j];
}
}
}
}
// Solve linear equation.
MAT_SIZE_T n = num_basis;
// This function assumes that jac_ and rhs_ contain the
// linear system to be solved. They are stored in orig_jac_
// and orig_rhs_, then the system is solved via LAPACK,
// overwriting the input data (jac_ and rhs_).
void TofDiscGalReorder::solveLinearSystem(const int cell)
{
MAT_SIZE_T n = basis_func_->numBasisFunc();
int num_tracer_to_compute = num_tracers_;
if (num_tracers_) {
if (tracerhead_by_cell_[cell] != NoTracerHead) {
@ -446,9 +517,9 @@ namespace Opm
}
}
MAT_SIZE_T nrhs = 1 + num_tracer_to_compute;
MAT_SIZE_T lda = num_basis;
std::vector<MAT_SIZE_T> piv(num_basis);
MAT_SIZE_T ldb = num_basis;
MAT_SIZE_T lda = n;
std::vector<MAT_SIZE_T> piv(n);
MAT_SIZE_T ldb = n;
MAT_SIZE_T info = 0;
orig_jac_ = jac_;
orig_rhs_ = rhs_;
@ -469,65 +540,6 @@ namespace Opm
}
OPM_THROW(std::runtime_error, "Lapack error: " << info << " encountered in cell " << cell);
}
// The solution ends up in rhs_, so we must copy it.
std::copy(rhs_.begin(), rhs_.begin() + num_basis, tof_coeff_ + num_basis*cell);
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
std::copy(rhs_.begin() + num_basis, rhs_.end(), tracer_coeff_ + num_tracers_*num_basis*cell);
}
// Apply limiter.
if (basis_func_->degree() > 0 && use_limiter_ && limiter_usage_ == DuringComputations) {
#ifdef EXTRA_VERBOSE
std::cout << "Cell: " << cell << " ";
std::cout << "v = ";
for (int dd = 0; dd < dim; ++dd) {
std::cout << velocity_[dd] << ' ';
}
std::cout << " grad tau = ";
for (int dd = 0; dd < dim; ++dd) {
std::cout << tof_coeff_[num_basis*cell + dd + 1] << ' ';
}
const double prod = std::inner_product(velocity_.begin(), velocity_.end(),
tof_coeff_ + num_basis*cell + 1, 0.0);
const double vv = std::inner_product(velocity_.begin(), velocity_.end(),
velocity_.begin(), 0.0);
const double gg = std::inner_product(tof_coeff_ + num_basis*cell + 1,
tof_coeff_ + num_basis*cell + num_basis,
tof_coeff_ + num_basis*cell + 1, 0.0);
std::cout << " prod = " << std::inner_product(velocity_.begin(), velocity_.end(),
tof_coeff_ + num_basis*cell + 1, 0.0);
std::cout << " normalized = " << prod/std::sqrt(vv*gg);
std::cout << " angle = " << std::acos(prod/std::sqrt(vv*gg))*360.0/(2.0*M_PI);
std::cout << std::endl;
#endif
applyLimiter(cell, tof_coeff_);
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
for (int tr = 0; tr < num_tracers_; ++tr) {
applyTracerLimiter(cell, tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis);
}
}
}
// Ensure that tracer averages sum to 1.
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
std::vector<double> tr_aver(num_tracers_);
double tr_sum = 0.0;
for (int tr = 0; tr < num_tracers_; ++tr) {
const double* local_basis = tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis;
tr_aver[tr] = basis_func_->functionAverage(local_basis);
tr_sum += tr_aver[tr];
}
if (tr_sum == 0.0) {
std::cout << "Tracer sum is zero in cell " << cell << std::endl;
} else {
for (int tr = 0; tr < num_tracers_; ++tr) {
const double increment = tr_aver[tr]/tr_sum - tr_aver[tr];
double* local_basis = tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis;
basis_func_->addConstant(increment, local_basis);
}
}
}
}

5
opm/core/tof/TofDiscGalReorder.hpp
View File

@ -121,6 +121,10 @@ namespace Opm
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
void cellContribs(const int cell);
void faceContribs(const int cell);
void solveLinearSystem(const int cell);
private:
// Disable copying and assignment.
TofDiscGalReorder(const TofDiscGalReorder&);
@ -146,6 +150,7 @@ namespace Opm
int num_tracers_;
enum { NoTracerHead = -1 };
std::vector<int> tracerhead_by_cell_;
bool tracers_ensure_unity_;
// Used by solveSingleCell().
std::vector<double> rhs_; // single-cell right-hand-sides
std::vector<double> jac_; // single-cell jacobian

261
opm/core/tof/TofReorder.cpp
View File

@ -42,8 +42,6 @@ namespace Opm
porevolume_(0),
source_(0),
tof_(0),
tracer_(0),
num_tracers_(0),
gauss_seidel_tol_(1e-3),
use_multidim_upwind_(use_multidim_upwind)
{
@ -71,7 +69,8 @@ namespace Opm
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src);
// OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src);
OPM_MESSAGE("Warning: sources do not sum to zero: " << cum_src);
}
#endif
tof.resize(grid_.number_of_cells);
@ -80,17 +79,11 @@ namespace Opm
if (use_multidim_upwind_) {
face_tof_.resize(grid_.number_of_faces);
std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
face_part_tof_.resize(grid_.face_nodepos[grid_.number_of_faces]);
std::fill(face_part_tof_.begin(), face_part_tof_.end(), 0.0);
}
num_tracers_ = 0;
num_multicell_ = 0;
max_size_multicell_ = 0;
max_iter_multicell_ = 0;
reorderAndTransport(grid_, darcyflux);
if (num_multicell_ > 0) {
std::cout << num_multicell_ << " multicell blocks with max size "
<< max_size_multicell_ << " cells in upto "
<< max_iter_multicell_ << " iterations." << std::endl;
}
compute_tracer_ = false;
executeSolve();
}
@ -117,43 +110,72 @@ namespace Opm
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
const int num_cells = grid_.number_of_cells;
#ifndef NDEBUG
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
const double cum_src = std::accumulate(source, source + num_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + num_cells)*1e-2) {
OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src);
}
#endif
tof.resize(grid_.number_of_cells);
tof.resize(num_cells);
std::fill(tof.begin(), tof.end(), 0.0);
tof_ = &tof[0];
// Find the tracer heads (injectors).
num_tracers_ = tracerheads.size();
tracer.resize(grid_.number_of_cells*num_tracers_);
std::fill(tracer.begin(), tracer.end(), 0.0);
if (num_tracers_ > 0) {
tracerhead_by_cell_.clear();
tracerhead_by_cell_.resize(grid_.number_of_cells, NoTracerHead);
if (use_multidim_upwind_) {
face_tof_.resize(grid_.number_of_faces);
std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
face_part_tof_.resize(grid_.face_nodepos[grid_.number_of_faces]);
std::fill(face_part_tof_.begin(), face_part_tof_.end(), 0.0);
}
for (int tr = 0; tr < num_tracers_; ++tr) {
// Execute solve for tof
compute_tracer_ = false;
executeSolve();
// Find the tracer heads (injectors).
const int num_tracers = tracerheads.size();
tracer.resize(num_cells*num_tracers);
std::fill(tracer.begin(), tracer.end(), 0.0);
if (num_tracers > 0) {
tracerhead_by_cell_.clear();
tracerhead_by_cell_.resize(num_cells, NoTracerHead);
}
for (int tr = 0; tr < num_tracers; ++tr) {
for (int i = 0; i < tracerheads[tr].size(); ++i) {
const int cell = tracerheads[tr][i];
tracer[cell*num_tracers_ + tr] = 1.0;
tracer[num_cells * tr + cell] = 1.0;
tracerhead_by_cell_[cell] = tr;
}
}
tracer_ = &tracer[0];
if (use_multidim_upwind_) {
face_tof_.resize(grid_.number_of_faces);
std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
OPM_THROW(std::runtime_error, "Multidimensional upwind not yet implemented for tracer.");
// Execute solve for tracers.
std::vector<double> fake_pv(num_cells, 0.0);
porevolume_ = fake_pv.data();
for (int tr = 0; tr < num_tracers; ++tr) {
tof_ = tracer.data() + tr * num_cells;
compute_tracer_ = true;
executeSolve();
}
// Write output tracer data (transposing the computed data).
std::vector<double> computed = tracer;
for (int cell = 0; cell < num_cells; ++cell) {
for (int tr = 0; tr < num_tracers; ++tr) {
tracer[num_tracers * cell + tr] = computed[num_cells * tr + cell];
}
}
}
void TofReorder::executeSolve()
{
num_multicell_ = 0;
max_size_multicell_ = 0;
max_iter_multicell_ = 0;
reorderAndTransport(grid_, darcyflux);
reorderAndTransport(grid_, darcyflux_);
if (num_multicell_ > 0) {
std::cout << num_multicell_ << " multicell blocks with max size "
<< max_size_multicell_ << " cells in upto "
@ -163,7 +185,6 @@ namespace Opm
void TofReorder::solveSingleCell(const int cell)
{
if (use_multidim_upwind_) {
@ -176,10 +197,9 @@ namespace Opm
// to the downwind_flux (note sign change resulting from
// different sign conventions: pos. source is injection,
// pos. flux is outflow).
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
for (int tr = 0; tr < num_tracers_; ++tr) {
tracer_[num_tracers_*cell + tr] = 0.0;
}
if (compute_tracer_ && tracerhead_by_cell_[cell] != NoTracerHead) {
// This is a tracer head cell, already has solution.
return;
}
double upwind_term = 0.0;
double downwind_flux = std::max(-source_[cell], 0.0);
@ -202,11 +222,6 @@ namespace Opm
// face.
if (other != -1) {
upwind_term += flux*tof_[other];
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
for (int tr = 0; tr < num_tracers_; ++tr) {
tracer_[num_tracers_*cell + tr] += flux*tracer_[num_tracers_*other + tr];
}
}
}
} else {
downwind_flux += flux;
@ -215,14 +230,6 @@ namespace Opm
// Compute tof.
tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux;
// Compute tracers (if any).
// Do not change tracer solution in source cells.
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
for (int tr = 0; tr < num_tracers_; ++tr) {
tracer_[num_tracers_*cell + tr] *= -1.0/downwind_flux;
}
}
}
@ -260,7 +267,11 @@ namespace Opm
}
// Compute tof for cell.
if (compute_tracer_ && tracerhead_by_cell_[cell] != NoTracerHead) {
// Do nothing to the value in this cell, since we are at a tracer head.
} else {
tof_[cell] = (porevolume_[cell] - upwind_term - downwind_term_face)/downwind_term_cell_factor;
}
// Compute tof for downwind faces.
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
@ -270,6 +281,20 @@ namespace Opm
double fterm, cterm_factor;
multidimUpwindTerms(f, cell, fterm, cterm_factor);
face_tof_[f] = fterm + cterm_factor*tof_[cell];
// Combine locally computed (for each adjacent vertex) terms, with uniform weighting.
const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[f];
const int* face_nodes_end = grid_.face_nodes + grid_.face_nodepos[f + 1];
const int num_terms = face_nodes_end - face_nodes_beg;
assert(num_terms == 2 || grid_.dimensions != 2);
for (const int* fn_iter = face_nodes_beg; fn_iter < face_nodes_end; ++fn_iter) {
double loc_face_term = 0.0;
double loc_cell_term_factor = 0.0;
const int node_pos = fn_iter - grid_.face_nodes;
localMultidimUpwindTerms(f, cell, node_pos,
loc_face_term, loc_cell_term_factor);
face_part_tof_[node_pos] = loc_face_term + loc_cell_term_factor * tof_[cell];
}
}
}
}
@ -303,7 +328,10 @@ namespace Opm
// Assumes that face_tof_[f] is known for all upstream faces f of upwind_cell.
// Assumes that face_part_tof_[node_pos] is known for all inflow
// faces to 'upwind_cell' sharing vertices with 'face'. The index
// 'node_pos' is the same as the one used for the grid face-node
// connectivity.
// Assumes that darcyflux_[face] is != 0.0.
// This function returns factors to compute the tof for 'face':
// tof(face) = face_term + cell_term_factor*tof(upwind_cell).
@ -314,71 +342,112 @@ namespace Opm
double& face_term,
double& cell_term_factor) const
{
// Implements multidim upwind according to
// Implements multidim upwind inspired by
// "Multidimensional upstream weighting for multiphase transport on general grids"
// by Keilegavlen, Kozdon, Mallison.
// However, that article does not give a 3d extension other than noting that using
// multidimensional upwinding in the XY-plane and not in the Z-direction may be
// a good idea. We have here attempted some generalization, by looking at all
// face-neighbours across edges as upwind candidates, and giving them all uniform weight.
// This will over-weight the immediate upstream cell value in an extruded 2d grid with
// one layer (top and bottom no-flow faces will enter the computation) compared to the
// original 2d case. Improvements are welcome.
// Note: Modified algorithm to consider faces that share even a single vertex with
// the input face. This reduces the problem of non-edge-conformal grids, but does not
// eliminate it entirely.
// a good idea. We have here attempted some generalization, by treating each face-part
// (association of a face and a vertex) as possibly influencing all downwind face-parts
// of the neighbouring cell that share the same vertex.
// The current implementation aims to reproduce 2d results for extruded 3d grids.
// Identify the adjacent faces of the upwind cell.
// Combine locally computed (for each adjacent vertex) terms, with uniform weighting.
const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[face];
const int* face_nodes_end = grid_.face_nodes + grid_.face_nodepos[face + 1];
assert(face_nodes_end - face_nodes_beg == 2 || grid_.dimensions != 2);
adj_faces_.clear();
const int num_terms = face_nodes_end - face_nodes_beg;
assert(num_terms == 2 || grid_.dimensions != 2);
face_term = 0.0;
cell_term_factor = 0.0;
for (const int* fn_iter = face_nodes_beg; fn_iter < face_nodes_end; ++fn_iter) {
double loc_face_term = 0.0;
double loc_cell_term_factor = 0.0;
localMultidimUpwindTerms(face, upwind_cell, fn_iter - grid_.face_nodes,
loc_face_term, loc_cell_term_factor);
face_term += loc_face_term;
cell_term_factor += loc_cell_term_factor;
}
face_term /= double(num_terms);
cell_term_factor /= double(num_terms);
}
namespace {
double weightFunc(const double w)
{
// SPU
// return 0.0;
// TMU
return w > 0.0 ? std::min(w, 1.0) : 0.0;
// SMU
// return w > 0.0 ? w/(1.0 + w) : 0.0;
}
}
void TofReorder::localMultidimUpwindTerms(const int face,
const int upwind_cell,
const int node_pos,
double& face_term,
double& cell_term_factor) const
{
// Loop over all faces adjacent to the given cell and the
// vertex in position node_pos.
// If that part's influx is positive, we store it, and also its associated
// node position.
std::vector<double> influx;
std::vector<int> node_pos_influx;
influx.reserve(5);
node_pos_influx.reserve(5);
const int node = grid_.face_nodes[node_pos];
for (int hf = grid_.cell_facepos[upwind_cell]; hf < grid_.cell_facepos[upwind_cell + 1]; ++hf) {
const int f = grid_.cell_faces[hf];
if (f != face) {
// Find out if the face 'f' is adjacent to vertex 'node'.
const int* f_nodes_beg = grid_.face_nodes + grid_.face_nodepos[f];
const int* f_nodes_end = grid_.face_nodes + grid_.face_nodepos[f + 1];
// Find out how many vertices they have in common.
// Using simple linear searches since sets are small.
int num_common = 0;
for (const int* f_iter = f_nodes_beg; f_iter < f_nodes_end; ++f_iter) {
num_common += std::count(face_nodes_beg, face_nodes_end, *f_iter);
const int* pos = std::find(f_nodes_beg, f_nodes_end, node);
const int node_pos2 = pos - grid_.face_nodes;
const bool is_adj = (pos != f_nodes_end);
if (is_adj) {
const int num_parts = f_nodes_end - f_nodes_beg;
const double influx_sign = (grid_.face_cells[2*f] == upwind_cell) ? -1.0 : 1.0;
const double part_influx = influx_sign * darcyflux_[f] / double(num_parts);
if (part_influx > 0.0) {
influx.push_back(part_influx);
node_pos_influx.push_back(node_pos2);
}
// Before: neighbours over an edge (3d) or vertex (2d).
// Now: neighbours across a vertex.
// if (num_common == grid_.dimensions - 1) {
if (num_common > 0) {
adj_faces_.push_back(f);
}
}
}
// Indentify adjacent faces with inflows, compute omega_star, omega,
// add up contributions.
const int num_adj = adj_faces_.size();
// The assertion below only holds if the grid is edge-conformal.
// No longer testing, since method no longer requires it.
// assert(num_adj == face_nodes_end - face_nodes_beg);
const double flux_face = std::fabs(darcyflux_[face]);
// Now we may compute the weighting of the upwind terms.
const int num_parts = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
const double outflux_sign = (grid_.face_cells[2*face] == upwind_cell) ? 1.0 : -1.0;
const double part_outflux = outflux_sign * darcyflux_[face] / double(num_parts);
const double sum_influx = std::accumulate(influx.begin(), influx.end(), 0.0);
const double w_factor = weightFunc(sum_influx / part_outflux);
const int num_influx = influx.size();
std::vector<double> w(num_influx);
face_term = 0.0;
cell_term_factor = 0.0;
for (int ii = 0; ii < num_adj; ++ii) {
const int f = adj_faces_[ii];
const double influx_f = (grid_.face_cells[2*f] == upwind_cell) ? -darcyflux_[f] : darcyflux_[f];
const double omega_star = influx_f/flux_face;
// SPU
// const double omega = 0.0;
// TMU
// const double omega = omega_star > 0.0 ? std::min(omega_star, 1.0) : 0.0;
// SMU
const double omega = omega_star > 0.0 ? omega_star/(1.0 + omega_star) : 0.0;
face_term += omega*face_tof_[f];
cell_term_factor += (1.0 - omega);
for (int ii = 0; ii < num_influx; ++ii) {
w[ii] = (influx[ii] / sum_influx) * w_factor;
face_term += w[ii] * face_part_tof_[node_pos_influx[ii]];
}
face_term /= double(num_adj);
cell_term_factor /= double(num_adj);
const double sum_w = std::accumulate(w.begin(), w.end(), 0.0);
cell_term_factor = 1.0 - sum_w;
const double tol = 1e-5;
if (cell_term_factor < -tol && cell_term_factor > 1.0 + tol) {
OPM_THROW(std::logic_error, "cell_term_factor outside [0,1]: " << cell_term_factor);
}
cell_term_factor = std::min(std::max(cell_term_factor, 0.0), 1.0);
assert(cell_term_factor >= 0.0);
assert(cell_term_factor <= 1.0);
}
} // namespace Opm

8
opm/core/tof/TofReorder.hpp
View File

@ -80,6 +80,7 @@ namespace Opm
std::vector<double>& tracer);
private:
void executeSolve();
virtual void solveSingleCell(const int cell);
void solveSingleCellMultidimUpwind(const int cell);
void assembleSingleCell(const int cell,
@ -90,6 +91,8 @@ namespace Opm
void multidimUpwindTerms(const int face, const int upwind_cell,
double& face_term, double& cell_term_factor) const;
void localMultidimUpwindTerms(const int face, const int upwind_cell, const int node_pos,
double& face_term, double& cell_term_factor) const;
private:
const UnstructuredGrid& grid_;
@ -97,8 +100,7 @@ namespace Opm
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
double* tof_;
double* tracer_;
int num_tracers_;
bool compute_tracer_;
enum { NoTracerHead = -1 };
std::vector<int> tracerhead_by_cell_;
// For solveMultiCell():
@ -109,7 +111,7 @@ namespace Opm
// For multidim upwinding:
bool use_multidim_upwind_;
std::vector<double> face_tof_; // For multidim upwind face tofs.
mutable std::vector<int> adj_faces_; // For multidim upwind logic.
std::vector<double> face_part_tof_; // For multidim upwind face tofs.
};
} // namespace Opm