OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/master' into add-phaseusage-to-interface

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2013-05-15 10:10:18 +02:00
commit 8c92a47b89
6 changed files with 410 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
examples/compute_tof_from_files.cpp
View File

@ -29,6 +29,7 @@
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/SparseTable.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
@ -120,6 +121,23 @@ main(int argc, char** argv)
}
}
const bool compute_tracer = param.getDefault("compute_tracer", false);
Opm::SparseTable<int> tracerheads;
if (compute_tracer) {
std::ifstream tr_stream(param.get<std::string>("tracerheads_filename").c_str());
int num_rows;
tr_stream >> num_rows;
for (int row = 0; row < num_rows; ++row) {
int row_size;
tr_stream >> row_size;
std::vector<int> rowdata(row_size);
for (int elem = 0; elem < row_size; ++elem) {
tr_stream >> rowdata[elem];
}
tracerheads.appendRow(rowdata.begin(), rowdata.end());
}
}
// Choice of tof solver.
bool use_dg = param.getDefault("use_dg", false);
bool use_multidim_upwind = false;
@ -130,10 +148,6 @@ main(int argc, char** argv)
} 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);
@ -161,11 +175,15 @@ main(int argc, char** argv)
std::vector<double> tof;
std::vector<double> tracer;
if (use_dg) {
dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof);
if (compute_tracer) {
dg_solver->solveTofTracer(&flux[0], &porevol[0], &src[0], tracerheads, tof, tracer);
} else {
dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof);
}
} else {
Opm::TofReorder tofsolver(grid, use_multidim_upwind);
if (compute_tracer) {
tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tof, tracer);
tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tracerheads, tof, tracer);
} else {
tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof);
}

16
opm/core/simulator/initState_impl.hpp
View File

