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

Merge pull request #15 from atgeirr/master

Browse Source 
Modifications to compressible transport solver, other minor mods.
This commit is contained in:
Bård Skaflestad 2012-08-17 03:49:30 -07:00
commit 21b33e40d7
7 changed files with 530 additions and 494 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
examples/sim_2p_incomp_reorder.cpp
View File

@ -47,12 +47,25 @@
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
} // anon namespace
// ----------------- Main program -----------------
@ -155,13 +168,6 @@ main(int argc, char** argv)
// Linear solver.
LinearSolverFactory linsolver(param);
//Warn if any parameters are unused.
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
if (output) {
@ -196,6 +202,7 @@ main(int argc, char** argv)
grav);
SimulatorTimer simtimer;
simtimer.init(param);
warnIfUnusedParams(param);
WellState well_state;
well_state.init(0, state);
rep = simulator.run(simtimer, state, well_state);
@ -229,7 +236,7 @@ main(int argc, char** argv)
<< "\n (number of steps: "
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
// Create new wells, well_satate
// Create new wells, well_state
WellsManager wells(*deck, *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
@ -248,6 +255,9 @@ main(int argc, char** argv)
bcs.c_bcs(),
linsolver,
grav);
if (epoch == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
// Update total timing report and remember step number.

12
examples/sim_wateroil.cpp
View File

@ -68,8 +68,6 @@
#include <vector>
#include <numeric>
#define TRANSPORT_SOLVER_FIXED 1
template <class State>
static void outputState(const UnstructuredGrid& grid,
@ -299,14 +297,13 @@ main(int argc, char** argv)
nl_press_res_tol, nl_press_change_tol, nl_press_maxiter,
grav, wells->c_wells());
// Reordering solver.
#if TRANSPORT_SOLVER_FIXED
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
if (use_gauss_seidel_gravity) {
reorder_model.initGravity(grav);
reorder_model.initGravity();
}
#endif // TRANSPORT_SOLVER_FIXED
// Column-based gravity segregation solver.
std::vector<std::vector<int> > columns;
if (use_column_solver) {
@ -414,7 +411,6 @@ main(int argc, char** argv)
// Solve transport.
transport_timer.start();
#if TRANSPORT_SOLVER_FIXED
double stepsize = simtimer.currentStepLength();
if (num_transport_substeps != 1) {
stepsize /= double(num_transport_substeps);
@ -427,11 +423,9 @@ main(int argc, char** argv)
&porevol[0], &porevol[0], &reorder_src[0], stepsize, state.saturation());
// Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
if (use_segregation_split) {
THROW("Segregation not implemented yet.");
// reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation());
reorder_model.solveGravity(columns, &state.pressure()[0], &porevol[0], stepsize, grav, state.saturation());
}
}
#endif // TRANSPORT_SOLVER_FIXED
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;

5
opm/core/simulator/SimulatorTwophase.cpp
View File

@ -44,6 +44,7 @@
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
@ -137,7 +138,7 @@ namespace Opm
// Write data in VTK format.
std::ostringstream vtkfilename;
vtkfilename << output_dir << "/vtk_files";
boost::filesystem::path fpath(vtkfilename.str().c_str());
boost::filesystem::path fpath(vtkfilename.str());
try {
create_directories(fpath);
}
@ -161,7 +162,7 @@ namespace Opm
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first;
boost::filesystem::path fpath(fname.str().c_str());
fpath = fname.str();
try {
create_directories(fpath);
}

53
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
View File

@ -254,6 +254,8 @@ namespace Opm
++update_count;
const int cell = cells[i];
const double old_s = saturation_[cell];
// solveSingleCell() requires saturation_[cell]
// to be s0.
saturation_[cell] = s0[i];
solveSingleCell(cell);
const double s_change = std::fabs(saturation_[cell] - old_s);
@ -305,8 +307,9 @@ namespace Opm
// Residual function r(s) for a single-cell implicit Euler gravity segregation
//
// r(s) = s - s0 + dt/pv*sum_{j adj i}( gravmod_ij * gf_ij ).
// [[ incompressible was: r(s) = s - s0 + dt/pv*sum_{j adj i}( gravmod_ij * gf_ij ) ]]
//
// r(s) = s - B*z0 + dt/pv*( influx + outflux*f(s) )
struct TransportModelCompressibleTwophase::GravityResidual
{
int cell;
@ -372,32 +375,52 @@ namespace Opm
void TransportModelCompressibleTwophase::initGravity(const double* grav)
void TransportModelCompressibleTwophase::initGravity()
{
// Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j)
// Set up transmissibilities.
std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
const int nf = grid_.number_of_faces;
const int dim = grid_.dimensions;
gravflux_.resize(nf);
trans_.resize(nf);
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &gravflux_[0]);
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &trans_[0]);
}
const double delta_rho = 0.0;// props_.density()[0] - props_.density()[1];
THROW("TransportModelCompressibleTwophase gravity solver not done yet."); // See line above...
void TransportModelCompressibleTwophase::initGravityDynamic(const double* grav)
{
// Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f)
// + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ]
// But b_w,i * rho_w,S = rho_w,i, which we conmpute with a call to props_.density().
// We assume that we already have stored T_ij in trans_.
// We also assume that the A_ matrices are updated from an earlier call to solve().
const int nc = grid_.number_of_cells;
const int nf = grid_.number_of_faces;
const int np = props_.numPhases();
ASSERT(np == 2);
const int dim = grid_.dimensions;
density_.resize(nc*np);
props_.density(grid_.number_of_cells, &A_[0], &density_[0]);
std::fill(gravflux_.begin(), gravflux_.end(), 0.0);
for (int f = 0; f < nf; ++f) {
const int* c = &grid_.face_cells[2*f];
double gdz = 0.0;
const double signs[2] = { 1.0, -1.0 };
if (c[0] != -1 && c[1] != -1) {
for (int d = 0; d < dim; ++d) {
gdz += grav[d]*(grid_.cell_centroids[dim*c[0] + d] - grid_.cell_centroids[dim*c[1] + d]);
for (int ci = 0; ci < 2; ++ci) {
double gdz = 0.0;
for (int d = 0; d < dim; ++d) {
gdz += grav[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]);
}
gravflux_[f] += signs[ci]*trans_[f]*gdz*(density_[2*c[ci]] - density_[2*c[ci] + 1]);
}
}
gravflux_[f] *= delta_rho*gdz;
}
}
void TransportModelCompressibleTwophase::solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
const double* gravflux)
@ -470,10 +493,13 @@ namespace Opm
void TransportModelCompressibleTwophase::solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure,
const double* porevolume0,
const double* porevolume,
const double dt,
const double* grav,
std::vector<double>& saturation)
{
// Assume that solve() has already been called, so that A_ is current.
initGravityDynamic(grav);
// Initialize mobilities.
const int nc = grid_.number_of_cells;
std::vector<int> cells(nc);
@ -493,7 +519,6 @@ namespace Opm
// Set up other variables.
porevolume0_ = porevolume0;
porevolume_ = porevolume;
dt_ = dt;
toWaterSat(saturation, saturation_);

96
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
View File

@ -33,91 +33,97 @@ namespace Opm
class TransportModelCompressibleTwophase : public TransportModelInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] tol Tolerance used in the solver.
/// \param[in] maxit Maximum number of non-linear iterations used.
TransportModelCompressibleTwophase(const UnstructuredGrid& grid,
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] tol Tolerance used in the solver.
/// \param[in] maxit Maximum number of non-linear iterations used.
TransportModelCompressibleTwophase(const UnstructuredGrid& grid,
const Opm::BlackoilPropertiesInterface& props,
const double tol,
const int maxit);
/// Solve for saturation at next timestep.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] pressure Array of cell pressures
/// \param[in] surfacevol0 Array of surface volumes at start of timestep
/// \param[in] porevolume0 Array of pore volumes at start of timestep.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
void solve(const double* darcyflux,
/// Solve for saturation at next timestep.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] pressure Array of cell pressures
/// \param[in] surfacevol0 Array of surface volumes at start of timestep
/// \param[in] porevolume0 Array of pore volumes at start of timestep.
/// \param[in] porevolume Array of pore volumes at end of timestep.
/// \param[in] source Transport source term.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
void solve(const double* darcyflux,
const double* pressure,
const double* surfacevol0,
const double* porevolume0,
const double* porevolume,
const double* source,
const double dt,
std::vector<double>& saturation);
const double* source,
const double dt,
std::vector<double>& saturation);
/// Initialise quantities needed by gravity solver.
/// \param[in] grav gravity vector
void initGravity(const double* grav);
void initGravity();
/// Solve for gravity segregation.
/// This uses a column-wise nonlinear Gauss-Seidel approach.
/// It assumes that the input columns contain cells in a single
/// vertical stack, that do not interact with other columns (for
/// gravity segregation.
/// \TODO: Implement this.
/// \param[in] columns Vector of cell-columns.
/// \param[in] porevolume0 Array of pore volumes at start of timestep.
/// \param[in] dt Time step.
/// \param[in] grav Gravity vector.
/// \param[in, out] saturation Phase saturations.
void solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure,
const double* porevolume0,
const double* porevolume,
const double dt,
const double* grav,
std::vector<double>& saturation);
private:
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
void solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
const double* gravflux);
int solveGravityColumn(const std::vector<int>& cells);
void initGravityDynamic(const double* grav);
private:
const UnstructuredGrid& grid_;
const BlackoilPropertiesInterface& props_;
const UnstructuredGrid& grid_;
const BlackoilPropertiesInterface& props_;
std::vector<int> allcells_;
std::vector<double> visc_;
std::vector<double> A_;
std::vector<double> smin_;
std::vector<double> smax_;
double tol_;
double maxit_;
std::vector<double> smin_;
std::vector<double> smax_;
double tol_;
double maxit_;
const double* darcyflux_; // one flux per grid face
const double* surfacevol0_; // one per phase per cell
const double* porevolume0_; // one volume per cell
const double* porevolume_; // one volume per cell
const double* source_; // one source per cell
double dt_;
const double* darcyflux_; // one flux per grid face
const double* surfacevol0_; // one per phase per cell
const double* porevolume0_; // one volume per cell
const double* porevolume_; // one volume per cell
const double* source_; // one source per cell
double dt_;
std::vector<double> saturation_; // P (= num. phases) per cell
std::vector<double> fractionalflow_; // = m[0]/(m[0] + m[1]) per cell
std::vector<double> fractionalflow_; // = m[0]/(m[0] + m[1]) per cell
// For gravity segregation.
std::vector<double> trans_;
std::vector<double> density_;
std::vector<double> gravflux_;
std::vector<double> mob_;
std::vector<double> s0_;
// Storing the upwind and downwind graphs for experiments.
std::vector<int> ia_upw_;
std::vector<int> ja_upw_;
std::vector<int> ia_downw_;
std::vector<int> ja_downw_;
// Storing the upwind and downwind graphs for experiments.
std::vector<int> ia_upw_;
std::vector<int> ja_upw_;
std::vector<int> ia_downw_;
std::vector<int> ja_downw_;
struct Residual;
double fracFlow(double s, int cell) const;
struct Residual;
double fracFlow(double s, int cell) const;
struct GravityResidual;
void mobility(double s, int cell, double* mob) const;

