My Project
Loading...
Searching...
No Matches
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 <opm/common/OpmLog/OpmLog.hpp>
28
29#include <cstddef>
30#include <map>
31#include <memory>
32#include <optional>
33#include <string>
34#include <vector>
35
36#include <opm/input/eclipse/Schedule/Group/Group.hpp>
37#include <opm/input/eclipse/Schedule/Group/GuideRate.hpp>
38#include <opm/input/eclipse/Schedule/Schedule.hpp>
39#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
40
41#include <opm/models/discretization/common/baseauxiliarymodule.hh>
42
43#include <opm/simulators/flow/countGlobalCells.hpp>
45#include <opm/simulators/flow/SubDomain.hpp>
46
47#include <opm/simulators/linalg/matrixblock.hh>
48
49#include <opm/simulators/wells/BlackoilWellModelGeneric.hpp>
50#include <opm/simulators/wells/BlackoilWellModelGuideRates.hpp>
51#include <opm/simulators/wells/GasLiftGroupInfo.hpp>
52#include <opm/simulators/wells/GasLiftSingleWell.hpp>
53#include <opm/simulators/wells/GasLiftSingleWellGeneric.hpp>
54#include <opm/simulators/wells/GasLiftWellState.hpp>
55#include <opm/simulators/wells/MultisegmentWell.hpp>
56#include <opm/simulators/wells/ParallelWBPCalculation.hpp>
57#include <opm/simulators/wells/ParallelWellInfo.hpp>
58#include <opm/simulators/wells/PerforationData.hpp>
61#include <opm/simulators/wells/StandardWell.hpp>
62#include <opm/simulators/wells/VFPInjProperties.hpp>
63#include <opm/simulators/wells/VFPProdProperties.hpp>
64#include <opm/simulators/wells/WGState.hpp>
65#include <opm/simulators/wells/WellGroupHelpers.hpp>
66#include <opm/simulators/wells/WellInterface.hpp>
67#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
68#include <opm/simulators/wells/WellState.hpp>
69
70#include <opm/simulators/timestepping/SimulatorReport.hpp>
71#include <opm/simulators/timestepping/gatherConvergenceReport.hpp>
72
73#include <dune/common/fmatrix.hh>
74#include <dune/istl/bcrsmatrix.hh>
75#include <dune/istl/matrixmatrix.hh>
76
77#include <opm/material/densead/Math.hpp>
78
79#include <opm/simulators/utils/DeferredLogger.hpp>
80
81namespace Opm::Properties {
82
83template<class TypeTag, class MyTypeTag>
85 using type = UndefinedProperty;
86};
87
88} // namespace Opm::Properties
89
90namespace Opm {
91
93 template<typename TypeTag>
94 class BlackoilWellModel : public BaseAuxiliaryModule<TypeTag>
96 {
97 public:
98 // --------- Types ---------
99 using ModelParameters = BlackoilModelParameters<TypeTag>;
100
101 using Grid = GetPropType<TypeTag, Properties::Grid>;
102 using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
103 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
104 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
105 using Indices = GetPropType<TypeTag, Properties::Indices>;
106 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
107 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
108 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
109 using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
110 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
111 using GasLiftSingleWell = typename WellInterface<TypeTag>::GasLiftSingleWell;
112 using GLiftOptWells = typename BlackoilWellModelGeneric::GLiftOptWells;
113 using GLiftProdWells = typename BlackoilWellModelGeneric::GLiftProdWells;
114 using GLiftWellStateMap =
115 typename BlackoilWellModelGeneric::GLiftWellStateMap;
116 using GLiftEclWells = typename GasLiftGroupInfo::GLiftEclWells;
117 using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups;
118 constexpr static std::size_t pressureVarIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
119 typedef typename BaseAuxiliaryModule<TypeTag>::NeighborSet NeighborSet;
120
121 static const int numEq = Indices::numEq;
122 static const int solventSaturationIdx = Indices::solventSaturationIdx;
123 static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
124 static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
125 static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
126 static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
127
128 // TODO: where we should put these types, WellInterface or Well Model?
129 // or there is some other strategy, like TypeTag
130 typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
131 typedef Dune::BlockVector<VectorBlockType> BVector;
132
133 typedef BlackOilPolymerModule<TypeTag> PolymerModule;
134 typedef BlackOilMICPModule<TypeTag> MICPModule;
135
136 // For the conversion between the surface volume rate and resrevoir voidage rate
137 using RateConverterType = RateConverter::
138 SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
139
140 // For computing average pressured used by gpmaint
141 using AverageRegionalPressureType = RegionAverageCalculator::
142 AverageRegionalPressure<FluidSystem, std::vector<int> >;
143
144 using Domain = SubDomain<Grid>;
145
146 BlackoilWellModel(Simulator& simulator);
147
148 void init();
149 void initWellContainer(const int reportStepIdx) override;
150
152 // <eWoms auxiliary module stuff>
154 unsigned numDofs() const override
155 // No extra dofs are inserted for wells. (we use a Schur complement.)
156 { return 0; }
157
158 void addNeighbors(std::vector<NeighborSet>& neighbors) const override;
159
160 void applyInitial() override
161 {}
162
163 void linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res) override;
164 void linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res);
165
166 void postSolve(GlobalEqVector& deltaX) override
167 {
168 recoverWellSolutionAndUpdateWellState(deltaX);
169 }
170
171 void postSolveDomain(GlobalEqVector& deltaX, const Domain& domain)
172 {
173 recoverWellSolutionAndUpdateWellStateDomain(deltaX, domain);
174 }
175
177 // </ eWoms auxiliary module stuff>
179
180 template <class Restarter>
181 void deserialize(Restarter& /* res */)
182 {
183 // TODO (?)
184 }
185
190 template <class Restarter>
191 void serialize(Restarter& /* res*/)
192 {
193 // TODO (?)
194 }
195
196 void beginEpisode()
197 {
198 OPM_TIMEBLOCK(beginEpsiode);
199 beginReportStep(simulator_.episodeIndex());
200 }
201
202 void beginTimeStep();
203
204 void beginIteration()
205 {
206 OPM_TIMEBLOCK(beginIteration);
207 assemble(simulator_.model().newtonMethod().numIterations(),
208 simulator_.timeStepSize());
209 }
210
211 void endIteration()
212 { }
213
214 void endTimeStep()
215 {
216 OPM_TIMEBLOCK(endTimeStep);
217 timeStepSucceeded(simulator_.time(), simulator_.timeStepSize());
218 }
219
220 void endEpisode()
221 {
222 endReportStep();
223 }
224
225 void computeTotalRatesForDof(RateVector& rate,
226 unsigned globalIdx) const;
227
228 template <class Context>
229 void computeTotalRatesForDof(RateVector& rate,
230 const Context& context,
231 unsigned spaceIdx,
232 unsigned timeIdx) const;
233
234
235 using WellInterfacePtr = std::shared_ptr<WellInterface<TypeTag> >;
236
237 using BlackoilWellModelGeneric::initFromRestartFile;
238 void initFromRestartFile(const RestartValue& restartValues)
239 {
240 initFromRestartFile(restartValues,
241 this->simulator_.vanguard().transferWTestState(),
242 grid().size(0),
243 param_.use_multisegment_well_);
244 }
245
246 using BlackoilWellModelGeneric::prepareDeserialize;
247 void prepareDeserialize(const int report_step)
248 {
249 prepareDeserialize(report_step, grid().size(0),
250 param_.use_multisegment_well_);
251 }
252
253 data::Wells wellData() const
254 {
255 auto wsrpt = this->wellState()
256 .report(simulator_.vanguard().globalCell().data(),
257 [this](const int well_index) -> bool
258 {
259 return this->wasDynamicallyShutThisTimeStep(well_index);
260 });
261
262 const auto& tracerRates = simulator_.problem().tracerModel().getWellTracerRates();
263 this->assignWellTracerRates(wsrpt, tracerRates);
264
265
266 BlackoilWellModelGuideRates(*this).assignWellGuideRates(wsrpt, this->reportStepIndex());
267 this->assignShutConnections(wsrpt, this->reportStepIndex());
268
269 return wsrpt;
270 }
271
272 data::WellBlockAveragePressures wellBlockAveragePressures() const
273 {
274 return this->computeWellBlockAveragePressures();
275 }
276
277 // subtract Binv(D)rw from r;
278 void apply( BVector& r) const;
279
280 // subtract B*inv(D)*C * x from A*x
281 void apply(const BVector& x, BVector& Ax) const;
282
283 // accumulate the contributions of all Wells in the WellContributions object
284 void getWellContributions(WellContributions& x) const;
285
286 // apply well model with scaling of alpha
287 void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
288
289 // Check if well equations is converged.
290 ConvergenceReport getWellConvergence(const std::vector<Scalar>& B_avg, const bool checkWellGroupControls = false) const;
291
292 // Check if well equations are converged locally.
293 ConvergenceReport getDomainWellConvergence(const Domain& domain,
294 const std::vector<Scalar>& B_avg,
295 DeferredLogger& local_deferredLogger) const;
296
297 const SimulatorReportSingle& lastReport() const;
298
299 void addWellContributions(SparseMatrixAdapter& jacobian) const;
300
301 // add source from wells to the reservoir matrix
302 void addReservoirSourceTerms(GlobalEqVector& residual,
303 std::vector<typename SparseMatrixAdapter::MatrixBlock*>& diagMatAddress) const;
304
305 // called at the beginning of a report step
306 void beginReportStep(const int time_step);
307
308 // it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation()
309 // makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions
310 // twice at the beginning of the time step
313 void calculateExplicitQuantities(DeferredLogger& deferred_logger) const;
314 // some preparation work, mostly related to group control and RESV,
315 // at the beginning of each time step (Not report step)
316 void prepareTimeStep(DeferredLogger& deferred_logger);
317 void initPrimaryVariablesEvaluation() const;
318 void initPrimaryVariablesEvaluationDomain(const Domain& domain) const;
319
320 std::pair<bool, bool>
321 updateWellControls(const bool mandatory_network_balance, DeferredLogger& deferred_logger, const bool relax_network_tolerance = false);
322
323 void updateAndCommunicate(const int reportStepIdx,
324 const int iterationIdx,
325 DeferredLogger& deferred_logger);
326
327 bool updateGroupControls(const Group& group,
328 DeferredLogger& deferred_logger,
329 const int reportStepIdx,
330 const int iterationIdx);
331
332 WellInterfacePtr getWell(const std::string& well_name) const;
333 bool hasWell(const std::string& well_name) const;
334
335 using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
336
337 void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const;
338
339 void addWellPressureEquationsStruct(PressureMatrix& jacobian) const;
340
341 void initGliftEclWellMap(GLiftEclWells &ecl_well_map);
342
344 const std::vector<WellInterfacePtr>& localNonshutWells() const
345 {
346 return well_container_;
347 }
348
349 // prototype for assemble function for ASPIN solveLocal()
350 // will try to merge back to assemble() when done prototyping
351 void assembleDomain(const int iterationIdx,
352 const double dt,
353 const Domain& domain);
354 void updateWellControlsDomain(DeferredLogger& deferred_logger, const Domain& domain);
355
356 void logPrimaryVars() const;
357 std::vector<double> getPrimaryVarsDomain(const Domain& domain) const;
358 void setPrimaryVarsDomain(const Domain& domain, const std::vector<double>& vars);
359
360 void setupDomains(const std::vector<Domain>& domains);
361
362 protected:
363 Simulator& simulator_;
364
365 // a vector of all the wells.
366 std::vector<WellInterfacePtr> well_container_{};
367
368 std::vector<bool> is_cell_perforated_{};
369
370 void initializeWellState(const int timeStepIdx);
371
372 // create the well container
373 void createWellContainer(const int report_step) override;
374
375 WellInterfacePtr
376 createWellPointer(const int wellID,
377 const int report_step) const;
378
379 template <typename WellType>
380 std::unique_ptr<WellType>
381 createTypedWellPointer(const int wellID,
382 const int time_step) const;
383
384 WellInterfacePtr createWellForWellTest(const std::string& well_name, const int report_step, DeferredLogger& deferred_logger) const;
385
386
387 const ModelParameters param_;
388 std::size_t global_num_cells_{};
389 // the number of the cells in the local grid
390 std::size_t local_num_cells_{};
391 double gravity_{};
392 std::vector<double> depth_{};
393 bool alternative_well_rate_init_{};
394
395 std::unique_ptr<RateConverterType> rateConverter_{};
396 std::map<std::string, std::unique_ptr<AverageRegionalPressureType>> regionalAveragePressureCalculator_{};
397
399 {
400 std::optional<typename std::vector<WellInterfacePtr>::size_type> openWellIdx_{};
401 std::size_t wbpCalcIdx_{};
402 };
403
404 std::vector<WBPCalcID> wbpCalcMap_{};
405
406 SimulatorReportSingle last_report_{};
407
408 // Pre-step network solve at static reservoir conditions (group and well states might be updated)
409 void doPreStepNetworkRebalance(DeferredLogger& deferred_logger);
410
411 // used to better efficiency of calcuation
412 mutable BVector scaleAddRes_{};
413
414 std::vector<Scalar> B_avg_{};
415
416 // Keep track of the domain of each well, if using subdomains.
417 std::map<std::string, int> well_domain_;
418
419 const Grid& grid() const
420 { return simulator_.vanguard().grid(); }
421
422 const EquilGrid& equilGrid() const
423 { return simulator_.vanguard().equilGrid(); }
424
425 const EclipseState& eclState() const
426 { return simulator_.vanguard().eclState(); }
427
428 // compute the well fluxes and assemble them in to the reservoir equations as source terms
429 // and in the well equations.
430 void assemble(const int iterationIdx,
431 const double dt);
432
433 // well controls and network pressures affect each other and are solved in an iterative manner.
434 // the function handles one iteration of updating well controls and network pressures.
435 // it is possible to decouple the update of well controls and network pressures further.
436 // the returned two booleans are {continue_due_to_network, well_group_control_changed}, respectively
437 std::pair<bool, bool> updateWellControlsAndNetworkIteration(const bool mandatory_network_balance,
438 const bool relax_network_tolerance,
439 const double dt,
440 DeferredLogger& local_deferredLogger);
441
442 bool updateWellControlsAndNetwork(const bool mandatory_network_balance,
443 const double dt,
444 DeferredLogger& local_deferredLogger);
445
454 void initializeLocalWellStructure(const int reportStepIdx,
455 const bool enableWellPIScaling);
456
460 void initializeGroupStructure(const int reportStepIdx);
461
462 // called at the end of a time step
463 void timeStepSucceeded(const double simulationTime, const double dt);
464
465 // called at the end of a report step
466 void endReportStep();
467
468 // using the solution x to recover the solution xw for wells and applying
469 // xw to update Well State
470 void recoverWellSolutionAndUpdateWellState(const BVector& x);
471
472 // using the solution x to recover the solution xw for wells and applying
473 // xw to update Well State
474 void recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, const Domain& domain);
475
476 // setting the well_solutions_ based on well_state.
477 void updatePrimaryVariables(DeferredLogger& deferred_logger);
478
479 void initializeWBPCalculationService();
480
481 data::WellBlockAveragePressures
482 computeWellBlockAveragePressures() const;
483
485 makeWellSourceEvaluatorFactory(const std::vector<Well>::size_type wellIdx) const;
486
487 void registerOpenWellsForWBPCalculation();
488
489 void updateAverageFormationFactor();
490
491 void computePotentials(const std::size_t widx,
492 const WellState& well_state_copy,
493 std::string& exc_msg,
494 ExceptionType::ExcEnum& exc_type,
495 DeferredLogger& deferred_logger) override;
496
497 const std::vector<double>& wellPerfEfficiencyFactors() const;
498
499 void calculateProductivityIndexValuesShutWells(const int reportStepIdx, DeferredLogger& deferred_logger) override;
500 void calculateProductivityIndexValues(DeferredLogger& deferred_logger) override;
501 void calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
502 DeferredLogger& deferred_logger);
503
504 // The number of components in the model.
505 int numComponents() const;
506
507 int reportStepIndex() const;
508
509 void assembleWellEq(const double dt, DeferredLogger& deferred_logger);
510 void assembleWellEqDomain(const double dt, const Domain& domain, DeferredLogger& deferred_logger);
511
512 void prepareWellsBeforeAssembling(const double dt, DeferredLogger& deferred_logger);
513
514 // TODO: finding a better naming
515 void assembleWellEqWithoutIteration(const double dt, DeferredLogger& deferred_logger);
516
517 bool maybeDoGasLiftOptimize(DeferredLogger& deferred_logger);
518
519 void gasLiftOptimizationStage1(DeferredLogger& deferred_logger,
520 GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
521 GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map);
522
523 // cannot be const since it accesses the non-const WellState
524 void gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
525 DeferredLogger& deferred_logger,
526 GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
527 GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
528 GLiftSyncGroups& groups_to_sync);
529
530 void extractLegacyCellPvtRegionIndex_();
531
532 void extractLegacyDepth_();
533
535 void updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const;
536
537 void wellTesting(const int timeStepIdx, const double simulationTime, DeferredLogger& deferred_logger);
538
539 void calcRates(const int fipnum,
540 const int pvtreg,
541 const std::vector<double>& production_rates,
542 std::vector<double>& resv_coeff) override;
543
544 void calcInjRates(const int fipnum,
545 const int pvtreg,
546 std::vector<double>& resv_coeff) override;
547
548 void computeWellTemperature();
549
550 int compressedIndexForInterior(int cartesian_cell_idx) const override {
551 return simulator_.vanguard().compressedIndexForInterior(cartesian_cell_idx);
552 }
553
554 private:
555 BlackoilWellModel(Simulator& simulator, const PhaseUsage& pu);
556 };
557
558
559} // namespace Opm
560
561#include "BlackoilWellModel_impl.hpp"
562#endif
Helper class for grid instantiation of ECL file-format using problems.
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:83
Class for handling the blackoil well model.
Definition BlackoilWellModel.hpp:96
void initializeGroupStructure(const int reportStepIdx)
Initialize group control modes/constraints and group solution state.
Definition BlackoilWellModel_impl.hpp:335
void serialize(Restarter &)
This method writes the complete state of the well to the harddisk.
Definition BlackoilWellModel.hpp:191
const std::vector< WellInterfacePtr > & localNonshutWells() const
Get list of local nonshut wells.
Definition BlackoilWellModel.hpp:344
void calculateExplicitQuantities(DeferredLogger &deferred_logger) const
Calculating the explict quantities used in the well calculation.
Definition BlackoilWellModel_impl.hpp:1833
void initializeLocalWellStructure(const int reportStepIdx, const bool enableWellPIScaling)
Update rank's notion of intersecting wells and their associate solution variables.
Definition BlackoilWellModel_impl.hpp:290
void updateWellTestState(const double &simulationTime, WellTestState &wellTestState) const
upate the wellTestState related to economic limits
Definition BlackoilWellModel_impl.hpp:2164
int compressedIndexForInterior(int cartesian_cell_idx) const override
get compressed index for interior cells (-1, otherwise
Definition BlackoilWellModel.hpp:550
Definition DeferredLogger.hpp:57
Definition GasLiftSingleWell.hpp:38
std::function< Evaluator()> EvaluatorFactory
Callback for constructing a source term evaluation function on the current MPI rank.
Definition ParallelWBPCalculation.hpp:62
Convert component rates at surface conditions to phase (voidage) rates at reservoir conditions.
Definition RateConverter.hpp:70
Computes hydrocarbon weighed average pressures over regions.
Definition RegionAverageCalculator.hpp:65
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
Solver parameters for the BlackoilModel.
Definition BlackoilModelParameters.hpp:484
Definition BlackoilWellModel.hpp:399
Definition BlackoilPhases.hpp:46
Definition BlackoilWellModel.hpp:84
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34
Representing a part of a grid, in a way suitable for performing local solves.
Definition SubDomain.hpp:62