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

Whitespace cleanup.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2012-08-17 12:40:15 +02:00
parent 8fa0140e66
commit 9a1bcaf496
3 changed files with 462 additions and 462 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -33,33 +33,33 @@ 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 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,
/// 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.
void initGravity();
@@ -69,11 +69,11 @@ 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] 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.
/// \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,
@@ -82,8 +82,8 @@ namespace Opm
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);
@@ -91,24 +91,24 @@ namespace Opm
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_;
@@ -116,14 +116,14 @@ namespace Opm
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;