@ -337,7 +337,7 @@ namespace Opm
}
}
state.setFirstSat(left_cells, props, State::MaxSat);
const double init_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
std::fill(state.pressure().begin(), state.pressure().end(), init_p);
} else if (segregation_testcase) {
// Warn against error-prone usage.
@ -351,7 +351,7 @@ namespace Opm
const double woc = param.get<double>("water_oil_contact");
initWaterOilContact(grid, props, woc, WaterAbove, state);
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
double dens[2] = { props.density()[1], props.density()[0] };
initHydrostaticPressure(grid, dens, woc, gravity, woc, ref_p, state);
} else if (param.has("water_oil_contact")) {
@ -366,7 +366,7 @@ namespace Opm
const double woc = param.get<double>("water_oil_contact");
initWaterOilContact(grid, props, woc, WaterBelow, state);
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
initHydrostaticPressure(grid, props.density(), woc, gravity, woc, ref_p, state);
} else if (param.has("init_saturation")) {
// Initialise water saturation to init_saturation parameter.
@ -376,7 +376,7 @@ namespace Opm
state.saturation()[2*cell + 1] = 1.0 - init_saturation;
}
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
const double rho = props.density()[0]*init_saturation + props.density()[1]*(1.0 - init_saturation);
const double dens[2] = { rho, rho };
const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
@ -384,7 +384,7 @@ namespace Opm
} else {
// Use default: water saturation is minimum everywhere.
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
const double rho = props.density()[1];
const double dens[2] = { rho, rho };
const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
@ -432,7 +432,7 @@ namespace Opm
}
}
state.setFirstSat(left_cells, props, State::MaxSat);
const double init_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
std::fill(state.pressure().begin(), state.pressure().end(), init_p);
} else if (param.has("water_oil_contact")) {
// Warn against error-prone usage.
@ -446,12 +446,12 @@ namespace Opm
const double woc = param.get<double>("water_oil_contact");
initWaterOilContact(grid, props, woc, WaterBelow, state);
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
initHydrostaticPressure(grid, props, woc, gravity, woc, ref_p, state);
} else {
// Use default: water saturation is minimum everywhere.
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
const double woc = -1e100;
initHydrostaticPressure(grid, props, woc, gravity, ref_z, ref_p, state);

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

@ -24,6 +24,7 @@
#include <opm/core/tof/DGBasis.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/SparseTable.hpp>
#include <opm/core/utility/VelocityInterpolation.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/linalg/blas_lapack.h>
@ -45,7 +46,8 @@ namespace Opm
limiter_method_(MinUpwindAverage),
limiter_usage_(DuringComputations),
coord_(grid.dimensions),
velocity_(grid.dimensions)
velocity_(grid.dimensions),
gauss_seidel_tol_(1e-3)
{
const int dg_degree = param.getDefault("dg_degree", 0);
const bool use_tensorial_basis = param.getDefault("use_tensorial_basis", false);
@ -97,9 +99,9 @@ namespace Opm
/// Solve for time-of-flight.
void TofDiscGalReorder::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof_coeff)
const double* porevolume,
const double* source,
std::vector<double>& tof_coeff)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
@ -123,6 +125,11 @@ namespace Opm
basis_nb_.resize(num_basis);
grad_basis_.resize(num_basis*grid_.dimensions);
velocity_interpolation_->setupFluxes(darcyflux);
num_tracers_ = 0;
num_multicell_ = 0;
max_size_multicell_ = 0;
max_iter_multicell_ = 0;
num_singlesolves_ = 0;
reorderAndTransport(grid_, darcyflux);
switch (limiter_usage_) {
case AsPostProcess:
@ -137,6 +144,109 @@ namespace Opm
default:
THROW("Unknown limiter usage choice: " << limiter_usage_);
}
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;
std::cout << "Average solves per cell (for all cells) was "
<< double(num_singlesolves_)/double(grid_.number_of_cells) << std::endl;
}
}
/// Solve for time-of-flight and a number of tracers.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] tracerheads Table containing one row per tracer, and each
/// row contains the source cells for that tracer.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
/// \param[out] tracer_coeff Array of tracer solution coefficients. N*K per cell,
/// where N is equal to tracerheads.size(). All K coefs
/// for a tracer are consecutive, and all tracers' coefs
/// for a cell come before those for the next cell.
void TofDiscGalReorder::solveTofTracer(const double* darcyflux,
const double* porevolume,
const double* source,
const SparseTable<int>& tracerheads,
std::vector<double>& tof_coeff,
std::vector<double>& tracer_coeff)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
#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) {
// THROW("Sources do not sum to zero: " << cum_src);
MESSAGE("Warning: sources do not sum to zero: " << cum_src);
}
#endif
const int num_basis = basis_func_->numBasisFunc();
num_tracers_ = tracerheads.size();
tof_coeff.resize(num_basis*grid_.number_of_cells);
std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0);
tof_coeff_ = &tof_coeff[0];
rhs_.resize(num_basis*(num_tracers_ + 1));
jac_.resize(num_basis*num_basis);
orig_jac_.resize(num_basis*num_basis);
basis_.resize(num_basis);
basis_nb_.resize(num_basis);
grad_basis_.resize(num_basis*grid_.dimensions);
velocity_interpolation_->setupFluxes(darcyflux);
// Set up tracer
tracer_coeff.resize(grid_.number_of_cells*num_tracers_*num_basis);
std::fill(tracer_coeff.begin(), tracer_coeff.end(), 0.0);
if (num_tracers_ > 0) {
tracerhead_by_cell_.clear();
tracerhead_by_cell_.resize(grid_.number_of_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];
basis_func_->addConstant(1.0, &tracer_coeff[cell*num_tracers_*num_basis + tr*num_basis]);
tracer_coeff[cell*num_tracers_ + tr] = 1.0;
tracerhead_by_cell_[cell] = tr;
}
}
tracer_coeff_ = &tracer_coeff[0];
num_multicell_ = 0;
max_size_multicell_ = 0;
max_iter_multicell_ = 0;
num_singlesolves_ = 0;
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_);
}
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;
std::cout << "Average solves per cell (for all cells) was "
<< double(num_singlesolves_)/double(grid_.number_of_cells) << std::endl;
}
}
@ -153,9 +263,17 @@ namespace Opm
// This is linear in c_i, so we do not need any nonlinear iterations.
// We assemble the jacobian and the right-hand side. The residual is
// equal to Res = Jac*c - rhs, and we compute rhs directly.
//
// For tracers, the equation is the same, except for the last
// term being zero (the one with \phi).
//
// The rhs_ vector contains a (Fortran ordering) matrix of all
// 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);
@ -170,6 +288,7 @@ namespace Opm
basis_func_->eval(cell, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
// Only adding to the tof rhs.
rhs_[j] += w * basis_[j] * porevolume_[cell] / grid_.cell_volumes[cell];
}
}
@ -193,6 +312,8 @@ namespace Opm
}
if (upstream_cell < 0) {
// This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
// For tracers, a cell with inflow should be marked as a tracer head cell,
// and not be modified.
continue;
}
// Do quadrature over the face to compute
@ -209,12 +330,23 @@ namespace Opm
quad.quadPtCoord(quad_pt, &coord_[0]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
basis_func_->eval(upstream_cell, &coord_[0], &basis_nb_[0]);
const double w = quad.quadPtWeight(quad_pt);
// Modify tof rhs
const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(),
tof_coeff_ + num_basis*upstream_cell, 0.0);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j];
}
// Modify tracer rhs
if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
for (int tr = 0; tr < num_tracers_; ++tr) {
const double* up_tr_co = tracer_coeff_ + num_tracers_*num_basis*upstream_cell + num_basis*tr;
const double tracer_up = std::inner_product(basis_nb_.begin(), basis_nb_.end(), up_tr_co, 0.0);
for (int j = 0; j < num_basis; ++j) {
rhs_[num_basis*(tr + 1) + j] -= w * tracer_up * normal_velocity * basis_[j];
}
}
}
}
}
@ -305,7 +437,13 @@ namespace Opm
// Solve linear equation.
MAT_SIZE_T n = num_basis;
MAT_SIZE_T nrhs = 1;
int num_tracer_to_compute = num_tracers_;
if (num_tracers_) {
if (tracerhead_by_cell_[cell] != NoTracerHead) {
num_tracer_to_compute = 0;
}
}
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;
@ -331,7 +469,10 @@ namespace Opm
}
// The solution ends up in rhs_, so we must copy it.
std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell);
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) {
@ -359,6 +500,31 @@ namespace Opm
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);
}
}
}
}
@ -367,10 +533,27 @@ namespace Opm
void TofDiscGalReorder::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
++num_multicell_;
max_size_multicell_ = std::max(max_size_multicell_, num_cells);
// std::cout << "Multiblock solve with " << num_cells << " cells." << std::endl;
// Using a Gauss-Seidel approach.
const int nb = basis_func_->numBasisFunc();
double max_delta = 1e100;
int num_iter = 0;
while (max_delta > gauss_seidel_tol_) {
max_delta = 0.0;
++num_iter;
for (int ci = 0; ci < num_cells; ++ci) {
const int cell = cells[ci];
const double tof_before = basis_func_->functionAverage(&tof_coeff_[nb*cell]);
solveSingleCell(cell);
const double tof_after = basis_func_->functionAverage(&tof_coeff_[nb*cell]);
max_delta = std::max(max_delta, std::fabs(tof_after - tof_before));
}
// std::cout << "Max delta = " << max_delta << std::endl;
}
max_iter_multicell_ = std::max(max_iter_multicell_, num_iter);
}
@ -462,7 +645,7 @@ namespace Opm
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;
// std::cout << "Trouble in cell " << cell << std::endl;
limiter = 0.0;
basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell);
}
@ -562,5 +745,47 @@ namespace Opm
void TofDiscGalReorder::applyTracerLimiter(const int cell, double* local_coeff)
{
// Evaluate the solution in all corners of all faces. Extract max and min.
const int dim = grid_.dimensions;
const int num_basis = basis_func_->numBasisFunc();
double min_cornerval = 1e100;
double max_cornerval = -1e100;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
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 tracer_corner = std::inner_product(basis_.begin(), basis_.end(),
local_coeff, 0.0);
min_cornerval = std::min(min_cornerval, tracer_corner);
max_cornerval = std::max(min_cornerval, tracer_corner);
}
}
const double average = basis_func_->functionAverage(local_coeff);
if (average < 0.0 || average > 1.0) {
// Adjust average. Flatten gradient.
std::fill(local_coeff, local_coeff + num_basis, 0.0);
if (average > 1.0) {
basis_func_->addConstant(1.0, local_coeff);
}
} else {
// Possibly adjust gradient.
double factor = 1.0;
if (min_cornerval < 0.0) {
factor = average/(average - min_cornerval);
}
if (max_cornerval > 1.0) {
factor = std::min(factor, (1.0 - average)/(max_cornerval - average));
}
if (factor != 1.0) {
basis_func_->multiplyGradient(factor, local_coeff);
}
}
}
} // namespace Opm

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

