Merge pull request #1614 from totto82/well_ppp

Output of well potential when asked for

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -314,6 +314,19 @@ namespace Opm {
             well->calculateReservoirRates(well_state_);
         }
         updateWellTestState(simulationTime, wellTestState_);
+
+        // calculate the well potentials for output
+        // TODO: when necessary
+        try
+        {
+            std::vector<double> well_potentials;
+            computeWellPotentials(well_potentials);
+        }
+        catch ( std::runtime_error& e )
+        {
+            const std::string msg = "A zero well potential is returned for output purposes. ";
+            OpmLog::warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
+        }
         previous_well_state_ = well_state_;
     }
 
@@ -897,15 +910,31 @@ namespace Opm {
         const int np = numPhases();
         well_potentials.resize(nw * np, 0.0);
 
+        const Opm::SummaryConfig summaryConfig = ebosSimulator_.vanguard().summaryConfig();
         for (const auto& well : well_container_) {
-            std::vector<double> potentials;
-            well->computeWellPotentials(ebosSimulator_, well_state_, potentials);
+            // Only compute the well potential when asked for
+            bool needed_for_output = ((summaryConfig.hasSummaryKey( "WWPI:" + well->name()) ||
+                                       summaryConfig.hasSummaryKey( "WOPI:" + well->name()) ||
+                                       summaryConfig.hasSummaryKey( "WGPI:" + well->name())) && well->wellType() == INJECTOR) ||
+                                    ((summaryConfig.hasSummaryKey( "WWPP:" + well->name()) ||
+                                                       summaryConfig.hasSummaryKey( "WOPP:" + well->name()) ||
+                                                       summaryConfig.hasSummaryKey( "WGPP:" + well->name())) && well->wellType() == PRODUCER);
 
-            // putting the sucessfully calculated potentials to the well_potentials
-            for (int p = 0; p < np; ++p) {
-                well_potentials[well->indexOfWell() * np + p] = std::abs(potentials[p]);
+            if (needed_for_output || wellCollection().requireWellPotentials())
+            {
+                std::vector<double> potentials;
+                well->computeWellPotentials(ebosSimulator_, well_state_, potentials);
+
+                // putting the sucessfully calculated potentials to the well_potentials
+                for (int p = 0; p < np; ++p) {
+                    well_potentials[well->indexOfWell() * np + p] = std::abs(potentials[p]);
+                }
             }
         } // end of for (int w = 0; w < nw; ++w)
+
+        // Store it in the well state
+        well_state_.wellPotentials() = well_potentials;
+
     }
 
 

opm/autodiff/MultisegmentWell_impl.hpp

@@ -559,9 +559,17 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     computeWellPotentials(const Simulator& /* ebosSimulator */,
                           const WellState& /* well_state */,
-                          std::vector<double>& /* well_potentials */)
+                          std::vector<double>& well_potentials)
     {
-        OPM_THROW(std::runtime_error, "well potential calculation for multisegment wells is not supported yet");
+        const std::string msg = std::string("Well potential calculation is not supported for multisegment wells \n")
+                + "A well potential of zero is returned for output purposes. \n"
+                + "If you need well potential to set the guide rate for group controled wells \n"
+                + "you will have to change the " + name() + " well to a standard well \n";
+
+        OpmLog::warning("WELL_POTENTIAL_NOT_IMPLEMENTED_FOR_MULTISEG_WELLS", msg);
+
+        const int np = number_of_phases_;
+        well_potentials.resize(np, 0.0);
     }
 
 

opm/autodiff/StandardWell_impl.hpp

@@ -439,6 +439,9 @@ namespace Opm
                    const double dt,
                    WellState& well_state)
     {
+
+        const Opm::SummaryConfig summaryConfig = ebosSimulator.vanguard().summaryConfig();
+
         const int np = number_of_phases_;
 
         // clear all entries
@@ -584,13 +587,19 @@ namespace Opm
             // Store the perforation pressure for later usage.
             well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf];
 
-            // Compute Productivity index
+            // Compute Productivity index if asked for
+             const auto& pu = phaseUsage();
             for (int p = 0; p < np; ++p) {
-                const unsigned int compIdx = flowPhaseToEbosCompIdx(p);
-                const double drawdown = well_state.perfPress()[first_perf_ + perf] - intQuants.fluidState().pressure(FluidSystem::oilPhaseIdx).value();
-                double productivity_index = cq_s[compIdx].value() / drawdown;
-                scaleProductivityIndex(perf, productivity_index);
-                well_state.productivityIndex()[np*index_of_well_ + p] += productivity_index;
+                if ( (pu.phase_pos[Water] == p && (summaryConfig.hasSummaryKey("WPIW:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name())))
+                        || (pu.phase_pos[Oil] == p && (summaryConfig.hasSummaryKey("WPIO:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name())))
+                        || (pu.phase_pos[Gas] == p && summaryConfig.hasSummaryKey("WPIG:" + name()))) {
+
+                    const unsigned int compIdx = flowPhaseToEbosCompIdx(p);
+                    const double drawdown = well_state.perfPress()[first_perf_ + perf] - intQuants.fluidState().pressure(FluidSystem::oilPhaseIdx).value();
+                    double productivity_index = cq_s[compIdx].value() / drawdown;
+                    scaleProductivityIndex(perf, productivity_index);
+                    well_state.productivityIndex()[np*index_of_well_ + p] += productivity_index;
+                }
             }
 
         }

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -145,6 +145,7 @@ namespace Opm
             perfRateSolvent_.clear();
             perfRateSolvent_.resize(nperf, 0.0);
             productivity_index_.resize(nw * np, 0.0);
+            well_potentials_.resize(nw * np, 0.0);
 
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
@@ -292,6 +293,7 @@ namespace Opm
             well_vaporized_oil_rates_.resize(nw, 0.0);
 
             productivity_index_.resize(nw * np, 0.0);
+            well_potentials_.resize(nw*np, 0.0);
 
             // Ensure that we start out with zero rates by default.
             perfphaserates_.clear();
@@ -518,6 +520,18 @@ namespace Opm
                     well.rates.set( rt::productivity_index_gas, this->productivity_index_[well_rate_index + pu.phase_pos[Gas]] );
                 }
 
+                if ( pu.phase_used[Water] ) {
+                    well.rates.set( rt::well_potential_water, this->well_potentials_[well_rate_index + pu.phase_pos[Water]] );
+                }
+
+                if ( pu.phase_used[Oil] ) {
+                    well.rates.set( rt::well_potential_oil, this->well_potentials_[well_rate_index + pu.phase_pos[Oil]] );
+                }
+
+                if ( pu.phase_used[Gas] ) {
+                    well.rates.set( rt::well_potential_gas, this->well_potentials_[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] );
 
@@ -807,6 +821,14 @@ namespace Opm
             return productivity_index_;
         }
 
+        std::vector<double>& wellPotentials() {
+            return well_potentials_;
+        }
+
+        const std::vector<double>& wellPotentials() const {
+            return well_potentials_;
+        }
+
     private:
         std::vector<double> perfphaserates_;
         std::vector<int> current_controls_;
@@ -841,6 +863,9 @@ namespace Opm
         // Productivity Index
         std::vector<double> productivity_index_;
 
+        // Well potentials
+        std::vector<double> well_potentials_;
+
     };
 
 } // namespace Opm