My Project
Loading...
Searching...
No Matches
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, 2016, 2017 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 <fmt/format.h>
28
29#include <ebos/eclproblem.hh>
30
31#include <opm/common/ErrorMacros.hpp>
32#include <opm/common/Exceptions.hpp>
33#include <opm/common/OpmLog/OpmLog.hpp>
34
35#include <opm/core/props/phaseUsageFromDeck.hpp>
36
37#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
38#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
39
40#include <opm/simulators/aquifers/AquiferGridUtils.hpp>
41#include <opm/simulators/aquifers/BlackoilAquiferModel.hpp>
42#include <opm/simulators/flow/BlackoilModelEbosNldd.hpp>
43#include <opm/simulators/flow/countGlobalCells.hpp>
44#include <opm/simulators/flow/NonlinearSolverEbos.hpp>
45#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
46#include <opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp>
47#include <opm/simulators/timestepping/ConvergenceReport.hpp>
48#include <opm/simulators/timestepping/SimulatorReport.hpp>
49#include <opm/simulators/timestepping/SimulatorTimer.hpp>
50#include <opm/simulators/utils/ComponentName.hpp>
51#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
52#include <opm/simulators/utils/ParallelCommunication.hpp>
53#include <opm/simulators/wells/BlackoilWellModel.hpp>
54
55#include <dune/common/timer.hh>
56
57#include <fmt/format.h>
58
59#include <algorithm>
60#include <cassert>
61#include <cmath>
62#include <iomanip>
63#include <ios>
64#include <limits>
65#include <memory>
66#include <sstream>
67#include <tuple>
68#include <utility>
69#include <vector>
70
71namespace Opm::Properties {
72
73namespace TTag {
78}
79template<class TypeTag>
80struct OutputDir<TypeTag, TTag::EclFlowProblem> {
81 static constexpr auto value = "";
82};
83template<class TypeTag>
84struct EnableDebuggingChecks<TypeTag, TTag::EclFlowProblem> {
85 static constexpr bool value = false;
86};
87// default in flow is to formulate the equations in surface volumes
88template<class TypeTag>
89struct BlackoilConserveSurfaceVolume<TypeTag, TTag::EclFlowProblem> {
90 static constexpr bool value = true;
91};
92template<class TypeTag>
93struct UseVolumetricResidual<TypeTag, TTag::EclFlowProblem> {
94 static constexpr bool value = false;
95};
96
97template<class TypeTag>
98struct EclAquiferModel<TypeTag, TTag::EclFlowProblem> {
100};
101
102// disable all extensions supported by black oil model. this should not really be
103// necessary but it makes things a bit more explicit
104template<class TypeTag>
105struct EnablePolymer<TypeTag, TTag::EclFlowProblem> {
106 static constexpr bool value = false;
107};
108template<class TypeTag>
109struct EnableSolvent<TypeTag, TTag::EclFlowProblem> {
110 static constexpr bool value = false;
111};
112template<class TypeTag>
113struct EnableTemperature<TypeTag, TTag::EclFlowProblem> {
114 static constexpr bool value = true;
115};
116template<class TypeTag>
117struct EnableEnergy<TypeTag, TTag::EclFlowProblem> {
118 static constexpr bool value = false;
119};
120template<class TypeTag>
121struct EnableFoam<TypeTag, TTag::EclFlowProblem> {
122 static constexpr bool value = false;
123};
124template<class TypeTag>
125struct EnableBrine<TypeTag, TTag::EclFlowProblem> {
126 static constexpr bool value = false;
127};
128template<class TypeTag>
129struct EnableSaltPrecipitation<TypeTag, TTag::EclFlowProblem> {
130 static constexpr bool value = false;
131};
132template<class TypeTag>
133struct EnableMICP<TypeTag, TTag::EclFlowProblem> {
134 static constexpr bool value = false;
135};
136
137template<class TypeTag>
138struct EclWellModel<TypeTag, TTag::EclFlowProblem> {
140};
141template<class TypeTag>
142struct LinearSolverSplice<TypeTag, TTag::EclFlowProblem> {
144};
145
146} // namespace Opm::Properties
147
148namespace Opm {
149
156 template <class TypeTag>
158 {
159 public:
160 // --------- Types and enums ---------
162
175
176 static constexpr int numEq = Indices::numEq;
177 static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
178 static constexpr int contiZfracEqIdx = Indices::contiZfracEqIdx;
179 static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
180 static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
181 static constexpr int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
182 static constexpr int contiFoamEqIdx = Indices::contiFoamEqIdx;
183 static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx;
184 static constexpr int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
185 static constexpr int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
186 static constexpr int contiUreaEqIdx = Indices::contiUreaEqIdx;
187 static constexpr int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
188 static constexpr int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
189 static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
190 static constexpr int zFractionIdx = Indices::zFractionIdx;
191 static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
192 static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
193 static constexpr int temperatureIdx = Indices::temperatureIdx;
194 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
195 static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx;
196 static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
197 static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
198 static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
199 static constexpr int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
200 static constexpr int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
201
202 using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
203 using MatrixBlockType = typename SparseMatrixAdapter::MatrixBlock;
204 using Mat = typename SparseMatrixAdapter::IstlMatrix;
205 using BVector = Dune::BlockVector<VectorBlockType>;
206
208
209 // --------- Public methods ---------
210
221 BlackoilModelEbos(Simulator& ebosSimulator,
222 const ModelParameters& param,
224 const bool terminal_output)
225 : ebosSimulator_(ebosSimulator)
226 , grid_(ebosSimulator_.vanguard().grid())
227 , phaseUsage_(phaseUsageFromDeck(eclState()))
228 , param_( param )
229 , well_model_ (well_model)
231 , current_relaxation_(1.0)
232 , dx_old_(ebosSimulator_.model().numGridDof())
233 {
234 // compute global sum of number of cells
235 global_nc_ = detail::countGlobalCells(grid_);
236 convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves.
237 // TODO: remember to fix!
238 if (param_.nonlinear_solver_ == "nldd") {
239 if (terminal_output) {
240 OpmLog::info("Using Non-Linear Domain Decomposition solver (nldd).");
241 }
242 nlddSolver_ = std::make_unique<BlackoilModelEbosNldd<TypeTag>>(*this);
243 } else if (param_.nonlinear_solver_ == "newton") {
244 if (terminal_output) {
245 OpmLog::info("Using Newton nonlinear solver.");
246 }
247 } else {
248 OPM_THROW(std::runtime_error, "Unknown nonlinear solver option: " + param_.nonlinear_solver_);
249 }
250 }
251
252
253 bool isParallel() const
254 { return grid_.comm().size() > 1; }
255
256
257 const EclipseState& eclState() const
258 { return ebosSimulator_.vanguard().eclState(); }
259
260
264 {
266 Dune::Timer perfTimer;
267 perfTimer.start();
268 // update the solution variables in ebos
269 if ( timer.lastStepFailed() ) {
270 ebosSimulator_.model().updateFailed();
271 } else {
272 ebosSimulator_.model().advanceTimeLevel();
273 }
274
275 // Set the timestep size, episode index, and non-linear iteration index
276 // for ebos explicitly. ebos needs to know the report step/episode index
277 // because of timing dependent data despite the fact that flow uses its
278 // own time stepper. (The length of the episode does not matter, though.)
279 ebosSimulator_.setTime(timer.simulationTimeElapsed());
280 ebosSimulator_.setTimeStepSize(timer.currentStepLength());
281 ebosSimulator_.model().newtonMethod().setIterationIndex(0);
282
283 ebosSimulator_.problem().beginTimeStep();
284
285 unsigned numDof = ebosSimulator_.model().numGridDof();
286 wasSwitched_.resize(numDof);
287 std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
288
289 if (param_.update_equations_scaling_) {
290 OpmLog::error("Equation scaling not supported");
291 //updateEquationsScaling();
292 }
293
294 if (nlddSolver_) {
295 nlddSolver_->prepareStep();
296 }
297
298 report.pre_post_time += perfTimer.stop();
299
300 return report;
301 }
302
303
304 void initialLinearization(SimulatorReportSingle& report,
305 const int iteration,
306 const int minIter,
308 {
309 // ----------- Set up reports and timer -----------
310 failureReport_ = SimulatorReportSingle();
311 Dune::Timer perfTimer;
312
313 perfTimer.start();
314 report.total_linearizations = 1;
315
316 // ----------- Assemble -----------
317 try {
319 report.assemble_time += perfTimer.stop();
320 }
321 catch (...) {
322 report.assemble_time += perfTimer.stop();
323 failureReport_ += report;
324 throw; // continue throwing the stick
325 }
326
327 // ----------- Check if converged -----------
328 std::vector<double> residual_norms;
329 perfTimer.reset();
330 perfTimer.start();
331 // the step is not considered converged until at least minIter iterations is done
332 {
333 auto convrep = getConvergence(timer, iteration, residual_norms);
334 report.converged = convrep.converged() && iteration > minIter;
335 ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
336 convergence_reports_.back().report.push_back(std::move(convrep));
337
338 // Throw if any NaN or too large residual found.
339 if (severity == ConvergenceReport::Severity::NotANumber) {
340 OPM_THROW(NumericalProblem, "NaN residual found!");
341 } else if (severity == ConvergenceReport::Severity::TooLarge) {
342 OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!");
343 }
344 }
345 report.update_time += perfTimer.stop();
346 residual_norms_history_.push_back(residual_norms);
347 }
348
349
357 template <class NonlinearSolverType>
361 {
362 if (iteration == 0) {
363 // For each iteration we store in a vector the norms of the residual of
364 // the mass balance for each active phase, the well flux and the well equations.
365 residual_norms_history_.clear();
366 current_relaxation_ = 1.0;
367 dx_old_ = 0.0;
368 convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}});
369 convergence_reports_.back().report.reserve(11);
370 }
371
372 if ((this->param_.nonlinear_solver_ != "nldd") ||
373 (iteration < this->param_.nldd_num_initial_newton_iter_))
374 {
375 return this->nonlinearIterationNewton(iteration, timer, nonlinear_solver);
376 }
377 else {
378 return this->nlddSolver_->nonlinearIterationNldd(iteration, timer, nonlinear_solver);
379 }
380 }
381
382
383 template <class NonlinearSolverType>
384 SimulatorReportSingle nonlinearIterationNewton(const int iteration,
387 {
388
389 // ----------- Set up reports and timer -----------
391 Dune::Timer perfTimer;
392
393 this->initialLinearization(report, iteration, nonlinear_solver.minIter(), timer);
394
395 // ----------- If not converged, solve linear system and do Newton update -----------
396 if (!report.converged) {
397 perfTimer.reset();
398 perfTimer.start();
399 report.total_newton_iterations = 1;
400
401 // Compute the nonlinear update.
402 unsigned nc = ebosSimulator_.model().numGridDof();
403 BVector x(nc);
404
405 // Solve the linear system.
406 linear_solve_setup_time_ = 0.0;
407 try {
408 // Apply the Schur complement of the well model to
409 // the reservoir linearized equations.
410 // Note that linearize may throw for MSwells.
411 wellModel().linearize(ebosSimulator().model().linearizer().jacobian(),
412 ebosSimulator().model().linearizer().residual());
413
414 // ---- Solve linear system ----
416
417 report.linear_solve_setup_time += linear_solve_setup_time_;
418 report.linear_solve_time += perfTimer.stop();
419 report.total_linear_iterations += linearIterationsLastSolve();
420 }
421 catch (...) {
422 report.linear_solve_setup_time += linear_solve_setup_time_;
423 report.linear_solve_time += perfTimer.stop();
424 report.total_linear_iterations += linearIterationsLastSolve();
425
426 failureReport_ += report;
427 throw; // re-throw up
428 }
429
430 perfTimer.reset();
431 perfTimer.start();
432
433 // handling well state update before oscillation treatment is a decision based
434 // on observation to avoid some big performance degeneration under some circumstances.
435 // there is no theorectical explanation which way is better for sure.
436 wellModel().postSolve(x);
437
438 if (param_.use_update_stabilization_) {
439 // Stabilize the nonlinear update.
440 bool isOscillate = false;
441 bool isStagnate = false;
442 nonlinear_solver.detectOscillations(residual_norms_history_, residual_norms_history_.size() - 1, isOscillate, isStagnate);
443 if (isOscillate) {
444 current_relaxation_ -= nonlinear_solver.relaxIncrement();
445 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
446 if (terminalOutputEnabled()) {
447 std::string msg = " Oscillating behavior detected: Relaxation set to "
448 + std::to_string(current_relaxation_);
449 OpmLog::info(msg);
450 }
451 }
452 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
453 }
454
455 // ---- Newton update ----
456 // Apply the update, with considering model-dependent limitations and
457 // chopping of the update.
459
460 report.update_time += perfTimer.stop();
461 }
462
463 return report;
464 }
465
466
471 {
473 Dune::Timer perfTimer;
474 perfTimer.start();
475 ebosSimulator_.problem().endTimeStep();
476 report.pre_post_time += perfTimer.stop();
477 return report;
478 }
479
482 const int iterationIdx)
483 {
484 // -------- Mass balance equations --------
485 ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
486 ebosSimulator_.problem().beginIteration();
487 ebosSimulator_.model().linearizer().linearizeDomain();
488 ebosSimulator_.problem().endIteration();
489 return wellModel().lastReport();
490 }
491
492 // compute the "relative" change of the solution between time steps
493 double relativeChange() const
494 {
495 Scalar resultDelta = 0.0;
496 Scalar resultDenom = 0.0;
497
498 const auto& elemMapper = ebosSimulator_.model().elementMapper();
499 const auto& gridView = ebosSimulator_.gridView();
500 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
501 unsigned globalElemIdx = elemMapper.index(elem);
502 const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
503
504 Scalar pressureNew;
505 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
506
507 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
508 Scalar oilSaturationNew = 1.0;
509 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
510 FluidSystem::numActivePhases() > 1 &&
511 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
512 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
513 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
514 }
515
516 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
517 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
518 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
519 assert(Indices::compositionSwitchIdx >= 0 );
520 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
521 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
522 }
523
524 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
525 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
526 }
527
528 const auto& priVarsOld = ebosSimulator_.model().solution(/*timeIdx=*/1)[globalElemIdx];
529
530 Scalar pressureOld;
531 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
532
533 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
534 Scalar oilSaturationOld = 1.0;
535
536 // NB fix me! adding pressures changes to satutation changes does not make sense
537 Scalar tmp = pressureNew - pressureOld;
538 resultDelta += tmp*tmp;
539 resultDenom += pressureNew*pressureNew;
540
541 if (FluidSystem::numActivePhases() > 1) {
542 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
543 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSwitchIdx];
544 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
545 }
546
547 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
548 {
549 assert(Indices::compositionSwitchIdx >= 0 );
550 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
551 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
552 }
553
554 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
555 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
556 }
557 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
558 Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
559 resultDelta += tmpSat*tmpSat;
560 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
561 assert(std::isfinite(resultDelta));
562 assert(std::isfinite(resultDenom));
563 }
564 }
565 }
566
567 resultDelta = gridView.comm().sum(resultDelta);
568 resultDenom = gridView.comm().sum(resultDenom);
569
570 if (resultDenom > 0.0)
571 return resultDelta/resultDenom;
572 return 0.0;
573 }
574
575
578 {
579 return ebosSimulator_.model().newtonMethod().linearSolver().iterations ();
580 }
581
582
583 // Obtain reference to linear solver setup time
584 double& linearSolveSetupTime()
585 {
586 return linear_solve_setup_time_;
587 }
588
589
592 void solveJacobianSystem(BVector& x)
593 {
594
595 auto& ebosJac = ebosSimulator_.model().linearizer().jacobian().istlMatrix();
596 auto& ebosResid = ebosSimulator_.model().linearizer().residual();
597 auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
598
599 const int numSolvers = ebosSolver.numAvailableSolvers();
600 if ((numSolvers > 1) && (ebosSolver.getSolveCount() % 100 == 0)) {
601
602 if ( terminal_output_ ) {
603 OpmLog::debug("\nRunning speed test for comparing available linear solvers.");
604 }
605
606 Dune::Timer perfTimer;
607 std::vector<double> times(numSolvers);
608 std::vector<double> setupTimes(numSolvers);
609
610 x = 0.0;
611 std::vector<BVector> x_trial(numSolvers, x);
612 for (int solver = 0; solver < numSolvers; ++solver) {
613 BVector x0(x);
614 ebosSolver.setActiveSolver(solver);
615 perfTimer.start();
616 ebosSolver.prepare(ebosJac, ebosResid);
617 setupTimes[solver] = perfTimer.stop();
618 perfTimer.reset();
619 ebosSolver.setResidual(ebosResid);
620 perfTimer.start();
621 ebosSolver.solve(x_trial[solver]);
622 times[solver] = perfTimer.stop();
623 perfTimer.reset();
624 if (terminal_output_) {
625 OpmLog::debug(fmt::format("Solver time {}: {}", solver, times[solver]));
626 }
627 }
628
629 int fastest_solver = std::min_element(times.begin(), times.end()) - times.begin();
630 // Use timing on rank 0 to determine fastest, must be consistent across ranks.
631 grid_.comm().broadcast(&fastest_solver, 1, 0);
632 linear_solve_setup_time_ = setupTimes[fastest_solver];
634 ebosSolver.setActiveSolver(fastest_solver);
635
636 } else {
637
638 // set initial guess
639 x = 0.0;
640
641 Dune::Timer perfTimer;
642 perfTimer.start();
643 ebosSolver.prepare(ebosJac, ebosResid);
644 linear_solve_setup_time_ = perfTimer.stop();
645 ebosSolver.setResidual(ebosResid);
646 // actually, the error needs to be calculated after setResidual in order to
647 // account for parallelization properly. since the residual of ECFV
648 // discretizations does not need to be synchronized across processes to be
649 // consistent, this is not relevant for OPM-flow...
650 ebosSolver.solve(x);
651 }
652 }
653
654
656 void updateSolution(const BVector& dx)
657 {
659 auto& ebosNewtonMethod = ebosSimulator_.model().newtonMethod();
660 SolutionVector& solution = ebosSimulator_.model().solution(/*timeIdx=*/0);
661
662 ebosNewtonMethod.update_(/*nextSolution=*/solution,
663 /*curSolution=*/solution,
664 /*update=*/dx,
665 /*resid=*/dx); // the update routines of the black
666 // oil model do not care about the
667 // residual
668
669 // if the solution is updated, the intensive quantities need to be recalculated
670 {
672 ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
673 ebosSimulator_.problem().eclWriter()->mutableEclOutputModule().invalidateLocalData();
674 }
675 }
676
679 {
680 return terminal_output_;
681 }
682
683 std::tuple<double,double> convergenceReduction(Parallel::Communication comm,
684 const double pvSumLocal,
685 const double numAquiferPvSumLocal,
686 std::vector< Scalar >& R_sum,
687 std::vector< Scalar >& maxCoeff,
688 std::vector< Scalar >& B_avg)
689 {
690 OPM_TIMEBLOCK(convergenceReduction);
691 // Compute total pore volume (use only owned entries)
692 double pvSum = pvSumLocal;
694
695 if( comm.size() > 1 )
696 {
697 // global reduction
698 std::vector< Scalar > sumBuffer;
699 std::vector< Scalar > maxBuffer;
700 const int numComp = B_avg.size();
701 sumBuffer.reserve( 2*numComp + 2 ); // +2 for (numAquifer)pvSum
702 maxBuffer.reserve( numComp );
703 for( int compIdx = 0; compIdx < numComp; ++compIdx )
704 {
705 sumBuffer.push_back( B_avg[ compIdx ] );
706 sumBuffer.push_back( R_sum[ compIdx ] );
707 maxBuffer.push_back( maxCoeff[ compIdx ] );
708 }
709
710 // Compute total pore volume
711 sumBuffer.push_back( pvSum );
712 sumBuffer.push_back( numAquiferPvSum );
713
714 // compute global sum
715 comm.sum( sumBuffer.data(), sumBuffer.size() );
716
717 // compute global max
718 comm.max( maxBuffer.data(), maxBuffer.size() );
719
720 // restore values to local variables
721 for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
722 {
723 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
724 ++buffIdx;
725
726 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
727 }
728
729 for( int compIdx = 0; compIdx < numComp; ++compIdx )
730 {
731 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
732 }
733
734 // restore global pore volume
735 pvSum = sumBuffer[sumBuffer.size()-2];
736 numAquiferPvSum = sumBuffer.back();
737 }
738
739 // return global pore volume
740 return {pvSum, numAquiferPvSum};
741 }
742
746 std::pair<double,double> localConvergenceData(std::vector<Scalar>& R_sum,
747 std::vector<Scalar>& maxCoeff,
748 std::vector<Scalar>& B_avg,
749 std::vector<int>& maxCoeffCell)
750 {
752 double pvSumLocal = 0.0;
753 double numAquiferPvSumLocal = 0.0;
754 const auto& ebosModel = ebosSimulator_.model();
755 const auto& ebosProblem = ebosSimulator_.problem();
756
757 const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
758
759 ElementContext elemCtx(ebosSimulator_);
760 const auto& gridView = ebosSimulator().gridView();
762 OPM_BEGIN_PARALLEL_TRY_CATCH();
763 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
764 elemCtx.updatePrimaryStencil(elem);
765 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
766
767 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
768 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
769 const auto& fs = intQuants.fluidState();
770
771 const auto pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) *
772 ebosModel.dofTotalVolume(cell_idx);
774
776 {
778 }
779
780 this->getMaxCoeff(cell_idx, intQuants, fs, ebosResid, pvValue,
782 }
783
784 OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm());
785
786 // compute local average in terms of global number of elements
787 const int bSize = B_avg.size();
788 for ( int i = 0; i<bSize; ++i )
789 {
790 B_avg[ i ] /= Scalar( global_nc_ );
791 }
792
794 }
795
796
799 double computeCnvErrorPv(const std::vector<Scalar>& B_avg, double dt)
800 {
802 double errorPV{};
803 const auto& ebosModel = ebosSimulator_.model();
804 const auto& ebosProblem = ebosSimulator_.problem();
805 const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
806 const auto& gridView = ebosSimulator().gridView();
807 ElementContext elemCtx(ebosSimulator_);
809
810 OPM_BEGIN_PARALLEL_TRY_CATCH();
811 for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder))
812 {
813 // Skip cells of numerical Aquifer
815 {
816 continue;
817 }
818 elemCtx.updatePrimaryStencil(elem);
819 // elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
820 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
821 const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx );
822 const auto& cellResidual = ebosResid[cell_idx];
823 bool cnvViolated = false;
824
825 for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx)
826 {
827 using std::abs;
829 cnvViolated = cnvViolated || (abs(CNV) > param_.tolerance_cnv_);
830 }
831
832 if (cnvViolated)
833 {
834 errorPV += pvValue;
835 }
836 }
837
838 OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::ComputeCnvError() failed: ", grid_.comm());
839
840 return grid_.comm().sum(errorPV);
841 }
842
843
844 void updateTUNING(const Tuning& tuning) {
845 param_.tolerance_mb_ = tuning.XXXMBE;
846 if ( terminal_output_ ) {
847 OpmLog::debug(fmt::format("Setting BlackoilModelEbos mass balance limit (XXXMBE) to {:.2e}", tuning.XXXMBE));
848 }
849 }
850
851
852 ConvergenceReport getReservoirConvergence(const double reportTime,
853 const double dt,
854 const int iteration,
855 std::vector<Scalar>& B_avg,
856 std::vector<Scalar>& residual_norms)
857 {
858 OPM_TIMEBLOCK(getReservoirConvergence);
859 using Vector = std::vector<Scalar>;
860
861 const int numComp = numEq;
862 Vector R_sum(numComp, 0.0 );
863 Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
864 std::vector<int> maxCoeffCell(numComp, -1);
865 const auto [ pvSumLocal, numAquiferPvSumLocal] = localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell);
866
867 // compute global sum and max of quantities
868 const auto [ pvSum, numAquiferPvSum ] =
869 convergenceReduction(grid_.comm(), pvSumLocal,
870 numAquiferPvSumLocal,
871 R_sum, maxCoeff, B_avg);
872
873 auto cnvErrorPvFraction = computeCnvErrorPv(B_avg, dt);
874 cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
875
876 const double tol_mb = param_.tolerance_mb_;
877 // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0.
878 // For each iteration, we need to determine whether to use the relaxed CNV tolerance.
879 // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0.
880 const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.min_strict_cnv_iter_;
881 const double tol_cnv = use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
882
883 // Finish computation
884 std::vector<Scalar> CNV(numComp);
885 std::vector<Scalar> mass_balance_residual(numComp);
886 for ( int compIdx = 0; compIdx < numComp; ++compIdx )
887 {
888 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
889 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
890 residual_norms.push_back(CNV[compIdx]);
891 }
892
893 // Create convergence report.
894 ConvergenceReport report{reportTime};
895 using CR = ConvergenceReport;
896 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
897 double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
898 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
899 CR::ReservoirFailure::Type::Cnv };
900 double tol[2] = { tol_mb, tol_cnv };
901 for (int ii : {0, 1}) {
902 if (std::isnan(res[ii])) {
903 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
904 if ( terminal_output_ ) {
905 OpmLog::debug("NaN residual for " + this->compNames_.name(compIdx) + " equation.");
906 }
907 } else if (res[ii] > maxResidualAllowed()) {
908 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
909 if ( terminal_output_ ) {
910 OpmLog::debug("Too large residual for " + this->compNames_.name(compIdx) + " equation.");
911 }
912 } else if (res[ii] < 0.0) {
913 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
914 if ( terminal_output_ ) {
915 OpmLog::debug("Negative residual for " + this->compNames_.name(compIdx) + " equation.");
916 }
917 } else if (res[ii] > tol[ii]) {
918 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
919 }
920 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
921 }
922 }
923
924 // Output of residuals.
925 if ( terminal_output_ )
926 {
927 // Only rank 0 does print to std::cout
928 if (iteration == 0) {
929 std::string msg = "Iter";
930 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
931 msg += " MB(";
932 msg += this->compNames_.name(compIdx)[0];
933 msg += ") ";
934 }
935 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
936 msg += " CNV(";
937 msg += this->compNames_.name(compIdx)[0];
938 msg += ") ";
939 }
940 OpmLog::debug(msg);
941 }
942 std::ostringstream ss;
943 const std::streamsize oprec = ss.precision(3);
944 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
945 ss << std::setw(4) << iteration;
946 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
947 ss << std::setw(11) << mass_balance_residual[compIdx];
948 }
949 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
950 ss << std::setw(11) << CNV[compIdx];
951 }
952 ss.precision(oprec);
953 ss.flags(oflags);
954 OpmLog::debug(ss.str());
955 }
956
957 return report;
958 }
959
966 const int iteration,
967 std::vector<double>& residual_norms)
968 {
970 // Get convergence reports for reservoir and wells.
971 std::vector<Scalar> B_avg(numEq, 0.0);
972 auto report = getReservoirConvergence(timer.simulationTimeElapsed(),
973 timer.currentStepLength(),
975 {
976 OPM_TIMEBLOCK(getWellConvergence);
977 report += wellModel().getWellConvergence(B_avg, /*checkWellGroupControls*/report.converged());
978 }
979 return report;
980 }
981
982
984 int numPhases() const
985 {
986 return phaseUsage_.num_phases;
987 }
988
990 template<class T>
991 std::vector<std::vector<double> >
992 computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
993 {
995 }
996
998 std::vector<std::vector<double> >
999 computeFluidInPlace(const std::vector<int>& /*fipnum*/) const
1000 {
1002 //assert(true)
1003 //return an empty vector
1004 std::vector<std::vector<double> > regionValues(0, std::vector<double>(0,0.0));
1005 return regionValues;
1006 }
1007
1008 const Simulator& ebosSimulator() const
1009 { return ebosSimulator_; }
1010
1011 Simulator& ebosSimulator()
1012 { return ebosSimulator_; }
1013
1016 { return failureReport_; }
1017
1020 {
1021 return nlddSolver_ ? nlddSolver_->localAccumulatedReports()
1023 }
1024
1025 const std::vector<StepReport>& stepReports() const
1026 {
1027 return convergence_reports_;
1028 }
1029
1030 protected:
1031 // --------- Data members ---------
1032
1033 Simulator& ebosSimulator_;
1034 const Grid& grid_;
1035 const PhaseUsage phaseUsage_;
1036 static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
1037 static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
1038 static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
1039 static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
1040 static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
1041 static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
1042 static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
1043 static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
1044
1045 ModelParameters param_;
1046 SimulatorReportSingle failureReport_;
1047
1048 // Well Model
1049 BlackoilWellModel<TypeTag>& well_model_;
1050
1054 long int global_nc_;
1055
1056 std::vector<std::vector<double>> residual_norms_history_;
1057 double current_relaxation_;
1058 BVector dx_old_;
1059
1060 std::vector<StepReport> convergence_reports_;
1061 ComponentName compNames_{};
1062
1063 std::unique_ptr<BlackoilModelEbosNldd<TypeTag>> nlddSolver_;
1064
1065 public:
1068 wellModel() { return well_model_; }
1069
1071 wellModel() const { return well_model_; }
1072
1073 void beginReportStep()
1074 {
1075 ebosSimulator_.problem().beginEpisode();
1076 }
1077
1078 void endReportStep()
1079 {
1080 ebosSimulator_.problem().endEpisode();
1081 }
1082
1083 template<class FluidState, class Residual>
1084 void getMaxCoeff(const unsigned cell_idx,
1085 const IntensiveQuantities& intQuants,
1086 const FluidState& fs,
1087 const Residual& ebosResid,
1088 const Scalar pvValue,
1089 std::vector<Scalar>& B_avg,
1090 std::vector<Scalar>& R_sum,
1091 std::vector<Scalar>& maxCoeff,
1092 std::vector<int>& maxCoeffCell)
1093 {
1094 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1095 {
1096 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1097 continue;
1098 }
1099
1100 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1101
1102 B_avg[compIdx] += 1.0 / fs.invB(phaseIdx).value();
1103 const auto R2 = ebosResid[cell_idx][compIdx];
1104
1105 R_sum[compIdx] += R2;
1106 const double Rval = std::abs(R2) / pvValue;
1107 if (Rval > maxCoeff[compIdx]) {
1108 maxCoeff[compIdx] = Rval;
1109 maxCoeffCell[compIdx] = cell_idx;
1110 }
1111 }
1112
1113 if constexpr (has_solvent_) {
1114 B_avg[contiSolventEqIdx] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
1115 const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
1116 R_sum[contiSolventEqIdx] += R2;
1117 maxCoeff[contiSolventEqIdx] = std::max(maxCoeff[contiSolventEqIdx],
1118 std::abs(R2) / pvValue);
1119 }
1120 if constexpr (has_extbo_) {
1121 B_avg[contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1122 const auto R2 = ebosResid[cell_idx][contiZfracEqIdx];
1123 R_sum[ contiZfracEqIdx ] += R2;
1124 maxCoeff[contiZfracEqIdx] = std::max(maxCoeff[contiZfracEqIdx],
1125 std::abs(R2) / pvValue);
1126 }
1127 if constexpr (has_polymer_) {
1128 B_avg[contiPolymerEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1129 const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
1130 R_sum[contiPolymerEqIdx] += R2;
1131 maxCoeff[contiPolymerEqIdx] = std::max(maxCoeff[contiPolymerEqIdx],
1132 std::abs(R2) / pvValue);
1133 }
1134 if constexpr (has_foam_) {
1135 B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1136 const auto R2 = ebosResid[cell_idx][contiFoamEqIdx];
1137 R_sum[contiFoamEqIdx] += R2;
1138 maxCoeff[contiFoamEqIdx] = std::max(maxCoeff[contiFoamEqIdx],
1139 std::abs(R2) / pvValue);
1140 }
1141 if constexpr (has_brine_) {
1142 B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1143 const auto R2 = ebosResid[cell_idx][contiBrineEqIdx];
1144 R_sum[contiBrineEqIdx] += R2;
1145 maxCoeff[contiBrineEqIdx] = std::max(maxCoeff[contiBrineEqIdx],
1146 std::abs(R2) / pvValue);
1147 }
1148
1149 if constexpr (has_polymermw_) {
1150 static_assert(has_polymer_);
1151
1152 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1153 // the residual of the polymer molecular equation is scaled down by a 100, since molecular weight
1154 // can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations
1155 // TODO: there should be a more general way to determine the scaling-down coefficient
1156 const auto R2 = ebosResid[cell_idx][contiPolymerMWEqIdx] / 100.;
1157 R_sum[contiPolymerMWEqIdx] += R2;
1158 maxCoeff[contiPolymerMWEqIdx] = std::max(maxCoeff[contiPolymerMWEqIdx],
1159 std::abs(R2) / pvValue);
1160 }
1161
1162 if constexpr (has_energy_) {
1163 B_avg[contiEnergyEqIdx] += 1.0 / (4.182e1); // converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water
1164 const auto R2 = ebosResid[cell_idx][contiEnergyEqIdx];
1165 R_sum[contiEnergyEqIdx] += R2;
1166 maxCoeff[contiEnergyEqIdx] = std::max(maxCoeff[contiEnergyEqIdx],
1167 std::abs(R2) / pvValue);
1168 }
1169
1170 if constexpr (has_micp_) {
1171 B_avg[contiMicrobialEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1172 const auto R1 = ebosResid[cell_idx][contiMicrobialEqIdx];
1173 R_sum[contiMicrobialEqIdx] += R1;
1174 maxCoeff[contiMicrobialEqIdx] = std::max(maxCoeff[contiMicrobialEqIdx],
1175 std::abs(R1) / pvValue);
1176 B_avg[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1177 const auto R2 = ebosResid[cell_idx][contiOxygenEqIdx];
1178 R_sum[contiOxygenEqIdx] += R2;
1179 maxCoeff[contiOxygenEqIdx] = std::max(maxCoeff[contiOxygenEqIdx],
1180 std::abs(R2) / pvValue);
1181 B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1182 const auto R3 = ebosResid[cell_idx][contiUreaEqIdx];
1183 R_sum[contiUreaEqIdx] += R3;
1184 maxCoeff[contiUreaEqIdx] = std::max(maxCoeff[contiUreaEqIdx],
1185 std::abs(R3) / pvValue);
1186 B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1187 const auto R4 = ebosResid[cell_idx][contiBiofilmEqIdx];
1188 R_sum[contiBiofilmEqIdx] += R4;
1189 maxCoeff[contiBiofilmEqIdx] = std::max(maxCoeff[contiBiofilmEqIdx],
1190 std::abs(R4) / pvValue);
1191 B_avg[contiCalciteEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1192 const auto R5 = ebosResid[cell_idx][contiCalciteEqIdx];
1193 R_sum[contiCalciteEqIdx] += R5;
1194 maxCoeff[contiCalciteEqIdx] = std::max(maxCoeff[contiCalciteEqIdx],
1195 std::abs(R5) / pvValue);
1196 }
1197 }
1198
1200 const ModelParameters& param() const
1201 {
1202 return param_;
1203 }
1204
1207 {
1208 return compNames_;
1209 }
1210
1211 private:
1212 double dpMaxRel() const { return param_.dp_max_rel_; }
1213 double dsMax() const { return param_.ds_max_; }
1214 double drMaxRel() const { return param_.dr_max_rel_; }
1215 double maxResidualAllowed() const { return param_.max_residual_allowed_; }
1216 double linear_solve_setup_time_;
1217
1218 public:
1219 std::vector<bool> wasSwitched_;
1220 };
1221} // namespace Opm
1222
1223#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
Definition AquiferInterface.hpp:35
A model implementation for three-phase black oil.
Definition BlackoilModelEbos.hpp:158
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition BlackoilModelEbos.hpp:1068
std::vector< std::vector< double > > computeFluidInPlace(const std::vector< int > &) const
Should not be called.
Definition BlackoilModelEbos.hpp:999
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition BlackoilModelEbos.hpp:577
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition BlackoilModelEbos.hpp:678
const SimulatorReportSingle & failureReport() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModelEbos.hpp:1015
ConvergenceReport 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:965
std::unique_ptr< BlackoilModelEbosNldd< TypeTag > > nlddSolver_
Non-linear DD solver.
Definition BlackoilModelEbos.hpp:1063
BlackoilModelEbos(Simulator &ebosSimulator, const ModelParameters &param, BlackoilWellModel< TypeTag > &well_model, const bool terminal_output)
Construct the model.
Definition BlackoilModelEbos.hpp:221
bool terminal_output_
Whether we print something to std::cout.
Definition BlackoilModelEbos.hpp:1052
SimulatorReportSingle afterStep(const SimulatorTimerInterface &)
Called once after each time step.
Definition BlackoilModelEbos.hpp:470
std::pair< double, double > localConvergenceData(std::vector< Scalar > &R_sum, std::vector< Scalar > &maxCoeff, std::vector< Scalar > &B_avg, std::vector< int > &maxCoeffCell)
Get reservoir quantities on this process needed for convergence calculations.
Definition BlackoilModelEbos.hpp:746
const ComponentName & compNames() const
Returns const reference to component names.
Definition BlackoilModelEbos.hpp:1206
SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Called once per nonlinear iteration.
Definition BlackoilModelEbos.hpp:358
SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &, const int iterationIdx)
Assemble the residual and Jacobian of the nonlinear system.
Definition BlackoilModelEbos.hpp:481
SimulatorReportSingle prepareStep(const SimulatorTimerInterface &timer)
Called once before each time step.
Definition BlackoilModelEbos.hpp:263
long int global_nc_
The number of cells of the global grid.
Definition BlackoilModelEbos.hpp:1054
double computeCnvErrorPv(const std::vector< Scalar > &B_avg, double dt)
Compute the total pore volume of cells violating CNV that are not part of a numerical aquifer.
Definition BlackoilModelEbos.hpp:799
int numPhases() const
The number of active fluid phases in the model.
Definition BlackoilModelEbos.hpp:984
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:992
SimulatorReportSingle localAccumulatedReports() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModelEbos.hpp:1019
const ModelParameters & param() const
Returns const reference to model parameters.
Definition BlackoilModelEbos.hpp:1200
void solveJacobianSystem(BVector &x)
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition BlackoilModelEbos.hpp:592
void updateSolution(const BVector &dx)
Apply an update to the primary variables.
Definition BlackoilModelEbos.hpp:656
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
Interface class for SimulatorTimer objects, to be improved.
Definition SimulatorTimerInterface.hpp:34
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition phaseUsageFromDeck.cpp:145
Definition AquiferGridUtils.hpp:35
Definition ISTLSolverEbos.hpp:69
Definition BlackoilModelEbos.hpp:74
Definition ISTLSolverEbos.hpp:63
Definition BlackoilModelParametersEbos.hpp:37
Definition NonlinearSolverEbos.hpp:41
Definition AdaptiveTimeSteppingEbos.hpp:51
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34