@ -35,6 +35,7 @@ namespace Opm
class VelocityInterpolationInterface;
class DGBasisInterface;
namespace parameter { class ParameterGroup; }
template <typename T> class SparseTable;
/// Implements a discontinuous Galerkin solver for
/// (single-phase) time-of-flight using reordering.
@ -83,7 +84,7 @@ namespace Opm
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// cell come before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
void solveTof(const double* darcyflux,
@ -91,6 +92,31 @@ namespace Opm
const double* source,
std::vector<double>& tof_coeff);
/// Solve for time-of-flight and a number of tracers.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] tracerheads Table containing one row per tracer, and each
/// row contains the source cells for that tracer.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
/// \param[out] tracer_coeff Array of tracer solution coefficients. N*K per cell,
/// where N is equal to tracerheads.size(). All K coefs
/// for a tracer are consecutive, and all tracers' coefs
/// for a cell come before those for the next cell.
void solveTofTracer(const double* darcyflux,
const double* porevolume,
const double* source,
const SparseTable<int>& tracerheads,
std::vector<double>& tof_coeff,
std::vector<double>& tracer_coeff);
private:
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
@ -115,16 +141,27 @@ namespace Opm
const double* source_; // one volumetric source term per cell
boost::shared_ptr<DGBasisInterface> basis_func_;
double* tof_coeff_;
std::vector<double> rhs_; // single-cell right-hand-side
// For tracers.
double* tracer_coeff_;
int num_tracers_;
enum { NoTracerHead = -1 };
std::vector<int> tracerhead_by_cell_;
// Used by solveSingleCell().
std::vector<double> rhs_; // single-cell right-hand-sides
std::vector<double> jac_; // single-cell jacobian
std::vector<double> orig_rhs_; // single-cell right-hand-side (copy)
std::vector<double> orig_rhs_; // single-cell right-hand-sides (copy)
std::vector<double> orig_jac_; // single-cell jacobian (copy)
// Below: storage for quantities needed by solveSingleCell().
std::vector<double> coord_;
mutable std::vector<double> basis_;
mutable std::vector<double> basis_nb_;
std::vector<double> grad_basis_;
std::vector<double> velocity_;
int num_singlesolves_;
// Used by solveMultiCell():
double gauss_seidel_tol_;
int num_multicell_;
int max_size_multicell_;
int max_iter_multicell_;
// Private methods
@ -137,6 +174,10 @@ namespace Opm
void applyLimiterAsSimultaneousPostProcess();
double totalFlux(const int cell) const;
double minCornerVal(const int cell, const int face) const;
// Apply a simple (restrict to [0,1]) limiter.
// Intended for tracers.
void applyTracerLimiter(const int cell, double* local_coeff);
};
} // namespace Opm

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

