update the reservoir volume in WellState

opm/autodiff/BlackoilOutputEbos.hpp

@@ -41,7 +41,6 @@
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
-#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 #include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
 
 #include <string>

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -180,6 +180,12 @@ namespace Opm {
     void
     BlackoilWellModel<TypeTag>::
     timeStepSucceeded() {
+        // TODO: when necessary
+        rateConverter_->template defineState<ElementContext>(ebosSimulator_);
+        for (const auto& well : well_container_) {
+            well->calculateReservoirRates(well_state_);
+        }
+
         previous_well_state_ = well_state_;
     }
 

opm/autodiff/RateConverter.hpp

@@ -640,6 +640,86 @@ namespace Opm {
             }
 
 
+            /**
+             * Converting surface volume rates to reservoir voidage rates
+             *
+             * \tparam Rates Type representing contiguous collection
+             * of surface-to-reservoir conversion coefficients.  Must
+             * support direct indexing through \code operator[]()
+             * \endcode.
+             *
+             *
+             * \param[in] r Fluid-in-place region of the well
+             * \param[in] pvtRegionIdx PVT region of the well
+             * \param[in] surface_rates surface voluem rates for
+             * all active phases
+             *
+             * \param[out] voidage_rates reservoir volume rates for
+             * all active phases
+             */
+            template <class Rates >
+            void
+            calcReservoirVoidageRates(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates,
+                                      Rates& voidage_rates) const
+            {
+                assert(voidage_rates.size() == surface_rates.size());
+
+                const auto& pu = phaseUsage_;
+                const auto& ra = attr_.attributes(r);
+                const double p = ra.pressure;
+                const double T = ra.temperature;
+
+                const int   iw = Details::PhasePos::water(pu);
+                const int   io = Details::PhasePos::oil  (pu);
+                const int   ig = Details::PhasePos::gas  (pu);
+
+                if (Details::PhaseUsed::water(pu)) {
+                    // q[w]_r = q[w]_s / bw
+
+                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p);
+
+                    voidage_rates[iw] = surface_rates[iw] / bw;
+                }
+
+                // Determinant of 'R' matrix
+                const double detR = 1.0 - (ra.rs * ra.rv);
+
+                if (Details::PhaseUsed::oil(pu)) {
+                    // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
+
+                    const double Rs = ra.rs;
+                    const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
+                    const double den = bo * detR;
+
+                    voidage_rates[io] = surface_rates[io];
+
+                    if (Details::PhaseUsed::gas(pu)) {
+                        const double Rv = ra.rv;
+                        voidage_rates[io] -= Rv * surface_rates[ig];
+                    }
+
+                    voidage_rates[io] /= den;
+                }
+
+                if (Details::PhaseUsed::gas(pu)) {
+                    // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
+
+                    const double Rv = ra.rv;
+                    const double bg  = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
+                    const double den = bg * detR;
+
+                    voidage_rates[ig] = surface_rates[ig];
+
+                    if (Details::PhaseUsed::oil(pu)) {
+                        const double Rs = ra.rs;
+                        voidage_rates[ig] -= Rs * surface_rates[io];
+                    }
+
+                    voidage_rates[ig] /= den;
+                }
+            }
+
+
             /**
              * Compute coefficients for surface-to-reservoir voidage
              * conversion for solvent.

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -288,7 +288,8 @@ public:
             // write simulation state at the report stage
             Dune::Timer perfTimer;
             perfTimer.start();
-            const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;            
+            const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
+
             output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model(), false, nextstep, report);
             report.output_write_time += perfTimer.stop();
 

opm/autodiff/StandardWell_impl.hpp

@@ -563,9 +563,9 @@ namespace Opm
 
         const EvalWell& bhp = getBhp();
 
-        // the solution gas rate and solution oil rate needs to be updated for well_state.
-        std::fill(well_state.wellVaporizedOilRates().begin(), well_state.wellVaporizedOilRates().end(), 0.0);
-        std::fill(well_state.wellDissolvedGasRates().begin(), well_state.wellDissolvedGasRates().end(), 0.0);
+        // the solution gas rate and solution oil rate needs to be reset to be zero for well_state.
+        well_state.wellVaporizedOilRates()[index_of_well_] = 0.;
+        well_state.wellDissolvedGasRates()[index_of_well_] = 0.;
 
         for (int perf = 0; perf < number_of_perforations_; ++perf) {
 

opm/autodiff/WellInterface.hpp

@@ -202,6 +202,9 @@ namespace Opm
         virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
                                                  const WellState& well_state) = 0; // should be const?
 
+        // updating the voidage rates in well_state when requested
+        void calculateReservoirRates(WellState& well_state) const;
+
     protected:
 
         // to indicate a invalid connection
@@ -318,7 +321,6 @@ namespace Opm
 
         double scalingFactor(const int comp_idx) const;
 
-
     };
 
 }

opm/autodiff/WellInterface_impl.hpp

@@ -844,4 +844,29 @@ namespace Opm
         return 1.0;
     }
 
+
+
+
+
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::calculateReservoirRates(WellState& well_state) const
+    {
+        const int fipreg = 0; // not considering the region for now
+        const int np = number_of_phases_;
+
+        std::vector<double> surface_rates(np, 0.0);
+        const int well_rate_index = np * index_of_well_;
+        for (int p =0; p < np; ++p) {
+            surface_rates[p] = well_state.wellRates()[well_rate_index + p];
+        }
+
+        std::vector<double> voidage_rates(np, 0.0);
+        rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegionIdx_, surface_rates, voidage_rates);
+
+        for (int p = 0; p < np; ++p) {
+            well_state.wellReservoirRates()[well_rate_index + p] = voidage_rates[p];
+        }
+    }
+
 }

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -296,7 +296,7 @@ namespace Opm
                 }
 
                 if ( pu.phase_used[Gas] ) {
-                    well.rates.set( rt::reservoir_oil, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Gas]] );
+                    well.rates.set( rt::reservoir_gas, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Gas]] );
                 }
 
                 well.rates.set( rt::dissolved_gas, this->well_dissolved_gas_rates_[w] );
@@ -549,7 +549,6 @@ namespace Opm
             return well_vaporized_oil_rates_;
         }
 
-
         const std::vector<double>& segRates() const
         {
             return segrates_;