My Project
BlackoilWellModel.hpp
1 /*
2  Copyright 2016 SINTEF ICT, Applied Mathematics.
3  Copyright 2016 - 2017 Statoil ASA.
4  Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
5  Copyright 2016 - 2018 IRIS AS
6 
7  This file is part of the Open Porous Media project (OPM).
8 
9  OPM is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 3 of the License, or
12  (at your option) any later version.
13 
14  OPM is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with OPM. If not, see <http://www.gnu.org/licenses/>.
21 */
22 
23 
24 #ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
25 #define OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
26 
27 #include <ebos/eclproblem.hh>
28 #include <opm/common/OpmLog/OpmLog.hpp>
29 
30 #include <cassert>
31 #include <map>
32 #include <memory>
33 #include <optional>
34 #include <set>
35 #include <string>
36 #include <tuple>
37 #include <unordered_map>
38 #include <vector>
39 
40 #include <stddef.h>
41 
42 #include <opm/input/eclipse/EclipseState/Runspec.hpp>
43 
44 #include <opm/input/eclipse/Schedule/Schedule.hpp>
45 #include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
46 #include <opm/input/eclipse/Schedule/Group/GuideRate.hpp>
47 #include <opm/input/eclipse/Schedule/Group/Group.hpp>
48 
49 #include <opm/simulators/timestepping/SimulatorReport.hpp>
50 #include <opm/simulators/flow/countGlobalCells.hpp>
51 #include <opm/simulators/wells/BlackoilWellModelGeneric.hpp>
52 #include <opm/simulators/wells/GasLiftSingleWell.hpp>
53 #include <opm/simulators/wells/GasLiftWellState.hpp>
54 #include <opm/simulators/wells/GasLiftSingleWellGeneric.hpp>
55 #include <opm/simulators/wells/GasLiftStage2.hpp>
56 #include <opm/simulators/wells/GasLiftGroupInfo.hpp>
57 #include <opm/simulators/wells/PerforationData.hpp>
58 #include <opm/simulators/wells/VFPInjProperties.hpp>
59 #include <opm/simulators/wells/VFPProdProperties.hpp>
60 #include <opm/simulators/wells/WellState.hpp>
61 #include <opm/simulators/wells/WGState.hpp>
64 #include <opm/simulators/wells/WellInterface.hpp>
65 #include <opm/simulators/wells/StandardWell.hpp>
66 #include <opm/simulators/wells/MultisegmentWell.hpp>
67 #include <opm/simulators/wells/WellGroupHelpers.hpp>
68 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
69 #include <opm/simulators/wells/ParallelWellInfo.hpp>
70 #include <opm/simulators/timestepping/gatherConvergenceReport.hpp>
71 #include <dune/common/fmatrix.hh>
72 #include <dune/istl/bcrsmatrix.hh>
73 #include <dune/istl/matrixmatrix.hh>
74 
75 #include <opm/material/densead/Math.hpp>
76 
77 #include <opm/simulators/utils/DeferredLogger.hpp>
78 
79 namespace Opm::Properties {
80 
81 template<class TypeTag, class MyTypeTag>
83  using type = UndefinedProperty;
84 };
85 
86 } // namespace Opm::Properties
87 
88 namespace Opm {
89 
91  template<typename TypeTag>
92  class BlackoilWellModel : public BaseAuxiliaryModule<TypeTag>
94  {
95  public:
96  // --------- Types ---------
98 
99  using Grid = GetPropType<TypeTag, Properties::Grid>;
100  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
101  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
102  using Indices = GetPropType<TypeTag, Properties::Indices>;
103  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
104  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
105  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
106  using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
107  using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
108  using GasLiftSingleWell = typename WellInterface<TypeTag>::GasLiftSingleWell;
109  using GLiftOptWells = typename BlackoilWellModelGeneric::GLiftOptWells;
110  using GLiftProdWells = typename BlackoilWellModelGeneric::GLiftProdWells;
111  using GLiftWellStateMap =
112  typename BlackoilWellModelGeneric::GLiftWellStateMap;
113  using GLiftEclWells = typename GasLiftGroupInfo::GLiftEclWells;
114  using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups;
115 
116  typedef typename BaseAuxiliaryModule<TypeTag>::NeighborSet NeighborSet;
117 
118  static const int numEq = Indices::numEq;
119  static const int solventSaturationIdx = Indices::solventSaturationIdx;
120  static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
121  static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
122  static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
123  static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
124 
125  // TODO: where we should put these types, WellInterface or Well Model?
126  // or there is some other strategy, like TypeTag
127  typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
128  typedef Dune::BlockVector<VectorBlockType> BVector;
129 
130  typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
131 
132  typedef BlackOilPolymerModule<TypeTag> PolymerModule;
133  typedef BlackOilMICPModule<TypeTag> MICPModule;
134 
135  // For the conversion between the surface volume rate and resrevoir voidage rate
136  using RateConverterType = RateConverter::
137  SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
138 
139  // For computing average pressured used by gpmaint
140  using AverageRegionalPressureType = RegionAverageCalculator::
141  AverageRegionalPressure<FluidSystem, std::vector<int> >;
142 
143  BlackoilWellModel(Simulator& ebosSimulator);
144 
145  void init();
146  void initWellContainer() override;
147 
149  // <eWoms auxiliary module stuff>
151  unsigned numDofs() const override
152  // No extra dofs are inserted for wells. (we use a Schur complement.)
153  { return 0; }
154 
155  void addNeighbors(std::vector<NeighborSet>& neighbors) const override;
156 
157  void applyInitial() override
158  {}
159 
160  void linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res) override;
161 
162  void postSolve(GlobalEqVector& deltaX) override
163  {
164  recoverWellSolutionAndUpdateWellState(deltaX);
165  }
166 
168  // </ eWoms auxiliary module stuff>
170 
171  template <class Restarter>
172  void deserialize(Restarter& /* res */)
173  {
174  // TODO (?)
175  }
176 
181  template <class Restarter>
182  void serialize(Restarter& /* res*/)
183  {
184  // TODO (?)
185  }
186 
187  void beginEpisode()
188  {
189  beginReportStep(ebosSimulator_.episodeIndex());
190  }
191 
192  void beginTimeStep();
193 
194  void beginIteration()
195  {
196  assemble(ebosSimulator_.model().newtonMethod().numIterations(),
197  ebosSimulator_.timeStepSize());
198  }
199 
200  void endIteration()
201  { }
202 
203  void endTimeStep()
204  {
205  timeStepSucceeded(ebosSimulator_.time(), ebosSimulator_.timeStepSize());
206  }
207 
208  void endEpisode()
209  {
210  endReportStep();
211  }
212 
213  template <class Context>
214  void computeTotalRatesForDof(RateVector& rate,
215  const Context& context,
216  unsigned spaceIdx,
217  unsigned timeIdx) const;
218 
219 
220  using WellInterfacePtr = std::shared_ptr<WellInterface<TypeTag> >;
221 
222  using BlackoilWellModelGeneric::initFromRestartFile;
223  void initFromRestartFile(const RestartValue& restartValues)
224  {
225  initFromRestartFile(restartValues,
226  this->ebosSimulator_.vanguard().transferWTestState(),
227  UgGridHelpers::numCells(grid()),
228  param_.use_multisegment_well_);
229  }
230 
231  data::Wells wellData() const
232  {
233  auto wsrpt = this->wellState()
234  .report(UgGridHelpers::globalCell(this->grid()),
235  [this](const int well_index) -> bool
236  {
237  return this->wasDynamicallyShutThisTimeStep(well_index);
238  });
239 
240  this->assignWellTracerRates(wsrpt);
241 
242  this->assignWellGuideRates(wsrpt, this->reportStepIndex());
243  this->assignShutConnections(wsrpt, this->reportStepIndex());
244 
245  return wsrpt;
246  }
247 
248  // subtract Binv(D)rw from r;
249  void apply( BVector& r) const;
250 
251  // subtract B*inv(D)*C * x from A*x
252  void apply(const BVector& x, BVector& Ax) const;
253 
254 #if HAVE_CUDA || HAVE_OPENCL
255  // accumulate the contributions of all Wells in the WellContributions object
256  void getWellContributions(WellContributions& x) const;
257 #endif
258 
259  // apply well model with scaling of alpha
260  void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
261 
262  // Check if well equations is converged.
263  ConvergenceReport getWellConvergence(const std::vector<Scalar>& B_avg, const bool checkGroupConvergence = false) const;
264 
265  const SimulatorReportSingle& lastReport() const;
266 
267  void addWellContributions(SparseMatrixAdapter& jacobian) const
268  {
269  for ( const auto& well: well_container_ ) {
270  well->addWellContributions(jacobian);
271  }
272  }
273 
274  // called at the beginning of a report step
275  void beginReportStep(const int time_step);
276 
277  void updatePerforationIntensiveQuantities();
278  // it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation()
279  // makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions
280  // twice at the beginning of the time step
283  void calculateExplicitQuantities(DeferredLogger& deferred_logger) const;
284  // some preparation work, mostly related to group control and RESV,
285  // at the beginning of each time step (Not report step)
286  void prepareTimeStep(DeferredLogger& deferred_logger);
287  void initPrimaryVariablesEvaluation() const;
288  void updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls);
289 
290  void updateAndCommunicate(const int reportStepIdx,
291  const int iterationIdx,
292  DeferredLogger& deferred_logger);
293 
294  WellInterfacePtr getWell(const std::string& well_name) const;
295  bool hasWell(const std::string& well_name) const;
296 
297  void initGliftEclWellMap(GLiftEclWells &ecl_well_map);
298 
300  const std::vector<WellInterfacePtr>& localNonshutWells()
301  {
302  return well_container_;
303  }
304 
305  protected:
306  Simulator& ebosSimulator_;
307 
308  // a vector of all the wells.
309  std::vector<WellInterfacePtr > well_container_{};
310 
311  std::vector<bool> is_cell_perforated_{};
312 
313  void initializeWellState(const int timeStepIdx,
314  const SummaryState& summaryState);
315 
316  // create the well container
317  void createWellContainer(const int time_step) override;
318 
319  WellInterfacePtr
320  createWellPointer(const int wellID,
321  const int time_step) const;
322 
323  template <typename WellType>
324  std::unique_ptr<WellType>
325  createTypedWellPointer(const int wellID,
326  const int time_step) const;
327 
328  WellInterfacePtr createWellForWellTest(const std::string& well_name, const int report_step, DeferredLogger& deferred_logger) const;
329 
330 
331  const ModelParameters param_;
332  size_t global_num_cells_{};
333  // the number of the cells in the local grid
334  size_t local_num_cells_{};
335  double gravity_{};
336  std::vector<double> depth_{};
337  bool alternative_well_rate_init_{};
338 
339  std::unique_ptr<RateConverterType> rateConverter_{};
340  std::unique_ptr<AverageRegionalPressureType> regionalAveragePressureCalculator_{};
341 
342 
343  SimulatorReportSingle last_report_{};
344 
345  // used to better efficiency of calcuation
346  mutable BVector scaleAddRes_{};
347 
348  std::vector<Scalar> B_avg_{};
349 
350  const Grid& grid() const
351  { return ebosSimulator_.vanguard().grid(); }
352 
353  const EclipseState& eclState() const
354  { return ebosSimulator_.vanguard().eclState(); }
355 
356  // compute the well fluxes and assemble them in to the reservoir equations as source terms
357  // and in the well equations.
358  void assemble(const int iterationIdx,
359  const double dt);
360 
361  // called at the end of a time step
362  void timeStepSucceeded(const double& simulationTime, const double dt);
363 
364  // called at the end of a report step
365  void endReportStep();
366 
367  // using the solution x to recover the solution xw for wells and applying
368  // xw to update Well State
369  void recoverWellSolutionAndUpdateWellState(const BVector& x);
370 
371  // setting the well_solutions_ based on well_state.
372  void updatePrimaryVariables(DeferredLogger& deferred_logger);
373 
374  void updateAverageFormationFactor();
375 
376  void computePotentials(const std::size_t widx,
377  const WellState& well_state_copy,
378  std::string& exc_msg,
379  ExceptionType::ExcEnum& exc_type,
380  DeferredLogger& deferred_logger) override;
381 
382  const std::vector<double>& wellPerfEfficiencyFactors() const;
383 
384  void calculateProductivityIndexValuesShutWells(const int reportStepIdx, DeferredLogger& deferred_logger) override;
385  void calculateProductivityIndexValues(DeferredLogger& deferred_logger) override;
386  void calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
387  DeferredLogger& deferred_logger);
388 
389  // The number of components in the model.
390  int numComponents() const;
391 
392  int reportStepIndex() const;
393 
394  void assembleWellEq(const double dt, DeferredLogger& deferred_logger);
395 
396  bool maybeDoGasLiftOptimize(DeferredLogger& deferred_logger);
397 
398  void gasLiftOptimizationStage1(DeferredLogger& deferred_logger,
399  GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
400  GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map);
401 
402  // cannot be const since it accesses the non-const WellState
403  void gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
404  DeferredLogger& deferred_logger,
405  GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
406  GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
407  GLiftSyncGroups& groups_to_sync);
408 
409  void extractLegacyCellPvtRegionIndex_();
410 
411  void extractLegacyDepth_();
412 
414  void updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const;
415 
416  void wellTesting(const int timeStepIdx, const double simulationTime, DeferredLogger& deferred_logger);
417 
418  void calcRates(const int fipnum,
419  const int pvtreg,
420  std::vector<double>& resv_coeff) override;
421 
422  void calcInjRates(const int fipnum,
423  const int pvtreg,
424  std::vector<double>& resv_coeff) override;
425 
426  void computeWellTemperature();
427 
428  void assignWellTracerRates(data::Wells& wsrpt) const;
429 
430  int compressedIndexForInterior(int cartesian_cell_idx) const override {
431  return ebosSimulator_.vanguard().compressedIndexForInterior(cartesian_cell_idx);
432  }
433 
434  private:
435  BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& pu);
436  };
437 
438 
439 } // namespace Opm
440 
441 #include "BlackoilWellModel_impl.hpp"
442 #endif
Facility for converting component rates at surface conditions to phase (voidage) rates at reservoir c...
Facility for converting component rates at surface conditions to phase (voidage) rates at reservoir c...
Class for handling the blackoil well model.
Definition: BlackoilWellModelGeneric.hpp:72
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:94
void serialize(Restarter &)
This method writes the complete state of the well to the harddisk.
Definition: BlackoilWellModel.hpp:182
const std::vector< WellInterfacePtr > & localNonshutWells()
Get list of local nonshut wells.
Definition: BlackoilWellModel.hpp:300
void calculateExplicitQuantities(DeferredLogger &deferred_logger) const
Calculating the explict quantities used in the well calculation.
Definition: BlackoilWellModel_impl.hpp:1273
void updateWellTestState(const double &simulationTime, WellTestState &wellTestState) const
upate the wellTestState related to economic limits
Definition: BlackoilWellModel_impl.hpp:1376
int compressedIndexForInterior(int cartesian_cell_idx) const override
get compressed index for interior cells (-1, otherwise
Definition: BlackoilWellModel.hpp:430
Definition: GasLiftSingleWell.hpp:39
Convert component rates at surface conditions to phase (voidage) rates at reservoir conditions.
Definition: RateConverter.hpp:68
Computes hydrocarbon weighed average pressures over regions.
Definition: RegionAverageCalculator.hpp:62
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParametersEbos.hpp:327
bool use_multisegment_well_
Whether to use MultisegmentWell to handle multisegment wells it is something temporary before the mul...
Definition: BlackoilModelParametersEbos.hpp:408
Definition: BlackoilPhases.hpp:46
Definition: BlackoilWellModel.hpp:82