All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
BlackoilModelEbos.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2014, 2015 Statoil ASA.
5  Copyright 2015 NTNU
6  Copyright 2015 IRIS AS
7 
8  This file is part of the Open Porous Media project (OPM).
9 
10  OPM is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  OPM is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with OPM. If not, see <http://www.gnu.org/licenses/>.
22 */
23 
24 #ifndef OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
25 #define OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
26 
27 #include <ebos/eclproblem.hh>
28 #include <ewoms/common/start.hh>
29 
30 #include <opm/autodiff/BlackoilModelParameters.hpp>
31 #include <opm/autodiff/BlackoilWellModel.hpp>
32 #include <opm/autodiff/AutoDiffBlock.hpp>
33 #include <opm/autodiff/AutoDiffHelpers.hpp>
34 #include <opm/autodiff/GridHelpers.hpp>
35 #include <opm/autodiff/WellHelpers.hpp>
36 #include <opm/autodiff/GeoProps.hpp>
37 #include <opm/autodiff/WellDensitySegmented.hpp>
38 #include <opm/autodiff/VFPProperties.hpp>
39 #include <opm/autodiff/VFPProdProperties.hpp>
40 #include <opm/autodiff/VFPInjProperties.hpp>
41 #include <opm/autodiff/BlackoilDetails.hpp>
42 #include <opm/autodiff/BlackoilModelEnums.hpp>
43 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
45 
46 #include <opm/core/grid.h>
47 #include <opm/core/simulator/SimulatorReport.hpp>
48 #include <opm/core/linalg/LinearSolverInterface.hpp>
49 #include <opm/core/linalg/ParallelIstlInformation.hpp>
50 #include <opm/core/props/phaseUsageFromDeck.hpp>
51 #include <opm/common/ErrorMacros.hpp>
52 #include <opm/common/Exceptions.hpp>
53 #include <opm/common/OpmLog/OpmLog.hpp>
54 #include <opm/parser/eclipse/Units/Units.hpp>
55 #include <opm/core/well_controls.h>
56 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
57 #include <opm/core/utility/parameters/ParameterGroup.hpp>
58 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
59 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
60 
61 #include <opm/autodiff/ISTLSolver.hpp>
62 #include <opm/common/data/SimulationDataContainer.hpp>
63 
64 #include <dune/istl/owneroverlapcopy.hh>
65 #include <dune/common/parallel/collectivecommunication.hh>
66 #include <dune/common/timer.hh>
67 #include <dune/common/unused.hh>
68 
69 #include <cassert>
70 #include <cmath>
71 #include <iostream>
72 #include <iomanip>
73 #include <limits>
74 #include <vector>
75 #include <algorithm>
76 //#include <fstream>
77 
78 
79 
80 namespace Ewoms {
81 namespace Properties {
82 NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
83 SET_BOOL_PROP(EclFlowProblem, DisableWells, true);
84 SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false);
85 SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility, true);
86 // default in flow is to formulate the equations in surface volumes
87 SET_BOOL_PROP(EclFlowProblem, BlackoilConserveSurfaceVolume, true);
88 SET_BOOL_PROP(EclFlowProblem, UseVolumetricResidual, false);
89 
90 
91 // SWATINIT is done by the flow part of flow_ebos. this can be removed once the legacy
92 // code for fluid and satfunc handling gets fully retired.
93 SET_BOOL_PROP(EclFlowProblem, EnableSwatinit, false);
94 }}
95 
96 namespace Opm {
103  template <class TypeTag>
105  {
106  public:
107  // --------- Types and enums ---------
108  typedef BlackoilState ReservoirState;
111 
112  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
113  typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
114  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
115  typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
116  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables ;
117  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
118  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
119  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
120  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
121 
122  typedef double Scalar;
123  static const int numEq = Indices::numEq;
124  static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
125  static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
126  static const int solventSaturationIdx = Indices::solventSaturationIdx;
127  static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
128 
129  typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
130  typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
131  typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
132  typedef Dune::BlockVector<VectorBlockType> BVector;
133 
135  //typedef typename SolutionVector :: value_type PrimaryVariables ;
136 
137  // For the conversion between the surface volume rate and resrevoir voidage rate
140 
141  typedef Opm::FIPData FIPDataType;
142 
143  // --------- Public methods ---------
144 
155  BlackoilModelEbos(Simulator& ebosSimulator,
156  const ModelParameters& param,
157  BlackoilWellModel<TypeTag>& well_model,
158  RateConverterType& rate_converter,
159  const NewtonIterationBlackoilInterface& linsolver,
160  const bool terminal_output
161  )
162  : ebosSimulator_(ebosSimulator)
163  , grid_(ebosSimulator_.gridManager().grid())
164  , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
165  , phaseUsage_(phaseUsageFromDeck(eclState()))
166  , vfp_properties_(
167  eclState().getTableManager().getVFPInjTables(),
168  eclState().getTableManager().getVFPProdTables())
169  , active_(detail::activePhases(phaseUsage_))
170  , has_disgas_(FluidSystem::enableDissolvedGas())
171  , has_vapoil_(FluidSystem::enableVaporizedOil())
172  , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
173  , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
174  , param_( param )
175  , well_model_ (well_model)
176  , terminal_output_ (terminal_output)
177  , rate_converter_(rate_converter)
178  , current_relaxation_(1.0)
179  , dx_old_(AutoDiffGrid::numCells(grid_))
180  , isBeginReportStep_(false)
181  {
182  // Wells are active if they are active wells on at least
183  // one process.
184  int wellsActive = localWellsActive() ? 1 : 0;
185  wellsActive = grid_.comm().max(wellsActive);
186  wellModel().setWellsActive( wellsActive );
187  // compute global sum of number of cells
188  global_nc_ = detail::countGlobalCells(grid_);
189  wellModel().setVFPProperties(&vfp_properties_);
190  if (!istlSolver_)
191  {
192  OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed");
193  }
194  }
195 
196  bool isParallel() const
197  { return grid_.comm().size() > 1; }
198 
199  const EclipseState& eclState() const
200  { return ebosSimulator_.gridManager().eclState(); }
201 
207  const ReservoirState& /*reservoir_state*/,
208  const WellState& /* well_state */)
209  {
210  if ( wellModel().wellCollection()->havingVREPGroups() ) {
211  updateRateConverter();
212  }
213 
214  // update the solution variables in ebos
215 
216  // if the last time step failed we need to update the curent solution
217  // and recalculate the Intesive Quantities.
218  if ( timer.lastStepFailed() ) {
219  ebosSimulator_.model().solution( 0 /* timeIdx */ ) = ebosSimulator_.model().solution( 1 /* timeIdx */ );
220  ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
221  } else {
222  // set the initial solution.
223  ebosSimulator_.model().solution( 1 /* timeIdx */ ) = ebosSimulator_.model().solution( 0 /* timeIdx */ );
224  }
225 
226  unsigned numDof = ebosSimulator_.model().numGridDof();
227  wasSwitched_.resize(numDof);
228  std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
229  }
230 
231 
241  template <class NonlinearSolverType>
242  SimulatorReport nonlinearIteration(const int iteration,
243  const SimulatorTimerInterface& timer,
244  NonlinearSolverType& nonlinear_solver,
245  ReservoirState& /*reservoir_state*/,
246  WellState& well_state)
247  {
248  SimulatorReport report;
249  failureReport_ = SimulatorReport();
250  Dune::Timer perfTimer;
251 
252  perfTimer.start();
253  if (iteration == 0) {
254  // For each iteration we store in a vector the norms of the residual of
255  // the mass balance for each active phase, the well flux and the well equations.
256  residual_norms_history_.clear();
257  current_relaxation_ = 1.0;
258  dx_old_ = 0.0;
259  }
260 
261  report.total_linearizations = 1;
262 
263  try {
264  report += assemble(timer, iteration, well_state);
265  report.assemble_time += perfTimer.stop();
266  }
267  catch (...) {
268  report.assemble_time += perfTimer.stop();
269  failureReport_ += report;
270  // todo (?): make the report an attribute of the class
271  throw; // continue throwing the stick
272  }
273 
274  std::vector<double> residual_norms;
275  perfTimer.reset();
276  perfTimer.start();
277  // the step is not considered converged until at least minIter iterations is done
278  report.converged = getConvergence(timer, iteration,residual_norms) && iteration > nonlinear_solver.minIter();
279 
280  // checking whether the group targets are converged
281  if (wellModel().wellCollection()->groupControlActive()) {
282  report.converged = report.converged && wellModel().wellCollection()->groupTargetConverged(well_state.wellRates());
283  }
284 
285  report.update_time += perfTimer.stop();
286  residual_norms_history_.push_back(residual_norms);
287  if (!report.converged) {
288  perfTimer.reset();
289  perfTimer.start();
290  report.total_newton_iterations = 1;
291 
292  // enable single precision for solvers when dt is smaller then 20 days
293  //residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
294 
295  // Compute the nonlinear update.
296  const int nc = AutoDiffGrid::numCells(grid_);
297  const int nw = numWells();
298  BVector x(nc);
299 
300  try {
302  report.linear_solve_time += perfTimer.stop();
303  report.total_linear_iterations += linearIterationsLastSolve();
304  }
305  catch (...) {
306  report.linear_solve_time += perfTimer.stop();
307  report.total_linear_iterations += linearIterationsLastSolve();
308 
309  failureReport_ += report;
310  throw; // re-throw up
311  }
312 
313  perfTimer.reset();
314  perfTimer.start();
315 
316  // handling well state update before oscillation treatment is a decision based
317  // on observation to avoid some big performance degeneration under some circumstances.
318  // there is no theorectical explanation which way is better for sure.
319 
320  if( nw > 0 )
321  {
322  wellModel().recoverWellSolutionAndUpdateWellState(x, well_state);
323  }
324 
325  if (param_.use_update_stabilization_) {
326  // Stabilize the nonlinear update.
327  bool isOscillate = false;
328  bool isStagnate = false;
329  nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
330  if (isOscillate) {
331  current_relaxation_ -= nonlinear_solver.relaxIncrement();
332  current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
333  if (terminalOutputEnabled()) {
334  std::string msg = " Oscillating behavior detected: Relaxation set to "
335  + std::to_string(current_relaxation_);
336  OpmLog::info(msg);
337  }
338  }
339  nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
340  }
341 
342  // Apply the update, with considering model-dependent limitations and
343  // chopping of the update.
344  updateState(x);
345 
346  report.update_time += perfTimer.stop();
347  }
348 
349  return report;
350  }
351 
352  void printIf(int c, double x, double y, double eps, std::string type) {
353  if (std::abs(x-y) > eps) {
354  std::cout << type << " " <<c << ": "<<x << " " << y << std::endl;
355  }
356  }
357 
358 
365  const ReservoirState& reservoir_state,
366  WellState& well_state)
367  {
368  DUNE_UNUSED_PARAMETER(timer);
369  DUNE_UNUSED_PARAMETER(reservoir_state);
370  DUNE_UNUSED_PARAMETER(well_state);
371  }
372 
377  SimulatorReport assemble(const SimulatorTimerInterface& timer,
378  const int iterationIdx,
379  WellState& well_state)
380  {
381  using namespace Opm::AutoDiffGrid;
382 
383  SimulatorReport report;
384 
385  // when having VREP group control, update the rate converter based on reservoir state
386  if ( wellModel().wellCollection()->havingVREPGroups() ) {
387  updateRateConverter();
388  }
389 
390  // -------- Mass balance equations --------
391  assembleMassBalanceEq(timer, iterationIdx);
392 
393  // -------- Well equations ----------
394  double dt = timer.currentStepLength();
395 
396  try
397  {
398  // assembles the well equations and applies the wells to
399  // the reservoir equations as a source term.
400  report = wellModel().assemble(ebosSimulator_, iterationIdx, dt, well_state);
401  }
402  catch ( const Dune::FMatrixError& e )
403  {
404  OPM_THROW(Opm::NumericalProblem,"Error encounted when solving well equations");
405  }
406 
407  return report;
408  }
409 
410  // compute the "relative" change of the solution between time steps
411  template <class Dummy>
412  double relativeChange(const Dummy&, const Dummy&) const
413  {
414  Scalar resultDelta = 0.0;
415  Scalar resultDenom = 0.0;
416 
417  const auto& elemMapper = ebosSimulator_.model().elementMapper();
418  const auto& gridView = ebosSimulator_.gridView();
419  auto elemIt = gridView.template begin</*codim=*/0>();
420  const auto& elemEndIt = gridView.template end</*codim=*/0>();
421  for (; elemIt != elemEndIt; ++elemIt) {
422  const auto& elem = *elemIt;
423  if (elem.partitionType() != Dune::InteriorEntity)
424  continue;
425 
426  unsigned globalElemIdx = elemMapper.index(elem);
427  const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
428 
429  Scalar pressureNew;
430  pressureNew = priVarsNew[Indices::pressureSwitchIdx];
431 
432  Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
433  Scalar oilSaturationNew = 1.0;
434  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
435  saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx];
436  oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
437  }
438 
439  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
440  saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
441  oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
442  }
443 
444  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
445  saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
446  }
447 
448  const auto& priVarsOld = ebosSimulator_.model().solution(/*timeIdx=*/1)[globalElemIdx];
449 
450  Scalar pressureOld;
451  pressureOld = priVarsOld[Indices::pressureSwitchIdx];
452 
453  Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
454  Scalar oilSaturationOld = 1.0;
455  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
456  saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
457  oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
458  }
459 
460  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
461  saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
462  oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
463  }
464 
465  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
466  saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
467  }
468 
469  Scalar tmp = pressureNew - pressureOld;
470  resultDelta += tmp*tmp;
471  resultDenom += pressureNew*pressureNew;
472 
473  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
474  Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
475  resultDelta += tmp*tmp;
476  resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
477  }
478  }
479 
480  resultDelta = gridView.comm().sum(resultDelta);
481  resultDenom = gridView.comm().sum(resultDenom);
482 
483  if (resultDenom > 0.0)
484  return resultDelta/resultDenom;
485  return 0.0;
486  }
487 
488 
490  int sizeNonLinear() const
491  {
492  const int nc = Opm::AutoDiffGrid::numCells(grid_);
493  const int nw = numWells();
494  return numComponents() * (nc + nw);
495  }
496 
499  {
500  return istlSolver().iterations();
501  }
502 
505  void solveJacobianSystem(BVector& x) const
506  {
507  const auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
508  auto& ebosResid = ebosSimulator_.model().linearizer().residual();
509 
510  // J = [A, B; C, D], where A is the reservoir equations, B and C the interaction of well
511  // with the reservoir and D is the wells itself.
512  // The full system is reduced to a number of cells X number of cells system via Schur complement
513  // A -= B^T D^-1 C
514  // Instead of modifying A, the Ax operator is modified. i.e Ax -= B^T D^-1 C x in the WellModelMatrixAdapter.
515  // The residual is modified similarly.
516  // r = [r, r_well], where r is the residual and r_well the well residual.
517  // r -= B^T * D^-1 r_well
518 
519  // apply well residual to the residual.
520  wellModel().apply(ebosResid);
521 
522  // set initial guess
523  x = 0.0;
524 
525  // Solve system.
526  if( isParallel() )
527  {
529  Operator opA(ebosJac, well_model_, istlSolver().parallelInformation() );
530  assert( opA.comm() );
531  istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
532  }
533  else
534  {
536  Operator opA(ebosJac, well_model_);
537  istlSolver().solve( opA, x, ebosResid );
538  }
539  }
540 
541  //=====================================================================
542  // Implementation for ISTL-matrix based operator
543  //=====================================================================
544 
550  template<class M, class X, class Y, class WellModel, bool overlapping >
551  class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
552  {
553  typedef Dune::AssembledLinearOperator<M,X,Y> BaseType;
554 
555  public:
556  typedef M matrix_type;
557  typedef X domain_type;
558  typedef Y range_type;
559  typedef typename X::field_type field_type;
560 
561 #if HAVE_MPI
562  typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
563 #else
564  typedef Dune::CollectiveCommunication< Grid > communication_type;
565 #endif
566 
567  enum {
569  category = overlapping ?
570  Dune::SolverCategory::overlapping :
571  Dune::SolverCategory::sequential
572  };
573 
575  WellModelMatrixAdapter (const M& A, const WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
576  : A_( A ), wellMod_( wellMod ), comm_()
577  {
578 #if HAVE_MPI
579  if( parallelInformation.type() == typeid(ParallelISTLInformation) )
580  {
581  const ParallelISTLInformation& info =
582  boost::any_cast<const ParallelISTLInformation&>( parallelInformation);
583  comm_.reset( new communication_type( info.communicator() ) );
584  }
585 #endif
586  }
587 
588  virtual void apply( const X& x, Y& y ) const
589  {
590  A_.mv( x, y );
591  // add well model modification to y
592  wellMod_.apply(x, y );
593 
594 #if HAVE_MPI
595  if( comm_ )
596  comm_->project( y );
597 #endif
598  }
599 
600  // y += \alpha * A * x
601  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
602  {
603  A_.usmv(alpha,x,y);
604  // add scaled well model modification to y
605  wellMod_.applyScaleAdd( alpha, x, y );
606 
607 #if HAVE_MPI
608  if( comm_ )
609  comm_->project( y );
610 #endif
611  }
612 
613  virtual const matrix_type& getmat() const { return A_; }
614 
615  communication_type* comm()
616  {
617  return comm_.operator->();
618  }
619 
620  protected:
621  const matrix_type& A_ ;
622  const WellModel& wellMod_;
623  std::unique_ptr< communication_type > comm_;
624  };
625 
630  void updateState(const BVector& dx)
631  {
632  using namespace Opm::AutoDiffGrid;
633 
634  const auto& ebosProblem = ebosSimulator_.problem();
635 
636  unsigned numSwitched = 0;
637  ElementContext elemCtx( ebosSimulator_ );
638  const auto& gridView = ebosSimulator_.gridView();
639  const auto& elemEndIt = gridView.template end</*codim=*/0>();
640  SolutionVector& solution = ebosSimulator_.model().solution( 0 /* timeIdx */ );
641 
642  for (auto elemIt = gridView.template begin</*codim=*/0>();
643  elemIt != elemEndIt;
644  ++elemIt)
645  {
646  const auto& elem = *elemIt;
647 
648  elemCtx.updatePrimaryStencil(elem);
649  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
650  const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
651 
652  PrimaryVariables& priVars = solution[ cell_idx ];
653 
654  const double& dp = dx[cell_idx][Indices::pressureSwitchIdx];
655  double& p = priVars[Indices::pressureSwitchIdx];
656  const double& dp_rel_max = dpMaxRel();
657  const int sign_dp = dp > 0 ? 1: -1;
658  p -= sign_dp * std::min(std::abs(dp), std::abs(p)*dp_rel_max);
659  p = std::max(p, 0.0);
660 
661  // Saturation updates.
662  const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0;
663  const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0;
664 
665  double dso = 0.0;
666  double dsg = 0.0;
667  double drs = 0.0;
668  double drv = 0.0;
669 
670  // determine the saturation delta values
671  if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
672  dsg = dxvar;
673  }
674  else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
675  drs = dxvar;
676  }
677  else {
678  assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv);
679  drv = dxvar;
680  dsg = 0.0;
681  }
682 
683  // solvent
684  const double dss = has_solvent_ ? dx[cell_idx][Indices::solventSaturationIdx] : 0.0;
685 
686  // polymer
687  const double dc = has_polymer_ ? dx[cell_idx][Indices::polymerConcentrationIdx] : 0.0;
688 
689  // oil
690  dso = - (dsw + dsg + dss);
691 
692  // compute a scaling factor for the saturation update so that the maximum
693  // allowed change of saturations between iterations is not exceeded
694  double maxVal = 0.0;
695  maxVal = std::max(std::abs(dsw),maxVal);
696  maxVal = std::max(std::abs(dsg),maxVal);
697  maxVal = std::max(std::abs(dso),maxVal);
698  maxVal = std::max(std::abs(dss),maxVal);
699 
700  double satScaleFactor = 1.0;
701  if (maxVal > dsMax()) {
702  satScaleFactor = dsMax()/maxVal;
703  }
704 
705  if (active_[Water]) {
706  double& sw = priVars[Indices::waterSaturationIdx];
707  sw -= satScaleFactor * dsw;
708  }
709 
710  if (active_[Gas]) {
711  if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
712  double& sg = priVars[Indices::compositionSwitchIdx];
713  sg -= satScaleFactor * dsg;
714  }
715  }
716 
717  if (has_solvent_) {
718  double& ss = priVars[Indices::solventSaturationIdx];
719  ss -= satScaleFactor * dss;
720  ss = std::min(std::max(ss, 0.0),1.0);
721  }
722  if (has_polymer_) {
723  double& c = priVars[Indices::polymerConcentrationIdx];
724  c -= satScaleFactor * dc;
725  c = std::max(c, 0.0);
726  }
727 
728  // Update rs and rv
729  if (active_[Gas] && active_[Oil] ) {
730  unsigned pvtRegionIdx = ebosSimulator_.problem().pvtRegionIndex(cell_idx);
731  const double drmaxrel = drMaxRel();
732  if (has_disgas_) {
733  if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
734  Scalar RsSat =
735  FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx, 300.0, p);
736 
737  double& rs = priVars[Indices::compositionSwitchIdx];
738  rs -= ((drs<0)?-1:1)*std::min(std::abs(drs), RsSat*drmaxrel);
739  rs = std::max(rs, 0.0);
740  }
741 
742  }
743  if (has_vapoil_) {
744  if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
745  Scalar RvSat =
746  FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx, 300.0, p);
747 
748  double& rv = priVars[Indices::compositionSwitchIdx];
749  rv -= ((drv<0)?-1:1)*std::min(std::abs(drv), RvSat*drmaxrel);
750  rv = std::max(rv, 0.0);
751  }
752  }
753  }
754 
755  // Add an epsilon to make it harder to switch back immediately after the primary variable was changed.
756  if (wasSwitched_[cell_idx])
757  wasSwitched_[cell_idx] = priVars.adaptPrimaryVariables(ebosProblem, cell_idx, 1e-5);
758  else
759  wasSwitched_[cell_idx] = priVars.adaptPrimaryVariables(ebosProblem, cell_idx);
760 
761  if (wasSwitched_[cell_idx])
762  ++numSwitched;
763  }
764 
765  // if the solution is updated the intensive Quantities need to be recalculated
766  ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
767 
768  }
769 
772  {
773  return terminal_output_;
774  }
775 
776  template <class CollectiveCommunication>
777  double convergenceReduction(const CollectiveCommunication& comm,
778  const double pvSumLocal,
779  std::vector< Scalar >& R_sum,
780  std::vector< Scalar >& maxCoeff,
781  std::vector< Scalar >& B_avg)
782  {
783  // Compute total pore volume (use only owned entries)
784  double pvSum = pvSumLocal;
785 
786  if( comm.size() > 1 )
787  {
788  // global reduction
789  std::vector< Scalar > sumBuffer;
790  std::vector< Scalar > maxBuffer;
791  const int numComp = B_avg.size();
792  sumBuffer.reserve( 2*numComp + 1 ); // +1 for pvSum
793  maxBuffer.reserve( numComp );
794  for( int compIdx = 0; compIdx < numComp; ++compIdx )
795  {
796  sumBuffer.push_back( B_avg[ compIdx ] );
797  sumBuffer.push_back( R_sum[ compIdx ] );
798  maxBuffer.push_back( maxCoeff[ compIdx ] );
799  }
800 
801  // Compute total pore volume
802  sumBuffer.push_back( pvSum );
803 
804  // compute global sum
805  comm.sum( sumBuffer.data(), sumBuffer.size() );
806 
807  // compute global max
808  comm.max( maxBuffer.data(), maxBuffer.size() );
809 
810  // restore values to local variables
811  for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
812  {
813  B_avg[ compIdx ] = sumBuffer[ buffIdx ];
814  ++buffIdx;
815 
816  R_sum[ compIdx ] = sumBuffer[ buffIdx ];
817  }
818 
819  for( int compIdx = 0; compIdx < numComp; ++compIdx )
820  {
821  maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
822  }
823 
824  // restore global pore volume
825  pvSum = sumBuffer.back();
826  }
827 
828  // return global pore volume
829  return pvSum;
830  }
831 
837  bool getConvergence(const SimulatorTimerInterface& timer, const int iteration, std::vector<double>& residual_norms)
838  {
839  typedef std::vector< Scalar > Vector;
840 
841  const double dt = timer.currentStepLength();
842  const double tol_mb = param_.tolerance_mb_;
843  const double tol_cnv = param_.tolerance_cnv_;
844 
845  const int np = numPhases();
846  const int numComp = numComponents();
847 
848  Vector R_sum(numComp, 0.0 );
849  Vector B_avg(numComp, 0.0 );
850  Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
851 
852  const auto& ebosModel = ebosSimulator_.model();
853  const auto& ebosProblem = ebosSimulator_.problem();
854 
855  const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
856 
857  ElementContext elemCtx(ebosSimulator_);
858  const auto& gridView = ebosSimulator().gridView();
859  const auto& elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
860 
861  double pvSumLocal = 0.0;
862  for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
863  elemIt != elemEndIt;
864  ++elemIt)
865  {
866  const auto& elem = *elemIt;
867  elemCtx.updatePrimaryStencil(elem);
868  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
869  const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
870  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
871  const auto& fs = intQuants.fluidState();
872 
873  const double pvValue = ebosProblem.porosity(cell_idx) * ebosModel.dofTotalVolume( cell_idx );
874  pvSumLocal += pvValue;
875 
876  for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx )
877  {
878  const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx);
879  const int ebosCompIdx = flowPhaseToEbosCompIdx(phaseIdx);
880 
881  B_avg[ phaseIdx ] += 1.0 / fs.invB(ebosPhaseIdx).value();
882  const auto R2 = ebosResid[cell_idx][ebosCompIdx];
883 
884  R_sum[ phaseIdx ] += R2;
885  maxCoeff[ phaseIdx ] = std::max( maxCoeff[ phaseIdx ], std::abs( R2 ) / pvValue );
886  }
887 
888  if ( has_solvent_ ) {
889  B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
890  const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
891  R_sum[ contiSolventEqIdx ] += R2;
892  maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue );
893  }
894  if (has_polymer_ ) {
895  B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
896  const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
897  R_sum[ contiPolymerEqIdx ] += R2;
898  maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
899  }
900 
901  }
902 
903  // compute local average in terms of global number of elements
904  const int bSize = B_avg.size();
905  for ( int i = 0; i<bSize; ++i )
906  {
907  B_avg[ i ] /= Scalar( global_nc_ );
908  }
909 
910  // TODO: we remove the maxNormWell for now because the convergence of wells are on a individual well basis.
911  // Anyway, we need to provide some infromation to help debug the well iteration process.
912 
913 
914  // compute global sum and max of quantities
915  const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal,
916  R_sum, maxCoeff, B_avg);
917 
918  Vector CNV(numComp);
919  Vector mass_balance_residual(numComp);
920 
921  bool converged_MB = true;
922  bool converged_CNV = true;
923  // Finish computation
924  for ( int compIdx = 0; compIdx < numComp; ++compIdx )
925  {
926  CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
927  mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
928  converged_MB = converged_MB && (mass_balance_residual[compIdx] < tol_mb);
929  converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv);
930 
931  residual_norms.push_back(CNV[compIdx]);
932  }
933 
934  const bool converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg);
935 
936  bool converged = converged_MB && converged_Well;
937 
938  // do not care about the cell based residual in the last two Newton
939  // iterations
940  if (iteration < param_.max_strict_iter_)
941  converged = converged && converged_CNV;
942 
943  if ( terminal_output_ )
944  {
945  // Only rank 0 does print to std::cout
946  if (iteration == 0) {
947  std::string msg = "Iter";
948 
949  std::vector< std::string > key( numComp );
950  for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) {
951  const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
952  key[ phaseIdx ] = std::toupper( phaseName.front() );
953  }
954  if (has_solvent_) {
955  key[ solventSaturationIdx ] = "S";
956  }
957 
958  if (has_polymer_) {
959  key[ polymerConcentrationIdx ] = "P";
960  }
961 
962  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
963  msg += " MB(" + key[ compIdx ] + ") ";
964  }
965  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
966  msg += " CNV(" + key[ compIdx ] + ") ";
967  }
968  OpmLog::note(msg);
969  }
970  std::ostringstream ss;
971  const std::streamsize oprec = ss.precision(3);
972  const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
973  ss << std::setw(4) << iteration;
974  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
975  ss << std::setw(11) << mass_balance_residual[compIdx];
976  }
977  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
978  ss << std::setw(11) << CNV[compIdx];
979  }
980  ss.precision(oprec);
981  ss.flags(oflags);
982  OpmLog::note(ss.str());
983  }
984 
985  for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) {
986  const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
987 
988  if (std::isnan(mass_balance_residual[phaseIdx])
989  || std::isnan(CNV[phaseIdx])) {
990  OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName);
991  }
992  if (mass_balance_residual[phaseIdx] > maxResidualAllowed()
993  || CNV[phaseIdx] > maxResidualAllowed()) {
994  OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName);
995  }
996  }
997 
998  return converged;
999  }
1000 
1001 
1003  int numPhases() const
1004  {
1005  return phaseUsage_.num_phases;
1006  }
1007 
1008  int numComponents() const
1009  {
1010  if (numPhases() == 2) {
1011  return 2;
1012  }
1013  int numComp = FluidSystem::numComponents;
1014  if (has_solvent_)
1015  numComp ++;
1016 
1017  if (has_polymer_)
1018  numComp ++;
1019 
1020  return numComp;
1021  }
1022 
1024  template<class T>
1025  std::vector<std::vector<double> >
1026  computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
1027  {
1028  return computeFluidInPlace(fipnum);
1029  }
1030 
1031  std::vector<std::vector<double> >
1032  computeFluidInPlace(const std::vector<int>& fipnum) const
1033  {
1034  const auto& comm = grid_.comm();
1035  const auto& gridView = ebosSimulator().gridView();
1036  const int nc = gridView.size(/*codim=*/0);
1037  const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
1038  int ntFip = *std::max_element(fipnum.begin(), fipnum.end());
1039  ntFip = comm.max(ntFip);
1040 
1041  std::vector<double> tpv(ntFip, 0.0);
1042  std::vector<double> hcpv(ntFip, 0.0);
1043 
1044  std::vector<std::vector<double> > regionValues(ntFip, std::vector<double>(FIPDataType::fipValues,0.0));
1045 
1046  for (int i = 0; i<FIPDataType::fipValues; i++) {
1047  fip_.fip[i].resize(nc,0.0);
1048  }
1049 
1050  ElementContext elemCtx(ebosSimulator_);
1051  const auto& elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
1052  for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
1053  elemIt != elemEndIt;
1054  ++elemIt)
1055  {
1056  elemCtx.updatePrimaryStencil(*elemIt);
1057  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1058 
1059  const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
1060  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
1061  const auto& fs = intQuants.fluidState();
1062 
1063  const int regionIdx = fipnum[cellIdx] - 1;
1064  if (regionIdx < 0) {
1065  // the given cell is not attributed to any region
1066  continue;
1067  }
1068 
1069  // calculate the pore volume of the current cell. Note that the porosity
1070  // returned by the intensive quantities is defined as the ratio of pore
1071  // space to total cell volume and includes all pressure dependent (->
1072  // rock compressibility) and static modifiers (MULTPV, MULTREGP, NTG,
1073  // PORV, MINPV and friends). Also note that because of this, the porosity
1074  // returned by the intensive quantities can be outside of the physical
1075  // range [0, 1] in pathetic cases.
1076  const double pv =
1077  ebosSimulator_.model().dofTotalVolume(cellIdx)
1078  * intQuants.porosity().value();
1079 
1080  for (int phase = 0; phase < maxnp; ++phase) {
1081  const double b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value();
1082  const double s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value();
1083 
1084  fip_.fip[phase][cellIdx] = b * s * pv;
1085 
1086  if (active_[ phase ]) {
1087  regionValues[regionIdx][phase] += fip_.fip[phase][cellIdx];
1088  }
1089  }
1090 
1091  if (active_[ Oil ] && active_[ Gas ]) {
1092  // Account for gas dissolved in oil and vaporized oil
1093  fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][cellIdx] = fs.Rs().value() * fip_.fip[FIPDataType::FIP_LIQUID][cellIdx];
1094  fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][cellIdx] = fs.Rv().value() * fip_.fip[FIPDataType::FIP_VAPOUR][cellIdx];
1095 
1096  regionValues[regionIdx][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][cellIdx];
1097  regionValues[regionIdx][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][cellIdx];
1098  }
1099 
1100  const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
1101 
1102  tpv[regionIdx] += pv;
1103  hcpv[regionIdx] += pv * hydrocarbon;
1104  }
1105 
1106  // sum tpv (-> total pore volume of the regions) and hcpv (-> pore volume of the
1107  // the regions that is occupied by hydrocarbons)
1108  comm.sum(tpv.data(), tpv.size());
1109  comm.sum(hcpv.data(), hcpv.size());
1110 
1111  for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
1112  elemIt != elemEndIt;
1113  ++elemIt)
1114  {
1115  const auto& elem = *elemIt;
1116 
1117  elemCtx.updatePrimaryStencil(elem);
1118  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1119 
1120  unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
1121  const int regionIdx = fipnum[cellIdx] - 1;
1122  if (regionIdx < 0) {
1123  // the cell is not attributed to any region. ignore it!
1124  continue;
1125  }
1126 
1127  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
1128  const auto& fs = intQuants.fluidState();
1129 
1130  // calculate the pore volume of the current cell. Note that the
1131  // porosity returned by the intensive quantities is defined as the
1132  // ratio of pore space to total cell volume and includes all pressure
1133  // dependent (-> rock compressibility) and static modifiers (MULTPV,
1134  // MULTREGP, NTG, PORV, MINPV and friends). Also note that because of
1135  // this, the porosity returned by the intensive quantities can be
1136  // outside of the physical range [0, 1] in pathetic cases.
1137  const double pv =
1138  ebosSimulator_.model().dofTotalVolume(cellIdx)
1139  * intQuants.porosity().value();
1140 
1141  fip_.fip[FIPDataType::FIP_PV][cellIdx] = pv;
1142  const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
1143 
1144  //Compute hydrocarbon pore volume weighted average pressure.
1145  //If we have no hydrocarbon in region, use pore volume weighted average pressure instead
1146  if (hcpv[regionIdx] > 1e-10) {
1147  fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[regionIdx];
1148  } else {
1149  fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() / tpv[regionIdx];
1150  }
1151 
1152  regionValues[regionIdx][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][cellIdx];
1153  regionValues[regionIdx][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx];
1154  }
1155 
1156  // sum the results over all processes
1157  for(int regionIdx=0; regionIdx < ntFip; ++regionIdx) {
1158  comm.sum(regionValues[regionIdx].data(), regionValues[regionIdx].size());
1159  }
1160 
1161  return regionValues;
1162  }
1163 
1164  SimulationDataContainer getSimulatorData ( const SimulationDataContainer& localState) const
1165  {
1166  typedef std::vector<double> VectorType;
1167 
1168  const auto& ebosModel = ebosSimulator().model();
1169  const auto& phaseUsage = phaseUsage_;
1170 
1171  // extract everything which can possibly be written to disk
1172  const int numCells = ebosModel.numGridDof();
1173  const int num_phases = numPhases();
1174 
1175  SimulationDataContainer simData( numCells, 0, num_phases );
1176 
1177  //Get shorthands for water, oil, gas
1178  const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
1179  const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
1180  const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
1181 
1182  const int aqua_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Aqua ];
1183  const int liquid_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Liquid ];
1184  const int vapour_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Vapour ];
1185 
1186  VectorType zero;
1187 
1188  VectorType& pressureOil = simData.pressure();
1189  VectorType& temperature = simData.temperature();
1190  VectorType& saturation = simData.saturation();
1191 
1192  // WATER
1193  if( aqua_active ) {
1194  simData.registerCellData( "1OVERBW", 1 );
1195  simData.registerCellData( "WAT_DEN", 1 );
1196  simData.registerCellData( "WAT_VISC", 1 );
1197  simData.registerCellData( "WATKR", 1 );
1198  }
1199 
1200  VectorType& bWater = aqua_active ? simData.getCellData( "1OVERBW" ) : zero;
1201  VectorType& rhoWater = aqua_active ? simData.getCellData( "WAT_DEN" ) : zero;
1202  VectorType& muWater = aqua_active ? simData.getCellData( "WAT_VISC" ) : zero;
1203  VectorType& krWater = aqua_active ? simData.getCellData( "WATKR" ) : zero;
1204 
1205  // OIL
1206  if( liquid_active ) {
1207  simData.registerCellData( "1OVERBO", 1 );
1208  simData.registerCellData( "OIL_DEN", 1 );
1209  simData.registerCellData( "OIL_VISC", 1 );
1210  simData.registerCellData( "OILKR", 1 );
1211  }
1212 
1213  VectorType& bOil = liquid_active ? simData.getCellData( "1OVERBO" ) : zero;
1214  VectorType& rhoOil = liquid_active ? simData.getCellData( "OIL_DEN" ) : zero;
1215  VectorType& muOil = liquid_active ? simData.getCellData( "OIL_VISC" ) : zero;
1216  VectorType& krOil = liquid_active ? simData.getCellData( "OILKR" ) : zero;
1217 
1218  // GAS
1219  if( vapour_active ) {
1220  simData.registerCellData( "1OVERBG", 1 );
1221  simData.registerCellData( "GAS_DEN", 1 );
1222  simData.registerCellData( "GAS_VISC", 1 );
1223  simData.registerCellData( "GASKR", 1 );
1224  }
1225 
1226  VectorType& bGas = vapour_active ? simData.getCellData( "1OVERBG" ) : zero;
1227  VectorType& rhoGas = vapour_active ? simData.getCellData( "GAS_DEN" ) : zero;
1228  VectorType& muGas = vapour_active ? simData.getCellData( "GAS_VISC" ) : zero;
1229  VectorType& krGas = vapour_active ? simData.getCellData( "GASKR" ) : zero;
1230 
1231  simData.registerCellData( BlackoilState::GASOILRATIO, 1 );
1232  simData.registerCellData( BlackoilState::RV, 1 );
1233  simData.registerCellData( "RSSAT", 1 );
1234  simData.registerCellData( "RVSAT", 1 );
1235 
1236  VectorType& Rs = simData.getCellData( BlackoilState::GASOILRATIO );
1237  VectorType& Rv = simData.getCellData( BlackoilState::RV );
1238  VectorType& RsSat = simData.getCellData( "RSSAT" );
1239  VectorType& RvSat = simData.getCellData( "RVSAT" );
1240 
1241  simData.registerCellData( "PBUB", 1 );
1242  simData.registerCellData( "PDEW", 1 );
1243 
1244  VectorType& Pb = simData.getCellData( "PBUB" );
1245  VectorType& Pd = simData.getCellData( "PDEW" );
1246 
1247  simData.registerCellData( "SOMAX", 1 );
1248  VectorType& somax = simData.getCellData( "SOMAX" );
1249 
1250  // Two components for hysteresis parameters
1251  // pcSwMdc/krnSwMdc, one for oil-water and one for gas-oil
1252  simData.registerCellData( "PCSWMDC_GO", 1 );
1253  simData.registerCellData( "KRNSWMDC_GO", 1 );
1254 
1255  simData.registerCellData( "PCSWMDC_OW", 1 );
1256  simData.registerCellData( "KRNSWMDC_OW", 1 );
1257 
1258  VectorType& pcSwMdc_go = simData.getCellData( "PCSWMDC_GO" );
1259  VectorType& krnSwMdc_go = simData.getCellData( "KRNSWMDC_GO" );
1260 
1261  VectorType& pcSwMdc_ow = simData.getCellData( "PCSWMDC_OW" );
1262  VectorType& krnSwMdc_ow = simData.getCellData( "KRNSWMDC_OW" );
1263 
1264  if (has_solvent_) {
1265  simData.registerCellData( "SSOL", 1 );
1266  }
1267  VectorType& ssol = has_solvent_ ? simData.getCellData( "SSOL" ) : zero;
1268 
1269  if (has_polymer_) {
1270  simData.registerCellData( "POLYMER", 1 );
1271  }
1272  VectorType& cpolymer = has_polymer_ ? simData.getCellData( "POLYMER" ) : zero;
1273 
1274  std::vector<int> failed_cells_pb;
1275  std::vector<int> failed_cells_pd;
1276  const auto& gridView = ebosSimulator().gridView();
1277  auto elemIt = gridView.template begin</*codim=*/ 0, Dune::Interior_Partition>();
1278  const auto& elemEndIt = gridView.template end</*codim=*/ 0, Dune::Interior_Partition>();
1279  ElementContext elemCtx(ebosSimulator());
1280 
1281  for (; elemIt != elemEndIt; ++elemIt) {
1282  const auto& elem = *elemIt;
1283 
1284  elemCtx.updatePrimaryStencil(elem);
1285  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1286 
1287  const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
1288  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
1289 
1290  const auto& fs = intQuants.fluidState();
1291 
1292  const int satIdx = cellIdx * num_phases;
1293 
1294  pressureOil[cellIdx] = fs.pressure(FluidSystem::oilPhaseIdx).value();
1295 
1296  temperature[cellIdx] = fs.temperature(FluidSystem::oilPhaseIdx).value();
1297 
1298  somax[cellIdx] = ebosSimulator().model().maxOilSaturation(cellIdx);
1299 
1300  const auto& matLawManager = ebosSimulator().problem().materialLawManager();
1301  if (matLawManager->enableHysteresis()) {
1302  matLawManager->oilWaterHysteresisParams(
1303  pcSwMdc_ow[cellIdx],
1304  krnSwMdc_ow[cellIdx],
1305  cellIdx);
1306  matLawManager->gasOilHysteresisParams(
1307  pcSwMdc_go[cellIdx],
1308  krnSwMdc_go[cellIdx],
1309  cellIdx);
1310  }
1311 
1312  if (aqua_active) {
1313  saturation[ satIdx + aqua_pos ] = fs.saturation(FluidSystem::waterPhaseIdx).value();
1314  bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value();
1315  rhoWater[cellIdx] = fs.density(FluidSystem::waterPhaseIdx).value();
1316  muWater[cellIdx] = fs.viscosity(FluidSystem::waterPhaseIdx).value();
1317  krWater[cellIdx] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
1318  }
1319  if (vapour_active) {
1320  saturation[ satIdx + vapour_pos ] = fs.saturation(FluidSystem::gasPhaseIdx).value();
1321  bGas[cellIdx] = fs.invB(FluidSystem::gasPhaseIdx).value();
1322  rhoGas[cellIdx] = fs.density(FluidSystem::gasPhaseIdx).value();
1323  muGas[cellIdx] = fs.viscosity(FluidSystem::gasPhaseIdx).value();
1324  krGas[cellIdx] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
1325  Rs[cellIdx] = fs.Rs().value();
1326  Rv[cellIdx] = fs.Rv().value();
1327  RsSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
1328  FluidSystem::oilPhaseIdx,
1329  intQuants.pvtRegionIndex(),
1330  /*maxOilSaturation=*/1.0).value();
1331  RvSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
1332  FluidSystem::gasPhaseIdx,
1333  intQuants.pvtRegionIndex(),
1334  /*maxOilSaturation=*/1.0).value();
1335  try {
1336  Pb[cellIdx] = FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()).value();
1337  }
1338  catch (const NumericalProblem& e) {
1339  const auto globalIdx = ebosSimulator_.gridManager().grid().globalCell()[cellIdx];
1340  failed_cells_pb.push_back(globalIdx);
1341  }
1342  try {
1343  Pd[cellIdx] = FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()).value();
1344  }
1345  catch (const NumericalProblem& e) {
1346  const auto globalIdx = ebosSimulator_.gridManager().grid().globalCell()[cellIdx];
1347  failed_cells_pd.push_back(globalIdx);
1348  }
1349  }
1350  if( liquid_active )
1351  {
1352  saturation[ satIdx + liquid_pos ] = fs.saturation(FluidSystem::oilPhaseIdx).value();
1353  bOil[cellIdx] = fs.invB(FluidSystem::oilPhaseIdx).value();
1354  rhoOil[cellIdx] = fs.density(FluidSystem::oilPhaseIdx).value();
1355  muOil[cellIdx] = fs.viscosity(FluidSystem::oilPhaseIdx).value();
1356  krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
1357  }
1358 
1359  if (has_solvent_)
1360  {
1361  ssol[cellIdx] = intQuants.solventSaturation().value();
1362  }
1363 
1364  if (has_polymer_)
1365  {
1366  cpolymer[cellIdx] = intQuants.polymerConcentration().value();
1367  }
1368 
1369  // hack to make the intial output of rs and rv Ecl compatible.
1370  // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step
1371  // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite
1372  // rs and rv with the values passed by the localState.
1373  // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values.
1374  if (ebosSimulator_.episodeIndex() < 0 && vapour_active && liquid_active ) {
1375 
1376  Rs[cellIdx] = localState.getCellData( BlackoilState::GASOILRATIO )[cellIdx];
1377  Rv[cellIdx] = localState.getCellData( BlackoilState::RV)[cellIdx];
1378 
1379  // copy the fluidstate and set the new rs and rv values
1380  auto fs_updated = fs;
1381  auto rs_eval = fs_updated.Rs();
1382  rs_eval.setValue( Rs[cellIdx] );
1383  fs_updated.setRs(rs_eval);
1384  auto rv_eval = fs_updated.Rv();
1385  rv_eval.setValue( Rv[cellIdx] );
1386  fs_updated.setRv(rv_eval);
1387 
1388  //re-compute the volume factors, viscosities and densities.
1389  rhoOil[cellIdx] = FluidSystem::density(fs_updated,
1390  FluidSystem::oilPhaseIdx,
1391  intQuants.pvtRegionIndex()).value();
1392  rhoGas[cellIdx] = FluidSystem::density(fs_updated,
1393  FluidSystem::gasPhaseIdx,
1394  intQuants.pvtRegionIndex()).value();
1395 
1396  bOil[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated,
1397  FluidSystem::oilPhaseIdx,
1398  intQuants.pvtRegionIndex()).value();
1399  bGas[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated,
1400  FluidSystem::gasPhaseIdx,
1401  intQuants.pvtRegionIndex()).value();
1402 
1403  muOil[cellIdx] = FluidSystem::viscosity(fs_updated,
1404  FluidSystem::oilPhaseIdx,
1405  intQuants.pvtRegionIndex()).value();
1406  muGas[cellIdx] = FluidSystem::viscosity(fs_updated,
1407  FluidSystem::gasPhaseIdx,
1408  intQuants.pvtRegionIndex()).value();
1409 
1410  }
1411  }
1412 
1413  const size_t max_num_cells_faillog = 20;
1414 
1415  int pb_size = failed_cells_pb.size(), pd_size = failed_cells_pd.size();
1416  std::vector<int> displ_pb, displ_pd, recv_len_pb, recv_len_pd;
1417  const auto& comm = grid_.comm();
1418 
1419  if ( comm.rank() == 0 )
1420  {
1421  displ_pb.resize(comm.size()+1, 0);
1422  displ_pd.resize(comm.size()+1, 0);
1423  recv_len_pb.resize(comm.size());
1424  recv_len_pd.resize(comm.size());
1425  }
1426 
1427  comm.gather(&pb_size, recv_len_pb.data(), 1, 0);
1428  comm.gather(&pd_size, recv_len_pd.data(), 1, 0);
1429  std::partial_sum(recv_len_pb.begin(), recv_len_pb.end(), displ_pb.begin()+1);
1430  std::partial_sum(recv_len_pd.begin(), recv_len_pd.end(), displ_pd.begin()+1);
1431  std::vector<int> global_failed_cells_pb, global_failed_cells_pd;
1432 
1433  if ( comm.rank() == 0 )
1434  {
1435  global_failed_cells_pb.resize(displ_pb.back());
1436  global_failed_cells_pd.resize(displ_pd.back());
1437  }
1438 
1439  comm.gatherv(failed_cells_pb.data(), static_cast<int>(failed_cells_pb.size()),
1440  global_failed_cells_pb.data(), recv_len_pb.data(),
1441  displ_pb.data(), 0);
1442  comm.gatherv(failed_cells_pd.data(), static_cast<int>(failed_cells_pd.size()),
1443  global_failed_cells_pd.data(), recv_len_pd.data(),
1444  displ_pd.data(), 0);
1445  std::sort(global_failed_cells_pb.begin(), global_failed_cells_pb.end());
1446  std::sort(global_failed_cells_pd.begin(), global_failed_cells_pd.end());
1447 
1448  if (global_failed_cells_pb.size() > 0) {
1449  std::stringstream errlog;
1450  errlog << "Finding the bubble point pressure failed for " << global_failed_cells_pb.size() << " cells [";
1451  errlog << global_failed_cells_pb[0];
1452  const size_t max_elems = std::min(max_num_cells_faillog, failed_cells_pb.size());
1453  for (size_t i = 1; i < max_elems; ++i) {
1454  errlog << ", " << global_failed_cells_pb[i];
1455  }
1456  if (global_failed_cells_pb.size() > max_num_cells_faillog) {
1457  errlog << ", ...";
1458  }
1459  errlog << "]";
1460  OpmLog::warning("Bubble point numerical problem", errlog.str());
1461  }
1462  if (global_failed_cells_pd.size() > 0) {
1463  std::stringstream errlog;
1464  errlog << "Finding the dew point pressure failed for " << global_failed_cells_pd.size() << " cells [";
1465  errlog << global_failed_cells_pd[0];
1466  const size_t max_elems = std::min(max_num_cells_faillog, global_failed_cells_pd.size());
1467  for (size_t i = 1; i < max_elems; ++i) {
1468  errlog << ", " << global_failed_cells_pd[i];
1469  }
1470  if (global_failed_cells_pd.size() > max_num_cells_faillog) {
1471  errlog << ", ...";
1472  }
1473  errlog << "]";
1474  OpmLog::warning("Dew point numerical problem", errlog.str());
1475  }
1476 
1477  return simData;
1478  }
1479 
1480  const FIPDataType& getFIPData() const {
1481  return fip_;
1482  }
1483 
1484  const Simulator& ebosSimulator() const
1485  { return ebosSimulator_; }
1486 
1488  const SimulatorReport& failureReport() const
1489  { return failureReport_; }
1490 
1491  protected:
1492  const ISTLSolverType& istlSolver() const
1493  {
1494  assert( istlSolver_ );
1495  return *istlSolver_;
1496  }
1497 
1498  // --------- Data members ---------
1499 
1500  Simulator& ebosSimulator_;
1501  const Grid& grid_;
1502  const ISTLSolverType* istlSolver_;
1503  const PhaseUsage phaseUsage_;
1504  VFPProperties vfp_properties_;
1505  // For each canonical phase -> true if active
1506  const std::vector<bool> active_;
1507  // Size = # active phases. Maps active -> canonical phase indices.
1508  const std::vector<int> cells_; // All grid cells
1509  const bool has_disgas_;
1510  const bool has_vapoil_;
1511  const bool has_solvent_;
1512  const bool has_polymer_;
1513 
1514  ModelParameters param_;
1515  SimulatorReport failureReport_;
1516 
1517  // Well Model
1518  BlackoilWellModel<TypeTag>& well_model_;
1519 
1523  long int global_nc_;
1524 
1525  // rate converter between the surface volume rates and reservoir voidage rates
1526  RateConverterType& rate_converter_;
1527 
1528  std::vector<std::vector<double>> residual_norms_history_;
1529  double current_relaxation_;
1530  BVector dx_old_;
1531  mutable FIPDataType fip_;
1532 
1533  public:
1536  wellModel() { return well_model_; }
1537 
1539  wellModel() const { return well_model_; }
1540 
1541  int numWells() const { return well_model_.numWells(); }
1542 
1544  bool localWellsActive() const { return well_model_.localWellsActive(); }
1545 
1546 
1547 
1548 
1549  public:
1550  int flowPhaseToEbosCompIdx( const int phaseIdx ) const
1551  {
1552  const auto& pu = phaseUsage_;
1553  if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
1554  return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1555  if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
1556  return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1557  if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
1558  return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1559 
1560  // for other phases return the index
1561  return phaseIdx;
1562  }
1563 
1564  int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
1565  {
1566  const auto& pu = phaseUsage_;
1567  if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
1568  return FluidSystem::waterPhaseIdx;
1569  if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
1570  return FluidSystem::oilPhaseIdx;
1571  if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
1572  return FluidSystem::gasPhaseIdx;
1573 
1574  assert(phaseIdx < 3);
1575  // for other phases return the index
1576  return phaseIdx;
1577  }
1578 
1579  private:
1580 
1581  void updateRateConverter()
1582  {
1583  rate_converter_.defineState<ElementContext>(ebosSimulator_);
1584  }
1585 
1586 
1587  public:
1588  void beginReportStep()
1589  {
1590  isBeginReportStep_ = true;
1591  }
1592 
1593  void endReportStep()
1594  {
1595  ebosSimulator_.problem().endEpisode();
1596  }
1597 
1598  private:
1599  void assembleMassBalanceEq(const SimulatorTimerInterface& timer,
1600  const int iterationIdx)
1601  {
1602  ebosSimulator_.startNextEpisode( timer.currentStepLength() );
1603  ebosSimulator_.setEpisodeIndex( timer.reportStepNum() );
1604  ebosSimulator_.setTimeStepIndex( timer.reportStepNum() );
1605  ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
1606 
1607  static int prevEpisodeIdx = 10000;
1608 
1609  // notify ebos about the end of the previous episode and time step if applicable
1610  if (isBeginReportStep_) {
1611  isBeginReportStep_ = false;
1612  ebosSimulator_.problem().beginEpisode();
1613  }
1614 
1615  // doing the notifactions here is conceptually wrong and also causes the
1616  // endTimeStep() and endEpisode() methods to be not called for the
1617  // simulation's last time step and episode.
1618  if (ebosSimulator_.model().newtonMethod().numIterations() == 0
1619  && prevEpisodeIdx < timer.reportStepNum())
1620  {
1621  ebosSimulator_.problem().endTimeStep();
1622  }
1623 
1624  ebosSimulator_.setTimeStepSize( timer.currentStepLength() );
1625  if (ebosSimulator_.model().newtonMethod().numIterations() == 0)
1626  {
1627  ebosSimulator_.problem().beginTimeStep();
1628  }
1629 
1630  ebosSimulator_.problem().beginIteration();
1631  ebosSimulator_.model().linearizer().linearize();
1632  ebosSimulator_.problem().endIteration();
1633 
1634  prevEpisodeIdx = ebosSimulator_.episodeIndex();
1635 
1636  if (param_.update_equations_scaling_) {
1637  std::cout << "equation scaling not suported yet" << std::endl;
1638  //updateEquationsScaling();
1639  }
1640  }
1641 
1642  double dpMaxRel() const { return param_.dp_max_rel_; }
1643  double dsMax() const { return param_.ds_max_; }
1644  double drMaxRel() const { return param_.dr_max_rel_; }
1645  double maxResidualAllowed() const { return param_.max_residual_allowed_; }
1646 
1647  public:
1648  bool isBeginReportStep_;
1649  std::vector<bool> wasSwitched_;
1650  };
1651 } // namespace Opm
1652 
1653 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
bool localWellsActive() const
return true if wells are available on this process
Definition: BlackoilModelEbos.hpp:1544
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &, const WellState &)
Called once before each time step.
Definition: BlackoilModelEbos.hpp:206
bool getConvergence(const SimulatorTimerInterface &timer, const int iteration, std::vector< double > &residual_norms)
Compute convergence based on total mass balance (tol_mb) and maximum residual mass balance (tol_cnv)...
Definition: BlackoilModelEbos.hpp:837
Definition: BlackoilModelEnums.hpp:66
SimulatorReport assemble(const SimulatorTimerInterface &timer, const int iterationIdx, WellState &well_state)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilModelEbos.hpp:377
int max_strict_iter_
Maximum number of Newton iterations before we give up on the CNV convergence criterion.
Definition: BlackoilModelParameters.hpp:73
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:66
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition: BlackoilModelEbos.hpp:771
Adapter to turn a matrix into a linear operator.
Definition: BlackoilModelEbos.hpp:551
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition: BlackoilModelEbos.hpp:1536
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelEbos.hpp:1003
const SimulatorReport & failureReport() const
return the statistics if the nonlinearIteration() method failed
Definition: BlackoilModelEbos.hpp:1488
double tolerance_cnv_
Local convergence tolerance (max of local saturation errors).
Definition: BlackoilModelParameters.hpp:48
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition: BlackoilModelEbos.hpp:498
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolver.hpp:323
A model implementation for three-phase black oil.
Definition: BlackoilModelEbos.hpp:104
AutoDiffMatrix is a wrapper class that optimizes matrix operations.
Definition: AutoDiffMatrix.hpp:43
double tolerance_mb_
Relative mass balance tolerance (total mass balance error).
Definition: BlackoilModelParameters.hpp:46
void solve(Operator &opA, Vector &x, Vector &istlb, ScalarProd &sp, Precond &precond, Dune::InverseOperatorResult &result) const
Solve the system using the given preconditioner and scalar product.
Definition: ISTLSolver.hpp:487
void afterStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, WellState &well_state)
Called once after each time step.
Definition: BlackoilModelEbos.hpp:364
SimulatorReport nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver, ReservoirState &, WellState &well_state)
Called once per nonlinear iteration.
Definition: BlackoilModelEbos.hpp:242
void solveJacobianSystem(BVector &x) const
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilModelEbos.hpp:505
Facility for converting component rates at surface conditions to phase (voidage) rates at reservoir c...
virtual bool lastStepFailed() const =0
Return true if last time step failed.
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:31
std::vector< std::vector< double > > computeFluidInPlace(const T &, const std::vector< int > &fipnum) const
Wrapper required due to not following generic API.
Definition: BlackoilModelEbos.hpp:1026
virtual double currentStepLength() const =0
Current step length.
BlackoilModelEbos(Simulator &ebosSimulator, const ModelParameters &param, BlackoilWellModel< TypeTag > &well_model, RateConverterType &rate_converter, const NewtonIterationBlackoilInterface &linsolver, const bool terminal_output)
Construct the model.
Definition: BlackoilModelEbos.hpp:155
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
The solver category.
Definition: BlackoilModelEbos.hpp:569
void updateState(const BVector &dx)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilModelEbos.hpp:630
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolver.hpp:374
int sizeNonLinear() const
The size (number of unknowns) of the nonlinear system of equations.
Definition: BlackoilModelEbos.hpp:490
bool terminal_output_
Whether we print something to std::cout.
Definition: BlackoilModelEbos.hpp:1521
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
long int global_nc_
The number of cells of the global grid.
Definition: BlackoilModelEbos.hpp:1523
void defineState(const BlackoilState &state, const boost::any &info=boost::any())
Compute average hydrocarbon pressure and maximum dissolution and evaporation at average hydrocarbon p...
Definition: RateConverter.hpp:440
bool use_update_stabilization_
Try to detect oscillation or stagnation.
Definition: BlackoilModelParameters.hpp:82