Merge pull request #5490 from totto82/useCellTempRESV

Use cell temperature in perforated cell to compute reservoir rates

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -812,6 +812,7 @@ namespace Opm {
     initializeWellState(const int timeStepIdx)
     {
         std::vector<Scalar> cellPressures(this->local_num_cells_, 0.0);
+        std::vector<Scalar> cellTemperatures(this->local_num_cells_, 0.0);
         ElementContext elemCtx(simulator_);
 
         const auto& gridView = simulator_.vanguard().gridView();
@@ -831,10 +832,11 @@ namespace Opm {
             } else {
                 perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value();
             }
+            cellTemperatures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)] = fs.temperature(0).value();
         }
         OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::initializeWellState() failed: ", simulator_.vanguard().grid().comm());
 
-        this->wellState().init(cellPressures, this->schedule(), this->wells_ecl_,
+        this->wellState().init(cellPressures, cellTemperatures, this->schedule(), this->wells_ecl_,
                                this->local_parallel_well_info_, timeStepIdx,
                                &this->prevWellState(), this->well_perf_data_,
                                this->summaryState());

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -711,7 +711,7 @@ namespace Opm
             this->segments_.copyPhaseDensities(ws.pu, ws.segments);
         }
 
-        Base::calculateReservoirRates(well_state.well(this->index_of_well_));
+        Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
     }
 
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -706,7 +706,7 @@ namespace Opm
 
         const auto& summary_state = simulator.vanguard().summaryState();
         updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
-        Base::calculateReservoirRates(well_state.well(this->index_of_well_));
+        Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
     }
 
 
@@ -2587,7 +2587,7 @@ namespace Opm
             if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
                 // only handles single phase injection now
                 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
-                fs.setTemperature(this->well_ecl_.temperature());
+                fs.setTemperature(this->well_ecl_.inj_temperature());
                 typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
                 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
                 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();

opm/simulators/wells/WellInterfaceFluidSystem.cpp

