From 6ab748721836f04890f55cc41f5c2d0fbd09d3a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Jul 2012 15:34:42 +0200 Subject: [PATCH] Whitespace cleanup. --- .../TransportModelCompressibleTwophase.cpp | 460 +++++++++--------- 1 file changed, 230 insertions(+), 230 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 8d221f924..ce8ebb97e 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -39,38 +39,38 @@ namespace Opm TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid, - const Opm::BlackoilPropertiesInterface& 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::BlackoilPropertiesInterface& 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), 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) + ja_upw_(grid.number_of_faces, -1), + ia_downw_(grid.number_of_cells + 1, -1), + ja_downw_(grid.number_of_faces, -1) { - if (props.numPhases() != 2) { - THROW("Property object must have 2 phases"); - } + if (props.numPhases() != 2) { + THROW("Property object must have 2 phases"); + } int np = props.numPhases(); - int num_cells = props.numCells(); - visc_.resize(np*num_cells); + int num_cells = props.numCells(); + visc_.resize(np*num_cells); A_.resize(np*np*num_cells); - smin_.resize(np*num_cells); - smax_.resize(np*num_cells); - allcells_.resize(num_cells); - for (int i = 0; i < num_cells; ++i) { - allcells_[i] = i; - } - props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]); + smin_.resize(np*num_cells); + smax_.resize(np*num_cells); + allcells_.resize(num_cells); + for (int i = 0; i < num_cells; ++i) { + allcells_[i] = i; + } + props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]); } void TransportModelCompressibleTwophase::solve(const double* darcyflux, @@ -82,30 +82,30 @@ namespace Opm const double dt, std::vector& saturation) { - darcyflux_ = darcyflux; + darcyflux_ = darcyflux; surfacevol0_ = surfacevol0; porevolume0_ = porevolume0; porevolume_ = porevolume; - source_ = source; - dt_ = dt; + source_ = source; + dt_ = dt; toWaterSat(saturation, saturation_); props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); - 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]); - reorderAndTransport(grid_, darcyflux); + 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]); + reorderAndTransport(grid_, darcyflux); toBothSat(saturation_, saturation); } @@ -116,187 +116,187 @@ namespace Opm // r(s) = s - B*z0 + dt/pv*( influx + outflux*f(s) ) // // @@@ What about the source term - // + // // where influx is water influx, outflux is total outflux. // We need the formula influx = B_i sum_{j->i} b_j v_{ij} - B_i q_w. // outflux = B_i sum_{i->j} b_i v_{ij} - B_i q = sum_{i->j} v_{ij} - B_i q // Influxes are negative, outfluxes positive. struct TransportModelCompressibleTwophase::Residual { - int cell; + int cell; double B_cell; - double z0; - double influx; // B_i sum_j b_j min(v_ij, 0)*f(s_j) - B_i q_w - double outflux; // sum_j max(v_ij, 0) - B_i q + double z0; + double influx; // B_i sum_j b_j min(v_ij, 0)*f(s_j) - B_i q_w + double outflux; // sum_j max(v_ij, 0) - B_i q // @@@ TODO: figure out change to rock-comp. terms with fluid compr. double comp_term; // Now: used to be: q - sum_j v_ij - double dtpv; // dt/pv(i) - const TransportModelCompressibleTwophase& tm; - explicit Residual(const TransportModelCompressibleTwophase& tmodel, int cell_index) - : tm(tmodel) - { - cell = cell_index; + double dtpv; // dt/pv(i) + const TransportModelCompressibleTwophase& tm; + explicit Residual(const TransportModelCompressibleTwophase& tmodel, int cell_index) + : tm(tmodel) + { + cell = cell_index; const int np = tm.props_.numPhases(); - z0 = tm.surfacevol0_[np*cell + 0]; // I.e. water surface volume + z0 = tm.surfacevol0_[np*cell + 0]; // I.e. water surface volume B_cell = 1.0/tm.A_[np*np*cell + 0]; double src_flux = -tm.source_[cell]; bool src_is_inflow = src_flux < 0.0; - influx = src_is_inflow ? B_cell*src_flux : 0.0; - outflux = !src_is_inflow ? B_cell*src_flux : 0.0; + influx = src_is_inflow ? B_cell*src_flux : 0.0; + outflux = !src_is_inflow ? B_cell*src_flux : 0.0; comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell]; - dtpv = tm.dt_/tm.porevolume0_[cell]; - for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { - const 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) { + dtpv = tm.dt_/tm.porevolume0_[cell]; + for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { + const 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) { const double b_face = tm.A_[np*np*other + 0]; - influx += B_cell*b_face*flux*tm.fractionalflow_[other]; - } else { - outflux += flux; // Because B_cell*b_face = 1 for outflow faces - } - } - } - } - double operator()(double s) const - { - // return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx + s*comp_term); - return s - B_cell*z0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx) + s*comp_term; - } + influx += B_cell*b_face*flux*tm.fractionalflow_[other]; + } else { + outflux += flux; // Because B_cell*b_face = 1 for outflow faces + } + } + } + } + double operator()(double s) const + { + // return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx + s*comp_term); + return s - B_cell*z0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx) + s*comp_term; + } }; void TransportModelCompressibleTwophase::solveSingleCell(const int cell) { - Residual res(*this, cell); - int 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); + int iters_used; + saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); + fractionalflow_[cell] = fracFlow(saturation_[cell], cell); } void TransportModelCompressibleTwophase::solveMultiCell(const int num_cells, const int* cells) { - // 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; } double TransportModelCompressibleTwophase::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_[2*cell + 0]; - mob[1] /= visc_[2*cell + 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_[2*cell + 0]; + mob[1] /= visc_[2*cell + 1]; + return mob[0]/(mob[0] + mob[1]); } @@ -309,19 +309,19 @@ namespace Opm // struct TransportModelCompressibleTwophase::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 TransportModelCompressibleTwophase& tm; - explicit GravityResidual(const TransportModelCompressibleTwophase& tmodel, + const TransportModelCompressibleTwophase& tm; + explicit GravityResidual(const TransportModelCompressibleTwophase& 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) { @@ -334,40 +334,40 @@ namespace Opm nbcell[1] = cells[pos + 1]; gf[1] = gravflux[pos]; } - s0 = tm.saturation_[cell]; - dtpv = tm.dt_/tm.porevolume0_[cell]; - - } - double operator()(double s) const - { - double res = s - s0; + s0 = tm.saturation_[cell]; + dtpv = tm.dt_/tm.porevolume0_[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 TransportModelCompressibleTwophase::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]; } @@ -409,7 +409,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]); } @@ -420,17 +420,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 @@ -440,7 +440,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; @@ -456,12 +456,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; }