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 
22 #include <opm/autodiff/MSWellHelpers.hpp>
23 
24 namespace Opm
25 {
26 
27 
28  template <typename TypeTag>
29  MultisegmentWell<TypeTag>::
30  MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
31  : Base(well, time_step, wells, param)
32  , segment_perforations_(numberOfSegments())
33  , segment_inlets_(numberOfSegments())
34  , cell_perforation_depth_diffs_(number_of_perforations_, 0.0)
35  , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0)
36  , perforation_segment_depth_diffs_(number_of_perforations_, 0.0)
37  , segment_comp_initial_(numberOfSegments(), std::vector<double>(numComponents(), 0.0))
38  , segment_densities_(numberOfSegments(), 0.0)
39  , segment_viscosities_(numberOfSegments(), 0.0)
40  , segment_mass_rates_(numberOfSegments(), 0.0)
41  , segment_depth_diffs_(numberOfSegments(), 0.0)
42  {
43  // not handling solvent or polymer for now with multisegment well
44  if (has_solvent) {
45  OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet");
46  }
47 
48  if (has_polymer) {
49  OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet");
50  }
51  // since we decide to use the SegmentSet from the well parser. we can reuse a lot from it.
52  // for other facilities needed but not available from parser, we need to process them here
53 
54  // initialize the segment_perforations_
55  const CompletionSet& completion_set = well_ecl_->getCompletions(current_step_);
56  for (int perf = 0; perf < number_of_perforations_; ++perf) {
57  const Completion& completion = completion_set.get(perf);
58  const int segment_number = completion.getSegmentNumber();
59  const int segment_index = segmentNumberToIndex(segment_number);
60  segment_perforations_[segment_index].push_back(perf);
61  }
62 
63  // initialize the segment_inlets_
64  for (int seg = 0; seg < numberOfSegments(); ++seg) {
65  const Segment& segment = segmentSet()[seg];
66  const int segment_number = segment.segmentNumber();
67  const int outlet_segment_number = segment.outletSegment();
68  if (outlet_segment_number > 0) {
69  const int segment_index = segmentNumberToIndex(segment_number);
70  const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
71  segment_inlets_[outlet_segment_index].push_back(segment_index);
72  }
73  }
74 
75  // callcuate the depth difference between perforations and their segments
76  perf_depth_.resize(number_of_perforations_, 0.);
77  for (int seg = 0; seg < numberOfSegments(); ++seg) {
78  const double segment_depth = segmentSet()[seg].depth();
79  for (const int perf : segment_perforations_[seg]) {
80  perf_depth_[perf] = completion_set.get(perf).getCenterDepth();
81  perforation_segment_depth_diffs_[perf] = perf_depth_[perf] - segment_depth;
82  }
83  }
84 
85  // calculating the depth difference between the segment and its oulet_segments
86  // for the top segment, we will make its zero unless we find other purpose to use this value
87  for (int seg = 1; seg < numberOfSegments(); ++seg) {
88  const double segment_depth = segmentSet()[seg].depth();
89  const int outlet_segment_number = segmentSet()[seg].outletSegment();
90  const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)];
91  const double outlet_depth = outlet_segment.depth();
92  segment_depth_diffs_[seg] = segment_depth - outlet_depth;
93  }
94  }
95 
96 
97 
98 
99 
100  template <typename TypeTag>
101  void
102  MultisegmentWell<TypeTag>::
103  init(const PhaseUsage* phase_usage_arg,
104  const std::vector<bool>* active_arg,
105  const std::vector<double>& depth_arg,
106  const double gravity_arg,
107  const int num_cells)
108  {
109  Base::init(phase_usage_arg, active_arg,
110  depth_arg, gravity_arg, num_cells);
111 
112  // TODO: for StandardWell, we need to update the perf depth here using depth_arg.
113  // for MultisegmentWell, it is much more complicated.
114  // It can be specified directly, it can be calculated from the segment depth,
115  // it can also use the cell center, which is the same for StandardWell.
116  // For the last case, should we update the depth with the depth_arg? For the
117  // future, it can be a source of wrong result with Multisegment well.
118  // An indicator from the opm-parser should indicate what kind of depth we should use here.
119 
120  // \Note: we do not update the depth here. And it looks like for now, we only have the option to use
121  // specified perforation depth
122  initMatrixAndVectors(num_cells);
123 
124  // calcuate the depth difference between the perforations and the perforated grid block
125  for (int perf = 0; perf < number_of_perforations_; ++perf) {
126  const int cell_idx = well_cells_[perf];
127  cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - perf_depth_[perf];
128  }
129  }
130 
131 
132 
133 
134 
135  template <typename TypeTag>
136  void
137  MultisegmentWell<TypeTag>::
138  initMatrixAndVectors(const int num_cells) const
139  {
140  duneB_.setBuildMode( OffDiagMatWell::row_wise );
141  duneC_.setBuildMode( OffDiagMatWell::row_wise );
142  duneD_.setBuildMode( DiagMatWell::row_wise );
143 
144  // set the size and patterns for all the matrices and vectors
145  // [A C^T [x = [ res
146  // B D] x_well] res_well]
147 
148  // calculatiing the NNZ for duneD_
149  // NNZ = number_of_segments + 2 * (number_of_inlets / number_of_outlets)
150  {
151  int nnz_d = numberOfSegments();
152  for (const std::vector<int>& inlets : segment_inlets_) {
153  nnz_d += 2 * inlets.size();
154  }
155  duneD_.setSize(numberOfSegments(), numberOfSegments(), nnz_d);
156  }
157  duneB_.setSize(numberOfSegments(), num_cells, number_of_perforations_);
158  duneC_.setSize(numberOfSegments(), num_cells, number_of_perforations_);
159 
160  // we need to add the off diagonal ones
161  for (auto row = duneD_.createbegin(), end = duneD_.createend(); row != end; ++row) {
162  // the number of the row corrspnds to the segment now
163  const int seg = row.index();
164  // adding the item related to outlet relation
165  const Segment& segment = segmentSet()[seg];
166  const int outlet_segment_number = segment.outletSegment();
167  if (outlet_segment_number > 0) { // if there is a outlet_segment
168  const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
169  row.insert(outlet_segment_index);
170  }
171 
172  // Add nonzeros for diagonal
173  row.insert(seg);
174 
175  // insert the item related to its inlets
176  for (const int& inlet : segment_inlets_[seg]) {
177  row.insert(inlet);
178  }
179  }
180 
181  // make the C matrix
182  for (auto row = duneC_.createbegin(), end = duneC_.createend(); row != end; ++row) {
183  // the number of the row corresponds to the segment number now.
184  for (const int& perf : segment_perforations_[row.index()]) {
185  const int cell_idx = well_cells_[perf];
186  row.insert(cell_idx);
187  }
188  }
189 
190  // make the B^T matrix
191  for (auto row = duneB_.createbegin(), end = duneB_.createend(); row != end; ++row) {
192  // the number of the row corresponds to the segment number now.
193  for (const int& perf : segment_perforations_[row.index()]) {
194  const int cell_idx = well_cells_[perf];
195  row.insert(cell_idx);
196  }
197  }
198 
199  resWell_.resize( numberOfSegments() );
200 
201  primary_variables_.resize(numberOfSegments());
202  primary_variables_evaluation_.resize(numberOfSegments());
203  }
204 
205 
206 
207 
208 
209  template <typename TypeTag>
210  void
211  MultisegmentWell<TypeTag>::
212  initPrimaryVariablesEvaluation() const
213  {
214  for (int seg = 0; seg < numberOfSegments(); ++seg) {
215  for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
216  primary_variables_evaluation_[seg][eq_idx] = 0.0;
217  primary_variables_evaluation_[seg][eq_idx].setValue(primary_variables_[seg][eq_idx]);
218  primary_variables_evaluation_[seg][eq_idx].setDerivative(eq_idx + numEq, 1.0);
219  }
220  }
221  }
222 
223 
224 
225 
226 
227  template <typename TypeTag>
228  void
229  MultisegmentWell<TypeTag>::
230  assembleWellEq(Simulator& ebosSimulator,
231  const double dt,
232  WellState& well_state,
233  bool only_wells)
234  {
235 
236  const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
237  if (use_inner_iterations) {
238  iterateWellEquations(ebosSimulator, dt, well_state);
239  }
240 
241  assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells);
242  }
243 
244 
245 
246 
247 
248  template <typename TypeTag>
249  void
251  updateWellStateWithTarget(const int current,
252  WellState& well_state) const
253  {
254  // Updating well state bas on well control
255  // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
256  const double target = well_controls_iget_target(well_controls_, current);
257  const double* distr = well_controls_iget_distr(well_controls_, current);
258  switch (well_controls_iget_type(well_controls_, current)) {
259  case BHP: {
260  well_state.bhp()[index_of_well_] = target;
261  const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
262  well_state.segPress()[top_segment_index] = well_state.bhp()[index_of_well_];
263  // TODO: similar to the way below to handle THP
264  // we should not something related to thp here when there is thp constraint
265  break;
266  }
267 
268  case THP: {
269  well_state.thp()[index_of_well_] = target;
270 
271  /* const Opm::PhaseUsage& pu = phaseUsage();
272  std::vector<double> rates(3, 0.0);
273  if (active()[ Water ]) {
274  rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ];
275  }
276  if (active()[ Oil ]) {
277  rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ];
278  }
279  if (active()[ Gas ]) {
280  rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
281  } */
282 
283  // const int table_id = well_controls_iget_vfp(well_controls_, current);
284  // const double& thp = well_controls_iget_target(well_controls_, current);
285  // const double& alq = well_controls_iget_alq(well_controls_, current);
286 
287  // TODO: implement calculateBhpFromThp function
288  // well_state.bhp()[index_of_well_] = calculateBhpFromThp(rates, current);
289  // also the top segment pressure
290  /* const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
291  well_state.segPress()[top_segment_index] = well_state.bhp()[index_of_well_]; */
292  break;
293  }
294 
295  case RESERVOIR_RATE: // intentional fall-through
296  case SURFACE_RATE:
297  // for the update of the rates, after we update the well rates, we can try to scale
298  // the segment rates and perforation rates with the same factor
299  // or the other way, we can use the same approach like the initialization of the well state,
300  // we devide the well rates to update the perforation rates, then we update the segment rates based
301  // on the perforation rates.
302  // the second way is safer, since if the well control is changed, the first way will not guarantee the consistence
303  // of between the segment rates and peforation rates
304 
305  // checking the number of the phases under control
306  int numPhasesWithTargetsUnderThisControl = 0;
307  for (int phase = 0; phase < number_of_phases_; ++phase) {
308  if (distr[phase] > 0.0) {
309  numPhasesWithTargetsUnderThisControl += 1;
310  }
311  }
312 
313  assert(numPhasesWithTargetsUnderThisControl > 0);
314 
315  if (well_type_ == INJECTOR) {
316  // assign target value as initial guess for injectors
317  // only handles single phase control at the moment
318  assert(numPhasesWithTargetsUnderThisControl == 1);
319 
320  for (int phase = 0; phase < number_of_phases_; ++phase) {
321  if (distr[phase] > 0.) {
322  well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target / distr[phase];
323  } else {
324  well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = 0.;
325  }
326  }
327 
328  initSegmentRatesWithWellRates(well_state);
329 
330  } else if (well_type_ == PRODUCER) {
331  // update the rates of phases under control based on the target,
332  // and also update rates of phases not under control to keep the rate ratio,
333  // assuming the mobility ratio does not change for the production wells
334  double original_rates_under_phase_control = 0.0;
335  for (int phase = 0; phase < number_of_phases_; ++phase) {
336  if (distr[phase] > 0.0) {
337  original_rates_under_phase_control += well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] * distr[phase];
338  }
339  }
340 
341  if (original_rates_under_phase_control != 0.0 ) {
342  double scaling_factor = target / original_rates_under_phase_control;
343 
344  for (int phase = 0; phase < number_of_phases_; ++phase) {
345  well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] *= scaling_factor;
346 
347  // scaling the segment rates with the same way with well rates
348  const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
349  for (int seg = 0; seg < numberOfSegments(); ++seg) {
350  well_state.segRates()[number_of_phases_ * (seg + top_segment_index) + phase] *= scaling_factor;
351  }
352  }
353  } else { // scaling factor is not well defined when original_rates_under_phase_control is zero
354  // separating targets equally between phases under control
355  const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl;
356  for (int phase = 0; phase < number_of_phases_; ++phase) {
357  if (distr[phase] > 0.0) {
358  well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target_rate_divided / distr[phase];
359  } else {
360  // this only happens for SURFACE_RATE control
361  well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target_rate_divided;
362  }
363  }
364  initSegmentRatesWithWellRates(well_state);
365  }
366  }
367 
368  break;
369  } // end of switch
370 
371  updatePrimaryVariables(well_state);
372  }
373 
374 
375 
376 
377 
378  template <typename TypeTag>
379  void
381  initSegmentRatesWithWellRates(WellState& well_state) const
382  {
383  for (int phase = 0; phase < number_of_phases_; ++phase) {
384  const double perf_phaserate = well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] / number_of_perforations_;
385  for (int perf = 0; perf < number_of_perforations_; ++perf) {
386  well_state.perfPhaseRates()[number_of_phases_ * (first_perf_ + perf) + phase] = perf_phaserate;
387  }
388  }
389 
390  const std::vector<double> perforation_rates(well_state.perfPhaseRates().begin() + number_of_phases_ * first_perf_,
391  well_state.perfPhaseRates().begin() +
392  number_of_phases_ * (first_perf_ + number_of_perforations_) );
393  std::vector<double> segment_rates;
394  WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, number_of_phases_,
395  0, segment_rates);
396  const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
397  std::copy(segment_rates.begin(), segment_rates.end(),
398  well_state.segRates().begin() + number_of_phases_ * top_segment_index );
399  // we need to check the top segment rates should be same with the well rates
400  }
401 
402 
403 
404 
405 
406  template <typename TypeTag>
407  typename MultisegmentWell<TypeTag>::ConvergenceReport
409  getWellConvergence(const std::vector<double>& B_avg) const
410  {
411  // assert((int(B_avg.size()) == numComponents()) || has_polymer);
412  assert( (int(B_avg.size()) == numComponents()) );
413 
414  // checking if any residual is NaN or too large. The two large one is only handled for the well flux
415  std::vector<std::vector<double>> abs_residual(numberOfSegments(), std::vector<double>(numWellEq, 0.0));
416  for (int seg = 0; seg < numberOfSegments(); ++seg) {
417  for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
418  abs_residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]);
419  }
420  }
421 
422  std::vector<double> maximum_residual(numWellEq, 0.0);
423 
424  ConvergenceReport report;
425  // TODO: the following is a little complicated, maybe can be simplified in some way?
426  for (int seg = 0; seg < numberOfSegments(); ++seg) {
427  for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
428  if (eq_idx < numComponents()) { // phase or component mass equations
429  const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
430  // TODO: the report can not handle the segment number yet.
431  if (std::isnan(flux_residual)) {
432  report.nan_residual_found = true;
433  const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx));
434  const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name};
435  report.nan_residual_wells.push_back(problem_well);
436  } else if (flux_residual > param_.max_residual_allowed_) {
437  report.too_large_residual_found = true;
438  const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx));
439  const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name};
440  report.nan_residual_wells.push_back(problem_well);
441  } else { // it is a normal residual
442  if (flux_residual > maximum_residual[eq_idx]) {
443  maximum_residual[eq_idx] = flux_residual;
444  }
445  }
446  } else { // pressure equation
447  // TODO: we should distinguish the rate control equations, bhp control equations
448  // and the oridnary pressure equations
449  const double pressure_residal = abs_residual[seg][eq_idx];
450  const std::string eq_name("Pressure");
451  if (std::isnan(pressure_residal)) {
452  report.nan_residual_found = true;
453  const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name};
454  report.nan_residual_wells.push_back(problem_well);
455  } else if (std::isinf(pressure_residal)) {
456  report.too_large_residual_found = true;
457  const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name};
458  report.nan_residual_wells.push_back(problem_well);
459  } else { // it is a normal residual
460  if (pressure_residal > maximum_residual[eq_idx]) {
461  maximum_residual[eq_idx] = pressure_residal;
462  }
463  }
464  }
465  }
466  }
467 
468  if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found
469  // check convergence for flux residuals
470  for ( int comp_idx = 0; comp_idx < numComponents(); ++comp_idx)
471  {
472  report.converged = report.converged && (maximum_residual[comp_idx] < param_.tolerance_wells_);
473  }
474 
475  report.converged = report.converged && (maximum_residual[SPres] < param_.tolerance_pressure_ms_wells_);
476  } else { // abnormal values found and no need to check the convergence
477  report.converged = false;
478  }
479 
480  return report;
481  }
482 
483 
484 
485 
486 
487  template <typename TypeTag>
488  void
490  apply(const BVector& x, BVector& Ax) const
491  {
492  BVectorWell Bx(duneB_.N());
493 
494  duneB_.mv(x, Bx);
495 
496  // invDBx = duneD^-1 * Bx_
497 #if HAVE_UMFPACK
498  const BVectorWell invDBx = mswellhelpers::invDXDirect(duneD_, Bx);
499 #else
500  const BVectorWell invDBx = mswellhelpers::invDX(duneD_, Bx);
501 #endif
502 
503  // Ax = Ax - duneC_^T * invDBx
504  duneC_.mmtv(invDBx,Ax);
505  }
506 
507 
508 
509 
510 
511  template <typename TypeTag>
512  void
514  apply(BVector& r) const
515  {
516  // invDrw_ = duneD^-1 * resWell_
517 #if HAVE_UMFPACK
518  const BVectorWell invDrw = mswellhelpers::invDXDirect(duneD_, resWell_);
519 #else
520  const BVectorWell invDrw = mswellhelpers::invDX(duneD_, resWell_);
521 #endif // HAVE_UMFPACK
522  // r = r - duneC_^T * invDrw
523  duneC_.mmtv(invDrw, r);
524  }
525 
526 
527 
528 
529 
530  template <typename TypeTag>
531  void
534  WellState& well_state) const
535  {
536  BVectorWell xw(1);
537  recoverSolutionWell(x, xw);
538  updateWellState(xw, false, well_state);
539  }
540 
541 
542 
543 
544 
545  template <typename TypeTag>
546  void
548  computeWellPotentials(const Simulator& /* ebosSimulator */,
549  const WellState& /* well_state */,
550  std::vector<double>& /* well_potentials */)
551  {
552  OPM_THROW(std::runtime_error, "well potential calculation for multisegment wells is not supported yet");
553  }
554 
555 
556 
557 
558 
559  template <typename TypeTag>
560  void
562  updatePrimaryVariables(const WellState& well_state) const
563  {
564  // TODO: to test using rate conversion coefficients to see if it will be better than
565  // this default one
566 
567  // the index of the top segment in the WellState
568  const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
569  const std::vector<double>& segment_rates = well_state.segRates();
570  const PhaseUsage& pu = phaseUsage();
571 
572  for (int seg = 0; seg < numberOfSegments(); ++seg) {
573  // calculate the total rate for each segment
574  double total_seg_rate = 0.0;
575  const int seg_index = top_segment_index + seg;
576  // the segment pressure
577  primary_variables_[seg][SPres] = well_state.segPress()[seg_index];
578  // TODO: under what kind of circustances, the following will be wrong?
579  // the definition of g makes the gas phase is always the last phase
580  for (int p = 0; p < number_of_phases_; p++) {
581  total_seg_rate += scalingFactor(p) * segment_rates[number_of_phases_ * seg_index + p];
582  }
583 
584  primary_variables_[seg][GTotal] = total_seg_rate;
585  if (std::abs(total_seg_rate) > 0.) {
586  if (active()[Water]) {
587  const int water_pos = pu.phase_pos[Water];
588  primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate;
589  }
590  if (active()[Gas]) {
591  const int gas_pos = pu.phase_pos[Gas];
592  primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate;
593  }
594  } else { // total_seg_rate == 0
595  if (well_type_ == INJECTOR) {
596  // only single phase injection handled
597  const double* distr = well_controls_get_current_distr(well_controls_);
598  if (active()[Water]) {
599  if (distr[pu.phase_pos[Water]] > 0.0) {
600  primary_variables_[seg][WFrac] = 1.0;
601  } else {
602  primary_variables_[seg][WFrac] = 0.0;
603  }
604  }
605 
606  if (active()[Gas]) {
607  if (distr[pu.phase_pos[Gas]] > 0.0) {
608  primary_variables_[seg][GFrac] = 1.0;
609  } else {
610  primary_variables_[seg][GFrac] = 0.0;
611  }
612  }
613  } else if (well_type_ == PRODUCER) { // producers
614  if (active()[Water]) {
615  primary_variables_[seg][WFrac] = 1.0 / number_of_phases_;
616  }
617 
618  if (active()[Gas]) {
619  primary_variables_[seg][GFrac] = 1.0 / number_of_phases_;
620  }
621  }
622  }
623  }
624  }
625 
626 
627 
628 
629 
630  template <typename TypeTag>
631  void
632  MultisegmentWell<TypeTag>::
633  recoverSolutionWell(const BVector& x, BVectorWell& xw) const
634  {
635  BVectorWell resWell = resWell_;
636  // resWell = resWell - B * x
637  duneB_.mmv(x, resWell);
638  // xw = D^-1 * resWell
639 #if HAVE_UMFPACK
640  xw = mswellhelpers::invDXDirect(duneD_, resWell);
641 #else
642  xw = mswellhelpers::invDX(duneD_, resWell);
643 #endif // HAVE_UMFPACK
644  }
645 
646 
647 
648 
649 
650  template <typename TypeTag>
651  void
652  MultisegmentWell<TypeTag>::
653  solveEqAndUpdateWellState(WellState& well_state)
654  {
655 #if HAVE_UMFPACK
656  // We assemble the well equations, then we check the convergence,
657  // which is why we do not put the assembleWellEq here.
658  const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
659 #else
660  const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_);
661 #endif
662 
663  updateWellState(dx_well, false, well_state);
664  }
665 
666 
667 
668 
669 
670  template <typename TypeTag>
671  void
672  MultisegmentWell<TypeTag>::
673  computePerfCellPressDiffs(const Simulator& ebosSimulator)
674  {
675  for (int perf = 0; perf < number_of_perforations_; ++perf) {
676 
677  std::vector<double> kr(number_of_phases_, 0.0);
678  std::vector<double> density(number_of_phases_, 0.0);
679 
680  const int cell_idx = well_cells_[perf];
681  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
682  const auto& fs = intQuants.fluidState();
683 
684  double sum_kr = 0.;
685 
686  const PhaseUsage& pu = phaseUsage();
687  if (pu.phase_used[BlackoilPhases::Aqua]) {
688  const int water_pos = pu.phase_pos[BlackoilPhases::Aqua];
689  kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
690  sum_kr += kr[water_pos];
691  density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
692  }
693 
694  if (pu.phase_used[BlackoilPhases::Liquid]) {
695  const int oil_pos = pu.phase_pos[BlackoilPhases::Liquid];
696  kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
697  sum_kr += kr[oil_pos];
698  density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
699  }
700 
701  if (pu.phase_used[BlackoilPhases::Vapour]) {
702  const int gas_pos = pu.phase_pos[BlackoilPhases::Vapour];
703  kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
704  sum_kr += kr[gas_pos];
705  density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
706  }
707 
708  assert(sum_kr != 0.);
709 
710  // calculate the average density
711  double average_density = 0.;
712  for (int p = 0; p < number_of_phases_; ++p) {
713  average_density += kr[p] * density[p];
714  }
715  average_density /= sum_kr;
716 
717  cell_perforation_pressure_diffs_[perf] = gravity_ * average_density * cell_perforation_depth_diffs_[perf];
718  }
719  }
720 
721 
722 
723 
724 
725  template <typename TypeTag>
726  void
727  MultisegmentWell<TypeTag>::
728  computeInitialComposition()
729  {
730  for (int seg = 0; seg < numberOfSegments(); ++seg) {
731  // TODO: probably it should be numWellEq -1 more accurately,
732  // while by meaning it should be num_comp
733  for (int comp_idx = 0; comp_idx < numComponents(); ++comp_idx) {
734  segment_comp_initial_[seg][comp_idx] = surfaceVolumeFraction(seg, comp_idx).value();
735  }
736  }
737  }
738 
739 
740 
741 
742 
743  template <typename TypeTag>
744  void
745  MultisegmentWell<TypeTag>::
746  updateWellState(const BVectorWell& dwells,
747  const bool inner_iteration,
748  WellState& well_state) const
749  {
750  const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
751 
752  const double relaxation_factor = (use_inner_iterations && inner_iteration) ? 0.2 : 1.0;
753 
754  const double dFLimit = param_.dwell_fraction_max_;
755  const double max_pressure_change = param_.max_pressure_change_ms_wells_;
756  const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
757 
758  for (int seg = 0; seg < numberOfSegments(); ++seg) {
759  if (active()[ Water ]) {
760  const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
761  const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit);
762  primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
763  }
764 
765  if (active()[ Gas ]) {
766  const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
767  const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit);
768  primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
769  }
770 
771  // handling the overshooting or undershooting of the fractions
772  processFractions(seg);
773 
774  // update the segment pressure
775  {
776  const int sign = dwells[seg][SPres] > 0.? 1 : -1;
777  const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change);
778  primary_variables_[seg][SPres] = old_primary_variables[seg][SPres] - dx_limited;
779  }
780 
781  // update the total rate // TODO: should we have a limitation of the total rate change?
782  {
783  primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal];
784  }
785 
786  }
787 
788  updateWellStateFromPrimaryVariables(well_state);
789  }
790 
791 
792 
793 
794 
795  template <typename TypeTag>
796  void
797  MultisegmentWell<TypeTag>::
798  calculateExplicitQuantities(const Simulator& ebosSimulator,
799  const WellState& /* well_state */)
800  {
801  computePerfCellPressDiffs(ebosSimulator);
802  computeInitialComposition();
803  }
804 
805 
806 
807 
808 
809  template <typename TypeTag>
810  const SegmentSet&
811  MultisegmentWell<TypeTag>::
812  segmentSet() const
813  {
814  return well_ecl_->getSegmentSet(current_step_);
815  }
816 
817 
818 
819 
820 
821  template <typename TypeTag>
822  int
825  {
826  return segmentSet().numberSegment();
827  }
828 
829 
830 
831 
832 
833  template <typename TypeTag>
834  int
836  numberOfPerforations() const
837  {
838  return segmentSet().number_of_perforations_;
839  }
840 
841 
842 
843 
844 
845  template <typename TypeTag>
846  WellSegment::CompPressureDropEnum
847  MultisegmentWell<TypeTag>::
848  compPressureDrop() const
849  {
850  return segmentSet().compPressureDrop();
851  }
852 
853 
854 
855 
856 
857  template <typename TypeTag>
858  WellSegment::MultiPhaseModelEnum
859  MultisegmentWell<TypeTag>::
860  multiphaseModel() const
861  {
862  return segmentSet().multiPhaseModel();
863  }
864 
865 
866 
867 
868 
869  template <typename TypeTag>
870  int
871  MultisegmentWell<TypeTag>::
872  segmentNumberToIndex(const int segment_number) const
873  {
874  return segmentSet().segmentNumberToIndex(segment_number);
875  }
876 
877 
878 
879 
880 
881  template <typename TypeTag>
882  typename MultisegmentWell<TypeTag>::EvalWell
883  MultisegmentWell<TypeTag>::
884  volumeFraction(const int seg, const int comp_idx) const
885  {
886  const PhaseUsage& pu = phaseUsage();
887 
888  if (active()[Water] && comp_idx == pu.phase_pos[Water]) {
889  return primary_variables_evaluation_[seg][WFrac];
890  }
891 
892  if (active()[Gas] && comp_idx == pu.phase_pos[Gas]) {
893  return primary_variables_evaluation_[seg][GFrac];
894  }
895 
896  // Oil fraction
897  EvalWell oil_fraction = 1.0;
898  if (active()[Water]) {
899  oil_fraction -= primary_variables_evaluation_[seg][WFrac];
900  }
901 
902  if (active()[Gas]) {
903  oil_fraction -= primary_variables_evaluation_[seg][GFrac];
904  }
905  /* if (has_solvent) {
906  oil_fraction -= primary_variables_evaluation_[seg][SFrac];
907  } */
908  return oil_fraction;
909  }
910 
911 
912 
913 
914 
915  template <typename TypeTag>
916  typename MultisegmentWell<TypeTag>::EvalWell
917  MultisegmentWell<TypeTag>::
918  volumeFractionScaled(const int seg, const int comp_idx) const
919  {
920  // For reservoir rate control, the distr in well control is used for the
921  // rate conversion coefficients. For the injection well, only the distr of the injection
922  // phase is not zero.
923  const double scale = scalingFactor(comp_idx);
924  if (scale > 0.) {
925  return volumeFraction(seg, comp_idx) / scale;
926  }
927 
928  return volumeFraction(seg, comp_idx);
929  }
930 
931 
932 
933 
934 
935  template <typename TypeTag>
936  typename MultisegmentWell<TypeTag>::EvalWell
937  MultisegmentWell<TypeTag>::
938  surfaceVolumeFraction(const int seg, const int comp_idx) const
939  {
940  EvalWell sum_volume_fraction_scaled = 0.;
941  const int num_comp = numComponents();
942  for (int idx = 0; idx < num_comp; ++idx) {
943  sum_volume_fraction_scaled += volumeFractionScaled(seg, idx);
944  }
945 
946  assert(sum_volume_fraction_scaled.value() != 0.);
947 
948  return volumeFractionScaled(seg, comp_idx) / sum_volume_fraction_scaled;
949  }
950 
951 
952 
953 
954 
955  template <typename TypeTag>
956  void
957  MultisegmentWell<TypeTag>::
958  computePerfRate(const IntensiveQuantities& int_quants,
959  const std::vector<EvalWell>& mob_perfcells,
960  const int seg,
961  const int perf,
962  const EvalWell& segment_pressure,
963  const bool& allow_cf,
964  std::vector<EvalWell>& cq_s) const
965  {
966  const int num_comp = numComponents();
967  std::vector<EvalWell> cmix_s(num_comp, 0.0);
968 
969  // the composition of the components inside wellbore
970  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
971  cmix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
972  }
973 
974  const auto& fs = int_quants.fluidState();
975 
976  const EvalWell pressure_cell = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
977  const EvalWell rs = extendEval(fs.Rs());
978  const EvalWell rv = extendEval(fs.Rv());
979 
980  // not using number_of_phases_ because of solvent
981  std::vector<EvalWell> b_perfcells(num_comp, 0.0);
982 
983  for (int phase = 0; phase < number_of_phases_; ++phase) {
984  const int phase_idx_ebos = flowPhaseToEbosPhaseIdx(phase);
985  b_perfcells[phase] = extendEval(fs.invB(phase_idx_ebos));
986  }
987 
988  // pressure difference between the segment and the perforation
989  const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf];
990  // pressure difference between the perforation and the grid cell
991  const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf];
992 
993  // Pressure drawdown (also used to determine direction of flow)
994  // TODO: not 100% sure about the sign of the seg_perf_press_diff
995  const EvalWell drawdown = (pressure_cell + cell_perf_press_diff) - (segment_pressure + perf_seg_press_diff);
996 
997  const Opm::PhaseUsage& pu = phaseUsage();
998 
999  // producing perforations
1000  if ( drawdown > 0.0) {
1001  // Do nothing is crossflow is not allowed
1002  if (!allow_cf && well_type_ == INJECTOR) {
1003  return;
1004  }
1005 
1006  // compute component volumetric rates at standard conditions
1007  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1008  const EvalWell cq_p = - well_index_[perf] * (mob_perfcells[comp_idx] * drawdown);
1009  cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
1010  }
1011 
1012  if (active()[Oil] && active()[Gas]) {
1013  const int oilpos = pu.phase_pos[Oil];
1014  const int gaspos = pu.phase_pos[Gas];
1015  const EvalWell cq_s_oil = cq_s[oilpos];
1016  const EvalWell cq_s_gas = cq_s[gaspos];
1017  cq_s[gaspos] += rs * cq_s_oil;
1018  cq_s[oilpos] += rv * cq_s_gas;
1019  }
1020  } else { // injecting perforations
1021  // Do nothing if crossflow is not allowed
1022  if (!allow_cf && well_type_ == PRODUCER) {
1023  return;
1024  }
1025 
1026  // for injecting perforations, we use total mobility
1027  EvalWell total_mob = mob_perfcells[0];
1028  for (int comp_idx = 1; comp_idx < num_comp; ++comp_idx) {
1029  total_mob += mob_perfcells[comp_idx];
1030  }
1031 
1032  // injection perforations total volume rates
1033  const EvalWell cqt_i = - well_index_[perf] * (total_mob * drawdown);
1034 
1035  // compute volume ratio between connection and at standard conditions
1036  EvalWell volume_ratio = 0.0;
1037  if (active()[Water]) {
1038  const int watpos = pu.phase_pos[Water];
1039  volume_ratio += cmix_s[watpos] / b_perfcells[watpos];
1040  }
1041 
1042  if (active()[Oil] && active()[Gas]) {
1043  const int oilpos = pu.phase_pos[Oil];
1044  const int gaspos = pu.phase_pos[Gas];
1045 
1046  // Incorporate RS/RV factors if both oil and gas active
1047  // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations
1048  // basically, for injecting perforations, the wellbore is the upstreaming side.
1049  const EvalWell d = 1.0 - rv * rs;
1050 
1051  if (d.value() == 0.0) {
1052  OPM_THROW(Opm::NumericalProblem, "Zero d value obtained for well " << name() << " during flux calcuation"
1053  << " with rs " << rs << " and rv " << rv);
1054  }
1055 
1056  const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d;
1057  volume_ratio += tmp_oil / b_perfcells[oilpos];
1058 
1059  const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d;
1060  volume_ratio += tmp_gas / b_perfcells[gaspos];
1061  } else { // not having gas and oil at the same time
1062  if (active()[Oil]) {
1063  const int oilpos = pu.phase_pos[Oil];
1064  volume_ratio += cmix_s[oilpos] / b_perfcells[oilpos];
1065  }
1066  if (active()[Gas]) {
1067  const int gaspos = pu.phase_pos[Gas];
1068  volume_ratio += cmix_s[gaspos] / b_perfcells[gaspos];
1069  }
1070  }
1071  // injecting connections total volumerates at standard conditions
1072  EvalWell cqt_is = cqt_i / volume_ratio;
1073  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1074  cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
1075  }
1076  } // end for injection perforations
1077  }
1078 
1079 
1080 
1081 
1082 
1083  template <typename TypeTag>
1084  typename MultisegmentWell<TypeTag>::EvalWell
1085  MultisegmentWell<TypeTag>::
1086  extendEval(const Eval& in) const
1087  {
1088  EvalWell out = 0.0;
1089  out.setValue(in.value());
1090  for(int eq_idx = 0; eq_idx < numEq;++eq_idx) {
1091  out.setDerivative(eq_idx, in.derivative(eq_idx));
1092  }
1093  return out;
1094  }
1095 
1096 
1097 
1098 
1099 
1100  template <typename TypeTag>
1101  void
1102  MultisegmentWell<TypeTag>::
1103  computeSegmentFluidProperties(const Simulator& ebosSimulator)
1104  {
1105  // TODO: the concept of phases and components are rather confusing in this function.
1106  // needs to be addressed sooner or later.
1107 
1108  // get the temperature for later use. It is only useful when we are not handling
1109  // thermal related simulation
1110  // basically, it is a single value for all the segments
1111 
1112  EvalWell temperature;
1113  // not sure how to handle the pvt region related to segment
1114  // for the current approach, we use the pvt region of the first perforated cell
1115  // although there are some text indicating using the pvt region of the lowest
1116  // perforated cell
1117  // TODO: later to investigate how to handle the pvt region
1118  int pvt_region_index;
1119  {
1120  // using the first perforated cell
1121  const int cell_idx = well_cells_[0];
1122  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1123  const auto& fs = intQuants.fluidState();
1124  temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1125  pvt_region_index = fs.pvtRegionIndex();
1126  }
1127 
1128  std::vector<double> surf_dens(number_of_phases_);
1129  // Surface density.
1130  // not using num_comp here is because solvent can be component
1131  for (int phase = 0; phase < number_of_phases_; ++phase) {
1132  surf_dens[phase] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx(phase), pvt_region_index );
1133  }
1134 
1135  const int num_comp = numComponents();
1136  const Opm::PhaseUsage& pu = phaseUsage();
1137  for (int seg = 0; seg < numberOfSegments(); ++seg) {
1138  // the compostion of the components inside wellbore under surface condition
1139  std::vector<EvalWell> mix_s(num_comp, 0.0);
1140  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1141  mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
1142  }
1143 
1144  std::vector<EvalWell> b(num_comp, 0.0);
1145  // it is the phase viscosities asked for
1146  std::vector<EvalWell> visc(number_of_phases_, 0.0);
1147  const EvalWell seg_pressure = getSegmentPressure(seg);
1148  if (pu.phase_used[BlackoilPhases::Aqua]) {
1149  // TODO: what is the difference between Water and BlackoilPhases::Aqua?
1150  const int water_pos = pu.phase_pos[BlackoilPhases::Aqua];
1151  b[water_pos] =
1152  FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1153  visc[water_pos] =
1154  FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure);
1155  }
1156 
1157  EvalWell rv(0.0);
1158  // gas phase
1159  if (pu.phase_used[BlackoilPhases::Vapour]) {
1160  const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
1161  if (pu.phase_used[BlackoilPhases::Liquid]) {
1162  const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
1163  const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure);
1164  if (mix_s[oilpos] > 0.0) {
1165  if (mix_s[gaspos] > 0.0) {
1166  rv = mix_s[oilpos] / mix_s[gaspos];
1167  }
1168 
1169  if (rv > rvmax) {
1170  rv = rvmax;
1171  }
1172  b[gaspos] =
1173  FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv);
1174  visc[gaspos] =
1175  FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv);
1176  } else { // no oil exists
1177  b[gaspos] =
1178  FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1179  visc[gaspos] =
1180  FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
1181  }
1182  } else { // no Liquid phase
1183  // it is the same with zero mix_s[Oil]
1184  b[gaspos] =
1185  FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1186  visc[gaspos] =
1187  FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
1188  }
1189  }
1190 
1191  EvalWell rs(0.0);
1192  // oil phase
1193  if (pu.phase_used[BlackoilPhases::Liquid]) {
1194  const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
1195  if (pu.phase_used[BlackoilPhases::Liquid]) {
1196  const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
1197  const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure);
1198  if (mix_s[gaspos] > 0.0) {
1199  if (mix_s[oilpos] > 0.0) {
1200  rs = mix_s[gaspos] / mix_s[oilpos];
1201  }
1202 
1203  if (rs > rsmax) {
1204  rs = rsmax;
1205  }
1206  b[oilpos] =
1207  FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs);
1208  visc[oilpos] =
1209  FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs);
1210  } else { // no oil exists
1211  b[oilpos] =
1212  FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1213  visc[oilpos] =
1214  FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
1215  }
1216  } else { // no Liquid phase
1217  // it is the same with zero mix_s[Oil]
1218  b[oilpos] =
1219  FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1220  visc[oilpos] =
1221  FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
1222  }
1223  }
1224 
1225  std::vector<EvalWell> mix(mix_s);
1226  if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
1227  const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
1228  const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
1229  if (rs != 0.0) { // rs > 0.0?
1230  mix[gaspos] = (mix_s[gaspos] - mix_s[oilpos] * rs) / (1. - rs * rv);
1231  }
1232  if (rv != 0.0) { // rv > 0.0?
1233  mix[oilpos] = (mix_s[oilpos] - mix_s[gaspos] * rv) / (1. - rs * rv);
1234  }
1235  }
1236 
1237  EvalWell volrat(0.0);
1238  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1239  volrat += mix[comp_idx] / b[comp_idx];
1240  }
1241 
1242  segment_viscosities_[seg] = 0.;
1243  // calculate the average viscosity
1244  for (int p = 0; p < number_of_phases_; ++p) {
1245  const EvalWell phase_fraction = mix[p] / b[p] / volrat;
1246  segment_viscosities_[seg] += visc[p] * phase_fraction;
1247  }
1248 
1249  EvalWell density(0.0);
1250  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1251  density += surf_dens[comp_idx] * mix_s[comp_idx];
1252  }
1253  segment_densities_[seg] = density / volrat;
1254 
1255  // calculate the mass rates
1256  segment_mass_rates_[seg] = 0.;
1257  for (int phase = 0; phase < number_of_phases_; ++phase) {
1258  const EvalWell rate = getSegmentRate(seg, phase);
1259  segment_mass_rates_[seg] += rate * surf_dens[phase];
1260  }
1261  }
1262  }
1263 
1264 
1265 
1266 
1267 
1268  template <typename TypeTag>
1269  typename MultisegmentWell<TypeTag>::EvalWell
1270  MultisegmentWell<TypeTag>::
1271  getSegmentPressure(const int seg) const
1272  {
1273  return primary_variables_evaluation_[seg][SPres];
1274  }
1275 
1276 
1277 
1278 
1279 
1280  template <typename TypeTag>
1281  typename MultisegmentWell<TypeTag>::EvalWell
1282  MultisegmentWell<TypeTag>::
1283  getSegmentRate(const int seg,
1284  const int comp_idx) const
1285  {
1286  return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx);
1287  }
1288 
1289 
1290 
1291 
1292 
1293  template <typename TypeTag>
1294  typename MultisegmentWell<TypeTag>::EvalWell
1295  MultisegmentWell<TypeTag>::
1296  getSegmentGTotal(const int seg) const
1297  {
1298  return primary_variables_evaluation_[seg][GTotal];
1299  }
1300 
1301 
1302 
1303 
1304 
1305  template <typename TypeTag>
1306  void
1307  MultisegmentWell<TypeTag>::
1308  getMobility(const Simulator& ebosSimulator,
1309  const int perf,
1310  std::vector<EvalWell>& mob) const
1311  {
1312  // TODO: most of this function, if not the whole function, can be moved to the base class
1313  const int np = number_of_phases_;
1314  const int cell_idx = well_cells_[perf];
1315  assert (int(mob.size()) == numComponents());
1316  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1317  const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1318 
1319  // either use mobility of the perforation cell or calcualte its own
1320  // based on passing the saturation table index
1321  const int satid = saturation_table_number_[perf] - 1;
1322  const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1323  if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
1324 
1325  for (int phase = 0; phase < np; ++phase) {
1326  int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
1327  mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx));
1328  }
1329  // if (has_solvent) {
1330  // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1331  // }
1332  } else {
1333 
1334  const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1335  Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
1336  MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1337 
1338  // reset the satnumvalue back to original
1339  materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1340 
1341  // compute the mobility
1342  for (int phase = 0; phase < np; ++phase) {
1343  int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
1344  mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx));
1345  }
1346 
1347  }
1348  }
1349 
1350 
1351 
1352 
1353 
1354  template <typename TypeTag>
1355  void
1356  MultisegmentWell<TypeTag>::
1357  assembleControlEq() const
1358  {
1359  EvalWell control_eq(0.0);
1360 
1361  switch (well_controls_get_current_type(well_controls_)) {
1362  case THP: // not handling this one for now
1363  {
1364  OPM_THROW(std::runtime_error, "Not handling THP control for Multisegment wells for now");
1365  }
1366  case BHP:
1367  {
1368  const double target_bhp = well_controls_get_current_target(well_controls_);
1369  control_eq = getSegmentPressure(0) - target_bhp;
1370  break;
1371  }
1372  case SURFACE_RATE:
1373  {
1374  // finding if it is a single phase control or combined phase control
1375  int number_phases_under_control = 0;
1376  const double* distr = well_controls_get_current_distr(well_controls_);
1377  for (int phase = 0; phase < number_of_phases_; ++phase) {
1378  if (distr[phase] > 0.0) {
1379  ++number_phases_under_control;
1380  }
1381  }
1382  assert(number_phases_under_control > 0);
1383 
1384  const std::vector<double> g = {1.0, 1.0, 0.01};
1385  const double target_rate = well_controls_get_current_target(well_controls_);
1386  // TODO: the two situations below should be able to be merged to be handled as one situation
1387 
1388  if (number_phases_under_control == 1) { // single phase control
1389  for (int phase = 0; phase < number_of_phases_; ++phase) {
1390  if (distr[phase] > 0.) { // under the control of this phase
1391  control_eq = getSegmentGTotal(0) * volumeFraction(0, phase) - g[phase] * target_rate;
1392  break;
1393  }
1394  }
1395  } else { // multiphase rate control
1396  EvalWell rate_for_control(0.0);
1397  const EvalWell G_total = getSegmentGTotal(0);
1398  for (int phase = 0; phase < number_of_phases_; ++phase) {
1399  if (distr[phase] > 0.) {
1400  rate_for_control += G_total * volumeFractionScaled(0, phase);
1401  }
1402  }
1403  // TODO: maybe the following equation can be scaled a little bit for gas phase
1404  control_eq = rate_for_control - target_rate;
1405  }
1406  break;
1407  }
1408  case RESERVOIR_RATE:
1409  {
1410  EvalWell rate_for_control(0.0); // reservoir rate
1411  const double* distr = well_controls_get_current_distr(well_controls_);
1412  for (int phase = 0; phase < number_of_phases_; ++phase) {
1413  if (distr[phase] > 0.) {
1414  rate_for_control += getSegmentGTotal(0) * volumeFraction(0, phase);
1415  }
1416  }
1417  const double target_rate = well_controls_get_current_target(well_controls_);
1418  control_eq = rate_for_control - target_rate;
1419  break;
1420  }
1421  default:
1422  OPM_THROW(std::runtime_error, "Unknown well control control types for well " << name());
1423  }
1424 
1425 
1426  // using control_eq to update the matrix and residuals
1427 
1428  resWell_[0][SPres] = control_eq.value();
1429  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1430  duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + numEq);
1431  }
1432  }
1433 
1434 
1435 
1436 
1437 
1438  template <typename TypeTag>
1439  void
1440  MultisegmentWell<TypeTag>::
1441  assemblePressureEq(const int seg) const
1442  {
1443  assert(seg != 0); // not top segment
1444 
1445  // for top segment, the well control equation will be used.
1446  EvalWell pressure_equation = getSegmentPressure(seg);
1447 
1448  // we need to handle the pressure difference between the two segments
1449  // we only consider the hydrostatic pressure loss first
1450  pressure_equation -= getHydroPressureLoss(seg);
1451 
1452  if (frictionalPressureLossConsidered()) {
1453  pressure_equation -= getFrictionPressureLoss(seg);
1454  }
1455 
1456  resWell_[seg][SPres] = pressure_equation.value();
1457  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1458  duneD_[seg][seg][SPres][pv_idx] = pressure_equation.derivative(pv_idx + numEq);
1459  }
1460 
1461  // contribution from the outlet segment
1462  const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
1463  const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
1464 
1465  resWell_[seg][SPres] -= outlet_pressure.value();
1466  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1467  duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
1468  }
1469 
1470  if (accelerationalPressureLossConsidered()) {
1471  handleAccelerationPressureLoss(seg);
1472  }
1473  }
1474 
1475 
1476 
1477 
1478 
1479  template <typename TypeTag>
1480  typename MultisegmentWell<TypeTag>::EvalWell
1481  MultisegmentWell<TypeTag>::
1482  getHydroPressureLoss(const int seg) const
1483  {
1484  return segment_densities_[seg] * gravity_ * segment_depth_diffs_[seg];
1485  }
1486 
1487 
1488 
1489 
1490 
1491  template <typename TypeTag>
1492  typename MultisegmentWell<TypeTag>::EvalWell
1493  MultisegmentWell<TypeTag>::
1494  getFrictionPressureLoss(const int seg) const
1495  {
1496  const EvalWell mass_rate = segment_mass_rates_[seg];
1497  const EvalWell density = segment_densities_[seg];
1498  const EvalWell visc = segment_viscosities_[seg];
1499  const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
1500  const double length = segmentSet()[seg].totalLength() - segmentSet()[outlet_segment_index].totalLength();
1501  assert(length > 0.);
1502  const double roughness = segmentSet()[seg].roughness();
1503  const double area = segmentSet()[seg].crossArea();
1504  const double diameter = segmentSet()[seg].internalDiameter();
1505 
1506  const double sign = mass_rate < 0. ? 1.0 : - 1.0;
1507 
1508  return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
1509  }
1510 
1511 
1512 
1513 
1514 
1515  template <typename TypeTag>
1516  void
1517  MultisegmentWell<TypeTag>::
1518  handleAccelerationPressureLoss(const int seg) const
1519  {
1520  // TODO: this pressure loss is not significant enough to be well tested yet.
1521  // handle the out velcocity head
1522  const double area = segmentSet()[seg].crossArea();
1523  const EvalWell mass_rate = segment_mass_rates_[seg];
1524  const EvalWell density = segment_densities_[seg];
1525  const EvalWell out_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density);
1526 
1527  resWell_[seg][SPres] -= out_velocity_head.value();
1528  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1529  duneD_[seg][seg][SPres][pv_idx] -= out_velocity_head.derivative(pv_idx + numEq);
1530  }
1531 
1532  // calcuate the maximum cross-area among the segment and its inlet segments
1533  double max_area = area;
1534  for (const int inlet : segment_inlets_[seg]) {
1535  const double inlet_area = segmentSet()[inlet].crossArea();
1536  if (inlet_area > max_area) {
1537  max_area = inlet_area;
1538  }
1539  }
1540 
1541  // handling the velocity head of intlet segments
1542  for (const int inlet : segment_inlets_[seg]) {
1543  const EvalWell density = segment_densities_[inlet];
1544  const EvalWell mass_rate = segment_mass_rates_[inlet];
1545  const EvalWell inlet_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density);
1546  resWell_[seg][SPres] += inlet_velocity_head.value();
1547  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1548  duneD_[seg][inlet][SPres][pv_idx] += inlet_velocity_head.derivative(pv_idx + numEq);
1549  }
1550  }
1551  }
1552 
1553 
1554 
1555 
1556 
1557  template <typename TypeTag>
1558  void
1559  MultisegmentWell<TypeTag>::
1560  processFractions(const int seg) const
1561  {
1562  const PhaseUsage& pu = phaseUsage();
1563 
1564  std::vector<double> fractions(number_of_phases_, 0.0);
1565 
1566  assert( active()[Oil] );
1567  const int oil_pos = pu.phase_pos[Oil];
1568  fractions[oil_pos] = 1.0;
1569 
1570  if ( active()[Water] ) {
1571  const int water_pos = pu.phase_pos[Water];
1572  fractions[water_pos] = primary_variables_[seg][WFrac];
1573  fractions[oil_pos] -= fractions[water_pos];
1574  }
1575 
1576  if ( active()[Gas] ) {
1577  const int gas_pos = pu.phase_pos[Gas];
1578  fractions[gas_pos] = primary_variables_[seg][GFrac];
1579  fractions[oil_pos] -= fractions[gas_pos];
1580  }
1581 
1582  if (active()[Water]) {
1583  const int water_pos = pu.phase_pos[Water];
1584  if (fractions[water_pos] < 0.0) {
1585  if ( active()[Gas] ) {
1586  fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]);
1587  }
1588  fractions[oil_pos] /= (1.0 - fractions[water_pos]);
1589  fractions[water_pos] = 0.0;
1590  }
1591  }
1592 
1593  if (active()[Gas]) {
1594  const int gas_pos = pu.phase_pos[Gas];
1595  if (fractions[gas_pos] < 0.0) {
1596  if ( active()[Water] ) {
1597  fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]);
1598  }
1599  fractions[oil_pos] /= (1.0 - fractions[gas_pos]);
1600  fractions[gas_pos] = 0.0;
1601  }
1602  }
1603 
1604  if (fractions[oil_pos] < 0.0) {
1605  if ( active()[Water] ) {
1606  fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]);
1607  }
1608  if ( active()[Gas] ) {
1609  fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]);
1610  }
1611  fractions[oil_pos] = 0.0;
1612  }
1613 
1614  if ( active()[Water] ) {
1615  primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]];
1616  }
1617 
1618  if ( active()[Gas] ) {
1619  primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
1620  }
1621  }
1622 
1623 
1624 
1625 
1626 
1627  template <typename TypeTag>
1628  void
1629  MultisegmentWell<TypeTag>::
1630  updateWellStateFromPrimaryVariables(WellState& well_state) const
1631  {
1632  const PhaseUsage& pu = phaseUsage();
1633  assert( active()[Oil] );
1634  const int oil_pos = pu.phase_pos[Oil];
1635 
1636  for (int seg = 0; seg < numberOfSegments(); ++seg) {
1637  std::vector<double> fractions(number_of_phases_, 0.0);
1638  fractions[oil_pos] = 1.0;
1639 
1640  if ( active()[Water] ) {
1641  const int water_pos = pu.phase_pos[Water];
1642  fractions[water_pos] = primary_variables_[seg][WFrac];
1643  fractions[oil_pos] -= fractions[water_pos];
1644  }
1645 
1646  if ( active()[Gas] ) {
1647  const int gas_pos = pu.phase_pos[Gas];
1648  fractions[gas_pos] = primary_variables_[seg][GFrac];
1649  fractions[oil_pos] -= fractions[gas_pos];
1650  }
1651 
1652  // convert the fractions to be Q_p / G_total to calculate the phase rates
1653  for (int p = 0; p < number_of_phases_; ++p) {
1654  const double scale = scalingFactor(p);
1655  // for injection wells, there should only one non-zero scaling factor
1656  if (scale > 0.) {
1657  fractions[p] /= scale;
1658  } else {
1659  // this should only happens to injection wells
1660  fractions[p] = 0.;
1661  }
1662  }
1663 
1664  // calculate the phase rates based on the primary variables
1665  const double g_total = primary_variables_[seg][GTotal];
1666  const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
1667  for (int p = 0; p < number_of_phases_; ++p) {
1668  const double phase_rate = g_total * fractions[p];
1669  well_state.segRates()[(seg + top_segment_index) * number_of_phases_ + p] = phase_rate;
1670  if (seg == 0) { // top segment
1671  well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = phase_rate;
1672  }
1673  }
1674 
1675  // update the segment pressure
1676  well_state.segPress()[seg + top_segment_index] = primary_variables_[seg][SPres];
1677  if (seg == 0) { // top segment
1678  well_state.bhp()[index_of_well_] = well_state.segPress()[seg + top_segment_index];
1679  }
1680  }
1681  }
1682 
1683 
1684 
1685 
1686 
1687  template <typename TypeTag>
1688  double
1689  MultisegmentWell<TypeTag>::
1690  scalingFactor(const int comp_idx) const
1691  {
1692  const double* distr = well_controls_get_current_distr(well_controls_);
1693 
1694  if (well_controls_get_current_type(well_controls_) == RESERVOIR_RATE) {
1695  // if (has_solvent && phaseIdx == contiSolventEqIdx )
1696  // OPM_THROW(std::runtime_error, "RESERVOIR_RATE control in combination with solvent is not implemented");
1697  return distr[comp_idx];
1698  }
1699 
1700  const PhaseUsage& pu = phaseUsage();
1701 
1702  if (active()[Water] && pu.phase_pos[Water] == comp_idx)
1703  return 1.0;
1704  if (active()[Oil] && pu.phase_pos[Oil] == comp_idx)
1705  return 1.0;
1706  if (active()[Gas] && pu.phase_pos[Gas] == comp_idx)
1707  return 0.01;
1708  // if (has_solvent && phaseIdx == contiSolventEqIdx )
1709  // return 0.01;
1710 
1711  // we should not come this far
1712  assert(false);
1713  return 1.0;
1714  }
1715 
1716 
1717 
1718 
1719 
1720  template <typename TypeTag>
1721  bool
1722  MultisegmentWell<TypeTag>::
1723  frictionalPressureLossConsidered() const
1724  {
1725  // HF- and HFA needs to consider frictional pressure loss
1726  return (segmentSet().compPressureDrop() != WellSegment::H__);
1727  }
1728 
1729 
1730 
1731 
1732 
1733  template <typename TypeTag>
1734  bool
1735  MultisegmentWell<TypeTag>::
1736  accelerationalPressureLossConsidered() const
1737  {
1738  return (segmentSet().compPressureDrop() == WellSegment::HFA);
1739  }
1740 
1741 
1742 
1743 
1744 
1745  template<typename TypeTag>
1746  void
1747  MultisegmentWell<TypeTag>::
1748  iterateWellEquations(Simulator& ebosSimulator,
1749  const double dt,
1750  WellState& well_state)
1751  {
1752  // basically, it only iterate through the equations.
1753  // we update the primary variables
1754  // if converged, we can update the well_state.
1755  // the function updateWellState() should have a flag to show
1756  // if we will update the well state.
1757  const int max_iter_number = param_.max_inner_iter_ms_wells_;
1758  int it = 0;
1759  for (; it < max_iter_number; ++it) {
1760 
1761  assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, true);
1762 
1763 #if HAVE_UMFPACK
1764  const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
1765 #else
1766  const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_);
1767 #endif
1768 
1769  // TODO: use these small values for now, not intend to reach the convergence
1770  // in this stage, but, should we?
1771  // We should try to avoid hard-code values in the code.
1772  // If we want to use the real one, we need to find a way to get them.
1773  // const std::vector<double> B {0.8, 0.8, 0.008};
1774  const std::vector<double> B {0.5, 0.5, 0.005};
1775 
1776  const ConvergenceReport report = getWellConvergence(B);
1777 
1778  if (report.converged) {
1779  break;
1780  }
1781 
1782  updateWellState(dx_well, true, well_state);
1783 
1784  initPrimaryVariablesEvaluation();
1785  }
1786  // TODO: maybe we should not use these values if they are not converged.
1787  }
1788 
1789 
1790 
1791 
1792 
1793  template<typename TypeTag>
1794  void
1795  MultisegmentWell<TypeTag>::
1796  assembleWellEqWithoutIteration(Simulator& ebosSimulator,
1797  const double dt,
1798  WellState& well_state,
1799  bool only_wells)
1800  {
1801  // calculate the fluid properties needed.
1802  computeSegmentFluidProperties(ebosSimulator);
1803 
1804  // clear all entries
1805  if (!only_wells) {
1806  duneB_ = 0.0;
1807  duneC_ = 0.0;
1808  }
1809 
1810  duneD_ = 0.0;
1811  resWell_ = 0.0;
1812 
1813  // for the black oil cases, there will be four equations,
1814  // the first three of them are the mass balance equations, the last one is the pressure equations.
1815  //
1816  // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
1817 
1818  auto& ebosJac = ebosSimulator.model().linearizer().matrix();
1819  auto& ebosResid = ebosSimulator.model().linearizer().residual();
1820 
1821  const bool allow_cf = getAllowCrossFlow();
1822 
1823  const int nseg = numberOfSegments();
1824  const int num_comp = numComponents();
1825 
1826  for (int seg = 0; seg < nseg; ++seg) {
1827  // calculating the accumulation term // TODO: without considering the efficiencty factor for now
1828  // volume of the segment
1829  {
1830  const double volume = segmentSet()[seg].volume();
1831  // for each component
1832  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1833  EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx])
1834  + getSegmentRate(seg, comp_idx);
1835 
1836  resWell_[seg][comp_idx] += accumulation_term.value();
1837  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1838  duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq);
1839  }
1840  }
1841  }
1842 
1843  // considering the contributions from the inlet segments
1844  {
1845  for (const int inlet : segment_inlets_[seg]) {
1846  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1847  const EvalWell inlet_rate = getSegmentRate(inlet, comp_idx);
1848  resWell_[seg][comp_idx] -= inlet_rate.value();
1849  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1850  duneD_[seg][inlet][comp_idx][pv_idx] -= inlet_rate.derivative(pv_idx + numEq);
1851  }
1852  }
1853  }
1854  }
1855 
1856  // calculating the perforation rate for each perforation that belongs to this segment
1857  const EvalWell seg_pressure = getSegmentPressure(seg);
1858  for (const int perf : segment_perforations_[seg]) {
1859  const int cell_idx = well_cells_[perf];
1860  const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1861  std::vector<EvalWell> mob(num_comp, 0.0);
1862  getMobility(ebosSimulator, perf, mob);
1863  std::vector<EvalWell> cq_s(num_comp, 0.0);
1864  computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s);
1865 
1866  for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1867  // the cq_s entering mass balance equations need to consider the efficiency factors.
1868  const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_;
1869 
1870  if (!only_wells) {
1871  // subtract sum of component fluxes in the reservoir equation.
1872  // need to consider the efficiency factor
1873  // TODO: the name of the function flowPhaseToEbosCompIdx is prolematic, since the input
1874  // is a component index :D
1875  ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value();
1876  }
1877 
1878  // subtract sum of phase fluxes in the well equations.
1879  resWell_[seg][comp_idx] -= cq_s_effective.value();
1880 
1881  // assemble the jacobians
1882  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1883  if (!only_wells) {
1884  // also need to consider the efficiency factor when manipulating the jacobians.
1885  duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix
1886  }
1887  // the index name for the D should be eq_idx / pv_idx
1888  duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq);
1889  }
1890 
1891  for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) {
1892  if (!only_wells) {
1893  // also need to consider the efficiency factor when manipulating the jacobians.
1894  ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(comp_idx)][pv_idx] -= cq_s_effective.derivative(pv_idx);
1895  duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
1896  }
1897  }
1898  }
1899  // TODO: we should save the perforation pressure and preforation rates?
1900  // we do not use it in the simulation for now, while we might need them if
1901  // we handle the pressure in SEG mode.
1902  }
1903 
1904  // the fourth dequation, the pressure drop equation
1905  if (seg == 0) { // top segment, pressure equation is the control equation
1906  assembleControlEq();
1907  } else {
1908  assemblePressureEq(seg);
1909  }
1910  }
1911  }
1912 
1913 }
virtual void recoverWellSolutionAndUpdateWellState(const BVector &x, WellState &well_state) const
using the solution x to recover the solution xw for wells and applying xw to update Well State ...
Definition: MultisegmentWell_impl.hpp:533
virtual ConvergenceReport getWellConvergence(const std::vector< double > &B_avg) const
check whether the well equations get converged for this well
Definition: MultisegmentWell_impl.hpp:409
Definition: AutoDiff.hpp:297
a struct to collect information about the convergence checking
Definition: WellInterface.hpp:117
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
int numberOfSegments() const
number of segments for this well int number_of_segments_;
Definition: MultisegmentWell_impl.hpp:824
virtual void apply(const BVector &x, BVector &Ax) const
Ax = Ax - C D^-1 B x.
Definition: MultisegmentWell_impl.hpp:490
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
virtual void computeWellPotentials(const Simulator &ebosSimulator, const WellState &well_state, std::vector< double > &well_potentials)
computing the well potentials for group control
Definition: MultisegmentWell_impl.hpp:548
virtual void updateWellStateWithTarget(const int current, WellState &well_state) const
updating the well state based the control mode specified with current
Definition: MultisegmentWell_impl.hpp:251
Definition: MultisegmentWell.hpp:32
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:719