From 8a52dd743bb30bbc70d8d0c779059409a836d4ba Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 1 Apr 2016 18:19:15 +0200 Subject: [PATCH] making StandardWells a simple class. --- opm/autodiff/BlackoilModelBase.hpp | 40 ++++-- opm/autodiff/BlackoilModelBase_impl.hpp | 125 ++++++++++++++---- .../BlackoilMultiSegmentModel_impl.hpp | 10 +- .../BlackoilPolymerModel_impl.hpp | 4 +- 4 files changed, 136 insertions(+), 43 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index c6a7107c1..099aeda4e 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -269,26 +269,42 @@ namespace Opm { ADB mob; // Phase mobility (per cell) }; - struct WellOps { - WellOps(const Wells* wells); - Eigen::SparseMatrix w2p; // well -> perf (scatter) - Eigen::SparseMatrix p2w; // perf -> well (gather) - std::vector well_cells; // the set of perforated cells - }; + class StandardWells { + protected: + struct WellOps { + WellOps(const Wells* wells); + Eigen::SparseMatrix w2p; // well -> perf (scatter) + Eigen::SparseMatrix p2w; // perf -> well (gather) + std::vector well_cells; // the set of perforated cells + }; - struct StandardWells { - // keeping the underline, later they will be private members + public: StandardWells(const Wells* wells); + const Wells& wells() const; + // return true if wells are available in the reservoir bool wellsActive() const; + bool& wellsActive(); // return true if wells are available on this process bool localWellsActive() const; + + const WellOps& wellOps() const; + + //Density of each well perforation + V& wellPerforationDensities(); + const V& wellPerforationDensities() const; + + // Diff to bhp for each well perforation. + V& wellPerforationPressureDiffs(); + const V& wellPerforationPressureDiffs() const; + + protected: bool wells_active_; - const Wells* wells_; - const WellOps wops_; - V well_perforation_densities_; //Density of each well perforation - V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. + const Wells* wells_; + const WellOps wops_; + V well_perforation_densities_; + V well_perforation_pressure_diffs_; }; // --------- Data members --------- diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 4bae2594c..62276eb70 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -194,9 +194,11 @@ 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 = std_wells_.wells_ ? std_wells_.wells_->number_of_wells : 0; + // int local_number_of_wells = std_wells_.wells_ ? std_wells_.wells_->number_of_wells : 0; + int local_number_of_wells = stdWells().localWellsActive() ? stdWells().wells().number_of_wells : 0; int global_number_of_wells = info.communicator().sum(local_number_of_wells); - std_wells_.wells_active_ = ( std_wells_.wells_ && global_number_of_wells > 0 ); + // std_wells_.wells_active_ = ( std_wells_.wells_ && global_number_of_wells > 0 ); + stdWells().wellsActive() = ( stdWells().localWellsActive() && global_number_of_wells > 0 ); // Compute the global number of cells std::vector v( Opm::AutoDiffGrid::numCells(grid_), 1); global_nc_ = 0; @@ -204,7 +206,7 @@ namespace detail { }else #endif { - std_wells_.wells_active_ = ( std_wells_.wells_ && std_wells_.wells_->number_of_wells > 0 ); + stdWells().wellsActive() = ( stdWells().localWellsActive() && stdWells().wells().number_of_wells > 0 ); global_nc_ = Opm::AutoDiffGrid::numCells(grid_); } } @@ -419,6 +421,7 @@ namespace detail { template BlackoilModelBase:: + StandardWells:: WellOps::WellOps(const Wells* wells) : w2p(), p2w(), @@ -461,6 +464,8 @@ namespace detail { StandardWells::StandardWells(const Wells* wells_arg) : wells_(wells_arg) , wops_(wells_arg) + , well_perforation_densities_(V()) + , well_perforation_pressure_diffs_(V()) { } @@ -493,6 +498,18 @@ namespace detail { + template + bool& + BlackoilModelBase:: + StandardWells::wellsActive() + { + return wells_active_; + } + + + + + template bool BlackoilModelBase:: @@ -505,6 +522,66 @@ namespace detail { + template + const typename BlackoilModelBase::StandardWells::WellOps& + BlackoilModelBase:: + StandardWells::wellOps() const + { + return wops_; + } + + + + + + template + V& + BlackoilModelBase:: + StandardWells::wellPerforationDensities() + { + return well_perforation_densities_; + } + + + + + + template + const V& + BlackoilModelBase:: + StandardWells::wellPerforationDensities() const + { + return well_perforation_densities_; + } + + + + + + template + V& + BlackoilModelBase:: + StandardWells::wellPerforationPressureDiffs() + { + return well_perforation_pressure_diffs_; + } + + + + + + template + const V& + BlackoilModelBase:: + StandardWells::wellPerforationPressureDiffs() const + { + return well_perforation_pressure_diffs_; + } + + + + + template int BlackoilModelBase::numWellVars() const @@ -865,7 +942,7 @@ namespace detail { } } - const std::vector& well_cells = stdWells().wops_.well_cells; + const std::vector& well_cells = stdWells().wellOps().well_cells; // Use cell values for the temperature as the wells don't knows its temperature yet. const ADB perf_temp = subset(state.temperature, well_cells); @@ -959,8 +1036,8 @@ namespace detail { stdWells().wells(), perf_depth, cd, grav); // 4. Store the results - stdWells().well_perforation_densities_ = Eigen::Map(cd.data(), nperf); - stdWells().well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); + stdWells().wellPerforationDensities() = Eigen::Map(cd.data(), nperf); + stdWells().wellPerforationPressureDiffs() = Eigen::Map(cdp.data(), nperf); } @@ -1159,7 +1236,7 @@ namespace detail { const int nc = Opm::AutoDiffGrid::numCells(grid_); const int np = asImpl().numPhases(); for (int phase = 0; phase < np; ++phase) { - residual_.material_balance_eq[phase] -= superset(cq_s[phase], stdWells().wops_.well_cells, nc); + residual_.material_balance_eq[phase] -= superset(cq_s[phase], stdWells().wellOps().well_cells, nc); } } @@ -1182,7 +1259,7 @@ namespace detail { return; } else { const int np = asImpl().numPhases(); - const std::vector& well_cells = stdWells().wops_.well_cells; + const std::vector& well_cells = stdWells().wellOps().well_cells; mob_perfcells.resize(np, ADB::null()); b_perfcells.resize(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { @@ -1212,17 +1289,17 @@ namespace detail { const int nperf = stdWells().wells().well_connpos[nw]; const Opm::PhaseUsage& pu = fluid_.phaseUsage(); V Tw = Eigen::Map(stdWells().wells().WI, nperf); - const std::vector& well_cells = stdWells().wops_.well_cells; + const std::vector& well_cells = stdWells().wellOps().well_cells; // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = std_wells_.well_perforation_pressure_diffs_; + const V& cdp = stdWells().wellPerforationPressureDiffs(); // 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); const ADB& rs_perfcells = subset(state.rs, well_cells); // Perforation pressure - const ADB perfpressure = (stdWells().wops_.w2p * state.bhp) + cdp; + const ADB perfpressure = (stdWells().wellOps().w2p * state.bhp) + cdp; // Pressure drawdown (also used to determine direction of flow) const ADB drawdown = p_perfcells - perfpressure; @@ -1242,8 +1319,8 @@ namespace detail { } // Handle cross flow - const V numInjectingPerforations = (stdWells().wops_.p2w * ADB::constant(selectInjectingPerforations)).value(); - const V numProducingPerforations = (stdWells().wops_.p2w * ADB::constant(selectProducingPerforations)).value(); + const V numInjectingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); + const V numProducingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); for (int w = 0; w < nw; ++w) { if (!stdWells().wells().allow_cf[w]) { for (int perf = stdWells().wells().well_connpos[w] ; perf < stdWells().wells().well_connpos[w+1]; ++perf) { @@ -1294,7 +1371,7 @@ namespace detail { std::vector wbq(np, ADB::null()); ADB wbqt = ADB::constant(V::Zero(nw)); for (int phase = 0; phase < np; ++phase) { - const ADB& q_ps = stdWells().wops_.p2w * cq_ps[phase]; + const ADB& q_ps = stdWells().wellOps().p2w * cq_ps[phase]; const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); const int pos = pu.phase_pos[phase]; @@ -1306,7 +1383,7 @@ namespace detail { std::vector cmix_s(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { const int pos = pu.phase_pos[phase]; - cmix_s[phase] = stdWells().wops_.w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); + cmix_s[phase] = stdWells().wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); } // compute volume ratio between connection at standard conditions @@ -1371,8 +1448,8 @@ namespace detail { xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); // Update the perforation pressures. - const V& cdp = std_wells_.well_perforation_pressure_diffs_; - const V perfpressure = (stdWells().wops_.w2p * state.bhp.value().matrix()).array() + cdp; + const V& cdp = stdWells().wellPerforationPressureDiffs(); + const V perfpressure = (stdWells().wellOps().w2p * state.bhp.value().matrix()).array() + cdp; xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); } @@ -1395,7 +1472,7 @@ namespace detail { const int nw = stdWells().wells().number_of_wells; ADB qs = state.qs; for (int phase = 0; phase < np; ++phase) { - qs -= superset(stdWells().wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); + qs -= superset(stdWells().wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); } @@ -1656,14 +1733,14 @@ namespace detail { if (well_type == INJECTOR) { double dp = detail::computeHydrostaticCorrection( stdWells().wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), - stdWells().well_perforation_densities_, gravity); + stdWells().wellPerforationDensities(), gravity); xw.bhp()[w] = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; } else if (well_type == PRODUCER) { double dp = detail::computeHydrostaticCorrection( stdWells().wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), - std_wells_.well_perforation_densities_, gravity); + std_wells_.wellPerforationDensities(), gravity); xw.bhp()[w] = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; } @@ -1895,7 +1972,7 @@ namespace detail { case THP: { const int perf = stdWells().wells().well_connpos[w]; - rho_v[w] = std_wells_.well_perforation_densities_[perf]; + rho_v[w] = stdWells().wellPerforationDensities()[perf]; const int table_id = well_controls_iget_vfp(wc, current); const double target = well_controls_iget_target(wc, current); @@ -1958,7 +2035,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(stdWells().wells(), vfp_ref_depth_v, std_wells_.well_perforation_densities_, gravity); + const ADB::V dp_v = detail::computeHydrostaticCorrection(stdWells().wells(), vfp_ref_depth_v, stdWells().wellPerforationDensities(), 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); @@ -2399,14 +2476,14 @@ namespace detail { if (well_type == INJECTOR) { double dp = detail::computeHydrostaticCorrection( stdWells().wells(), w, vfp_properties_.getInj()->getTable(table_id)->getDatumDepth(), - std_wells_.well_perforation_densities_, gravity); + stdWells().wellPerforationDensities(), 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( stdWells().wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(), - std_wells_.well_perforation_densities_, gravity); + stdWells().wellPerforationDensities(), gravity); well_state.thp()[w] = vfp_properties_.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq); } diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index d6a4ca8cc..9afcbd9b2 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -370,7 +370,7 @@ namespace Opm { std::vector& well_cells = wops_ms_.well_cells; - stdWells().well_perforation_densities_ = V::Zero(nperf_total); + stdWells().wellPerforationDensities() = V::Zero(nperf_total); const V perf_press = Eigen::Map(xw.perfPress().data(), nperf_total); @@ -471,8 +471,8 @@ namespace Opm { stdWells().wells(), perf_cell_depth, cd, grav); // 4. Store the results - stdWells().well_perforation_densities_ = Eigen::Map(cd.data(), nperf_total); // This one is not useful for segmented wells at all - stdWells().well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf_total); + stdWells().wellPerforationDensities() = Eigen::Map(cd.data(), nperf_total); // This one is not useful for segmented wells at all + stdWells().wellPerforationPressureDiffs() = Eigen::Map(cdp.data(), nperf_total); if ( !wops_ms_.has_multisegment_wells ) { well_perforation_cell_densities_ = V::Zero(nperf_total); @@ -698,7 +698,7 @@ namespace Opm { // Compute drawdown. ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_, - ADB::constant(stdWells().well_perforation_pressure_diffs_)); + ADB::constant( stdWells().wellPerforationPressureDiffs() )); const V h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, V::Zero(nperf)); // Special handling for when we are called from solveWellEq(). @@ -862,7 +862,7 @@ namespace Opm { // we need th concept of preforation pressures xw.perfPress().resize(nperf_total, -1.e100); - const V& cdp = stdWells().well_perforation_pressure_diffs_; + const V& cdp = stdWells().wellPerforationPressureDiffs(); int start_segment = 0; int start_perforation = 0; for (int i = 0; i < nw; ++i) { diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 72036351c..cea2b8ba6 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -731,8 +731,8 @@ namespace Opm { ADB b_perfcells = subset(rq_[water_pos].b, well_cells); const ADB& p_perfcells = subset(state.pressure, well_cells); - const V& cdp = stdWells().well_perforation_pressure_diffs_; - const ADB perfpressure = (stdWells().wops_.w2p * state.bhp) + cdp; + const V& cdp = stdWells().wellPerforationPressureDiffs(); + const ADB perfpressure = (stdWells().wellOps().w2p * state.bhp) + cdp; // Pressure drawdown (also used to determine direction of flow) const ADB drawdown = p_perfcells - perfpressure;