adding a few member variables for MultisegmentWells

and reduce the number of parameters of the related functions by using
the new introduced member variables.

opm/autodiff/BlackoilMultiSegmentModel.hpp

@@ -136,6 +136,7 @@ namespace Opm {
         using Base::cells_;
         using Base::param_;
         using Base::linsolver_;
+        using Base::phaseCondition_;
 
         MultisegmentWells ms_wells_;
 

opm/autodiff/BlackoilMultiSegmentModel_impl.hpp

@@ -71,7 +71,7 @@ namespace Opm {
                   const std::vector<WellMultiSegmentConstPtr>& wells_multisegment)
         : Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver,
                eclState, has_disgas, has_vapoil, terminal_output)
-        , ms_wells_(wells_multisegment, fluid.numPhases())
+        , ms_wells_(wells_multisegment, fluid_, active_, phaseCondition_)
     {
     }
 
@@ -454,7 +454,7 @@ namespace Opm {
             // Compute initial accumulation contributions
             // and well connection pressures.
             asImpl().computeAccum(state0, 0);
-            msWells().computeSegmentFluidProperties(state0, phaseCondition(), active_, fluid_);
+            msWells().computeSegmentFluidProperties(state0);
             const int np = numPhases();
             assert(np == int(msWells().segmentCompSurfVolumeInitial().size()));
             for (int phase = 0; phase < np; ++phase) {
@@ -482,7 +482,7 @@ namespace Opm {
         }
 
         // asImpl().computeSegmentFluidProperties(state);
-        msWells().computeSegmentFluidProperties(state, phaseCondition(), active_, fluid_);
+        msWells().computeSegmentFluidProperties(state);
 
         // asImpl().computeSegmentPressuresDelta(state);
         const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
@@ -504,14 +504,13 @@ namespace Opm {
         const int np = numPhases();
         const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
         const V perf_press_diffs = stdWells().wellPerforationPressureDiffs();
-        msWells().computeWellFlux(state, fluid_.phaseUsage(), active_,
+        msWells().computeWellFlux(state, perf_press_diffs, compi,
-                                  perf_press_diffs, compi,
                                   mob_perfcells, b_perfcells, aliveWells, cq_s);
         asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
         msWells().addWellFluxEq(cq_s, state, residual_);
         asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
         // asImpl().addWellControlEq(state, well_state, aliveWells);
-        msWells().addWellControlEq(state, well_state, aliveWells, active_, residual_);
+        msWells().addWellControlEq(state, well_state, aliveWells, residual_);
     }
 
 
@@ -663,14 +662,13 @@ namespace Opm {
             const int nw = wellsMultiSegment().size();
             const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
             const V perf_press_diffs = stdWells().wellPerforationPressureDiffs();
-            msWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_,
+            msWells().computeWellFlux(wellSolutionState, perf_press_diffs, compi,
-                                      perf_press_diffs, compi,
                                       mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
 
             updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
             msWells().addWellFluxEq(cq_s, wellSolutionState, residual_);
             // addWellControlEq(wellSolutionState, well_state, aliveWells);
-            msWells().addWellControlEq(wellSolutionState, well_state, aliveWells, active_, residual_);
+            msWells().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_);
             converged = Base::getWellConvergence(it);
 
             if (converged) {

opm/autodiff/MultisegmentWells.cpp

@@ -140,10 +140,16 @@ namespace Opm {
 
 
     MultisegmentWells::
-    MultisegmentWells(const std::vector<WellMultiSegmentConstPtr>& wells_ms, const int np)
+    MultisegmentWells(const std::vector<WellMultiSegmentConstPtr>& wells_ms,
+                      const BlackoilPropsAdInterface& fluid_arg,
+                      const std::vector<bool>& active_arg,
+                      const std::vector<PhasePresence>& pc_arg)
       : wells_multisegment_(wells_ms)
       , wops_ms_(wells_ms)
-      , num_phases_(np)
+      , num_phases_(wells_ms.empty()? 0 : wells_ms[0]->numberOfPhases())
+      , fluid_(fluid_arg)
+      , active_(active_arg)
+      , phase_condition_(pc_arg)
       , well_segment_perforation_pressure_diffs_(ADB::null())
       , well_segment_densities_(ADB::null())
       , well_segment_pressures_delta_(ADB::null())

opm/autodiff/MultisegmentWells.hpp

@@ -80,7 +80,10 @@ namespace Opm {
             // ---------  Public methods  ---------
             // TODO: using a vector of WellMultiSegmentConstPtr for now
             // TODO: it should use const Wells or something else later.
-            MultisegmentWells(const std::vector<WellMultiSegmentConstPtr>& wells_multisegment, const int np);
+            MultisegmentWells(const std::vector<WellMultiSegmentConstPtr>& wells_multisegment,
+                              const BlackoilPropsAdInterface& fluid_arg,
+                              const std::vector<bool>& active_arg,
+                              const std::vector<PhasePresence>& pc_arg);
 
             const std::vector<WellMultiSegmentConstPtr>& wells() const;
             const MultisegmentWellOps& wellOps() const;
@@ -121,8 +124,6 @@ namespace Opm {
             template <class SolutionState>
             void
             computeWellFlux(const SolutionState& state,
-                            const Opm::PhaseUsage& pu,
-                            const std::vector<bool>& active,
                             const Vector& well_perforation_pressure_diffs,
                             const DataBlock& compi,
                             const std::vector<ADB>& mob_perfcells,
@@ -135,10 +136,7 @@ namespace Opm {
             // And the surface volume of the components in the segments by dt
             template <class SolutionState>
             void
-            computeSegmentFluidProperties(const SolutionState& state,
+            computeSegmentFluidProperties(const SolutionState& state);
-                                          const std::vector<PhasePresence>& pc,
-                                          const std::vector<bool>& active,
-                                          const BlackoilPropsAdInterface& fluid);
 
             void
             computeSegmentPressuresDelta(const double grav);
@@ -154,7 +152,6 @@ namespace Opm {
             addWellControlEq(const SolutionState& state,
                              const WellState& xw,
                              const Vector& aliveWells,
-                             const std::vector<bool>& active,
                              LinearisedBlackoilResidual& residual);
 
             template <class WellState>
@@ -169,6 +166,9 @@ namespace Opm {
         const int num_phases_;
         int nseg_total_;
         int nperf_total_;
+        const BlackoilPropsAdInterface& fluid_;
+        const std::vector<bool>& active_;
+        const std::vector<PhasePresence>& phase_condition_;
 
         // Pressure correction due to the different depth of the perforation
         // and the cell center of the grid block

opm/autodiff/MultisegmentWells_impl.hpp

@@ -122,8 +122,6 @@ namespace Opm
     void
     MultisegmentWells::
     computeWellFlux(const SolutionState& state,
-                    const Opm::PhaseUsage& pu,
-                    const std::vector<bool>& active,
                     const Vector& well_perforation_pressure_diffs,
                     const DataBlock& compi,
                     const std::vector<ADB>& mob_perfcells,
@@ -141,6 +139,8 @@ namespace Opm
         const int nseg = nseg_total_;
         const int nperf = nperf_total_;
 
+        const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+
         cq_s.resize(np, ADB::null());
 
         {
@@ -201,7 +201,7 @@ namespace Opm
                 cq_ps[phase] = b_perfcells[phase] * cq_p;
             }
 
-            if (active[Oil] && active[Gas]) {
+            if (active_[Oil] && active_[Gas]) {
                 const int oilpos = pu.phase_pos[Oil];
                 const int gaspos = pu.phase_pos[Gas];
                 const ADB cq_psOil = cq_ps[oilpos];
@@ -285,11 +285,11 @@ namespace Opm
 
             for (int phase = 0; phase < np; ++phase) {
                 ADB tmp = cmix_s[phase];
-                if (phase == Oil && active[Gas]) {
+                if (phase == Oil && active_[Gas]) {
                     const int gaspos = pu.phase_pos[Gas];
                     tmp = tmp - rv_perfcells * cmix_s[gaspos] / d;
                 }
-                if (phase == Gas && active[Oil]) {
+                if (phase == Gas && active_[Oil]) {
                     const int oilpos = pu.phase_pos[Oil];
                     tmp = tmp - rs_perfcells * cmix_s[oilpos] / d;
                 }
@@ -313,10 +313,7 @@ namespace Opm
     template <class SolutionState>
     void
     MultisegmentWells::
-    computeSegmentFluidProperties(const SolutionState& state,
+    computeSegmentFluidProperties(const SolutionState& state)
-                                  const std::vector<PhasePresence>& pc,
-                                  const std::vector<bool>& active,
-                                  const BlackoilPropsAdInterface& fluid)
     {
         const int np = numPhases();
         const int nw = wells().size();
@@ -360,37 +357,37 @@ namespace Opm
         // Compute PVT properties for segments.
         std::vector<PhasePresence> segment_cond(nseg_total);
         for (int s = 0; s < nseg_total; ++s) {
-            segment_cond[s] = pc[segment_cells[s]];
+            segment_cond[s] = phase_condition_[segment_cells[s]];
         }
         std::vector<ADB> b_seg(np, ADB::null());
         // Viscosities for different phases
         std::vector<ADB> mu_seg(np, ADB::null());
         ADB rsmax_seg = ADB::null();
         ADB rvmax_seg = ADB::null();
-        const PhaseUsage& pu = fluid.phaseUsage();
+        const PhaseUsage& pu = fluid_.phaseUsage();
         if (pu.phase_used[Water]) {
-            b_seg[pu.phase_pos[Water]] = fluid.bWat(segment_press, segment_temp, segment_cells);
+            b_seg[pu.phase_pos[Water]] = fluid_.bWat(segment_press, segment_temp, segment_cells);
-            mu_seg[pu.phase_pos[Water]] = fluid.muWat(segment_press, segment_temp, segment_cells);
+            mu_seg[pu.phase_pos[Water]] = fluid_.muWat(segment_press, segment_temp, segment_cells);
         }
-        assert(active[Oil]);
+        assert(active_[Oil]);
         const ADB segment_so = subset(state.saturation[pu.phase_pos[Oil]], segment_cells);
         if (pu.phase_used[Oil]) {
             const ADB segment_rs = subset(state.rs, segment_cells);
-            b_seg[pu.phase_pos[Oil]] = fluid.bOil(segment_press, segment_temp, segment_rs,
+            b_seg[pu.phase_pos[Oil]] = fluid_.bOil(segment_press, segment_temp, segment_rs,
                                                    segment_cond, segment_cells);
             // rsmax_seg = fluidRsSat(segment_press, segment_so, segment_cells);
-            rsmax_seg = fluid.rsSat(segment_press, segment_so, segment_cells);
+            rsmax_seg = fluid_.rsSat(segment_press, segment_so, segment_cells);
-            mu_seg[pu.phase_pos[Oil]] = fluid.muOil(segment_press, segment_temp, segment_rs,
+            mu_seg[pu.phase_pos[Oil]] = fluid_.muOil(segment_press, segment_temp, segment_rs,
                                                      segment_cond, segment_cells);
         }
-        assert(active[Gas]);
+        assert(active_[Gas]);
         if (pu.phase_used[Gas]) {
             const ADB segment_rv = subset(state.rv, segment_cells);
-            b_seg[pu.phase_pos[Gas]] = fluid.bGas(segment_press, segment_temp, segment_rv,
+            b_seg[pu.phase_pos[Gas]] = fluid_.bGas(segment_press, segment_temp, segment_rv,
                                                    segment_cond, segment_cells);
             // rvmax_seg = fluidRvSat(segment_press, segment_so, segment_cells);
-            rvmax_seg = fluid.rvSat(segment_press, segment_so, segment_cells);
+            rvmax_seg = fluid_.rvSat(segment_press, segment_so, segment_cells);
-            mu_seg[pu.phase_pos[Gas]] = fluid.muGas(segment_press, segment_temp, segment_rv,
+            mu_seg[pu.phase_pos[Gas]] = fluid_.muGas(segment_press, segment_temp, segment_rv,
                                                    segment_cond, segment_cells);
         }
 
@@ -444,7 +441,7 @@ namespace Opm
         ADB big_values = ADB::constant(Vector::Constant(nseg_total, 1.e100));
         ADB mix_gas_oil = non_zero_mix_oilpos.select(mix[gaspos] / mix[oilpos], big_values);
         ADB mix_oil_gas = non_zero_mix_gaspos.select(mix[oilpos] / mix[gaspos], big_values);
-        if (active[Oil]) {
+        if (active_[Oil]) {
             Vector selectorUnderRsmax = Vector::Zero(nseg_total);
             Vector selectorAboveRsmax = Vector::Zero(nseg_total);
             for (int s = 0; s < nseg_total; ++s) {
@@ -456,7 +453,7 @@ namespace Opm
             }
             rs = non_zero_mix_oilpos.select(selectorAboveRsmax * rsmax_seg + selectorUnderRsmax * mix_gas_oil, rs);
         }
-        if (active[Gas]) {
+        if (active_[Gas]) {
             Vector selectorUnderRvmax = Vector::Zero(nseg_total);
             Vector selectorAboveRvmax = Vector::Zero(nseg_total);
             for (int s = 0; s < nseg_total; ++s) {
@@ -474,7 +471,7 @@ namespace Opm
         for (int phase = 0; phase < np; ++phase) {
             x[phase] = mix[phase];
         }
-        if (active[Gas] && active[Oil]) {
+        if (active_[Gas] && active_[Oil]) {
             x[gaspos] = (mix[gaspos] - mix[oilpos] * rs) / (Vector::Ones(nseg_total) - rs * rv);
             x[oilpos] = (mix[oilpos] - mix[gaspos] * rv) / (Vector::Ones(nseg_total) - rs * rv);
         }
@@ -488,7 +485,7 @@ namespace Opm
         // Compute segment densities.
         ADB dens = ADB::constant(Vector::Zero(nseg_total));
         for (int phase = 0; phase < np; ++phase) {
-            const Vector surface_density = fluid.surfaceDensity(phase, segment_cells);
+            const Vector surface_density = fluid_.surfaceDensity(phase, segment_cells);
             dens += surface_density * mix[phase];
         }
         well_segment_densities_ = dens / volrat;
@@ -504,7 +501,7 @@ namespace Opm
         segment_mass_flow_rates_ = ADB::constant(Vector::Zero(nseg_total));
         for (int phase = 0; phase < np; ++phase) {
             // TODO: how to remove one repeated surfaceDensity()
-            const Vector surface_density = fluid.surfaceDensity(phase, segment_cells);
+            const Vector surface_density = fluid_.surfaceDensity(phase, segment_cells);
             segment_mass_flow_rates_ += surface_density * segqs[phase];
         }
 
@@ -580,7 +577,6 @@ namespace Opm
     addWellControlEq(const SolutionState& state,
                      const WellState& xw,
                      const Vector& aliveWells,
-                     const std::vector<bool>& active,
                      LinearisedBlackoilResidual& residual)
     {
         // the name of the function is a a little misleading.
@@ -596,13 +592,13 @@ namespace Opm
         ADB liquid = ADB::constant(Vector::Zero(nseg_total));
         ADB vapour = ADB::constant(Vector::Zero(nseg_total));
 
-        if (active[Water]) {
+        if (active_[Water]) {
             aqua += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Aqua * nseg_total));
         }
-        if (active[Oil]) {
+        if (active_[Oil]) {
             liquid += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Liquid * nseg_total));
         }
-        if (active[Gas]) {
+        if (active_[Gas]) {
             vapour += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Vapour * nseg_total));
         }