@ -21,6 +21,7 @@
#include <opm/core/tof/TofReorder.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/SparseTable.hpp>
#include <algorithm>
#include <numeric>
#include <cmath>
@ -41,6 +42,7 @@ namespace Opm
tof_(0),
tracer_(0),
num_tracers_(0),
gauss_seidel_tol_(1e-3),
use_multidim_upwind_(use_multidim_upwind)
{
}
@ -56,9 +58,9 @@ namespace Opm
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values.
void TofReorder::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof)
const double* porevolume,
const double* source,
std::vector<double>& tof)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
@ -78,26 +80,35 @@ namespace Opm
std::fill(face_tof_.begin(), face_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;
}
}
/// Solve for time-of-flight and a number of tracers.
/// One tracer will be used for each inflow flux specified in
/// the source parameter.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] tracerheads Table containing one row per tracer, and each
/// row contains the source cells for that tracer.
/// \param[out] tof Array of time-of-flight values (1 per cell).
/// \param[out] tracer Array of tracer values (N per cell, where N is
/// the number of cells c for which source[c] > 0.0).
/// \param[out] tracer Array of tracer values. N per cell, where N is
/// equalt to tracerheads.size().
void TofReorder::solveTofTracer(const double* darcyflux,
const double* porevolume,
const double* source,
const SparseTable<int>& tracerheads,
std::vector<double>& tof,
std::vector<double>& tracer)
{
@ -114,26 +125,38 @@ namespace Opm
tof.resize(grid_.number_of_cells);
std::fill(tof.begin(), tof.end(), 0.0);
tof_ = &tof[0];
// Find the tracer heads (injectors).
std::vector<int> tracerheads;
for (int c = 0; c < grid_.number_of_cells; ++c) {
if (source[c] > 0.0) {
tracerheads.push_back(c);
}
}
num_tracers_ = tracerheads.size();
tracer.resize(grid_.number_of_cells*num_tracers_);
std::fill(tracer.begin(), tracer.end(), 0.0);
for (int tr = 0; tr < num_tracers_; ++tr) {
tracer[tracerheads[tr]*num_tracers_ + tr] = 1.0;
if (num_tracers_ > 0) {
tracerhead_by_cell_.clear();
tracerhead_by_cell_.resize(grid_.number_of_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;
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);
THROW("Multidimensional upwind not yet implemented for tracer.");
}
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;
}
}
@ -151,6 +174,11 @@ 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;
}
}
double upwind_term = 0.0;
double downwind_flux = std::max(-source_[cell], 0.0);
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
@ -172,8 +200,10 @@ namespace Opm
// face.
if (other != -1) {
upwind_term += flux*tof_[other];
for (int tr = 0; tr < num_tracers_; ++tr) {
tracer_[num_tracers_*cell + tr] += flux*tracer_[num_tracers_*other + tr];
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 {
@ -186,7 +216,7 @@ namespace Opm
// Compute tracers (if any).
// Do not change tracer solution in source cells.
if (source_[cell] <= 0.0) {
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;
}
@ -219,7 +249,7 @@ namespace Opm
// Add flux to upwind_term or downwind_term_[face|cell_factor].
if (flux < 0.0) {
upwind_term += flux*face_tof_[f];
} else {
} else if (flux > 0.0) {
double fterm, cterm_factor;
multidimUpwindTerms(f, cell, fterm, cterm_factor);
downwind_term_face += fterm*flux;
@ -228,7 +258,7 @@ namespace Opm
}
// Compute tof for cell.
tof_[cell] = (porevolume_[cell] - upwind_term - downwind_term_face)/downwind_term_cell_factor; // }
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) {
@ -247,10 +277,25 @@ namespace Opm
void TofReorder::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
++num_multicell_;
max_size_multicell_ = std::max(max_size_multicell_, num_cells);
// std::cout << "Multiblock solve with " << num_cells << " cells." << std::endl;
// Using a Gauss-Seidel approach.
double max_delta = 1e100;
int num_iter = 0;
while (max_delta > gauss_seidel_tol_) {
max_delta = 0.0;
++num_iter;
for (int ci = 0; ci < num_cells; ++ci) {
const int cell = cells[ci];
const double tof_before = tof_[cell];
solveSingleCell(cell);
max_delta = std::max(max_delta, std::fabs(tof_[cell] - tof_before));
}
// std::cout << "Max delta = " << max_delta << std::endl;
}
max_iter_multicell_ = std::max(max_iter_multicell_, num_iter);
}
@ -277,6 +322,9 @@ namespace Opm
// 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.
// Identify the adjacent faces of the upwind cell.
const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[face];
@ -294,11 +342,11 @@ namespace Opm
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);
}
if (num_common == grid_.dimensions - 1) {
// Neighbours over an edge (3d) or vertex (2d).
// 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);
} else {
ASSERT(num_common == 0);
}
}
}
@ -306,7 +354,9 @@ namespace Opm
// Indentify adjacent faces with inflows, compute omega_star, omega,
// add up contributions.
const int num_adj = adj_faces_.size();
ASSERT(num_adj == face_nodes_end - face_nodes_beg);
// 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]);
face_term = 0.0;
cell_term_factor = 0.0;

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

