From 055f2d5f193b60a01faa9170d372acf704a86307 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Sep 2012 14:04:52 +0200 Subject: [PATCH] Whitespace cleanup and adding a check for #phases <= 3. --- opm/core/utility/miscUtilities.cpp | 340 ++++++++++++++--------------- 1 file changed, 170 insertions(+), 170 deletions(-) diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 0b8d2010..e1658b83 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -40,14 +40,14 @@ namespace Opm /// @param[out] porevol the pore volume by cell. void computePorevolume(const UnstructuredGrid& grid, const double* porosity, - std::vector& porevol) + std::vector& porevol) { - int num_cells = grid.number_of_cells; - porevol.resize(num_cells); - std::transform(porosity, porosity + num_cells, - grid.cell_volumes, - porevol.begin(), - std::multiplies()); + int num_cells = grid.number_of_cells; + porevol.resize(num_cells); + std::transform(porosity, porosity + num_cells, + grid.cell_volumes, + porevol.begin(), + std::multiplies()); } @@ -98,20 +98,20 @@ namespace Opm /// For each phase p, we compute /// sat_vol_p = sum_i s_p_i pv_i void computeSaturatedVol(const std::vector& pv, - const std::vector& s, - double* sat_vol) + const std::vector& s, + double* sat_vol) { - const int num_cells = pv.size(); - const int np = s.size()/pv.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and pv vectors do not match."); - } - std::fill(sat_vol, sat_vol + np, 0.0); - for (int c = 0; c < num_cells; ++c) { - for (int p = 0; p < np; ++p) { - sat_vol[p] += pv[c]*s[np*c + p]; - } - } + const int num_cells = pv.size(); + const int np = s.size()/pv.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and pv vectors do not match."); + } + std::fill(sat_vol, sat_vol + np, 0.0); + for (int c = 0; c < num_cells; ++c) { + for (int p = 0; p < np; ++p) { + sat_vol[p] += pv[c]*s[np*c + p]; + } + } } @@ -123,28 +123,28 @@ namespace Opm /// For each phase p, we compute /// aver_sat_p = (sum_i s_p_i pv_i) / (sum_i pv_i). void computeAverageSat(const std::vector& pv, - const std::vector& s, - double* aver_sat) + const std::vector& s, + double* aver_sat) { - const int num_cells = pv.size(); - const int np = s.size()/pv.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and pv vectors do not match."); - } - double tot_pv = 0.0; - // Note that we abuse the output array to accumulate the - // saturated pore volumes. - std::fill(aver_sat, aver_sat + np, 0.0); - for (int c = 0; c < num_cells; ++c) { - tot_pv += pv[c]; - for (int p = 0; p < np; ++p) { - aver_sat[p] += pv[c]*s[np*c + p]; - } - } - // Must divide by pore volumes to get saturations. - for (int p = 0; p < np; ++p) { - aver_sat[p] /= tot_pv; - } + const int num_cells = pv.size(); + const int np = s.size()/pv.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and pv vectors do not match."); + } + double tot_pv = 0.0; + // Note that we abuse the output array to accumulate the + // saturated pore volumes. + std::fill(aver_sat, aver_sat + np, 0.0); + for (int c = 0; c < num_cells; ++c) { + tot_pv += pv[c]; + for (int p = 0; p < np; ++p) { + aver_sat[p] += pv[c]*s[np*c + p]; + } + } + // Must divide by pore volumes to get saturations. + for (int p = 0; p < np; ++p) { + aver_sat[p] /= tot_pv; + } } @@ -161,38 +161,38 @@ namespace Opm /// where P = s.size()/src.size(). /// @param[out] produced must also point to a valid array with P elements. void computeInjectedProduced(const IncompPropertiesInterface& props, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced) + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced) { - const int num_cells = src.size(); - const int np = s.size()/src.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and src vectors do not match."); - } - std::fill(injected, injected + np, 0.0); - std::fill(produced, produced + np, 0.0); - const double* visc = props.viscosity(); - std::vector mob(np); - for (int c = 0; c < num_cells; ++c) { - if (src[c] > 0.0) { - injected[0] += src[c]*dt; - } else if (src[c] < 0.0) { - const double flux = -src[c]*dt; - const double* sat = &s[np*c]; - props.relperm(1, sat, &c, &mob[0], 0); - double totmob = 0.0; - for (int p = 0; p < np; ++p) { - mob[p] /= visc[p]; - totmob += mob[p]; - } - for (int p = 0; p < np; ++p) { - produced[p] += (mob[p]/totmob)*flux; - } - } - } + const int num_cells = src.size(); + const int np = s.size()/src.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and src vectors do not match."); + } + std::fill(injected, injected + np, 0.0); + std::fill(produced, produced + np, 0.0); + const double* visc = props.viscosity(); + std::vector mob(np); + for (int c = 0; c < num_cells; ++c) { + if (src[c] > 0.0) { + injected[0] += src[c]*dt; + } else if (src[c] < 0.0) { + const double flux = -src[c]*dt; + const double* sat = &s[np*c]; + props.relperm(1, sat, &c, &mob[0], 0); + double totmob = 0.0; + for (int p = 0; p < np; ++p) { + mob[p] /= visc[p]; + totmob += mob[p]; + } + for (int p = 0; p < np; ++p) { + produced[p] += (mob[p]/totmob)*flux; + } + } + } } @@ -203,9 +203,9 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::IncompPropertiesInterface& props, - const std::vector& cells, - const std::vector& s, - std::vector& totmob) + const std::vector& cells, + const std::vector& s, + std::vector& totmob) { std::vector pmobc; @@ -231,10 +231,10 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::IncompPropertiesInterface& props, - const std::vector& cells, - const std::vector& s, - std::vector& totmob, - std::vector& omega) + const std::vector& cells, + const std::vector& s, + std::vector& totmob, + std::vector& omega) { std::vector pmobc; @@ -331,33 +331,33 @@ namespace Opm /// (+) positive inflow of first phase (water) /// (-) negative total outflow of both phases void computeTransportSource(const UnstructuredGrid& grid, - const std::vector& src, - const std::vector& faceflux, - const double inflow_frac, + const std::vector& src, + const std::vector& faceflux, + const double inflow_frac, const Wells* wells, const std::vector& well_perfrates, - std::vector& transport_src) + std::vector& transport_src) { - int nc = grid.number_of_cells; - transport_src.resize(nc); + int nc = grid.number_of_cells; + transport_src.resize(nc); // Source term and boundary contributions. - for (int c = 0; c < nc; ++c) { - transport_src[c] = 0.0; - transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c]; - for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) { - int f = grid.cell_faces[hf]; - const int* f2c = &grid.face_cells[2*f]; - double bdy_influx = 0.0; - if (f2c[0] == c && f2c[1] == -1) { - bdy_influx = -faceflux[f]; - } else if (f2c[0] == -1 && f2c[1] == c) { - bdy_influx = faceflux[f]; - } - if (bdy_influx != 0.0) { - transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx; - } - } - } + for (int c = 0; c < nc; ++c) { + transport_src[c] = 0.0; + transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c]; + for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) { + int f = grid.cell_faces[hf]; + const int* f2c = &grid.face_cells[2*f]; + double bdy_influx = 0.0; + if (f2c[0] == c && f2c[1] == -1) { + bdy_influx = -faceflux[f]; + } else if (f2c[0] == -1 && f2c[1] == c) { + bdy_influx = faceflux[f]; + } + if (bdy_influx != 0.0) { + transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx; + } + } + } // Well contributions. if (wells) { @@ -392,52 +392,52 @@ namespace Opm /// @param[in] face_flux signed per-face fluxes /// @param[out] cell_velocity the estimated velocities. void estimateCellVelocity(const UnstructuredGrid& grid, - const std::vector& face_flux, - std::vector& cell_velocity) + const std::vector& face_flux, + std::vector& cell_velocity) { - const int dim = grid.dimensions; - cell_velocity.clear(); - cell_velocity.resize(grid.number_of_cells*dim, 0.0); - for (int face = 0; face < grid.number_of_faces; ++face) { - int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] }; - const double* fc = &grid.face_centroids[face*dim]; - double flux = face_flux[face]; - for (int i = 0; i < 2; ++i) { - if (c[i] >= 0) { - const double* cc = &grid.cell_centroids[c[i]*dim]; - for (int d = 0; d < dim; ++d) { - double v_contrib = fc[d] - cc[d]; - v_contrib *= flux/grid.cell_volumes[c[i]]; - cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib; - } - } - } - } + const int dim = grid.dimensions; + cell_velocity.clear(); + cell_velocity.resize(grid.number_of_cells*dim, 0.0); + for (int face = 0; face < grid.number_of_faces; ++face) { + int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] }; + const double* fc = &grid.face_centroids[face*dim]; + double flux = face_flux[face]; + for (int i = 0; i < 2; ++i) { + if (c[i] >= 0) { + const double* cc = &grid.cell_centroids[c[i]*dim]; + for (int d = 0; d < dim; ++d) { + double v_contrib = fc[d] - cc[d]; + v_contrib *= flux/grid.cell_volumes[c[i]]; + cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib; + } + } + } + } } /// Extract a vector of water saturations from a vector of /// interleaved water and oil saturations. void toWaterSat(const std::vector& sboth, - std::vector& sw) + std::vector& sw) { - int num = sboth.size()/2; - sw.resize(num); - for (int i = 0; i < num; ++i) { - sw[i] = sboth[2*i]; - } + int num = sboth.size()/2; + sw.resize(num); + for (int i = 0; i < num; ++i) { + sw[i] = sboth[2*i]; + } } /// Make a a vector of interleaved water and oil saturations from /// a vector of water saturations. void toBothSat(const std::vector& sw, - std::vector& sboth) + std::vector& sboth) { - int num = sw.size(); - sboth.resize(2*num); - for (int i = 0; i < num; ++i) { - sboth[2*i] = sw[i]; - sboth[2*i + 1] = 1.0 - sw[i]; - } + int num = sw.size(); + sboth.resize(2*num); + for (int i = 0; i < num; ++i) { + sboth[2*i] = sw[i]; + sboth[2*i + 1] = 1.0 - sw[i]; + } } @@ -450,30 +450,30 @@ namespace Opm if (np != 2) { THROW("wellsToSrc() requires a 2 phase case."); } - src.resize(num_cells); - for (int w = 0; w < wells.number_of_wells; ++w) { + src.resize(num_cells); + for (int w = 0; w < wells.number_of_wells; ++w) { const int cur = wells.ctrls[w]->current; - if (wells.ctrls[w]->num != 1) { - MESSAGE("In wellsToSrc(): well has more than one control, all but current control will be ignored."); - } - if (wells.ctrls[w]->type[cur] != RESERVOIR_RATE) { - THROW("In wellsToSrc(): well is something other than RESERVOIR_RATE."); - } - if (wells.well_connpos[w+1] - wells.well_connpos[w] != 1) { - THROW("In wellsToSrc(): well has multiple perforations."); - } + if (wells.ctrls[w]->num != 1) { + MESSAGE("In wellsToSrc(): well has more than one control, all but current control will be ignored."); + } + if (wells.ctrls[w]->type[cur] != RESERVOIR_RATE) { + THROW("In wellsToSrc(): well is something other than RESERVOIR_RATE."); + } + if (wells.well_connpos[w+1] - wells.well_connpos[w] != 1) { + THROW("In wellsToSrc(): well has multiple perforations."); + } for (int p = 0; p < np; ++p) { if (wells.ctrls[w]->distr[np*cur + p] != 1.0) { THROW("In wellsToSrc(): well not controlled on total rate."); } } - double flow = wells.ctrls[w]->target[cur]; + double flow = wells.ctrls[w]->target[cur]; if (wells.type[w] == INJECTOR) { flow *= wells.comp_frac[np*w + 0]; // Obtaining water rate for inflow source. } - const double cell = wells.well_cells[wells.well_connpos[w]]; - src[cell] = flow; - } + const double cell = wells.well_cells[wells.well_connpos[w]]; + src[cell] = flow; + } } @@ -594,9 +594,10 @@ namespace Opm int nw = well_bhp.size(); ASSERT(nw == wells.number_of_wells); int np = props.numPhases(); - //if (props.numPhases() != 2) { - // THROW("WellReport for now assumes two phase flow."); - //} + const int max_np = 3; + if (np > max_np) { + THROW("WellReport for now assumes #phases <= " << max_np); + } const double* visc = props.viscosity(); std::vector data_now; data_now.reserve(1 + 3*nw); @@ -606,8 +607,8 @@ namespace Opm double well_rate_total = 0.0; double well_rate_water = 0.0; for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { - const double perf_rate = unit::convert::to(well_perfrates[perf], - unit::cubic(unit::meter)/unit::day); + const double perf_rate = unit::convert::to(well_perfrates[perf], + unit::cubic(unit::meter)/unit::day); well_rate_total += perf_rate; if (perf_rate > 0.0) { // Injection. @@ -615,13 +616,13 @@ namespace Opm } else { // Production. const int cell = wells.well_cells[perf]; - double mob[np]; + double mob[max_np]; props.relperm(1, &saturation[2*cell], &cell, mob, 0); - double tmob=0; - for(int i=0; i < np; ++i){ + double tmob = 0; + for(int i = 0; i < np; ++i) { mob[i] /= visc[i]; tmob += mob[i]; - } + } const double fracflow = mob[0]/tmob; well_rate_water += perf_rate*fracflow; } @@ -652,9 +653,10 @@ namespace Opm int nw = well_bhp.size(); ASSERT(nw == wells.number_of_wells); int np = props.numPhases(); - //if (props.numPhases() != 2) { - // THROW("WellReport for now assumes two phase flow."); - //} + const int max_np = 3; + if (np > max_np) { + THROW("WellReport for now assumes #phases <= " << max_np); + } std::vector data_now; data_now.reserve(1 + 3*nw); data_now.push_back(time/unit::day); @@ -665,25 +667,23 @@ namespace Opm for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { const double perf_rate = well_perfrates[perf]*(unit::day/unit::second); well_rate_total += perf_rate; - if (perf_rate > 0.0) { + if (perf_rate > 0.0) { // Injection. well_rate_water += perf_rate*wells.comp_frac[0]; } else { // Production. const int cell = wells.well_cells[perf]; - double mob[np]; + double mob[max_np]; props.relperm(1, &s[2*cell], &cell, mob, 0); - double visc[np]; + double visc[max_np]; props.viscosity(1, &p[cell], &z[2*cell], &cell, visc, 0); - double tmob=0; - for(int i=0; i < np; ++i){ + double tmob = 0; + for(int i = 0; i < np; ++i) { mob[i] /= visc[i]; tmob += mob[i]; - } + } const double fracflow = mob[0]/(tmob); well_rate_water += perf_rate*fracflow; - //const double fracflow = mob[0]/(tmob); - //well_rate_water += perf_rate*fracflow; } } data_now.push_back(well_rate_total);