@@ -64,12 +64,12 @@ WellInterfaceFluidSystem(const Well& well,
 template <typename FluidSystem>
 void
 WellInterfaceFluidSystem<FluidSystem>::
-calculateReservoirRates(SingleWellState<Scalar>& ws) const
+calculateReservoirRates(const bool co2store, SingleWellState<Scalar>& ws) const
 {
     const int np = this->number_of_phases_;
     const auto& pu = this->phaseUsage();
     // Calculate reservoir rates from average pressure and temperature
-    if ( !pu.has_energy || this->wellEcl().isProducer()) {
+    if ( !(co2store || pu.has_energy) || this->wellEcl().isProducer()) {
         const int fipreg = 0; // not considering the region for now
         this->rateConverter_
             .calcReservoirVoidageRates(fipreg,
@@ -99,7 +99,8 @@ calculateReservoirRates(SingleWellState<Scalar>& ws) const
         }
         return;
     }
-    // For injectors in a thermal case we convert using the well bhp and temperature
+    // For injectors in a co2 storage case or a thermal case
+    // we convert using the well bhp and temperature
     // Assume pure phases in the injector
     const Scalar saltConc = 0.0;
     Scalar rsMax = 0.0;

opm/simulators/wells/WellInterfaceFluidSystem.hpp

@@ -81,7 +81,7 @@ protected:
                              const std::vector<PerforationData<Scalar>>& perf_data);
 
     // updating the voidage rates in well_state when requested
-    void calculateReservoirRates(SingleWellState<Scalar>& ws) const;
+    void calculateReservoirRates(const bool co2store, SingleWellState<Scalar>& ws) const;
 
     bool checkIndividualConstraints(SingleWellState<Scalar>& ws,
                                     const SummaryState& summaryState,

opm/simulators/wells/WellState.cpp

@@ -149,6 +149,7 @@ serializationTestObject(const ParallelWellInfo<Scalar>& pinfo)
 
 template<class Scalar>
 void WellState<Scalar>::base_init(const std::vector<Scalar>& cellPressures,
+                                  const std::vector<Scalar>& cellTemperatures,
                                   const std::vector<Well>& wells_ecl,
                                   const std::vector<std::reference_wrapper<ParallelWellInfo<Scalar>>>& parallel_well_info,
                                   const std::vector<std::vector<PerforationData<Scalar>>>& well_perf_data,
@@ -164,7 +165,7 @@ void WellState<Scalar>::base_init(const std::vector<Scalar>& cellPressures,
             const Well& well = wells_ecl[w];
 
             // Initialize bhp(), thp(), wellRates(), temperature().
-            initSingleWell(cellPressures, well, well_perf_data[w], parallel_well_info[w], summary_state);
+            initSingleWell(cellPressures, cellTemperatures, well, well_perf_data[w], parallel_well_info[w], summary_state);
         }
     }
 }
@@ -202,12 +203,12 @@ template<class Scalar>
 void WellState<Scalar>::initSingleInjector(const Well& well,
                                            const ParallelWellInfo<Scalar>& well_info,
                                            Scalar pressure_first_connection,
+                                           Scalar temperature_first_connection,
                                            const std::vector<PerforationData<Scalar>>& well_perf_data,
                                            const SummaryState& summary_state)
 {
     const auto& pu = this->phase_usage_;
-    const Scalar temp = well.temperature();
+    const Scalar temp = well.hasInjTemperature() ? well.inj_temperature() : temperature_first_connection;
-
     auto& ws = this->wells_.add(well.name(), SingleWellState<Scalar>{well.name(),
                                                                      well_info,
                                                                      false,
@@ -228,18 +229,25 @@ void WellState<Scalar>::initSingleInjector(const Well& well,
 
 template<class Scalar>
 void WellState<Scalar>::initSingleWell(const std::vector<Scalar>& cellPressures,
+                                       const std::vector<Scalar>& cellTemperatures,
                                        const Well& well,
                                        const std::vector<PerforationData<Scalar>>& well_perf_data,
                                        const ParallelWellInfo<Scalar>& well_info,
                                        const SummaryState& summary_state)
 {
     Scalar pressure_first_connection = -1;
-    if (!well_perf_data.empty())
+    if (!well_perf_data.empty()) {
         pressure_first_connection = cellPressures[well_perf_data[0].cell_index];
+    }
     pressure_first_connection = well_info.broadcastFirstPerforationValue(pressure_first_connection);
 
     if (well.isInjector()) {
-        this->initSingleInjector(well, well_info, pressure_first_connection,
+        Scalar temperature_first_connection = -1;
+        if (!well_perf_data.empty()) {
+            temperature_first_connection = cellTemperatures[well_perf_data[0].cell_index];
+        }
+        temperature_first_connection = well_info.broadcastFirstPerforationValue(temperature_first_connection);
+        this->initSingleInjector(well, well_info, pressure_first_connection, temperature_first_connection,
                                  well_perf_data, summary_state);
     } else {
         this->initSingleProducer(well, well_info, pressure_first_connection,
@@ -249,6 +257,7 @@ void WellState<Scalar>::initSingleWell(const std::vector<Scalar>& cellPressures,
 
 template<class Scalar>
 void WellState<Scalar>::init(const std::vector<Scalar>& cellPressures,
+                             const std::vector<Scalar>& cellTemperatures,
                              const Schedule& schedule,
                              const std::vector<Well>& wells_ecl,
                              const std::vector<std::reference_wrapper<ParallelWellInfo<Scalar>>>& parallel_well_info,
@@ -258,7 +267,7 @@ void WellState<Scalar>::init(const std::vector<Scalar>& cellPressures,
                              const SummaryState& summary_state)
 {
     // call init on base class
-    this->base_init(cellPressures, wells_ecl, parallel_well_info,
+    this->base_init(cellPressures, cellTemperatures, wells_ecl, parallel_well_info,
                     well_perf_data, summary_state);
     this->global_well_info = std::make_optional<GlobalWellInfo>(schedule,
                                                                 report_step,
@@ -428,7 +437,7 @@ void WellState<Scalar>::resize(const std::vector<Well>& wells_ecl,
                                const SummaryState& summary_state)
 {
     const std::vector<Scalar> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
-    init(tmp, schedule, wells_ecl, parallel_well_info, 0, nullptr, well_perf_data, summary_state);
+    init(tmp, tmp, schedule, wells_ecl, parallel_well_info, 0, nullptr, well_perf_data, summary_state);
 
     if (handle_ms_well) {
         initWellStateMSWell(wells_ecl, nullptr);

opm/simulators/wells/WellState.hpp

@@ -97,6 +97,7 @@ public:
     /// to give useful initial values to the bhp(), wellRates()
     /// and perfPhaseRatesORG() fields, depending on controls
     void init(const std::vector<Scalar>& cellPressures,
+              const std::vector<Scalar>& cellTemperatures,
               const Schedule& schedule,
               const std::vector<Well>& wells_ecl,
               const std::vector<std::reference_wrapper<ParallelWellInfo<Scalar>>>& parallel_well_info,
@@ -377,12 +378,14 @@ private:
     /// perfRates() field is filled with zero, and perfPress()
     /// with -1e100.
     void base_init(const std::vector<Scalar>& cellPressures,
+                   const std::vector<Scalar>& cellTemperatures,
                    const std::vector<Well>& wells_ecl,
                    const std::vector<std::reference_wrapper<ParallelWellInfo<Scalar>>>& parallel_well_info,
                    const std::vector<std::vector<PerforationData<Scalar>>>& well_perf_data,
                    const SummaryState& summary_state);
 
     void initSingleWell(const std::vector<Scalar>& cellPressures,
+                        const std::vector<Scalar>& cellTemperatures,
                         const Well& well,
                         const std::vector<PerforationData<Scalar>>& well_perf_data,
                         const ParallelWellInfo<Scalar>& well_info,
@@ -397,6 +400,7 @@ private:
     void initSingleInjector(const Well& well,
                             const ParallelWellInfo<Scalar>& well_info,
                             Scalar pressure_first_connection,
+                            Scalar temperature_first_connection,
                             const std::vector<PerforationData<Scalar>>& well_perf_data,
                             const SummaryState& summary_state);
 };

tests/test_wellstate.cpp

@@ -158,6 +158,12 @@ namespace {
             std::vector<double>(setup.grid.c_grid()->number_of_cells,
                                 100.0*Opm::unit::barsa);
 
+        const auto& unit_sytem = setup.es.getDeckUnitSystem();
+        const double temp = unit_sytem.to_si(Opm::UnitSystem::measure::temperature, 25);
+        const auto ctemp =
+            std::vector<double>(setup.grid.c_grid()->number_of_cells,
+                                temp);
+
         auto wells = setup.sched.getWells(timeStep);
         pinfos.resize(wells.size());
         std::vector<std::reference_wrapper<Opm::ParallelWellInfo<double>>> ppinfos;
@@ -171,7 +177,7 @@ namespace {
             ++pw;
         }
 
-        state.init(cpress, setup.sched,
+        state.init(cpress, ctemp, setup.sched,
                    wells, ppinfos,
                    timeStep, nullptr, setup.well_perf_data, setup.st);