@ -24,12 +24,14 @@
#include <vector>
#include <map>
#include <ostream>
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
template <typename T> class SparseTable;
/// Implements a first-order finite volume solver for
/// (single-phase) time-of-flight using reordering.
@ -60,25 +62,30 @@ namespace Opm
std::vector<double>& tof);
/// Solve for time-of-flight and a number of tracers.
/// One tracer will be used for each inflow flux specified in
/// the source parameter.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] tracerheads Table containing one row per tracer, and each
/// row contains the source cells for that tracer.
/// \param[out] tof Array of time-of-flight values (1 per cell).
/// \param[out] tracer Array of tracer values (N per cell, where N is
/// the number of cells c for which source[c] > 0.0).
/// \param[out] tracer Array of tracer values. N per cell, where N is
/// equalt to tracerheads.size().
void solveTofTracer(const double* darcyflux,
const double* porevolume,
const double* source,
const SparseTable<int>& tracerheads,
std::vector<double>& tof,
std::vector<double>& tracer);
private:
virtual void solveSingleCell(const int cell);
void solveSingleCellMultidimUpwind(const int cell);
void assembleSingleCell(const int cell,
std::vector<int>& local_column,
std::vector<double>& local_coefficient,
double& rhs);
virtual void solveMultiCell(const int num_cells, const int* cells);
void multidimUpwindTerms(const int face, const int upwind_cell,
@ -92,6 +99,14 @@ namespace Opm
double* tof_;
double* tracer_;
int num_tracers_;
enum { NoTracerHead = -1 };
std::vector<int> tracerhead_by_cell_;
// For solveMultiCell():
double gauss_seidel_tol_;
int num_multicell_;
int max_size_multicell_;
int max_iter_multicell_;
// 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.