744
opm/core/transport/reorder/TransportModelTwophase.cpp
View File

@ -41,68 +41,68 @@ namespace Opm
TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const double tol,
const int maxit)
: grid_(grid),
props_(props),
tol_(tol),
maxit_(maxit),
darcyflux_(0),
source_(0),
dt_(0.0),
saturation_(grid.number_of_cells, -1.0),
fractionalflow_(grid.number_of_cells, -1.0),
mob_(2*grid.number_of_cells, -1.0)
const Opm::IncompPropertiesInterface& props,
const double tol,
const int maxit)
: grid_(grid),
props_(props),
tol_(tol),
maxit_(maxit),
darcyflux_(0),
source_(0),
dt_(0.0),
saturation_(grid.number_of_cells, -1.0),
fractionalflow_(grid.number_of_cells, -1.0),
mob_(2*grid.number_of_cells, -1.0)
#ifdef EXPERIMENT_GAUSS_SEIDEL
, ia_upw_(grid.number_of_cells + 1, -1),
ja_upw_(grid.number_of_faces, -1),
ia_downw_(grid.number_of_cells + 1, -1),
ja_downw_(grid.number_of_faces, -1)
, ia_upw_(grid.number_of_cells + 1, -1),
ja_upw_(grid.number_of_faces, -1),
ia_downw_(grid.number_of_cells + 1, -1),
ja_downw_(grid.number_of_faces, -1)
#endif
{
if (props.numPhases() != 2) {
THROW("Property object must have 2 phases");
}
visc_ = props.viscosity();
int num_cells = props.numCells();
smin_.resize(props.numPhases()*num_cells);
smax_.resize(props.numPhases()*num_cells);
std::vector<int> cells(num_cells);
for (int i = 0; i < num_cells; ++i) {
cells[i] = i;
}
props.satRange(props.numCells(), &cells[0], &smin_[0], &smax_[0]);
if (props.numPhases() != 2) {
THROW("Property object must have 2 phases");
}
visc_ = props.viscosity();
int num_cells = props.numCells();
smin_.resize(props.numPhases()*num_cells);
smax_.resize(props.numPhases()*num_cells);
std::vector<int> cells(num_cells);
for (int i = 0; i < num_cells; ++i) {
cells[i] = i;
}
props.satRange(props.numCells(), &cells[0], &smin_[0], &smax_[0]);
}
void TransportModelTwophase::solve(const double* darcyflux,
const double* porevolume,
const double* source,
const double dt,
std::vector<double>& saturation)
const double* source,
const double dt,
std::vector<double>& saturation)
{
darcyflux_ = darcyflux;
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
dt_ = dt;
source_ = source;
dt_ = dt;
toWaterSat(saturation, saturation_);
#ifdef EXPERIMENT_GAUSS_SEIDEL
std::vector<int> seq(grid_.number_of_cells);
std::vector<int> comp(grid_.number_of_cells + 1);
int ncomp;
compute_sequence_graph(&grid_, darcyflux_,
&seq[0], &comp[0], &ncomp,
&ia_upw_[0], &ja_upw_[0]);
const int nf = grid_.number_of_faces;
std::vector<double> neg_darcyflux(nf);
std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate<double>());
compute_sequence_graph(&grid_, &neg_darcyflux[0],
&seq[0], &comp[0], &ncomp,
&ia_downw_[0], &ja_downw_[0]);
std::vector<int> seq(grid_.number_of_cells);
std::vector<int> comp(grid_.number_of_cells + 1);
int ncomp;
compute_sequence_graph(&grid_, darcyflux_,
&seq[0], &comp[0], &ncomp,
&ia_upw_[0], &ja_upw_[0]);
const int nf = grid_.number_of_faces;
std::vector<double> neg_darcyflux(nf);
std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate<double>());
compute_sequence_graph(&grid_, &neg_darcyflux[0],
&seq[0], &comp[0], &ncomp,
&ia_downw_[0], &ja_downw_[0]);
#endif
reorderAndTransport(grid_, darcyflux);
reorderAndTransport(grid_, darcyflux);
toBothSat(saturation_, saturation);
}
@ -114,330 +114,330 @@ namespace Opm
// Influxes are negative, outfluxes positive.
struct TransportModelTwophase::Residual
{
int cell;
double s0;
double influx; // sum_j min(v_ij, 0)*f(s_j) + q_w
double outflux; // sum_j max(v_ij, 0) - q
int cell;
double s0;
double influx; // sum_j min(v_ij, 0)*f(s_j) + q_w
double outflux; // sum_j max(v_ij, 0) - q
double comp_term; // q - sum_j v_ij
double dtpv; // dt/pv(i)
const TransportModelTwophase& tm;
explicit Residual(const TransportModelTwophase& tmodel, int cell_index)
: tm(tmodel)
{
cell = cell_index;
s0 = tm.saturation_[cell];
double dtpv; // dt/pv(i)
const TransportModelTwophase& tm;
explicit Residual(const TransportModelTwophase& tmodel, int cell_index)
: tm(tmodel)
{
cell = cell_index;
s0 = tm.saturation_[cell];
double src_flux = -tm.source_[cell];
bool src_is_inflow = src_flux < 0.0;
influx = src_is_inflow ? src_flux : 0.0;
outflux = !src_is_inflow ? src_flux : 0.0;
influx = src_is_inflow ? src_flux : 0.0;
outflux = !src_is_inflow ? src_flux : 0.0;
comp_term = tm.source_[cell]; // Note: this assumes that all source flux is water.
dtpv = tm.dt_/tm.porevolume_[cell];
dtpv = tm.dt_/tm.porevolume_[cell];
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
int f = tm.grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == tm.grid_.face_cells[2*f]) {
flux = tm.darcyflux_[f];
other = tm.grid_.face_cells[2*f+1];
} else {
flux =-tm.darcyflux_[f];
other = tm.grid_.face_cells[2*f];
}
// Add flux to influx or outflux, if interior.
if (other != -1) {
if (flux < 0.0) {
influx += flux*tm.fractionalflow_[other];
} else {
outflux += flux;
}
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
int f = tm.grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == tm.grid_.face_cells[2*f]) {
flux = tm.darcyflux_[f];
other = tm.grid_.face_cells[2*f+1];
} else {
flux =-tm.darcyflux_[f];
other = tm.grid_.face_cells[2*f];
}
// Add flux to influx or outflux, if interior.
if (other != -1) {
if (flux < 0.0) {
influx += flux*tm.fractionalflow_[other];
} else {
outflux += flux;
}
comp_term -= flux;
}
}
}
double operator()(double s) const
{
return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx + s*comp_term);
}
}
}
}
double operator()(double s) const
{
return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx + s*comp_term);
}
};
void TransportModelTwophase::solveSingleCell(const int cell)
{
Residual res(*this, cell);
// const double r0 = res(saturation_[cell]);
// if (std::fabs(r0) < tol_) {
// return;
// }
int iters_used;
// saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used);
fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
Residual res(*this, cell);
// const double r0 = res(saturation_[cell]);
// if (std::fabs(r0) < tol_) {
// return;
// }
int iters_used;
// saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used);
fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
}
// namespace {
// class TofComputer
// {
// public:
// TofComputer(const int num_cells,
// const int* ia,
// const int* ja,
// const int startcell,
// std::vector<int>& tof)
// : ia_(ia),
// ja_(ja)
// {
// tof.clear();
// tof.resize(num_cells, num_cells);
// tof[startcell] = 0;
// tof_ = &tof[0];
// visitTof(startcell);
// }
// class TofComputer
// {
// public:
// TofComputer(const int num_cells,
// const int* ia,
// const int* ja,
// const int startcell,
// std::vector<int>& tof)
// : ia_(ia),
// ja_(ja)
// {
// tof.clear();
// tof.resize(num_cells, num_cells);
// tof[startcell] = 0;
// tof_ = &tof[0];
// visitTof(startcell);
// }
// private:
// const int* ia_;
// const int* ja_;
// int* tof_;
// private:
// const int* ia_;
// const int* ja_;
// int* tof_;
// void visitTof(const int cell)
// {
// for (int j = ia_[cell]; j < ia_[cell+1]; ++j) {
// const int nb_cell = ja_[j];
// if (tof_[nb_cell] > tof_[cell] + 1) {
// tof_[nb_cell] = tof_[cell] + 1;
// visitTof(nb_cell);
// }
// }
// }
// void visitTof(const int cell)
// {
// for (int j = ia_[cell]; j < ia_[cell+1]; ++j) {
// const int nb_cell = ja_[j];
// if (tof_[nb_cell] > tof_[cell] + 1) {
// tof_[nb_cell] = tof_[cell] + 1;
// visitTof(nb_cell);
// }
// }
// }
// };
// };
// } // anon namespace
void TransportModelTwophase::solveMultiCell(const int num_cells, const int* cells)
{
// std::ofstream os("dump");
// std::copy(cells, cells + num_cells, std::ostream_iterator<double>(os, "\n"));
// std::ofstream os("dump");
// std::copy(cells, cells + num_cells, std::ostream_iterator<double>(os, "\n"));
// Experiment: try a breath-first search to build a more suitable ordering.
// Verdict: failed to improve #iterations.
// {
// std::vector<int> pos(grid_.number_of_cells, -1);
// for (int i = 0; i < num_cells; ++i) {
// const int cell = cells[i];
// pos[cell] = i;
// }
// std::vector<int> done_pos(num_cells, 0);
// std::vector<int> upstream_pos;
// std::vector<int> new_pos;
// upstream_pos.push_back(0);
// done_pos[0] = 1;
// int current = 0;
// while (int(new_pos.size()) < num_cells) {
// const int i = upstream_pos[current++];
// new_pos.push_back(i);
// const int cell = cells[i];
// for (int j = ia_[cell]; j < ia_[cell+1]; ++j) {
// const int opos = pos[ja_[j]];
// if (!done_pos[opos]) {
// upstream_pos.push_back(opos);
// done_pos[opos] = 1;
// }
// }
// }
// std::reverse(new_pos.begin(), new_pos.end());
// std::copy(new_pos.begin(), new_pos.end(), const_cast<int*>(cells));
// }
// Experiment: try a breath-first search to build a more suitable ordering.
// Verdict: failed to improve #iterations.
// {
// std::vector<int> pos(grid_.number_of_cells, -1);
// for (int i = 0; i < num_cells; ++i) {
// const int cell = cells[i];
// pos[cell] = i;
// }
// std::vector<int> done_pos(num_cells, 0);
// std::vector<int> upstream_pos;
// std::vector<int> new_pos;
// upstream_pos.push_back(0);
// done_pos[0] = 1;
// int current = 0;
// while (int(new_pos.size()) < num_cells) {
// const int i = upstream_pos[current++];
// new_pos.push_back(i);
// const int cell = cells[i];
// for (int j = ia_[cell]; j < ia_[cell+1]; ++j) {
// const int opos = pos[ja_[j]];
// if (!done_pos[opos]) {
// upstream_pos.push_back(opos);
// done_pos[opos] = 1;
// }
// }
// }
// std::reverse(new_pos.begin(), new_pos.end());
// std::copy(new_pos.begin(), new_pos.end(), const_cast<int*>(cells));
// }
// Experiment: try a random ordering.
// Verdict: amazingly, reduced #iterations by approx. 25%!
// int* c = const_cast<int*>(cells);
// std::random_shuffle(c, c + num_cells);
// Experiment: try a random ordering.
// Verdict: amazingly, reduced #iterations by approx. 25%!
// int* c = const_cast<int*>(cells);
// std::random_shuffle(c, c + num_cells);
// Experiment: compute topological tof from first cell.
// Verdict: maybe useful, not tried to exploit it yet.
// std::vector<int> tof;
// TofComputer comp(grid_.number_of_cells, &ia_[0], &ja_[0], cells[0], tof);
// std::ofstream tofdump("tofdump");
// std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tofdump, "\n"));
// Experiment: compute topological tof from first cell.
// Verdict: maybe useful, not tried to exploit it yet.
// std::vector<int> tof;
// TofComputer comp(grid_.number_of_cells, &ia_[0], &ja_[0], cells[0], tof);
// std::ofstream tofdump("tofdump");
// std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tofdump, "\n"));
// Experiment: implement a metric measuring badness of ordering
// as average distance in (cyclic) ordering from
// upstream neighbours.
// Verdict: does not seem to predict #iterations very well, if at all.
// std::vector<int> pos(grid_.number_of_cells, -1);
// for (int i = 0; i < num_cells; ++i) {
// const int cell = cells[i];
// pos[cell] = i;
// }
// double diffsum = 0;
// for (int i = 0; i < num_cells; ++i) {
// const int cell = cells[i];
// int num_upstream = 0;
// int loc_diffsum = 0;
// for (int j = ia_[cell]; j < ia_[cell+1]; ++j) {
// const int opos = pos[ja_[j]];
// if (opos == -1) {
// std::cout << "Hmmm" << std::endl;
// continue;
// }
// ++num_upstream;
// const int diff = (i - opos + num_cells) % num_cells;
// loc_diffsum += diff;
// }
// diffsum += double(loc_diffsum)/double(num_upstream);
// }
// std::cout << "Average distance from upstream neighbours: " << diffsum/double(num_cells)
// << std::endl;
// Experiment: implement a metric measuring badness of ordering
// as average distance in (cyclic) ordering from
// upstream neighbours.
// Verdict: does not seem to predict #iterations very well, if at all.
// std::vector<int> pos(grid_.number_of_cells, -1);
// for (int i = 0; i < num_cells; ++i) {
// const int cell = cells[i];
// pos[cell] = i;
// }
// double diffsum = 0;
// for (int i = 0; i < num_cells; ++i) {
// const int cell = cells[i];
// int num_upstream = 0;
// int loc_diffsum = 0;
// for (int j = ia_[cell]; j < ia_[cell+1]; ++j) {
// const int opos = pos[ja_[j]];
// if (opos == -1) {
// std::cout << "Hmmm" << std::endl;
// continue;
// }
// ++num_upstream;
// const int diff = (i - opos + num_cells) % num_cells;
// loc_diffsum += diff;
// }
// diffsum += double(loc_diffsum)/double(num_upstream);
// }
// std::cout << "Average distance from upstream neighbours: " << diffsum/double(num_cells)
// << std::endl;
#ifdef EXPERIMENT_GAUSS_SEIDEL
// Experiment: when a cell changes more than the tolerance,
// mark all downwind cells as needing updates. After
// computing a single update in each cell, use marks
// to guide further updating. Clear mark in cell when
// its solution gets updated.
// Verdict: this is a good one! Approx. halved total time.
std::vector<int> needs_update(num_cells, 1);
// This one also needs the mapping from all cells to
// the strongly connected subset to filter out connections
std::vector<int> pos(grid_.number_of_cells, -1);
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
pos[cell] = i;
}
// Experiment: when a cell changes more than the tolerance,
// mark all downwind cells as needing updates. After
// computing a single update in each cell, use marks
// to guide further updating. Clear mark in cell when
// its solution gets updated.
// Verdict: this is a good one! Approx. halved total time.
std::vector<int> needs_update(num_cells, 1);
// This one also needs the mapping from all cells to
// the strongly connected subset to filter out connections
std::vector<int> pos(grid_.number_of_cells, -1);
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
pos[cell] = i;
}
// Note: partially copied from below.
const double tol = 1e-9;
const int max_iters = 300;
// Must store s0 before we start.
std::vector<double> s0(num_cells);
// Must set initial fractional flows before we start.
// Also, we compute the # of upstream neighbours.
// std::vector<int> num_upstream(num_cells);
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
s0[i] = saturation_[cell];
// num_upstream[i] = ia_upw_[cell + 1] - ia_upw_[cell];
}
// Solve once in each cell.
// std::vector<int> fully_marked_stack;
// fully_marked_stack.reserve(num_cells);
int num_iters = 0;
int update_count = 0; // Change name/meaning to cells_updated?
do {
update_count = 0; // Must reset count for every iteration.
for (int i = 0; i < num_cells; ++i) {
// while (!fully_marked_stack.empty()) {
// // std::cout << "# fully marked cells = " << fully_marked_stack.size() << std::endl;
// const int fully_marked_ci = fully_marked_stack.back();
// fully_marked_stack.pop_back();
// ++update_count;
// const int cell = cells[fully_marked_ci];
// const double old_s = saturation_[cell];
// saturation_[cell] = s0[fully_marked_ci];
// solveSingleCell(cell);
// const double s_change = std::fabs(saturation_[cell] - old_s);
// if (s_change > tol) {
// // Mark downwind cells.
// for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) {
// const int downwind_cell = ja_downw_[j];
// int ci = pos[downwind_cell];
// ++needs_update[ci];
// if (needs_update[ci] == num_upstream[ci]) {
// fully_marked_stack.push_back(ci);
// }
// }
// }
// // Unmark this cell.
// needs_update[fully_marked_ci] = 0;
// }
if (!needs_update[i]) {
continue;
}
++update_count;
const int cell = cells[i];
const double old_s = saturation_[cell];
saturation_[cell] = s0[i];
solveSingleCell(cell);
const double s_change = std::fabs(saturation_[cell] - old_s);
if (s_change > tol) {
// Mark downwind cells.
for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) {
const int downwind_cell = ja_downw_[j];
int ci = pos[downwind_cell];
// Note: partially copied from below.
const double tol = 1e-9;
const int max_iters = 300;
// Must store s0 before we start.
std::vector<double> s0(num_cells);
// Must set initial fractional flows before we start.
// Also, we compute the # of upstream neighbours.
// std::vector<int> num_upstream(num_cells);
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
s0[i] = saturation_[cell];
// num_upstream[i] = ia_upw_[cell + 1] - ia_upw_[cell];
}
// Solve once in each cell.
// std::vector<int> fully_marked_stack;
// fully_marked_stack.reserve(num_cells);
int num_iters = 0;
int update_count = 0; // Change name/meaning to cells_updated?
do {
update_count = 0; // Must reset count for every iteration.
for (int i = 0; i < num_cells; ++i) {
// while (!fully_marked_stack.empty()) {
// // std::cout << "# fully marked cells = " << fully_marked_stack.size() << std::endl;
// const int fully_marked_ci = fully_marked_stack.back();
// fully_marked_stack.pop_back();
// ++update_count;
// const int cell = cells[fully_marked_ci];
// const double old_s = saturation_[cell];
// saturation_[cell] = s0[fully_marked_ci];
// solveSingleCell(cell);
// const double s_change = std::fabs(saturation_[cell] - old_s);
// if (s_change > tol) {
// // Mark downwind cells.
// for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) {
// const int downwind_cell = ja_downw_[j];
// int ci = pos[downwind_cell];
// ++needs_update[ci];
// if (needs_update[ci] == num_upstream[ci]) {
// fully_marked_stack.push_back(ci);
// }
// }
// }
// // Unmark this cell.
// needs_update[fully_marked_ci] = 0;
// }
if (!needs_update[i]) {
continue;
}
++update_count;
const int cell = cells[i];
const double old_s = saturation_[cell];
saturation_[cell] = s0[i];
solveSingleCell(cell);
const double s_change = std::fabs(saturation_[cell] - old_s);
if (s_change > tol) {
// Mark downwind cells.
for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) {
const int downwind_cell = ja_downw_[j];
int ci = pos[downwind_cell];
if (ci != -1) {
needs_update[ci] = 1;
}
// ++needs_update[ci];
// if (needs_update[ci] == num_upstream[ci]) {
// fully_marked_stack.push_back(ci);
// }
}
}
// Unmark this cell.
needs_update[i] = 0;
}
// std::cout << "Iter = " << num_iters << " update_count = " << update_count
// << " # marked cells = "
// << std::accumulate(needs_update.begin(), needs_update.end(), 0) << std::endl;
} while (update_count > 0 && ++num_iters < max_iters);
// ++needs_update[ci];
// if (needs_update[ci] == num_upstream[ci]) {
// fully_marked_stack.push_back(ci);
// }
}
}
// Unmark this cell.
needs_update[i] = 0;
}
// std::cout << "Iter = " << num_iters << " update_count = " << update_count
// << " # marked cells = "
// << std::accumulate(needs_update.begin(), needs_update.end(), 0) << std::endl;
} while (update_count > 0 && ++num_iters < max_iters);
// Done with iterations, check if we succeeded.
if (update_count > 0) {
THROW("In solveMultiCell(), we did not converge after "
<< num_iters << " iterations. Remaining update count = " << update_count);
}
std::cout << "Solved " << num_cells << " cell multicell problem in "
<< num_iters << " iterations." << std::endl;
// Done with iterations, check if we succeeded.
if (update_count > 0) {
THROW("In solveMultiCell(), we did not converge after "
<< num_iters << " iterations. Remaining update count = " << update_count);
}
std::cout << "Solved " << num_cells << " cell multicell problem in "
<< num_iters << " iterations." << std::endl;
#else
double max_s_change = 0.0;
const double tol = 1e-9;
const int max_iters = 300;
int num_iters = 0;
// Must store s0 before we start.
std::vector<double> s0(num_cells);
// Must set initial fractional flows before we start.
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
s0[i] = saturation_[cell];
}
do {
max_s_change = 0.0;
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
const double old_s = saturation_[cell];
saturation_[cell] = s0[i];
solveSingleCell(cell);
double s_change = std::fabs(saturation_[cell] - old_s);
// std::cout << "cell = " << cell << " delta s = " << s_change << std::endl;
if (max_s_change < s_change) {
max_s_change = s_change;
}
}
// std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change
// << " in cell " << max_change_cell << std::endl;
} while (max_s_change > tol && ++num_iters < max_iters);
if (max_s_change > tol) {
THROW("In solveMultiCell(), we did not converge after "
<< num_iters << " iterations. Delta s = " << max_s_change);
}
std::cout << "Solved " << num_cells << " cell multicell problem in "
<< num_iters << " iterations." << std::endl;
double max_s_change = 0.0;
const double tol = 1e-9;
const int max_iters = 300;
int num_iters = 0;
// Must store s0 before we start.
std::vector<double> s0(num_cells);
// Must set initial fractional flows before we start.
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
s0[i] = saturation_[cell];
}
do {
max_s_change = 0.0;
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
const double old_s = saturation_[cell];
saturation_[cell] = s0[i];
solveSingleCell(cell);
double s_change = std::fabs(saturation_[cell] - old_s);
// std::cout << "cell = " << cell << " delta s = " << s_change << std::endl;
if (max_s_change < s_change) {
max_s_change = s_change;
}
}
// std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change
// << " in cell " << max_change_cell << std::endl;
} while (max_s_change > tol && ++num_iters < max_iters);
if (max_s_change > tol) {
THROW("In solveMultiCell(), we did not converge after "
<< num_iters << " iterations. Delta s = " << max_s_change);
}
std::cout << "Solved " << num_cells << " cell multicell problem in "
<< num_iters << " iterations." << std::endl;
#endif // EXPERIMENT_GAUSS_SEIDEL
}
double TransportModelTwophase::fracFlow(double s, int cell) const
{
double sat[2] = { s, 1.0 - s };
double mob[2];
props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[0];
mob[1] /= visc_[1];
return mob[0]/(mob[0] + mob[1]);
double sat[2] = { s, 1.0 - s };
double mob[2];
props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[0];
mob[1] /= visc_[1];
return mob[0]/(mob[0] + mob[1]);
}
@ -450,19 +450,19 @@ namespace Opm
//
struct TransportModelTwophase::GravityResidual
{
int cell;
int cell;
int nbcell[2];
double s0;
double dtpv; // dt/pv(i)
double s0;
double dtpv; // dt/pv(i)
double gf[2];
const TransportModelTwophase& tm;
explicit GravityResidual(const TransportModelTwophase& tmodel,
const TransportModelTwophase& tm;
explicit GravityResidual(const TransportModelTwophase& tmodel,
const std::vector<int>& cells,
const int pos,
const double* gravflux) // Always oriented towards next in column. Size = colsize - 1.
: tm(tmodel)
{
cell = cells[pos];
: tm(tmodel)
{
cell = cells[pos];
nbcell[0] = -1;
gf[0] = 0.0;
if (pos > 0) {
@ -475,40 +475,40 @@ namespace Opm
nbcell[1] = cells[pos + 1];
gf[1] = gravflux[pos];
}
s0 = tm.saturation_[cell];
dtpv = tm.dt_/tm.porevolume_[cell];
}
double operator()(double s) const
{
double res = s - s0;
s0 = tm.saturation_[cell];
dtpv = tm.dt_/tm.porevolume_[cell];
}
double operator()(double s) const
{
double res = s - s0;
double mobcell[2];
tm.mobility(s, cell, mobcell);
for (int nb = 0; nb < 2; ++nb) {
if (nbcell[nb] != -1) {
if (nbcell[nb] != -1) {
double m[2];
if (gf[nb] < 0.0) {
m[0] = mobcell[0];
m[1] = tm.mob_[2*nbcell[nb] + 1];
} else {
} else {
m[0] = tm.mob_[2*nbcell[nb]];
m[1] = mobcell[1];
}
if (m[0] + m[1] > 0.0) {
}
if (m[0] + m[1] > 0.0) {
res += -dtpv*gf[nb]*m[0]*m[1]/(m[0] + m[1]);
}
}
}
}
return res;
}
}
};
void TransportModelTwophase::mobility(double s, int cell, double* mob) const
{
double sat[2] = { s, 1.0 - s };
props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[0];
mob[1] /= visc_[1];
double sat[2] = { s, 1.0 - s };
props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[0];
mob[1] /= visc_[1];
}
@ -548,7 +548,7 @@ namespace Opm
saturation_[cell] = RootFinder::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
}
saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]);
mobility(saturation_[cell], cell, &mob_[2*cell]);
mobility(saturation_[cell], cell, &mob_[2*cell]);
}
@ -559,17 +559,17 @@ namespace Opm
const int nc = cells.size();
std::vector<double> col_gravflux(nc - 1);
for (int ci = 0; ci < nc - 1; ++ci) {
const int cell = cells[ci];
const int next_cell = cells[ci + 1];
for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) {
const int face = grid_.cell_faces[j];
const int c1 = grid_.face_cells[2*face + 0];
const int cell = cells[ci];
const int next_cell = cells[ci + 1];
for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) {
const int face = grid_.cell_faces[j];
const int c1 = grid_.face_cells[2*face + 0];
const int c2 = grid_.face_cells[2*face + 1];
if (c1 == next_cell || c2 == next_cell) {
if (c1 == next_cell || c2 == next_cell) {
const double gf = gravflux_[face];
col_gravflux[ci] = (c1 == cell) ? gf : -gf;
}
}
}
}
}
// Store initial saturation s0
@ -579,7 +579,7 @@ namespace Opm
}
// Solve single cell problems, repeating if necessary.
double max_s_change = 0.0;
double max_s_change = 0.0;
int num_iters = 0;
do {
max_s_change = 0.0;
@ -595,12 +595,12 @@ namespace Opm
std::fabs(saturation_[cells[ci2]] - old_s[1])));
}
// std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change << std::endl;
} while (max_s_change > tol_ && ++num_iters < maxit_);
} while (max_s_change > tol_ && ++num_iters < maxit_);
if (max_s_change > tol_) {
THROW("In solveGravityColumn(), we did not converge after "
<< num_iters << " iterations. Delta s = " << max_s_change);
}
if (max_s_change > tol_) {
THROW("In solveGravityColumn(), we did not converge after "
<< num_iters << " iterations. Delta s = " << max_s_change);
}
return num_iters + 1;
}

