Merge pull request #1405 from GitPaean/populating_reservoir_rated_output

[WIP] Populating reservoir rate related output

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

@@ -582,7 +582,13 @@ namespace Opm {
              *
              *
              * \param[out] coeff Surface-to-reservoir conversion
-             * coefficients for all active phases.
+             * coefficients that can be used to compute total reservoir
+             * volumes from surface volumes with the formula
+             *               q_{rT} = \sum_p coeff[p] q_{sp}.
+             * However, individual phase reservoir volumes cannot be calculated from
+             * these coefficients (i.e. q_{rp} is not equal to coeff[p] q_{sp})
+             * since they can depend on more than one surface volume rate when
+             * we have dissolved gas or vaporized oil.
              */
             template <class Coeff>
             void
@@ -640,6 +646,88 @@ 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());
+
+                std::fill(voidage_rates.begin(), voidage_rates.end(), 0.0);
+
+                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.hpp

@@ -281,7 +281,8 @@ namespace Opm
         void computePerfRate(const IntensiveQuantities& intQuants,
                              const std::vector<EvalWell>& mob_perfcells_dense,
                              const double Tw, const EvalWell& bhp, const double& cdp,
-                             const bool& allow_cf, std::vector<EvalWell>& cq_s) const;
+                             const bool& allow_cf, std::vector<EvalWell>& cq_s,
+                             double& perf_dis_gas_rate, double& perf_vap_oil_rate) const;
 
         // TODO: maybe we should provide a light version of computePerfRate, which does not include the
         // calculation of the derivatives

opm/autodiff/StandardWell_impl.hpp

@@ -395,7 +395,8 @@ namespace Opm
     computePerfRate(const IntensiveQuantities& intQuants,
                     const std::vector<EvalWell>& mob_perfcells_dense,
                     const double Tw, const EvalWell& bhp, const double& cdp,
-                    const bool& allow_cf, std::vector<EvalWell>& cq_s) const
+                    const bool& allow_cf, std::vector<EvalWell>& cq_s,
+                    double& perf_dis_gas_rate, double& perf_vap_oil_rate) const
     {
         std::vector<EvalWell> cmix_s(num_components_,0.0);
         for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
@@ -440,8 +441,17 @@ namespace Opm
                 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
                 const EvalWell cq_sOil = cq_s[oilCompIdx];
                 const EvalWell cq_sGas = cq_s[gasCompIdx];
-                cq_s[gasCompIdx] += rs * cq_sOil;
+                const EvalWell dis_gas = rs * cq_sOil;
-                cq_s[oilCompIdx] += rv * cq_sGas;
+                const EvalWell vap_oil = rv * cq_sGas;
+
+                cq_s[gasCompIdx] += dis_gas;
+                cq_s[oilCompIdx] += vap_oil;
+
+                // recording the perforation solution gas rate and solution oil rates
+                if (well_type_ == PRODUCER) {
+                    perf_dis_gas_rate = dis_gas.value();
+                    perf_vap_oil_rate = vap_oil.value();
+                }
             }
 
         } else {
@@ -506,6 +516,29 @@ namespace Opm
             for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
                 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is; // * b_perfcells_dense[phase];
             }
+
+            // calculating the perforation solution gas rate and solution oil rates
+            if (well_type_ == PRODUCER) {
+                if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                    const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
+                    const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+                    // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
+                    // s means standard condition, r means reservoir condition
+                    // q_os = q_or * b_o + rv * q_gr * b_g
+                    // q_gs = q_gr * g_g + rs * q_or * b_o
+                    // d = 1.0 - rs * rv
+                    // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
+                    // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
+
+                    const double d = 1.0 - rv.value() * rs.value();
+                    // vaporized oil into gas
+                    // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
+                    perf_vap_oil_rate = rv.value() * (cq_s[gasCompIdx].value() - rs.value() * cq_s[oilCompIdx].value()) / d;
+                    // dissolved of gas in oil
+                    // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
+                    perf_dis_gas_rate = rs.value() * (cq_s[oilCompIdx].value() - rv.value() * cq_s[gasCompIdx].value()) / d;
+                }
+            }
         }
     }
 
@@ -541,6 +574,10 @@ namespace Opm
 
         const EvalWell& bhp = getBhp();
 
+        // 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) {
 
             const int cell_idx = well_cells_[perf];
