My Project
Loading...
Searching...
No Matches
MultisegmentWell_impl.hpp
1/*
2 Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3 Copyright 2017 Statoil ASA.
4
5 This file is part of the Open Porous Media project (OPM).
6
7 OPM is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11
12 OPM is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with OPM. If not, see <http://www.gnu.org/licenses/>.
19*/
20
21// Improve IDE experience
22#ifndef OPM_MULTISEGMENTWELL_IMPL_HEADER_INCLUDED
23#define OPM_MULTISEGMENTWELL_IMPL_HEADER_INCLUDED
24
25#ifndef OPM_MULTISEGMENTWELL_HEADER_INCLUDED
26#include <config.h>
27#include <opm/simulators/wells/MultisegmentWell.hpp>
28#endif
29
30#include <opm/common/Exceptions.hpp>
31#include <opm/common/OpmLog/OpmLog.hpp>
32
33#include <opm/input/eclipse/Schedule/MSW/Segment.hpp>
34#include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
35#include <opm/input/eclipse/Schedule/MSW/WellSegments.hpp>
36#include <opm/input/eclipse/Schedule/Well/Connection.hpp>
37#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
38
39#include <opm/input/eclipse/Units/Units.hpp>
40
41#include <opm/material/densead/EvaluationFormat.hpp>
42
43#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
44#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
45#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
46
47#include <algorithm>
48#include <cstddef>
49#include <string>
50
51#if COMPILE_GPU_BRIDGE && (HAVE_CUDA || HAVE_OPENCL)
52#include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
53#endif
54
55namespace Opm
56{
57
58
59 template <typename TypeTag>
60 MultisegmentWell<TypeTag>::
61 MultisegmentWell(const Well& well,
63 const int time_step,
64 const ModelParameters& param,
65 const RateConverterType& rate_converter,
66 const int pvtRegionIdx,
67 const int num_components,
68 const int num_phases,
69 const int index_of_well,
70 const std::vector<PerforationData<Scalar>>& perf_data)
71 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
72 , MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices>&>(*this), pw_info)
73 , regularize_(false)
74 , segment_fluid_initial_(this->numberOfSegments(), std::vector<Scalar>(this->num_components_, 0.0))
75 {
76 // not handling solvent or polymer for now with multisegment well
77 if constexpr (has_solvent) {
78 OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet");
79 }
80
81 if constexpr (has_polymer) {
82 OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet");
83 }
84
85 if constexpr (Base::has_energy) {
86 OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet");
87 }
88
89 if constexpr (Base::has_foam) {
90 OPM_THROW(std::runtime_error, "foam is not supported by multisegment well yet");
91 }
92
93 if constexpr (Base::has_brine) {
94 OPM_THROW(std::runtime_error, "brine is not supported by multisegment well yet");
95 }
96
97 if constexpr (Base::has_watVapor) {
98 OPM_THROW(std::runtime_error, "water evaporation is not supported by multisegment well yet");
99 }
100
101 if(this->rsRvInj() > 0) {
102 OPM_THROW(std::runtime_error,
103 "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."
104 " \n See (WCONINJE item 10 / WCONHIST item 8)");
105 }
106
107 this->thp_update_iterations = true;
108 }
109
110
111
112
113
114 template <typename TypeTag>
115 void
116 MultisegmentWell<TypeTag>::
117 init(const PhaseUsage* phase_usage_arg,
118 const std::vector<Scalar>& depth_arg,
119 const Scalar gravity_arg,
120 const std::vector< Scalar >& B_avg,
121 const bool changed_to_open_this_step)
122 {
123 Base::init(phase_usage_arg, depth_arg, gravity_arg, B_avg, changed_to_open_this_step);
124
125 // TODO: for StandardWell, we need to update the perf depth here using depth_arg.
126 // for MultisegmentWell, it is much more complicated.
127 // It can be specified directly, it can be calculated from the segment depth,
128 // it can also use the cell center, which is the same for StandardWell.
129 // For the last case, should we update the depth with the depth_arg? For the
130 // future, it can be a source of wrong result with Multisegment well.
131 // An indicator from the opm-parser should indicate what kind of depth we should use here.
132
133 // \Note: we do not update the depth here. And it looks like for now, we only have the option to use
134 // specified perforation depth
135 this->initMatrixAndVectors();
136
137 // calculate the depth difference between the perforations and the perforated grid block
138 for (int local_perf_index = 0; local_perf_index < this->number_of_local_perforations_; ++local_perf_index) {
139 // This variable loops over the number_of_local_perforations_ of *this* process, hence it is *local*.
140 const int cell_idx = this->well_cells_[local_perf_index];
141 // Here we need to access the perf_depth_ at the global perforation index though!
142 this->cell_perforation_depth_diffs_[local_perf_index] = depth_arg[cell_idx] - this->perf_depth_[this->pw_info_.localToGlobal(local_perf_index)];
143 }
144 }
145
146
147
148
149
150 template <typename TypeTag>
151 void
152 MultisegmentWell<TypeTag>::
153 updatePrimaryVariables(const Simulator& simulator,
154 const WellState<Scalar>& well_state,
155 DeferredLogger& deferred_logger)
156 {
157 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
158 this->primary_variables_.update(well_state, stop_or_zero_rate_target);
159 }
160
161
162
163
164
165
166 template <typename TypeTag>
167 void
169 updateWellStateWithTarget(const Simulator& simulator,
170 const GroupState<Scalar>& group_state,
171 WellState<Scalar>& well_state,
172 DeferredLogger& deferred_logger) const
173 {
174 Base::updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
175 // scale segment rates based on the wellRates
176 // and segment pressure based on bhp
177 this->scaleSegmentRatesWithWellRates(this->segments_.inlets(),
178 this->segments_.perforations(),
179 well_state);
180 this->scaleSegmentPressuresWithBhp(well_state);
181 }
182
183
184
185
186
187 template <typename TypeTag>
190 getWellConvergence(const Simulator& /* simulator */,
191 const WellState<Scalar>& well_state,
192 const std::vector<Scalar>& B_avg,
193 DeferredLogger& deferred_logger,
194 const bool relax_tolerance) const
195 {
196 return this->MSWEval::getWellConvergence(well_state,
197 B_avg,
198 deferred_logger,
199 this->param_.max_residual_allowed_,
200 this->param_.tolerance_wells_,
201 this->param_.relaxed_tolerance_flow_well_,
202 this->param_.tolerance_pressure_ms_wells_,
203 this->param_.relaxed_tolerance_pressure_ms_well_,
204 relax_tolerance,
205 this->wellIsStopped());
206
207 }
208
209
210
211
212
213 template <typename TypeTag>
214 void
216 apply(const BVector& x, BVector& Ax) const
217 {
218 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
219 return;
220 }
221
222 if (this->param_.matrix_add_well_contributions_) {
223 // Contributions are already in the matrix itself
224 return;
225 }
226
227 this->linSys_.apply(x, Ax);
228 }
229
230
231
232
233
234 template <typename TypeTag>
235 void
237 apply(BVector& r) const
238 {
239 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
240 return;
241 }
242
243 this->linSys_.apply(r);
244 }
245
246
247
248 template <typename TypeTag>
249 void
251 recoverWellSolutionAndUpdateWellState(const Simulator& simulator,
252 const BVector& x,
253 WellState<Scalar>& well_state,
254 DeferredLogger& deferred_logger)
255 {
256 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
257 return;
258 }
259
260 BVectorWell xw(1);
261 this->linSys_.recoverSolutionWell(x, xw);
262 updateWellState(simulator, xw, well_state, deferred_logger);
263 }
264
265
266
267
268
269 template <typename TypeTag>
270 void
272 computeWellPotentials(const Simulator& simulator,
273 const WellState<Scalar>& well_state,
274 std::vector<Scalar>& well_potentials,
275 DeferredLogger& deferred_logger)
276 {
277 const auto [compute_potential, bhp_controlled_well] =
278 this->WellInterfaceGeneric<Scalar>::computeWellPotentials(well_potentials, well_state);
279
280 if (!compute_potential) {
281 return;
282 }
283
284 debug_cost_counter_ = 0;
285 bool converged_implicit = false;
286 if (this->param_.local_well_solver_control_switching_) {
287 converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger);
288 if (!converged_implicit) {
289 deferred_logger.debug("Implicit potential calculations failed for well "
290 + this->name() + ", reverting to original aproach.");
291 }
292 }
293 if (!converged_implicit) {
294 // does the well have a THP related constraint?
295 const auto& summaryState = simulator.vanguard().summaryState();
296 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
297 computeWellRatesAtBhpLimit(simulator, well_potentials, deferred_logger);
298 } else {
299 well_potentials = computeWellPotentialWithTHP(
300 well_state, simulator, deferred_logger);
301 }
302 }
303 deferred_logger.debug("Cost in iterations of finding well potential for well "
304 + this->name() + ": " + std::to_string(debug_cost_counter_));
305
306 this->checkNegativeWellPotentials(well_potentials,
307 this->param_.check_well_operability_,
308 deferred_logger);
309 }
310
311
312
313
314 template<typename TypeTag>
315 void
318 std::vector<Scalar>& well_flux,
319 DeferredLogger& deferred_logger) const
320 {
321 if (this->well_ecl_.isInjector()) {
322 const auto controls = this->well_ecl_.injectionControls(simulator.vanguard().summaryState());
323 computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger);
324 } else {
325 const auto controls = this->well_ecl_.productionControls(simulator.vanguard().summaryState());
326 computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger);
327 }
328 }
329
330 template<typename TypeTag>
331 void
332 MultisegmentWell<TypeTag>::
333 computeWellRatesWithBhp(const Simulator& simulator,
334 const Scalar& bhp,
335 std::vector<Scalar>& well_flux,
336 DeferredLogger& deferred_logger) const
337 {
338 const int np = this->number_of_phases_;
339
340 well_flux.resize(np, 0.0);
341 const bool allow_cf = this->getAllowCrossFlow();
342 const int nseg = this->numberOfSegments();
343 const WellState<Scalar>& well_state = simulator.problem().wellModel().wellState();
344 const auto& ws = well_state.well(this->indexOfWell());
345 auto segments_copy = ws.segments;
346 segments_copy.scale_pressure(bhp);
347 const auto& segment_pressure = segments_copy.pressure;
348 for (int seg = 0; seg < nseg; ++seg) {
349 for (const int perf : this->segments_.perforations()[seg]) {
350 const int local_perf_index = this->pw_info_.globalToLocal(perf);
351 if (local_perf_index < 0) // then the perforation is not on this process
352 continue;
353 const int cell_idx = this->well_cells_[local_perf_index];
354 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
355 // flux for each perforation
356 std::vector<Scalar> mob(this->num_components_, 0.);
357 getMobility(simulator, local_perf_index, mob, deferred_logger);
358 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
359 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
360 const std::vector<Scalar> Tw = this->wellIndex(local_perf_index, intQuants, trans_mult, wellstate_nupcol);
361 const Scalar seg_pressure = segment_pressure[seg];
362 std::vector<Scalar> cq_s(this->num_components_, 0.);
363 Scalar perf_press = 0.0;
364 PerforationRates<Scalar> perf_rates;
365 computePerfRate(intQuants, mob, Tw, seg, perf, seg_pressure,
366 allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
367
368 for(int p = 0; p < np; ++p) {
369 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
370 }
371 }
372 }
373 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
374 }
375
376
377 template<typename TypeTag>
378 void
379 MultisegmentWell<TypeTag>::
380 computeWellRatesWithBhpIterations(const Simulator& simulator,
381 const Scalar& bhp,
382 std::vector<Scalar>& well_flux,
383 DeferredLogger& deferred_logger) const
384 {
385 OPM_TIMEFUNCTION();
386 // creating a copy of the well itself, to avoid messing up the explicit information
387 // during this copy, the only information not copied properly is the well controls
388 MultisegmentWell<TypeTag> well_copy(*this);
389 well_copy.resetDampening();
390
391 well_copy.debug_cost_counter_ = 0;
392
393 // store a copy of the well state, we don't want to update the real well state
394 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
395 const auto& group_state = simulator.problem().wellModel().groupState();
396 auto& ws = well_state_copy.well(this->index_of_well_);
397
398 // Get the current controls.
399 const auto& summary_state = simulator.vanguard().summaryState();
400 auto inj_controls = well_copy.well_ecl_.isInjector()
401 ? well_copy.well_ecl_.injectionControls(summary_state)
402 : Well::InjectionControls(0);
403 auto prod_controls = well_copy.well_ecl_.isProducer()
404 ? well_copy.well_ecl_.productionControls(summary_state) :
405 Well::ProductionControls(0);
406
407 // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
408 if (well_copy.well_ecl_.isInjector()) {
409 inj_controls.bhp_limit = bhp;
410 ws.injection_cmode = Well::InjectorCMode::BHP;
411 } else {
412 prod_controls.bhp_limit = bhp;
413 ws.production_cmode = Well::ProducerCMode::BHP;
414 }
415 ws.bhp = bhp;
416 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
417
418 // initialized the well rates with the potentials i.e. the well rates based on bhp
419 const int np = this->number_of_phases_;
420 bool trivial = true;
421 for (int phase = 0; phase < np; ++phase){
422 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
423 }
424 if (!trivial) {
425 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
426 for (int phase = 0; phase < np; ++phase) {
427 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
428 }
429 }
430 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
431 this->segments_.perforations(),
432 well_state_copy);
433
434 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
435 const double dt = simulator.timeStepSize();
436 // iterate to get a solution at the given bhp.
437 well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
438 deferred_logger);
439
440 // compute the potential and store in the flux vector.
441 well_flux.clear();
442 well_flux.resize(np, 0.0);
443 for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
444 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
445 well_flux[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
446 }
447 debug_cost_counter_ += well_copy.debug_cost_counter_;
448 }
449
450
451
452 template<typename TypeTag>
453 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
454 MultisegmentWell<TypeTag>::
455 computeWellPotentialWithTHP(const WellState<Scalar>& well_state,
456 const Simulator& simulator,
457 DeferredLogger& deferred_logger) const
458 {
459 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
460 const auto& summary_state = simulator.vanguard().summaryState();
461
462 const auto& well = this->well_ecl_;
463 if (well.isInjector()) {
464 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
465 if (bhp_at_thp_limit) {
466 const auto& controls = well.injectionControls(summary_state);
467 const Scalar bhp = std::min(*bhp_at_thp_limit,
468 static_cast<Scalar>(controls.bhp_limit));
469 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
470 deferred_logger.debug("Converged thp based potential calculation for well "
471 + this->name() + ", at bhp = " + std::to_string(bhp));
472 } else {
473 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
474 "Failed in getting converged thp based potential calculation for well "
475 + this->name() + ". Instead the bhp based value is used");
476 const auto& controls = well.injectionControls(summary_state);
477 const Scalar bhp = controls.bhp_limit;
478 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
479 }
480 } else {
481 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
482 well_state, simulator, summary_state, deferred_logger);
483 if (bhp_at_thp_limit) {
484 const auto& controls = well.productionControls(summary_state);
485 const Scalar bhp = std::max(*bhp_at_thp_limit,
486 static_cast<Scalar>(controls.bhp_limit));
487 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
488 deferred_logger.debug("Converged thp based potential calculation for well "
489 + this->name() + ", at bhp = " + std::to_string(bhp));
490 } else {
491 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
492 "Failed in getting converged thp based potential calculation for well "
493 + this->name() + ". Instead the bhp based value is used");
494 const auto& controls = well.productionControls(summary_state);
495 const Scalar bhp = controls.bhp_limit;
496 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
497 }
498 }
499
500 return potentials;
501 }
502
503 template<typename TypeTag>
504 bool
505 MultisegmentWell<TypeTag>::
506 computeWellPotentialsImplicit(const Simulator& simulator,
507 const WellState<Scalar>& well_state,
508 std::vector<Scalar>& well_potentials,
509 DeferredLogger& deferred_logger) const
510 {
511 // Create a copy of the well.
512 // TODO: check if we can avoid taking multiple copies. Call from updateWellPotentials
513 // is allready a copy, but not from other calls.
514 MultisegmentWell<TypeTag> well_copy(*this);
515 well_copy.debug_cost_counter_ = 0;
516
517 // store a copy of the well state, we don't want to update the real well state
518 WellState<Scalar> well_state_copy = well_state;
519 const auto& group_state = simulator.problem().wellModel().groupState();
520 auto& ws = well_state_copy.well(this->index_of_well_);
521
522 // get current controls
523 const auto& summary_state = simulator.vanguard().summaryState();
524 auto inj_controls = well_copy.well_ecl_.isInjector()
525 ? well_copy.well_ecl_.injectionControls(summary_state)
526 : Well::InjectionControls(0);
527 auto prod_controls = well_copy.well_ecl_.isProducer()
528 ? well_copy.well_ecl_.productionControls(summary_state)
529 : Well::ProductionControls(0);
530
531 // prepare/modify well state and control
532 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
533
534 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
535
536 // initialize rates from previous potentials
537 const int np = this->number_of_phases_;
538 bool trivial = true;
539 for (int phase = 0; phase < np; ++phase){
540 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
541 }
542 if (!trivial) {
543 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
544 for (int phase = 0; phase < np; ++phase) {
545 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
546 }
547 }
548 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
549 this->segments_.perforations(),
550 well_state_copy);
551
552 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
553 const double dt = simulator.timeStepSize();
554 // solve equations
555 bool converged = false;
556 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
557 converged = well_copy.solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
558 } else {
559 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
560 }
561
562 // fetch potentials (sign is updated on the outside).
563 well_potentials.clear();
564 well_potentials.resize(np, 0.0);
565 for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
566 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
567 well_potentials[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
568 }
569 debug_cost_counter_ += well_copy.debug_cost_counter_;
570 return converged;
571 }
572
573 template <typename TypeTag>
574 void
575 MultisegmentWell<TypeTag>::
576 solveEqAndUpdateWellState(const Simulator& simulator,
577 WellState<Scalar>& well_state,
578 DeferredLogger& deferred_logger)
579 {
580 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
581
582 // We assemble the well equations, then we check the convergence,
583 // which is why we do not put the assembleWellEq here.
584 try{
585 const BVectorWell dx_well = this->linSys_.solve();
586
587 updateWellState(simulator, dx_well, well_state, deferred_logger);
588 }
589 catch(const NumericalProblem& exp) {
590 // Add information about the well and log to deferred logger
591 // (Logging done inside of solve() method will only be seen if
592 // this is the process with rank zero)
593 deferred_logger.problem("In MultisegmentWell::solveEqAndUpdateWellState for well "
594 + this->name() +": "+exp.what());
595 throw;
596 }
597 }
598
599
600
601
602
603 template <typename TypeTag>
604 void
605 MultisegmentWell<TypeTag>::
606 computePerfCellPressDiffs(const Simulator& simulator)
607 {
608 // We call this function on every process for the number_of_local_perforations_ on that process
609 // Each process updates the pressure for his perforations
610 for (int local_perf_index = 0; local_perf_index < this->number_of_local_perforations_; ++local_perf_index) {
611 // This variable loops over the number_of_local_perforations_ of *this* process, hence it is *local*.
612
613 std::vector<Scalar> kr(this->number_of_phases_, 0.0);
614 std::vector<Scalar> density(this->number_of_phases_, 0.0);
615
616 const int cell_idx = this->well_cells_[local_perf_index];
617 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
618 const auto& fs = intQuants.fluidState();
619
620 Scalar sum_kr = 0.;
621
622 const PhaseUsage& pu = this->phaseUsage();
623 if (pu.phase_used[Water]) {
624 const int water_pos = pu.phase_pos[Water];
625 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
626 sum_kr += kr[water_pos];
627 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
628 }
629
630 if (pu.phase_used[Oil]) {
631 const int oil_pos = pu.phase_pos[Oil];
632 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
633 sum_kr += kr[oil_pos];
634 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
635 }
636
637 if (pu.phase_used[Gas]) {
638 const int gas_pos = pu.phase_pos[Gas];
639 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
640 sum_kr += kr[gas_pos];
641 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
642 }
643
644 assert(sum_kr != 0.);
645
646 // calculate the average density
647 Scalar average_density = 0.;
648 for (int p = 0; p < this->number_of_phases_; ++p) {
649 average_density += kr[p] * density[p];
650 }
651 average_density /= sum_kr;
652
653 this->cell_perforation_pressure_diffs_[local_perf_index] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[local_perf_index];
654 }
655 }
656
657
658
659
660
661 template <typename TypeTag>
662 void
663 MultisegmentWell<TypeTag>::
664 computeInitialSegmentFluids(const Simulator& simulator)
665 {
666 for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
667 // TODO: trying to reduce the times for the surfaceVolumeFraction calculation
668 const Scalar surface_volume = getSegmentSurfaceVolume(simulator, seg).value();
669 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
670 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
671 }
672 }
673 }
674
675
676
677
678
679 template <typename TypeTag>
680 void
681 MultisegmentWell<TypeTag>::
682 updateWellState(const Simulator& simulator,
683 const BVectorWell& dwells,
684 WellState<Scalar>& well_state,
685 DeferredLogger& deferred_logger,
686 const Scalar relaxation_factor)
687 {
688 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
689
690 const Scalar dFLimit = this->param_.dwell_fraction_max_;
691 const Scalar max_pressure_change = this->param_.max_pressure_change_ms_wells_;
692 const bool stop_or_zero_rate_target =
693 this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
694 this->primary_variables_.updateNewton(dwells,
695 relaxation_factor,
696 dFLimit,
697 stop_or_zero_rate_target,
698 max_pressure_change);
699
700 const auto& summary_state = simulator.vanguard().summaryState();
701 this->primary_variables_.copyToWellState(*this, getRefDensity(),
702 well_state,
703 summary_state,
704 deferred_logger);
705
706 {
707 auto& ws = well_state.well(this->index_of_well_);
708 this->segments_.copyPhaseDensities(ws.pu, ws.segments);
709 }
710
711 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
712 }
713
714
715
716
717
718 template <typename TypeTag>
719 void
720 MultisegmentWell<TypeTag>::
721 calculateExplicitQuantities(const Simulator& simulator,
722 const WellState<Scalar>& well_state,
723 DeferredLogger& deferred_logger)
724 {
725 updatePrimaryVariables(simulator, well_state, deferred_logger);
726 computePerfCellPressDiffs(simulator);
727 computeInitialSegmentFluids(simulator);
728 }
729
730
731
732
733
734 template<typename TypeTag>
735 void
736 MultisegmentWell<TypeTag>::
737 updateProductivityIndex(const Simulator& simulator,
738 const WellProdIndexCalculator<Scalar>& wellPICalc,
739 WellState<Scalar>& well_state,
740 DeferredLogger& deferred_logger) const
741 {
742 auto fluidState = [&simulator, this](const int local_perf_index)
743 {
744 const auto cell_idx = this->well_cells_[local_perf_index];
745 return simulator.model()
746 .intensiveQuantities(cell_idx, /*timeIdx=*/ 0).fluidState();
747 };
748
749 const int np = this->number_of_phases_;
750 auto setToZero = [np](Scalar* x) -> void
751 {
752 std::fill_n(x, np, 0.0);
753 };
754
755 auto addVector = [np](const Scalar* src, Scalar* dest) -> void
756 {
757 std::transform(src, src + np, dest, dest, std::plus<>{});
758 };
759
760 auto& ws = well_state.well(this->index_of_well_);
761 auto& perf_data = ws.perf_data;
762 auto* connPI = perf_data.prod_index.data();
763 auto* wellPI = ws.productivity_index.data();
764
765 setToZero(wellPI);
766
767 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
768 auto subsetPerfID = 0;
769
770 for ( const auto& perf : *this->perf_data_){
771 auto allPerfID = perf.ecl_index;
772
773 auto connPICalc = [&wellPICalc, allPerfID](const Scalar mobility) -> Scalar
774 {
775 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
776 };
777
778 std::vector<Scalar> mob(this->num_components_, 0.0);
779 // The subsetPerfID loops over 0 .. this->perf_data_->size().
780 // *(this->perf_data_) contains info about the local processes only,
781 // hence subsetPerfID is a local perf id and we can call getMobility
782 // as well as fluidState directly with that.
783 getMobility(simulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
784
785 const auto& fs = fluidState(subsetPerfID);
786 setToZero(connPI);
787
788 if (this->isInjector()) {
789 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
790 mob, connPI, deferred_logger);
791 }
792 else { // Production or zero flow rate
793 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
794 }
795
796 addVector(connPI, wellPI);
797
798 ++subsetPerfID;
799 connPI += np;
800 }
801
802 // Sum with communication in case of distributed well.
803 const auto& comm = this->parallel_well_info_.communication();
804 if (comm.size() > 1) {
805 comm.sum(wellPI, np);
806 }
807
808 assert (static_cast<int>(subsetPerfID) == this->number_of_local_perforations_ &&
809 "Internal logic error in processing connections for PI/II");
810 }
811
812
813
814
815
816 template<typename TypeTag>
817 typename MultisegmentWell<TypeTag>::Scalar
818 MultisegmentWell<TypeTag>::
819 connectionDensity(const int globalConnIdx,
820 [[maybe_unused]] const int openConnIdx) const
821 {
822 // Simple approximation: Mixture density at reservoir connection is
823 // mixture density at connection's segment.
824
825 const auto segNum = this->wellEcl()
826 .getConnections()[globalConnIdx].segment();
827
828 const auto segIdx = this->wellEcl()
829 .getSegments().segmentNumberToIndex(segNum);
830
831 return this->segments_.density(segIdx).value();
832 }
833
834
835
836
837
838 template<typename TypeTag>
839 void
840 MultisegmentWell<TypeTag>::
841 addWellContributions(SparseMatrixAdapter& jacobian) const
842 {
843 this->linSys_.extract(jacobian);
844 }
845
846
847 template<typename TypeTag>
848 void
849 MultisegmentWell<TypeTag>::
850 addWellPressureEquations(PressureMatrix& jacobian,
851 const BVector& weights,
852 const int pressureVarIndex,
853 const bool use_well_weights,
854 const WellState<Scalar>& well_state) const
855 {
856 // Add the pressure contribution to the cpr system for the well
857 this->linSys_.extractCPRPressureMatrix(jacobian,
858 weights,
859 pressureVarIndex,
860 use_well_weights,
861 *this,
862 this->SPres,
863 well_state);
864 }
865
866
867 template<typename TypeTag>
868 template<class Value>
869 void
870 MultisegmentWell<TypeTag>::
871 computePerfRate(const Value& pressure_cell,
872 const Value& rs,
873 const Value& rv,
874 const std::vector<Value>& b_perfcells,
875 const std::vector<Value>& mob_perfcells,
876 const std::vector<Scalar>& Tw,
877 const int perf,
878 const Value& segment_pressure,
879 const Value& segment_density,
880 const bool& allow_cf,
881 const std::vector<Value>& cmix_s,
882 std::vector<Value>& cq_s,
883 Value& perf_press,
884 PerforationRates<Scalar>& perf_rates,
885 DeferredLogger& deferred_logger) const
886 {
887 const int local_perf_index = this->pw_info_.globalToLocal(perf);
888 if (local_perf_index < 0) // then the perforation is not on this process
889 return;
890
891 // pressure difference between the segment and the perforation
892 const Value perf_seg_press_diff = this->gravity() * segment_density *
893 this->segments_.local_perforation_depth_diff(local_perf_index);
894 // pressure difference between the perforation and the grid cell
895 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index];
896
897 // perforation pressure is the wellbore pressure corrected to perforation depth
898 // (positive sign due to convention in segments_.local_perforation_depth_diff() )
899 perf_press = segment_pressure + perf_seg_press_diff;
900
901 // cell pressure corrected to perforation depth
902 const Value cell_press_at_perf = pressure_cell - cell_perf_press_diff;
903
904 // Pressure drawdown (also used to determine direction of flow)
905 const Value drawdown = cell_press_at_perf - perf_press;
906
907 // producing perforations
908 if (drawdown > 0.0) {
909 // Do nothing if crossflow is not allowed
910 if (!allow_cf && this->isInjector()) {
911 return;
912 }
913
914 // compute component volumetric rates at standard conditions
915 for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
916 const Value cq_p = - Tw[comp_idx] * (mob_perfcells[comp_idx] * drawdown);
917 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
918 }
919
920 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
921 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
922 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
923 const Value cq_s_oil = cq_s[oilCompIdx];
924 const Value cq_s_gas = cq_s[gasCompIdx];
925 cq_s[gasCompIdx] += rs * cq_s_oil;
926 cq_s[oilCompIdx] += rv * cq_s_gas;
927 }
928 } else { // injecting perforations
929 // Do nothing if crossflow is not allowed
930 if (!allow_cf && this->isProducer()) {
931 return;
932 }
933
934 // for injecting perforations, we use total mobility
935 Value total_mob = mob_perfcells[0];
936 for (int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
937 total_mob += mob_perfcells[comp_idx];
938 }
939
940 // compute volume ratio between connection and at standard conditions
941 Value volume_ratio = 0.0;
942 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
943 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
944 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
945 }
946
947 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
948 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
949 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
950
951 // Incorporate RS/RV factors if both oil and gas active
952 // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations
953 // basically, for injecting perforations, the wellbore is the upstreaming side.
954 const Value d = 1.0 - rv * rs;
955
956 if (getValue(d) == 0.0) {
957 OPM_DEFLOG_PROBLEM(NumericalProblem,
958 fmt::format("Zero d value obtained for well {} "
959 "during flux calculation with rs {} and rv {}",
960 this->name(), rs, rv),
961 deferred_logger);
962 }
963
964 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
965 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
966
967 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
968 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
969 } else { // not having gas and oil at the same time
970 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
971 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
972 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
973 }
974 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
975 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
976 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
977 }
978 }
979 // injecting connections total volumerates at standard conditions
980 for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
981 const Value cqt_i = - Tw[componentIdx] * (total_mob * drawdown);
982 Value cqt_is = cqt_i / volume_ratio;
983 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
984 }
985 } // end for injection perforations
986
987 // calculating the perforation solution gas rate and solution oil rates
988 if (this->isProducer()) {
989 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
990 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
991 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
992 // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
993 // s means standard condition, r means reservoir condition
994 // q_os = q_or * b_o + rv * q_gr * b_g
995 // q_gs = q_gr * g_g + rs * q_or * b_o
996 // d = 1.0 - rs * rv
997 // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
998 // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
999
1000 const Scalar d = 1.0 - getValue(rv) * getValue(rs);
1001 // vaporized oil into gas
1002 // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
1003 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
1004 // dissolved of gas in oil
1005 // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
1006 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
1007 }
1008 }
1009 }
1010
1011 template <typename TypeTag>
1012 template<class Value>
1013 void
1014 MultisegmentWell<TypeTag>::
1015 computePerfRate(const IntensiveQuantities& int_quants,
1016 const std::vector<Value>& mob_perfcells,
1017 const std::vector<Scalar>& Tw,
1018 const int seg,
1019 const int perf,
1020 const Value& segment_pressure,
1021 const bool& allow_cf,
1022 std::vector<Value>& cq_s,
1023 Value& perf_press,
1024 PerforationRates<Scalar>& perf_rates,
1025 DeferredLogger& deferred_logger) const
1026
1027 {
1028 auto obtain = [this](const Eval& value)
1029 {
1030 if constexpr (std::is_same_v<Value, Scalar>) {
1031 static_cast<void>(this); // suppress clang warning
1032 return getValue(value);
1033 } else {
1034 return this->extendEval(value);
1035 }
1036 };
1037 auto obtainN = [](const auto& value)
1038 {
1039 if constexpr (std::is_same_v<Value, Scalar>) {
1040 return getValue(value);
1041 } else {
1042 return value;
1043 }
1044 };
1045 const auto& fs = int_quants.fluidState();
1046
1047 const Value pressure_cell = obtain(this->getPerfCellPressure(fs));
1048 const Value rs = obtain(fs.Rs());
1049 const Value rv = obtain(fs.Rv());
1050
1051 // not using number_of_phases_ because of solvent
1052 std::vector<Value> b_perfcells(this->num_components_, 0.0);
1053
1054 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1055 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1056 continue;
1057 }
1058
1059 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1060 b_perfcells[compIdx] = obtain(fs.invB(phaseIdx));
1061 }
1062
1063 std::vector<Value> cmix_s(this->numComponents(), 0.0);
1064 for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
1065 cmix_s[comp_idx] = obtainN(this->primary_variables_.surfaceVolumeFraction(seg, comp_idx));
1066 }
1067
1068 this->computePerfRate(pressure_cell,
1069 rs,
1070 rv,
1071 b_perfcells,
1072 mob_perfcells,
1073 Tw,
1074 perf,
1075 segment_pressure,
1076 obtainN(this->segments_.density(seg)),
1077 allow_cf,
1078 cmix_s,
1079 cq_s,
1080 perf_press,
1081 perf_rates,
1082 deferred_logger);
1083 }
1084
1085 template <typename TypeTag>
1086 void
1087 MultisegmentWell<TypeTag>::
1088 computeSegmentFluidProperties(const Simulator& simulator, DeferredLogger& deferred_logger)
1089 {
1090 // TODO: the concept of phases and components are rather confusing in this function.
1091 // needs to be addressed sooner or later.
1092
1093 // get the temperature for later use. It is only useful when we are not handling
1094 // thermal related simulation
1095 // basically, it is a single value for all the segments
1096
1097 EvalWell temperature;
1098 EvalWell saltConcentration;
1099 // not sure how to handle the pvt region related to segment
1100 // for the current approach, we use the pvt region of the first perforated cell
1101 // although there are some text indicating using the pvt region of the lowest
1102 // perforated cell
1103 // TODO: later to investigate how to handle the pvt region
1104 int pvt_region_index;
1105 {
1106 // using the first perforated cell, so we look for global index 0
1107 const int cell_idx = this->well_cells_[0];
1108 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
1109 const auto& fs = intQuants.fluidState();
1110
1111 // The following broadcast calls are neccessary to ensure that processes that do *not* contain
1112 // the first perforation get the correct temperature, saltConcentration and pvt_region_index
1113 auto fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1114 fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature);
1115 temperature.setValue(fsTemperature);
1116
1117 auto fsSaltConcentration = fs.saltConcentration();
1118 fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration);
1119 saltConcentration = this->extendEval(fsSaltConcentration);
1120
1121 pvt_region_index = fs.pvtRegionIndex();
1122 pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index);
1123 }
1124
1125 this->segments_.computeFluidProperties(temperature,
1126 saltConcentration,
1127 this->primary_variables_,
1128 pvt_region_index,
1129 deferred_logger);
1130 }
1131
1132 template <typename TypeTag>
1133 template<class Value>
1134 void
1135 MultisegmentWell<TypeTag>::
1136 getMobility(const Simulator& simulator,
1137 const int local_perf_index,
1138 std::vector<Value>& mob,
1139 DeferredLogger& deferred_logger) const
1140 {
1141 auto obtain = [this](const Eval& value)
1142 {
1143 if constexpr (std::is_same_v<Value, Scalar>) {
1144 static_cast<void>(this); // suppress clang warning
1145 return getValue(value);
1146 } else {
1147 return this->extendEval(value);
1148 }
1149 };
1150
1151 WellInterface<TypeTag>::getMobility(simulator, local_perf_index, mob, obtain, deferred_logger);
1152
1153 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
1154 const auto perf_ecl_index = this->perforationData()[local_perf_index].ecl_index;
1155 const Connection& con = this->well_ecl_.getConnections()[perf_ecl_index];
1156 const int seg = this->segmentNumberToIndex(con.segment());
1157 // from the reference results, it looks like MSW uses segment pressure instead of BHP here
1158 // Note: this is against the documented definition.
1159 // we can change this depending on what we want
1160 const Scalar segment_pres = this->primary_variables_.getSegmentPressure(seg).value();
1161 const Scalar perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
1162 * this->segments_.local_perforation_depth_diff(local_perf_index);
1163 const Scalar perf_press = segment_pres + perf_seg_press_diff;
1164 const Scalar multiplier = this->getInjMult(local_perf_index, segment_pres, perf_press, deferred_logger);
1165 for (std::size_t i = 0; i < mob.size(); ++i) {
1166 mob[i] *= multiplier;
1167 }
1168 }
1169 }
1170
1171
1172
1173 template<typename TypeTag>
1174 typename MultisegmentWell<TypeTag>::Scalar
1175 MultisegmentWell<TypeTag>::
1176 getRefDensity() const
1177 {
1178 return this->segments_.getRefDensity();
1179 }
1180
1181 template<typename TypeTag>
1182 void
1183 MultisegmentWell<TypeTag>::
1184 checkOperabilityUnderBHPLimit(const WellState<Scalar>& /*well_state*/,
1185 const Simulator& simulator,
1186 DeferredLogger& deferred_logger)
1187 {
1188 const auto& summaryState = simulator.vanguard().summaryState();
1189 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1190 // Crude but works: default is one atmosphere.
1191 // TODO: a better way to detect whether the BHP is defaulted or not
1192 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1193 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1194 // if the BHP limit is not defaulted or the well does not have a THP limit
1195 // we need to check the BHP limit
1196 Scalar total_ipr_mass_rate = 0.0;
1197 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1198 {
1199 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1200 continue;
1201 }
1202
1203 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1204 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1205
1206 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1207 total_ipr_mass_rate += ipr_rate * rho;
1208 }
1209 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1210 this->operability_status_.operable_under_only_bhp_limit = false;
1211 }
1212
1213 // checking whether running under BHP limit will violate THP limit
1214 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1215 // option 1: calculate well rates based on the BHP limit.
1216 // option 2: stick with the above IPR curve
1217 // we use IPR here
1218 std::vector<Scalar> well_rates_bhp_limit;
1219 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1220
1221 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1222 const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
1223 bhp_limit,
1224 this->getRefDensity(),
1225 this->wellEcl().alq_value(summaryState),
1226 thp_limit,
1227 deferred_logger);
1228 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1229 this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1230 }
1231 }
1232 } else {
1233 // defaulted BHP and there is a THP constraint
1234 // default BHP limit is about 1 atm.
1235 // when applied the hydrostatic pressure correction dp,
1236 // most likely we get a negative value (bhp + dp)to search in the VFP table,
1237 // which is not desirable.
1238 // we assume we can operate under defaulted BHP limit and will violate the THP limit
1239 // when operating under defaulted BHP limit.
1240 this->operability_status_.operable_under_only_bhp_limit = true;
1241 this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1242 }
1243 }
1244
1245
1246
1247 template<typename TypeTag>
1248 void
1249 MultisegmentWell<TypeTag>::
1250 updateIPR(const Simulator& simulator, DeferredLogger& deferred_logger) const
1251 {
1252 // TODO: not handling solvent related here for now
1253
1254 // initialize all the values to be zero to begin with
1255 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1256 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1257
1258 const int nseg = this->numberOfSegments();
1259 std::vector<Scalar> seg_dp(nseg, 0.0);
1260 for (int seg = 0; seg < nseg; ++seg) {
1261 // calculating the perforation rate for each perforation that belongs to this segment
1262 const Scalar dp = this->getSegmentDp(seg,
1263 this->segments_.density(seg).value(),
1264 seg_dp);
1265 seg_dp[seg] = dp;
1266 for (const int perf : this->segments_.perforations()[seg]) {
1267 const int local_perf_index = this->pw_info_.globalToLocal(perf);
1268 if (local_perf_index < 0) // then the perforation is not on this process
1269 continue;
1270 std::vector<Scalar> mob(this->num_components_, 0.0);
1271
1272 // TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
1273 getMobility(simulator, local_perf_index, mob, deferred_logger);
1274
1275 const int cell_idx = this->well_cells_[local_perf_index];
1276 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1277 const auto& fs = int_quantities.fluidState();
1278 // pressure difference between the segment and the perforation
1279 const Scalar perf_seg_press_diff = this->segments_.getPressureDiffSegLocalPerf(seg, local_perf_index);
1280 // pressure difference between the perforation and the grid cell
1281 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index];
1282 const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
1283
1284 // calculating the b for the connection
1285 std::vector<Scalar> b_perf(this->num_components_);
1286 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1287 if (!FluidSystem::phaseIsActive(phase)) {
1288 continue;
1289 }
1290 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1291 b_perf[comp_idx] = fs.invB(phase).value();
1292 }
1293
1294 // the pressure difference between the connection and BHP
1295 const Scalar h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
1296 const Scalar pressure_diff = pressure_cell - h_perf;
1297
1298 // do not take into consideration the crossflow here.
1299 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1300 deferred_logger.debug("CROSSFLOW_IPR",
1301 "cross flow found when updateIPR for well " + this->name());
1302 }
1303
1304 // the well index associated with the connection
1305 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
1306 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1307 const std::vector<Scalar> tw_perf = this->wellIndex(local_perf_index, int_quantities, trans_mult, wellstate_nupcol);
1308 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
1309 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
1310 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1311 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
1312 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1313 ipr_b_perf[comp_idx] += tw_mob;
1314 }
1315
1316 // we need to handle the rs and rv when both oil and gas are present
1317 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1318 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1319 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1320 const Scalar rs = (fs.Rs()).value();
1321 const Scalar rv = (fs.Rv()).value();
1322
1323 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1324 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1325
1326 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1327 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1328
1329 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1330 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1331
1332 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1333 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1334 }
1335
1336 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1337 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1338 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1339 }
1340 }
1341 }
1342 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1343 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
1344 }
1345
1346 template<typename TypeTag>
1347 void
1348 MultisegmentWell<TypeTag>::
1349 updateIPRImplicit(const Simulator& simulator,
1350 WellState<Scalar>& well_state,
1351 DeferredLogger& deferred_logger)
1352 {
1353 // Compute IPR based on *converged* well-equation:
1354 // For a component rate r the derivative dr/dbhp is obtained by
1355 // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial bhp_target)
1356 // where Eq(x)=0 is the well equation setup with bhp control and primary variables x
1357
1358 // We shouldn't have zero rates at this stage, but check
1359 bool zero_rates;
1360 auto rates = well_state.well(this->index_of_well_).surface_rates;
1361 zero_rates = true;
1362 for (std::size_t p = 0; p < rates.size(); ++p) {
1363 zero_rates &= rates[p] == 0.0;
1364 }
1365 auto& ws = well_state.well(this->index_of_well_);
1366 if (zero_rates) {
1367 const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
1368 deferred_logger.debug(msg);
1369 /*
1370 // could revert to standard approach here:
1371 updateIPR(simulator, deferred_logger);
1372 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
1373 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
1374 ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
1375 ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
1376 }
1377 return;
1378 */
1379 }
1380 const auto& group_state = simulator.problem().wellModel().groupState();
1381
1382 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
1383 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
1384 //WellState well_state_copy = well_state;
1385 auto inj_controls = Well::InjectionControls(0);
1386 auto prod_controls = Well::ProductionControls(0);
1387 prod_controls.addControl(Well::ProducerCMode::BHP);
1388 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
1389
1390 // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
1391 const auto cmode = ws.production_cmode;
1392 ws.production_cmode = Well::ProducerCMode::BHP;
1393 const double dt = simulator.timeStepSize();
1394 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1395
1396 BVectorWell rhs(this->numberOfSegments());
1397 rhs = 0.0;
1398 rhs[0][SPres] = -1.0;
1399
1400 const BVectorWell x_well = this->linSys_.solve(rhs);
1401 constexpr int num_eq = MSWEval::numWellEq;
1402 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
1403 const EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
1404 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
1405 for (size_t pvIdx = 0; pvIdx < num_eq; ++pvIdx) {
1406 // well primary variable derivatives in EvalWell start at position Indices::numEq
1407 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
1408 }
1409 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
1410 }
1411 // reset cmode
1412 ws.production_cmode = cmode;
1413 }
1414
1415 template<typename TypeTag>
1416 void
1417 MultisegmentWell<TypeTag>::
1418 checkOperabilityUnderTHPLimit(const Simulator& simulator,
1419 const WellState<Scalar>& well_state,
1420 DeferredLogger& deferred_logger)
1421 {
1422 const auto& summaryState = simulator.vanguard().summaryState();
1423 const auto obtain_bhp = this->isProducer()
1424 ? computeBhpAtThpLimitProd(
1425 well_state, simulator, summaryState, deferred_logger)
1426 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1427
1428 if (obtain_bhp) {
1429 this->operability_status_.can_obtain_bhp_with_thp_limit = true;
1430
1431 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1432 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1433
1434 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1435 if (this->isProducer() && *obtain_bhp < thp_limit) {
1436 const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1437 + " bars is SMALLER than thp limit "
1438 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1439 + " bars as a producer for well " + this->name();
1440 deferred_logger.debug(msg);
1441 }
1442 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1443 const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1444 + " bars is LARGER than thp limit "
1445 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1446 + " bars as a injector for well " + this->name();
1447 deferred_logger.debug(msg);
1448 }
1449 } else {
1450 // Shutting wells that can not find bhp value from thp
1451 // when under THP control
1452 this->operability_status_.can_obtain_bhp_with_thp_limit = false;
1453 this->operability_status_.obey_bhp_limit_with_thp_limit = false;
1454 if (!this->wellIsStopped()) {
1455 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1456 deferred_logger.debug(" could not find bhp value at thp limit "
1457 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1458 + " bar for well " + this->name() + ", the well might need to be closed ");
1459 }
1460 }
1461 }
1462
1463
1464
1465
1466
1467 template<typename TypeTag>
1468 bool
1469 MultisegmentWell<TypeTag>::
1470 iterateWellEqWithControl(const Simulator& simulator,
1471 const double dt,
1472 const Well::InjectionControls& inj_controls,
1473 const Well::ProductionControls& prod_controls,
1474 WellState<Scalar>& well_state,
1475 const GroupState<Scalar>& group_state,
1476 DeferredLogger& deferred_logger)
1477 {
1478 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
1479
1480 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1481
1482 {
1483 // getWellFiniteResiduals returns false for nan/inf residuals
1484 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1485 if(!isFinite)
1486 return false;
1487 }
1488
1489 std::vector<std::vector<Scalar> > residual_history;
1490 std::vector<Scalar> measure_history;
1491 int it = 0;
1492 // relaxation factor
1493 Scalar relaxation_factor = 1.;
1494 bool converged = false;
1495 bool relax_convergence = false;
1496 this->regularize_ = false;
1497 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1498
1499 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1500
1501 BVectorWell dx_well;
1502 try{
1503 dx_well = this->linSys_.solve();
1504 }
1505 catch(const NumericalProblem& exp) {
1506 // Add information about the well and log to deferred logger
1507 // (Logging done inside of solve() method will only be seen if
1508 // this is the process with rank zero)
1509 deferred_logger.problem("In MultisegmentWell::iterateWellEqWithControl for well "
1510 + this->name() +": "+exp.what());
1511 throw;
1512 }
1513
1514 if (it > this->param_.strict_inner_iter_wells_) {
1515 relax_convergence = true;
1516 this->regularize_ = true;
1517 }
1518
1519 const auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1520 if (report.converged()) {
1521 converged = true;
1522 break;
1523 }
1524
1525 {
1526 // getFinteWellResiduals returns false for nan/inf residuals
1527 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1528 if (!isFinite)
1529 return false;
1530
1531 residual_history.push_back(residuals);
1532 measure_history.push_back(this->getResidualMeasureValue(well_state,
1533 residual_history[it],
1534 this->param_.tolerance_wells_,
1535 this->param_.tolerance_pressure_ms_wells_,
1536 deferred_logger) );
1537 }
1538 bool min_relaxation_reached = this->update_relaxation_factor(measure_history, relaxation_factor, this->regularize_, deferred_logger);
1539 if (min_relaxation_reached || this->repeatedStagnation(measure_history, this->regularize_, deferred_logger)) {
1540 // try last attempt with relaxed tolerances
1541 const auto reportStag = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, true);
1542 if (reportStag.converged()) {
1543 converged = true;
1544 std::string message = fmt::format("Well stagnates/oscillates but {} manages to get converged with relaxed tolerances in {} inner iterations."
1545 ,this->name(), it);
1546 deferred_logger.debug(message);
1547 } else {
1548 converged = false;
1549 }
1550 break;
1551 }
1552 updateWellState(simulator, dx_well, well_state, deferred_logger, relaxation_factor);
1553 }
1554
1555 // TODO: we should decide whether to keep the updated well_state, or recover to use the old well_state
1556 if (converged) {
1557 std::ostringstream sstr;
1558 sstr << " Well " << this->name() << " converged in " << it << " inner iterations.";
1559 if (relax_convergence)
1560 sstr << " (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ << " iterations)";
1561 deferred_logger.debug(sstr.str());
1562 } else {
1563 std::ostringstream sstr;
1564 sstr << " Well " << this->name() << " did not converge in " << it << " inner iterations.";
1565#define EXTRA_DEBUG_MSW 0
1566#if EXTRA_DEBUG_MSW
1567 sstr << "***** Outputting the residual history for well " << this->name() << " during inner iterations:";
1568 for (int i = 0; i < it; ++i) {
1569 const auto& residual = residual_history[i];
1570 sstr << " residual at " << i << "th iteration ";
1571 for (const auto& res : residual) {
1572 sstr << " " << res;
1573 }
1574 sstr << " " << measure_history[i] << " \n";
1575 }
1576#endif
1577#undef EXTRA_DEBUG_MSW
1578 deferred_logger.debug(sstr.str());
1579 }
1580
1581 return converged;
1582 }
1583
1584
1585 template<typename TypeTag>
1586 bool
1587 MultisegmentWell<TypeTag>::
1588 iterateWellEqWithSwitching(const Simulator& simulator,
1589 const double dt,
1590 const Well::InjectionControls& inj_controls,
1591 const Well::ProductionControls& prod_controls,
1592 WellState<Scalar>& well_state,
1593 const GroupState<Scalar>& group_state,
1594 DeferredLogger& deferred_logger,
1595 const bool fixed_control /*false*/,
1596 const bool fixed_status /*false*/)
1597 {
1598 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1599
1600 {
1601 // getWellFiniteResiduals returns false for nan/inf residuals
1602 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1603 if(!isFinite)
1604 return false;
1605 }
1606
1607 std::vector<std::vector<Scalar> > residual_history;
1608 std::vector<Scalar> measure_history;
1609 int it = 0;
1610 // relaxation factor
1611 Scalar relaxation_factor = 1.;
1612 bool converged = false;
1613 bool relax_convergence = false;
1614 this->regularize_ = false;
1615 const auto& summary_state = simulator.vanguard().summaryState();
1616
1617 // Always take a few (more than one) iterations after a switch before allowing a new switch
1618 // The optimal number here is subject to further investigation, but it has been observerved
1619 // that unless this number is >1, we may get stuck in a cycle
1620 const int min_its_after_switch = 3;
1621 int its_since_last_switch = min_its_after_switch;
1622 int switch_count= 0;
1623 // if we fail to solve eqs, we reset status/operability before leaving
1624 const auto well_status_orig = this->wellStatus_;
1625 const auto operability_orig = this->operability_status_;
1626 // don't allow opening wells that are stopped from schedule or has a stopped well state
1627 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
1628 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
1629 // don't allow switcing for wells under zero rate target or requested fixed status and control
1630 const bool allow_switching = !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
1631 (!fixed_control || !fixed_status) && allow_open;
1632 bool final_check = false;
1633 // well needs to be set operable or else solving/updating of re-opened wells is skipped
1634 this->operability_status_.resetOperability();
1635 this->operability_status_.solvable = true;
1636
1637 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1638 ++its_since_last_switch;
1639 if (allow_switching && its_since_last_switch >= min_its_after_switch){
1640 const Scalar wqTotal = this->primary_variables_.getWQTotal().value();
1641 bool changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
1642 inj_controls, prod_controls, wqTotal,
1643 deferred_logger, fixed_control,
1644 fixed_status);
1645 if (changed) {
1646 its_since_last_switch = 0;
1647 ++switch_count;
1648 }
1649 if (!changed && final_check) {
1650 break;
1651 } else {
1652 final_check = false;
1653 }
1654 }
1655
1656 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls,
1657 well_state, group_state, deferred_logger);
1658
1659 const BVectorWell dx_well = this->linSys_.solve();
1660
1661 if (it > this->param_.strict_inner_iter_wells_) {
1662 relax_convergence = true;
1663 this->regularize_ = true;
1664 }
1665
1666 const auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1667 converged = report.converged();
1668 if (this->parallel_well_info_.communication().size() > 1 && converged != this->parallel_well_info_.communication().min(converged)) {
1669 OPM_THROW(std::runtime_error, fmt::format("Misalignment of the parallel simulation run in iterateWellEqWithSwitching - the well calculation succeeded on rank {} but failed on other ranks.", this->parallel_well_info_.communication().rank()));
1670 }
1671 if (converged) {
1672 // if equations are sufficiently linear they might converge in less than min_its_after_switch
1673 // in this case, make sure all constraints are satisfied before returning
1674 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
1675 final_check = true;
1676 its_since_last_switch = min_its_after_switch;
1677 } else {
1678 break;
1679 }
1680 }
1681
1682 // getFinteWellResiduals returns false for nan/inf residuals
1683 {
1684 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1685 if (!isFinite)
1686 return false;
1687
1688 residual_history.push_back(residuals);
1689 }
1690
1691 if (!converged) {
1692 measure_history.push_back(this->getResidualMeasureValue(well_state,
1693 residual_history[it],
1694 this->param_.tolerance_wells_,
1695 this->param_.tolerance_pressure_ms_wells_,
1696 deferred_logger));
1697 bool min_relaxation_reached = this->update_relaxation_factor(measure_history, relaxation_factor, this->regularize_, deferred_logger);
1698 if (min_relaxation_reached || this->repeatedStagnation(measure_history, this->regularize_, deferred_logger)) {
1699 converged = false;
1700 break;
1701 }
1702 }
1703 updateWellState(simulator, dx_well, well_state, deferred_logger, relaxation_factor);
1704 }
1705
1706 if (converged) {
1707 if (allow_switching){
1708 // update operability if status change
1709 const bool is_stopped = this->wellIsStopped();
1710 if (this->wellHasTHPConstraints(summary_state)){
1711 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
1712 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
1713 } else {
1714 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
1715 }
1716 }
1717 std::string message = fmt::format(" Well {} converged in {} inner iterations ("
1718 "{} control/status switches).", this->name(), it, switch_count);
1719 if (relax_convergence) {
1720 message.append(fmt::format(" (A relaxed tolerance was used after {} iterations)",
1721 this->param_.strict_inner_iter_wells_));
1722 }
1723 deferred_logger.debug(message);
1724 } else {
1725 this->wellStatus_ = well_status_orig;
1726 this->operability_status_ = operability_orig;
1727 const std::string message = fmt::format(" Well {} did not converge in {} inner iterations ("
1728 "{} control/status switches).", this->name(), it, switch_count);
1729 deferred_logger.debug(message);
1730 this->primary_variables_.outputLowLimitPressureSegments(deferred_logger);
1731 }
1732
1733 return converged;
1734 }
1735
1736
1737 template<typename TypeTag>
1738 void
1739 MultisegmentWell<TypeTag>::
1740 assembleWellEqWithoutIteration(const Simulator& simulator,
1741 const double dt,
1742 const Well::InjectionControls& inj_controls,
1743 const Well::ProductionControls& prod_controls,
1744 WellState<Scalar>& well_state,
1745 const GroupState<Scalar>& group_state,
1746 DeferredLogger& deferred_logger)
1747 {
1748 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1749
1750 // update the upwinding segments
1751 this->segments_.updateUpwindingSegments(this->primary_variables_);
1752
1753 // calculate the fluid properties needed.
1754 computeSegmentFluidProperties(simulator, deferred_logger);
1755
1756 // clear all entries
1757 this->linSys_.clear();
1758
1759 auto& ws = well_state.well(this->index_of_well_);
1760 ws.phase_mixing_rates.fill(0.0);
1761
1762 // for the black oil cases, there will be four equations,
1763 // the first three of them are the mass balance equations, the last one is the pressure equations.
1764 //
1765 // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
1766
1767 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1768
1769 const int nseg = this->numberOfSegments();
1770
1771 for (int seg = 0; seg < nseg; ++seg) {
1772 // calculating the perforation rate for each perforation that belongs to this segment
1773 const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg);
1774 auto& perf_data = ws.perf_data;
1775 auto& perf_rates = perf_data.phase_rates;
1776 auto& perf_press_state = perf_data.pressure;
1777 for (const int perf : this->segments_.perforations()[seg]) {
1778 const int local_perf_index = this->pw_info_.globalToLocal(perf);
1779 if (local_perf_index < 0) // then the perforation is not on this process
1780 continue;
1781 const int cell_idx = this->well_cells_[local_perf_index];
1782 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1783 std::vector<EvalWell> mob(this->num_components_, 0.0);
1784 getMobility(simulator, local_perf_index, mob, deferred_logger);
1785 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quants, cell_idx);
1786 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1787 const std::vector<Scalar> Tw = this->wellIndex(local_perf_index, int_quants, trans_mult, wellstate_nupcol);
1788 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1789 EvalWell perf_press;
1790 PerforationRates<Scalar> perfRates;
1791 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
1792 allow_cf, cq_s, perf_press, perfRates, deferred_logger);
1793
1794 // updating the solution gas rate and solution oil rate
1795 if (this->isProducer()) {
1796 ws.phase_mixing_rates[ws.dissolved_gas] += perfRates.dis_gas;
1797 ws.phase_mixing_rates[ws.vaporized_oil] += perfRates.vap_oil;
1798 perf_data.phase_mixing_rates[local_perf_index][ws.dissolved_gas] = perfRates.dis_gas;
1799 perf_data.phase_mixing_rates[local_perf_index][ws.vaporized_oil] = perfRates.vap_oil;
1800 }
1801
1802 // store the perf pressure and rates
1803 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1804 perf_rates[local_perf_index*this->number_of_phases_ + this->modelCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1805 }
1806 perf_press_state[local_perf_index] = perf_press.value();
1807
1808 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1809 // the cq_s entering mass balance equations need to consider the efficiency factors.
1810 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1811
1812 this->connectionRates_[local_perf_index][comp_idx] = Base::restrictEval(cq_s_effective);
1813
1814 MultisegmentWellAssemble(*this).
1815 assemblePerforationEq(seg, local_perf_index, comp_idx, cq_s_effective, this->linSys_);
1816 }
1817 }
1818 }
1819 // Accumulate dissolved gas and vaporized oil flow rates across all ranks sharing this well.
1820 {
1821 const auto& comm = this->parallel_well_info_.communication();
1822 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
1823 }
1824
1825 if (this->parallel_well_info_.communication().size() > 1) {
1826 // accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed)
1827 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
1828 }
1829 for (int seg = 0; seg < nseg; ++seg) {
1830 // calculating the accumulation term
1831 // TODO: without considering the efficiency factor for now
1832 {
1833 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(simulator, seg);
1834
1835 // Add a regularization_factor to increase the accumulation term
1836 // This will make the system less stiff and help convergence for
1837 // difficult cases
1838 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
1839 // for each component
1840 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1841 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx)
1842 - segment_fluid_initial_[seg][comp_idx]) / dt;
1843 MultisegmentWellAssemble(*this).
1844 assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_);
1845 }
1846 }
1847 // considering the contributions due to flowing out from the segment
1848 {
1849 const int seg_upwind = this->segments_.upwinding_segment(seg);
1850 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1851 const EvalWell segment_rate =
1852 this->primary_variables_.getSegmentRateUpwinding(seg,
1853 seg_upwind,
1854 comp_idx) *
1855 this->well_efficiency_factor_;
1856 MultisegmentWellAssemble(*this).
1857 assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
1858 }
1859 }
1860
1861 // considering the contributions from the inlet segments
1862 {
1863 for (const int inlet : this->segments_.inlets()[seg]) {
1864 const int inlet_upwind = this->segments_.upwinding_segment(inlet);
1865 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1866 const EvalWell inlet_rate =
1867 this->primary_variables_.getSegmentRateUpwinding(inlet,
1868 inlet_upwind,
1869 comp_idx) *
1870 this->well_efficiency_factor_;
1871 MultisegmentWellAssemble(*this).
1872 assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
1873 }
1874 }
1875 }
1876
1877 // the fourth equation, the pressure drop equation
1878 if (seg == 0) { // top segment, pressure equation is the control equation
1879 const auto& summaryState = simulator.vanguard().summaryState();
1880 const Schedule& schedule = simulator.vanguard().schedule();
1881 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1882 MultisegmentWellAssemble(*this).
1883 assembleControlEq(well_state,
1884 group_state,
1885 schedule,
1886 summaryState,
1887 inj_controls,
1888 prod_controls,
1889 getRefDensity(),
1890 this->primary_variables_,
1891 this->linSys_,
1892 stopped_or_zero_target,
1893 deferred_logger);
1894 } else {
1895 const UnitSystem& unit_system = simulator.vanguard().eclState().getDeckUnitSystem();
1896 const auto& summary_state = simulator.vanguard().summaryState();
1897 this->assemblePressureEq(seg, unit_system, well_state, summary_state, this->param_.use_average_density_ms_wells_, deferred_logger);
1898 }
1899 }
1900
1901 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1902 this->linSys_.createSolver();
1903 }
1904
1905
1906
1907
1908 template<typename TypeTag>
1909 bool
1910 MultisegmentWell<TypeTag>::
1911 openCrossFlowAvoidSingularity(const Simulator& simulator) const
1912 {
1913 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1914 }
1915
1916
1917 template<typename TypeTag>
1918 bool
1919 MultisegmentWell<TypeTag>::
1920 allDrawDownWrongDirection(const Simulator& simulator) const
1921 {
1922 bool all_drawdown_wrong_direction = true;
1923 const int nseg = this->numberOfSegments();
1924
1925 for (int seg = 0; seg < nseg; ++seg) {
1926 const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
1927 for (const int perf : this->segments_.perforations()[seg]) {
1928 const int local_perf_index = this->pw_info_.globalToLocal(perf);
1929 if (local_perf_index < 0) // then the perforation is not on this process
1930 continue;
1931
1932 const int cell_idx = this->well_cells_[local_perf_index];
1933 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1934 const auto& fs = intQuants.fluidState();
1935
1936 // pressure difference between the segment and the perforation
1937 const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegLocalPerf(seg, local_perf_index);
1938 // pressure difference between the perforation and the grid cell
1939 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index];
1940
1941 const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
1942 const Scalar perf_press = pressure_cell - cell_perf_press_diff;
1943 // Pressure drawdown (also used to determine direction of flow)
1944 // TODO: not 100% sure about the sign of the seg_perf_press_diff
1945 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1946
1947 // for now, if there is one perforation can produce/inject in the correct
1948 // direction, we consider this well can still produce/inject.
1949 // TODO: it can be more complicated than this to cause wrong-signed rates
1950 if ( (drawdown < 0. && this->isInjector()) ||
1951 (drawdown > 0. && this->isProducer()) ) {
1952 all_drawdown_wrong_direction = false;
1953 break;
1954 }
1955 }
1956 }
1957 const auto& comm = this->parallel_well_info_.communication();
1958 if (comm.size() > 1)
1959 {
1960 all_drawdown_wrong_direction =
1961 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1962 }
1963
1964 return all_drawdown_wrong_direction;
1965 }
1966
1967
1968
1969
1970 template<typename TypeTag>
1971 void
1972 MultisegmentWell<TypeTag>::
1973 updateWaterThroughput(const double /*dt*/, WellState<Scalar>& /*well_state*/) const
1974 {
1975 }
1976
1977
1978
1979
1980
1981 template<typename TypeTag>
1982 typename MultisegmentWell<TypeTag>::EvalWell
1983 MultisegmentWell<TypeTag>::
1984 getSegmentSurfaceVolume(const Simulator& simulator, const int seg_idx) const
1985 {
1986 EvalWell temperature;
1987 EvalWell saltConcentration;
1988 int pvt_region_index;
1989 {
1990 // using the pvt region of first perforated cell, so we look for global index 0
1991 // TODO: it should be a member of the WellInterface, initialized properly
1992 const int cell_idx = this->well_cells_[0];
1993 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
1994 const auto& fs = intQuants.fluidState();
1995
1996 // The following broadcast calls are neccessary to ensure that processes that do *not* contain
1997 // the first perforation get the correct temperature, saltConcentration and pvt_region_index
1998 auto fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1999 fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature);
2000 temperature.setValue(fsTemperature);
2001
2002 auto fsSaltConcentration = fs.saltConcentration();
2003 fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration);
2004 saltConcentration = this->extendEval(fsSaltConcentration);
2005
2006 pvt_region_index = fs.pvtRegionIndex();
2007 pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index);
2008 }
2009
2010 return this->segments_.getSurfaceVolume(temperature,
2011 saltConcentration,
2012 this->primary_variables_,
2013 pvt_region_index,
2014 seg_idx);
2015 }
2016
2017
2018 template<typename TypeTag>
2019 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2020 MultisegmentWell<TypeTag>::
2021 computeBhpAtThpLimitProd(const WellState<Scalar>& well_state,
2022 const Simulator& simulator,
2023 const SummaryState& summary_state,
2024 DeferredLogger& deferred_logger) const
2025 {
2026 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
2027 simulator,
2028 summary_state,
2029 this->getALQ(well_state),
2030 deferred_logger,
2031 /*iterate_if_no_solution */ true);
2032 }
2033
2034
2035
2036 template<typename TypeTag>
2037 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2038 MultisegmentWell<TypeTag>::
2039 computeBhpAtThpLimitProdWithAlq(const Simulator& simulator,
2040 const SummaryState& summary_state,
2041 const Scalar alq_value,
2042 DeferredLogger& deferred_logger,
2043 bool iterate_if_no_solution) const
2044 {
2045 OPM_TIMEFUNCTION();
2046 // Make the frates() function.
2047 auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) {
2048 // Not solving the well equations here, which means we are
2049 // calculating at the current Fg/Fw values of the
2050 // well. This does not matter unless the well is
2051 // crossflowing, and then it is likely still a good
2052 // approximation.
2053 std::vector<Scalar> rates(3);
2054 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2055 return rates;
2056 };
2057
2058 auto bhpAtLimit = WellBhpThpCalculator(*this).
2059 computeBhpAtThpLimitProd(frates,
2060 summary_state,
2061 this->maxPerfPress(simulator),
2062 this->getRefDensity(),
2063 alq_value,
2064 this->getTHPConstraint(summary_state),
2065 deferred_logger);
2066
2067 if (bhpAtLimit)
2068 return bhpAtLimit;
2069
2070 if (!iterate_if_no_solution)
2071 return std::nullopt;
2072
2073 auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) {
2074 // Solver the well iterations to see if we are
2075 // able to get a solution with an update
2076 // solution
2077 std::vector<Scalar> rates(3);
2078 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2079 return rates;
2080 };
2081
2082 return WellBhpThpCalculator(*this).
2083 computeBhpAtThpLimitProd(fratesIter,
2084 summary_state,
2085 this->maxPerfPress(simulator),
2086 this->getRefDensity(),
2087 alq_value,
2088 this->getTHPConstraint(summary_state),
2089 deferred_logger);
2090 }
2091
2092 template<typename TypeTag>
2093 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2094 MultisegmentWell<TypeTag>::
2095 computeBhpAtThpLimitInj(const Simulator& simulator,
2096 const SummaryState& summary_state,
2097 DeferredLogger& deferred_logger) const
2098 {
2099 // Make the frates() function.
2100 auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) {
2101 // Not solving the well equations here, which means we are
2102 // calculating at the current Fg/Fw values of the
2103 // well. This does not matter unless the well is
2104 // crossflowing, and then it is likely still a good
2105 // approximation.
2106 std::vector<Scalar> rates(3);
2107 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2108 return rates;
2109 };
2110
2111 auto bhpAtLimit = WellBhpThpCalculator(*this).
2112 computeBhpAtThpLimitInj(frates,
2113 summary_state,
2114 this->getRefDensity(),
2115 0.05,
2116 100,
2117 false,
2118 deferred_logger);
2119
2120 if (bhpAtLimit)
2121 return bhpAtLimit;
2122
2123 auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) {
2124 // Solver the well iterations to see if we are
2125 // able to get a solution with an update
2126 // solution
2127 std::vector<Scalar> rates(3);
2128 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2129 return rates;
2130 };
2131
2132 return WellBhpThpCalculator(*this).
2133 computeBhpAtThpLimitInj(fratesIter,
2134 summary_state,
2135 this->getRefDensity(),
2136 0.05,
2137 100,
2138 false,
2139 deferred_logger);
2140 }
2141
2142
2143
2144
2145
2146 template<typename TypeTag>
2147 typename MultisegmentWell<TypeTag>::Scalar
2148 MultisegmentWell<TypeTag>::
2149 maxPerfPress(const Simulator& simulator) const
2150 {
2151 Scalar max_pressure = 0.0;
2152 const int nseg = this->numberOfSegments();
2153 for (int seg = 0; seg < nseg; ++seg) {
2154 for (const int perf : this->segments_.perforations()[seg]) {
2155 const int local_perf_index = this->pw_info_.globalToLocal(perf);
2156 if (local_perf_index < 0) // then the perforation is not on this process
2157 continue;
2158
2159 const int cell_idx = this->well_cells_[local_perf_index];
2160 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2161 const auto& fs = int_quants.fluidState();
2162 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2163 max_pressure = std::max(max_pressure, pressure_cell);
2164 }
2165 }
2166 return max_pressure;
2167 }
2168
2169
2170
2171
2172
2173 template<typename TypeTag>
2174 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
2176 computeCurrentWellRates(const Simulator& simulator,
2177 DeferredLogger& deferred_logger) const
2178 {
2179 // Calculate the rates that follow from the current primary variables.
2180 std::vector<Scalar> well_q_s(this->num_components_, 0.0);
2181 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2182 const int nseg = this->numberOfSegments();
2183 for (int seg = 0; seg < nseg; ++seg) {
2184 // calculating the perforation rate for each perforation that belongs to this segment
2185 const Scalar seg_pressure = getValue(this->primary_variables_.getSegmentPressure(seg));
2186 for (const int perf : this->segments_.perforations()[seg]) {
2187 const int local_perf_index = this->pw_info_.globalToLocal(perf);
2188 if (local_perf_index < 0) // then the perforation is not on this process
2189 continue;
2190
2191 const int cell_idx = this->well_cells_[local_perf_index];
2192 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2193 std::vector<Scalar> mob(this->num_components_, 0.0);
2194 getMobility(simulator, local_perf_index, mob, deferred_logger);
2195 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quants, cell_idx);
2196 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2197 const std::vector<Scalar> Tw = this->wellIndex(local_perf_index, int_quants, trans_mult, wellstate_nupcol);
2198 std::vector<Scalar> cq_s(this->num_components_, 0.0);
2199 Scalar perf_press = 0.0;
2200 PerforationRates<Scalar> perf_rates;
2201 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
2202 allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
2203 for (int comp = 0; comp < this->num_components_; ++comp) {
2204 well_q_s[comp] += cq_s[comp];
2205 }
2206 }
2207 }
2208 const auto& comm = this->parallel_well_info_.communication();
2209 if (comm.size() > 1)
2210 {
2211 comm.sum(well_q_s.data(), well_q_s.size());
2212 }
2213 return well_q_s;
2214 }
2215
2216
2217 template <typename TypeTag>
2218 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
2220 getPrimaryVars() const
2221 {
2222 const int num_seg = this->numberOfSegments();
2223 constexpr int num_eq = MSWEval::numWellEq;
2224 std::vector<Scalar> retval(num_seg * num_eq);
2225 for (int ii = 0; ii < num_seg; ++ii) {
2226 const auto& pv = this->primary_variables_.value(ii);
2227 std::copy(pv.begin(), pv.end(), retval.begin() + ii * num_eq);
2228 }
2229 return retval;
2230 }
2231
2232
2233
2234
2235 template <typename TypeTag>
2236 int
2237 MultisegmentWell<TypeTag>::
2238 setPrimaryVars(typename std::vector<Scalar>::const_iterator it)
2239 {
2240 const int num_seg = this->numberOfSegments();
2241 constexpr int num_eq = MSWEval::numWellEq;
2242 std::array<Scalar, num_eq> tmp;
2243 for (int ii = 0; ii < num_seg; ++ii) {
2244 const auto start = it + ii * num_eq;
2245 std::copy(start, start + num_eq, tmp.begin());
2246 this->primary_variables_.setValue(ii, tmp);
2247 }
2248 return num_seg * num_eq;
2249 }
2250
2251} // namespace Opm
2252
2253#endif
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
Definition DeferredLogger.hpp:57
Definition GroupState.hpp:43
Definition MultisegmentWell.hpp:38
void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition MultisegmentWell_impl.hpp:216
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
Vanguard & vanguard()
Return a reference to the grid manager of simulation.
Definition simulator.hh:273
Definition WellInterfaceGeneric.hpp:53
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:66
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition phaseUsageFromDeck.cpp:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Definition PerforationData.hpp:41