From ded1424cc928599081a66a580d8dc07a30026fee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 17 Aug 2012 12:40:15 +0200 Subject: [PATCH] Whitespace cleanup. --- .../TransportModelCompressibleTwophase.hpp | 92 +-- .../reorder/TransportModelTwophase.cpp | 744 +++++++++--------- .../reorder/TransportModelTwophase.hpp | 88 +-- 3 files changed, 462 insertions(+), 462 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index c9f7e96a21..19c44b22dd 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -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& saturation); + const double* source, + const double dt, + std::vector& 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 >& columns, const double* pressure, const double* porevolume0, @@ -82,8 +82,8 @@ namespace Opm std::vector& 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& 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 allcells_; std::vector visc_; std::vector A_; - std::vector smin_; - std::vector smax_; - double tol_; - double maxit_; + std::vector smin_; + std::vector 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 saturation_; // P (= num. phases) per cell - std::vector fractionalflow_; // = m[0]/(m[0] + m[1]) per cell + std::vector fractionalflow_; // = m[0]/(m[0] + m[1]) per cell // For gravity segregation. std::vector trans_; std::vector density_; @@ -116,14 +116,14 @@ namespace Opm std::vector mob_; std::vector s0_; - // Storing the upwind and downwind graphs for experiments. - std::vector ia_upw_; - std::vector ja_upw_; - std::vector ia_downw_; - std::vector ja_downw_; + // Storing the upwind and downwind graphs for experiments. + std::vector ia_upw_; + std::vector ja_upw_; + std::vector ia_downw_; + std::vector 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; diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index ff6d2c05d9..b3630ea731 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -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 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 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& saturation) + const double* source, + const double dt, + std::vector& saturation) { - darcyflux_ = darcyflux; + darcyflux_ = darcyflux; porevolume_ = porevolume; - source_ = source; - dt_ = dt; + source_ = source; + dt_ = dt; toWaterSat(saturation, saturation_); #ifdef EXPERIMENT_GAUSS_SEIDEL - std::vector seq(grid_.number_of_cells); - std::vector 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 neg_darcyflux(nf); - std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate()); - compute_sequence_graph(&grid_, &neg_darcyflux[0], - &seq[0], &comp[0], &ncomp, - &ia_downw_[0], &ja_downw_[0]); + std::vector seq(grid_.number_of_cells); + std::vector 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 neg_darcyflux(nf); + std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate()); + 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& 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& 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(os, "\n")); + // std::ofstream os("dump"); + // std::copy(cells, cells + num_cells, std::ostream_iterator(os, "\n")); - // Experiment: try a breath-first search to build a more suitable ordering. - // Verdict: failed to improve #iterations. - // { - // std::vector pos(grid_.number_of_cells, -1); - // for (int i = 0; i < num_cells; ++i) { - // const int cell = cells[i]; - // pos[cell] = i; - // } - // std::vector done_pos(num_cells, 0); - // std::vector upstream_pos; - // std::vector 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(cells)); - // } + // Experiment: try a breath-first search to build a more suitable ordering. + // Verdict: failed to improve #iterations. + // { + // std::vector pos(grid_.number_of_cells, -1); + // for (int i = 0; i < num_cells; ++i) { + // const int cell = cells[i]; + // pos[cell] = i; + // } + // std::vector done_pos(num_cells, 0); + // std::vector upstream_pos; + // std::vector 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(cells)); + // } - // Experiment: try a random ordering. - // Verdict: amazingly, reduced #iterations by approx. 25%! - // int* c = const_cast(cells); - // std::random_shuffle(c, c + num_cells); + // Experiment: try a random ordering. + // Verdict: amazingly, reduced #iterations by approx. 25%! + // int* c = const_cast(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 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(tofdump, "\n")); + // Experiment: compute topological tof from first cell. + // Verdict: maybe useful, not tried to exploit it yet. + // std::vector 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(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 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 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 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 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 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 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 s0(num_cells); - // Must set initial fractional flows before we start. - // Also, we compute the # of upstream neighbours. - // std::vector 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 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 s0(num_cells); + // Must set initial fractional flows before we start. + // Also, we compute the # of upstream neighbours. + // std::vector 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 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 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 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& 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 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; } diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index ca302c18ad..74ed2bfe4f 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -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& saturation); + const double* source, + const double dt, + std::vector& 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 >& columns, const double* porevolume, const double dt, std::vector& 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& cells, const int pos, @@ -85,33 +85,33 @@ namespace Opm int solveGravityColumn(const std::vector& cells); private: - const UnstructuredGrid& grid_; - const IncompPropertiesInterface& props_; - const double* visc_; - std::vector smin_; - std::vector smax_; - double tol_; - double maxit_; + const UnstructuredGrid& grid_; + const IncompPropertiesInterface& props_; + const double* visc_; + std::vector smin_; + std::vector 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 saturation_; // one per cell, only water saturation! - std::vector fractionalflow_; // = m[0]/(m[0] + m[1]) per cell + std::vector fractionalflow_; // = m[0]/(m[0] + m[1]) per cell // For gravity segregation. std::vector gravflux_; std::vector mob_; std::vector s0_; - // Storing the upwind and downwind graphs for experiments. - std::vector ia_upw_; - std::vector ja_upw_; - std::vector ia_downw_; - std::vector ja_downw_; + // Storing the upwind and downwind graphs for experiments. + std::vector ia_upw_; + std::vector ja_upw_; + std::vector ia_downw_; + std::vector 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;