@@ -548,7 +585,16 @@ namespace Opm
             std::vector<EvalWell> cq_s(num_components_,0.0);
             std::vector<EvalWell> mob(num_components_, 0.0);
             getMobility(ebosSimulator, perf, mob);
-            computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
+            double perf_dis_gas_rate = 0.;
+            double perf_vap_oil_rate = 0.;
+            computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
+                            cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
+
+            // updating the solution gas rate and solution oil rate
+            if (well_type_ == PRODUCER) {
+                well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate;
+                well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate;
+            }
 
             for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
                 // the cq_s entering mass balance equations need to consider the efficiency factors.
@@ -1609,7 +1655,10 @@ namespace Opm
             std::vector<EvalWell> cq_s(num_components_, 0.0);
             std::vector<EvalWell> mob(num_components_, 0.0);
             getMobility(ebosSimulator, perf, mob);
-            computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
+            double perf_dis_gas_rate = 0.;
+            double perf_vap_oil_rate = 0.;
+            computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
+                            cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
 
             for(int p = 0; p < np; ++p) {
                 well_flux[ebosCompIdxToFlowCompIdx(p)] += cq_s[p].value();
@@ -1995,7 +2044,10 @@ namespace Opm
             const bool allow_cf = crossFlowAllowed(ebos_simulator);
             const EvalWell& bhp = getBhp();
             std::vector<EvalWell> cq_s(num_components_,0.0);
-            computePerfRate(int_quant, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
+            double perf_dis_gas_rate = 0.;
+            double perf_vap_oil_rate = 0.;
+            computePerfRate(int_quant, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
+                            cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
             // TODO: make area a member
             const double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf];
             const auto& material_law_manager = ebos_simulator.problem().materialLawManager();

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

@@ -80,11 +80,17 @@ namespace Opm
             }
 
             const int nw = wells->number_of_wells;
+
             if( nw == 0 ) return ;
 
             // Initialize perfphaserates_, which must be done here.
             const int np = wells->number_of_phases;
             const int nperf = wells->well_connpos[nw];
+
+            well_reservoir_rates_.resize(nw * np, 0.0);
+            well_dissolved_gas_rates_.resize(nw, 0.0);
+            well_vaporized_oil_rates_.resize(nw, 0.0);
+
             // Ensure that we start out with zero rates by default.
             perfphaserates_.clear();
             perfphaserates_.resize(nperf * np, 0.0);
@@ -279,6 +285,23 @@ namespace Opm
                 auto& well = res.at( wt.first );
                 well.control = this->currentControls()[ w ];
 
+                const int well_rate_index = w * pu.num_phases;
+
+                if ( pu.phase_used[Water] ) {
+                    well.rates.set( rt::reservoir_water, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Water]] );
+                }
+
+                if ( pu.phase_used[Oil] ) {
+                    well.rates.set( rt::reservoir_oil, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Oil]] );
+                }
+
+                if ( pu.phase_used[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] );
+                well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_rates_[w] );
+
                 int local_comp_index = 0;
                 for( auto& comp : well.completions ) {
                     const auto rates = this->perfPhaseRates().begin()
@@ -511,6 +534,21 @@ namespace Opm
             return solvent_well_rate;
         }
 
+        std::vector<double>& wellReservoirRates()
+        {
+            return well_reservoir_rates_;
+        }
+
+        std::vector<double>& wellDissolvedGasRates()
+        {
+            return well_dissolved_gas_rates_;
+        }
+
+        std::vector<double>& wellVaporizedOilRates()
+        {
+            return well_vaporized_oil_rates_;
+        }
+
         const std::vector<double>& segRates() const
         {
             return segrates_;
@@ -548,6 +586,18 @@ namespace Opm
         std::vector<int> current_controls_;
         std::vector<double> perfRateSolvent_;
 
+        // phase rates under reservoir condition for wells
+        // or voidage phase rates
+        std::vector<double> well_reservoir_rates_;
+
+        // dissolved gas rates or solution gas production rates
+        // should be zero for injection wells
+        std::vector<double> well_dissolved_gas_rates_;
+
+        // vaporized oil rates or solution oil producation rates
+        // should be zero for injection wells
+        std::vector<double> well_vaporized_oil_rates_;
+
         // marking whether the well is just added
         // for newly added well, the current initialized rates from WellState
         // will have very wrong compositions for production wells, will mostly cause