88
opm/core/transport/reorder/TransportModelTwophase.hpp
View File

@ -35,27 +35,27 @@ namespace Opm
class TransportModelTwophase : public TransportModelInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] tol Tolerance used in the solver.
/// \param[in] maxit Maximum number of non-linear iterations used.
TransportModelTwophase(const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const double tol,
const int maxit);
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] tol Tolerance used in the solver.
/// \param[in] maxit Maximum number of non-linear iterations used.
TransportModelTwophase(const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const double tol,
const int maxit);
/// Solve for saturation at next timestep.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
void solve(const double* darcyflux,
/// Solve for saturation at next timestep.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
void solve(const double* darcyflux,
const double* porevolume,
const double* source,
const double dt,
std::vector<double>& saturation);
const double* source,
const double dt,
std::vector<double>& saturation);
/// Initialise quantities needed by gravity solver.
/// \param[in] grav gravity vector
@ -66,18 +66,18 @@ namespace Opm
/// It assumes that the input columns contain cells in a single
/// vertical stack, that do not interact with other columns (for
/// gravity segregation.
/// \param[in] columns Vector of cell-columns.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
/// \param[in] columns Vector of cell-columns.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
void solveGravity(const std::vector<std::vector<int> >& columns,
const double* porevolume,
const double dt,
std::vector<double>& saturation);
private:
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
void solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
@ -85,33 +85,33 @@ namespace Opm
int solveGravityColumn(const std::vector<int>& cells);
private:
const UnstructuredGrid& grid_;
const IncompPropertiesInterface& props_;
const double* visc_;
std::vector<double> smin_;
std::vector<double> smax_;
double tol_;
double maxit_;
const UnstructuredGrid& grid_;
const IncompPropertiesInterface& props_;
const double* visc_;
std::vector<double> smin_;
std::vector<double> smax_;
double tol_;
double maxit_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one source per cell
double dt_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one source per cell
double dt_;
std::vector<double> saturation_; // one per cell, only water saturation!
std::vector<double> fractionalflow_; // = m[0]/(m[0] + m[1]) per cell
std::vector<double> fractionalflow_; // = m[0]/(m[0] + m[1]) per cell
// For gravity segregation.
std::vector<double> gravflux_;
std::vector<double> mob_;
std::vector<double> s0_;
// Storing the upwind and downwind graphs for experiments.
std::vector<int> ia_upw_;
std::vector<int> ja_upw_;
std::vector<int> ia_downw_;
std::vector<int> ja_downw_;
// Storing the upwind and downwind graphs for experiments.
std::vector<int> ia_upw_;
std::vector<int> ja_upw_;
std::vector<int> ia_downw_;
std::vector<int> ja_downw_;
struct Residual;
double fracFlow(double s, int cell) const;
struct Residual;
double fracFlow(double s, int cell) const;
struct GravityResidual;
void mobility(double s, int cell, double* mob) const;