Store well potential in well_state and pass it for output if asked for

A zero well potential is passed if the computation fails or
if it is a multisegmented well.

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -235,6 +235,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_;
     }
 
@@ -770,15 +783,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

@@ -556,9 +556,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/WellStateFullyImplicitBlackoil.hpp

@@ -129,6 +129,7 @@ namespace Opm
             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
@@ -265,6 +266,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();
@@ -496,6 +498,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] );
 
@@ -785,6 +799,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_;
@@ -820,6 +842,9 @@ namespace Opm
         // Productivity Index
         std::vector<double> productivity_index_;
 
+        // Well potentials
+        std::vector<double> well_potentials_;
+
     };
 
 } // namespace Opm