From 8367eaacb30d11689abdb5c1f2fdf861a066239e Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 31 Mar 2016 15:19:46 +0200 Subject: [PATCH] first attempt to introduce StandardWells struct. To replace the const Wells struct in BlackoilModelBase. Only testing with Flow for the moment. Will update other Flow siblings later. --- opm/autodiff/BlackoilModelBase.hpp | 24 +++- opm/autodiff/BlackoilModelBase_impl.hpp | 146 ++++++++++++++---------- 2 files changed, 103 insertions(+), 67 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 8b388ac10..449d83a63 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -276,13 +276,23 @@ namespace Opm { std::vector well_cells; // the set of perforated cells }; + struct StandardWells { + // keeping the underline, later they will be private members + StandardWells(const Wells* wells); + const Wells& wells() const; + bool wells_active_; + const Wells* wells_; + V well_perforation_densities_; //Density of each well perforation + V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. + }; + // --------- Data members --------- const Grid& grid_; const BlackoilPropsAdInterface& fluid_; const DerivedGeology& geo_; const RockCompressibility* rock_comp_props_; - const Wells* wells_; + StandardWells std_wells_; VFPProperties vfp_properties_; const NewtonIterationBlackoilInterface& linsolver_; // For each canonical phase -> true if active @@ -305,8 +315,6 @@ namespace Opm { V isRs_; V isRv_; V isSg_; - V well_perforation_densities_; //Density of each well perforation - V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. LinearisedBlackoilResidual residual_; @@ -339,11 +347,15 @@ namespace Opm { } // return true if wells are available in the reservoir - bool wellsActive() const { return wells_active_; } + bool wellsActive() const { return std_wells_.wells_active_; } // return true if wells are available on this process - bool localWellsActive() const { return wells_ ? (wells_->number_of_wells > 0 ) : false; } + bool localWellsActive() const { return std_wells_.wells_ ? (std_wells_.wells_->number_of_wells > 0 ) : false; } // return wells object - const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; } + // const Wells& wells () const { assert( bool(std_wells_.wells_ != 0) ); return *(std_wells_.wells_); } + + // return the StandardWells object + StandardWells& stdWells() { return std_wells_; } + const StandardWells& stdWells() const { return std_wells_; } int numWellVars() const; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index c9703992b..3f034e78a 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -159,14 +159,14 @@ namespace detail { , fluid_ (fluid) , geo_ (geo) , rock_comp_props_(rock_comp_props) - , wells_ (wells_arg) + , std_wells_ (wells_arg) , vfp_properties_(eclState->getTableManager()->getVFPInjTables(), eclState->getTableManager()->getVFPProdTables()) , linsolver_ (linsolver) , active_(detail::activePhases(fluid.phaseUsage())) , canph_ (detail::active2Canonical(fluid.phaseUsage())) , cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid))) , ops_ (grid, geo.nnc()) - , wops_ (wells_) + , wops_ (wells_arg) , has_disgas_(has_disgas) , has_vapoil_(has_vapoil) , param_( param ) @@ -195,9 +195,9 @@ namespace detail { // Only rank 0 does print to std::cout if terminal_output is enabled terminal_output_ = (info.communicator().rank()==0); } - int local_number_of_wells = wells_ ? wells_->number_of_wells : 0; + int local_number_of_wells = std_wells_.wells_ ? std_wells_.wells_->number_of_wells : 0; int global_number_of_wells = info.communicator().sum(local_number_of_wells); - wells_active_ = ( wells_ && global_number_of_wells > 0 ); + std_wells_.wells_active_ = ( std_wells_.wells_ && global_number_of_wells > 0 ); // Compute the global number of cells std::vector v( Opm::AutoDiffGrid::numCells(grid_), 1); global_nc_ = 0; @@ -205,7 +205,7 @@ namespace detail { }else #endif { - wells_active_ = ( wells_ && wells_->number_of_wells > 0 ); + std_wells_.wells_active_ = ( std_wells_.wells_ && std_wells_.wells_->number_of_wells > 0 ); global_nc_ = Opm::AutoDiffGrid::numCells(grid_); } } @@ -457,12 +457,36 @@ namespace detail { + template + BlackoilModelBase:: + StandardWells::StandardWells(const Wells* wells) + : wells_(wells) + { + } + + + + + + template + const Wells& + BlackoilModelBase:: + StandardWells::wells() const + { + assert(wells_ != 0); + return *(wells_); + } + + + + + template int BlackoilModelBase::numWellVars() const { // For each well, we have a bhp variable, and one flux per phase. - const int nw = localWellsActive() ? wells().number_of_wells : 0; + const int nw = localWellsActive() ? stdWells().wells().number_of_wells : 0; return (numPhases() + 1) * nw; } @@ -582,8 +606,8 @@ namespace detail { { // Need to reshuffle well rates, from phase running fastest // to wells running fastest. - const int nw = wells().number_of_wells; - const int np = wells().number_of_phases; + const int nw = stdWells().wells().number_of_wells; + const int np = stdWells().wells().number_of_phases; // The transpose() below switches the ordering. const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); @@ -803,15 +827,15 @@ namespace detail { std::vector& rvmax_perf, std::vector& surf_dens_perf) { - const int nperf = wells().well_connpos[wells().number_of_wells]; - const int nw = wells().number_of_wells; + const int nperf = stdWells().wells().well_connpos[stdWells().wells().number_of_wells]; + const int nw = stdWells().wells().number_of_wells; // Compute the average pressure in each well block const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); V avg_press = perf_press*0; for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + for (int perf = stdWells().wells().well_connpos[w]; perf < stdWells().wells().well_connpos[w+1]; ++perf) { + const double p_above = perf == stdWells().wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; const double p_avg = (perf_press[perf] + p_above)/2; avg_press[perf] = p_avg; } @@ -891,7 +915,7 @@ namespace detail { // 2. Compute densities std::vector cd = WellDensitySegmented::computeConnectionDensities( - wells(), xw, fluid_.phaseUsage(), + stdWells().wells(), xw, fluid_.phaseUsage(), b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); const int nperf = wells().well_connpos[wells().number_of_wells]; @@ -908,11 +932,11 @@ namespace detail { // 3. Compute pressure deltas std::vector cdp = WellDensitySegmented::computeConnectionPressureDelta( - wells(), perf_depth, cd, grav); + stdWells().wells(), perf_depth, cd, grav); // 4. Store the results - well_perforation_densities_ = Eigen::Map(cd.data(), nperf); - well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); + stdWells().well_perforation_densities_ = Eigen::Map(cd.data(), nperf); + stdWells().well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); } @@ -1159,15 +1183,15 @@ namespace detail { { if( ! localWellsActive() ) return ; - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; + const int np = stdWells().wells().number_of_phases; + const int nw = stdWells().wells().number_of_wells; + const int nperf = stdWells().wells().well_connpos[nw]; const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - V Tw = Eigen::Map(wells().WI, nperf); + V Tw = Eigen::Map(stdWells().wells().WI, nperf); const std::vector& well_cells = wops_.well_cells; // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = well_perforation_pressure_diffs_; + const V& cdp = std_wells_.well_perforation_pressure_diffs_; // Extract needed quantities for the perforation cells const ADB& p_perfcells = subset(state.pressure, well_cells); const ADB& rv_perfcells = subset(state.rv, well_cells); @@ -1197,15 +1221,15 @@ namespace detail { const V numInjectingPerforations = (wops_.p2w * ADB::constant(selectInjectingPerforations)).value(); const V numProducingPerforations = (wops_.p2w * ADB::constant(selectProducingPerforations)).value(); for (int w = 0; w < nw; ++w) { - if (!wells().allow_cf[w]) { - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + if (!stdWells().wells().allow_cf[w]) { + for (int perf = stdWells().wells().well_connpos[w] ; perf < stdWells().wells().well_connpos[w+1]; ++perf) { // Crossflow is not allowed; reverse flow is prevented. // At least one of the perforation must be open in order to have a meeningful // equation to solve. For the special case where all perforations have reverse flow, // and the target rate is non-zero all of the perforations are keept open. - if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { + if (stdWells().wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { selectProducingPerforations[perf] = 0.0; - } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ + } else if (stdWells().wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ selectInjectingPerforations[perf] = 0.0; } } @@ -1242,7 +1266,7 @@ namespace detail { // and the well injection rates. // compute avg. and total wellbore phase volumetric rates at standard conds - const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); + const DataBlock compi = Eigen::Map(stdWells().wells().comp_frac, nw, np); std::vector wbq(np, ADB::null()); ADB wbqt = ADB::constant(V::Zero(nw)); for (int phase = 0; phase < np; ++phase) { @@ -1313,9 +1337,9 @@ namespace detail { } // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; + const int np = stdWells().wells().number_of_phases; + const int nw = stdWells().wells().number_of_wells; + const int nperf = stdWells().wells().well_connpos[nw]; V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); for (int phase = 1; phase < np; ++phase) { cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); @@ -1323,7 +1347,7 @@ namespace detail { xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); // Update the perforation pressures. - const V& cdp = well_perforation_pressure_diffs_; + const V& cdp = std_wells_.well_perforation_pressure_diffs_; const V perfpressure = (wops_.w2p * state.bhp.value().matrix()).array() + cdp; xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); } @@ -1343,8 +1367,8 @@ namespace detail { return; } - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; + const int np = stdWells().wells().number_of_phases; + const int nw = stdWells().wells().number_of_wells; ADB qs = state.qs; for (int phase = 0; phase < np; ++phase) { qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); @@ -1501,10 +1525,10 @@ namespace detail { return false; } - const int nw = wells().number_of_wells; + const int nw = stdWells().wells().number_of_wells; //Loop over all wells for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; + const WellControls* wc = stdWells().wells().ctrls[w]; const int nwc = well_controls_get_num(wc); @@ -1530,12 +1554,12 @@ namespace detail { std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; // Find, for each well, if any constraints are broken. If so, // switch control to first broken constraint. - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; + const int np = stdWells().wells().number_of_phases; + const int nw = stdWells().wells().number_of_wells; const Opm::PhaseUsage& pu = fluid_.phaseUsage(); #pragma omp parallel for schedule(dynamic) for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; + const WellControls* wc = stdWells().wells().ctrls[w]; // The current control in the well state overrides // the current control set in the Wells struct, which // is instead treated as a default. @@ -1554,7 +1578,7 @@ namespace detail { } if (detail::constraintBroken( xw.bhp(), xw.thp(), xw.wellRates(), - w, np, wells().type[w], wc, ctrl_index)) { + w, np, stdWells().wells().type[w], wc, ctrl_index)) { // ctrl_index will be the index of the broken constraint after the loop. break; } @@ -1563,7 +1587,7 @@ namespace detail { // Constraint number ctrl_index was broken, switch to it. if (terminal_output_) { - std::cout << "Switching control mode for well " << wells().name[w] + std::cout << "Switching control mode for well " << stdWells().wells().name[w] << " from " << modestring[well_controls_iget_type(wc, current)] << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; } @@ -1603,19 +1627,19 @@ namespace detail { const double& alq = well_controls_iget_alq(wc, current); //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[w]; + const WellType& well_type = stdWells().wells().type[w]; if (well_type == INJECTOR) { double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), - well_perforation_densities_, gravity); + stdWells().wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), + stdWells().well_perforation_densities_, gravity); xw.bhp()[w] = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; } else if (well_type == PRODUCER) { double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), - well_perforation_densities_, gravity); + stdWells().wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), + std_wells_.well_perforation_densities_, gravity); xw.bhp()[w] = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; } @@ -1681,7 +1705,7 @@ namespace detail { WellState& well_state) { V aliveWells; - const int np = wells().number_of_phases; + const int np = stdWells().wells().number_of_phases; std::vector cq_s(np, ADB::null()); std::vector indices = variableWellStateIndices(); SolutionState state0 = state; @@ -1746,7 +1770,7 @@ namespace detail { if ( terminal_output_ ) { std::cout << "well converged iter: " << it << std::endl; } - const int nw = wells().number_of_wells; + const int nw = stdWells().wells().number_of_wells; { // We will set the bhp primary variable to the new ones, // but we do not change the derivatives here. @@ -1786,8 +1810,8 @@ namespace detail { { if( ! localWellsActive() ) return; - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; + const int np = stdWells().wells().number_of_phases; + const int nw = stdWells().wells().number_of_wells; ADB aqua = ADB::constant(ADB::V::Zero(nw)); ADB liquid = ADB::constant(ADB::V::Zero(nw)); @@ -1828,7 +1852,7 @@ namespace detail { //Run through all wells to calculate BHP/RATE targets //and gather info about current control for (int w = 0; w < nw; ++w) { - auto wc = wells().ctrls[w]; + auto wc = stdWells().wells().ctrls[w]; // The current control in the well state overrides // the current control set in the Wells struct, which @@ -1846,13 +1870,13 @@ namespace detail { case THP: { - const int perf = wells().well_connpos[w]; - rho_v[w] = well_perforation_densities_[perf]; + const int perf = stdWells().wells().well_connpos[w]; + rho_v[w] = std_wells_.well_perforation_densities_[perf]; const int table_id = well_controls_iget_vfp(wc, current); const double target = well_controls_iget_target(wc, current); - const WellType& well_type = wells().type[w]; + const WellType& well_type = stdWells().wells().type[w]; if (well_type == INJECTOR) { inj_table_id[w] = table_id; thp_inj_target_v[w] = target; @@ -1910,7 +1934,7 @@ namespace detail { //Perform hydrostatic correction to computed targets double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - const ADB::V dp_v = detail::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, well_perforation_densities_, gravity); + const ADB::V dp_v = detail::computeHydrostaticCorrection(stdWells().wells(), vfp_ref_depth_v, std_wells_.well_perforation_densities_, gravity); const ADB dp = ADB::constant(dp_v); const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw); const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw); @@ -2286,8 +2310,8 @@ namespace detail { if( localWellsActive() ) { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; + const int np = stdWells().wells().number_of_phases; + const int nw = stdWells().wells().number_of_wells; // Extract parts of dwells corresponding to each part. int varstart = 0; @@ -2323,7 +2347,7 @@ namespace detail { //Loop over all wells #pragma omp parallel for schedule(static) for (int w=0; wgetTable(table_id)->getDatumDepth(), - well_perforation_densities_, gravity); + stdWells().wells(), w, vfp_properties_.getInj()->getTable(table_id)->getDatumDepth(), + std_wells_.well_perforation_densities_, gravity); well_state.thp()[w] = vfp_properties_.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp); } else if (well_type == PRODUCER) { double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(), - well_perforation_densities_, gravity); + stdWells().wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(), + std_wells_.well_perforation_densities_, gravity); well_state.thp()[w] = vfp_properties_.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq); }