My Project
Loading...
Searching...
No Matches
BlackoilModel.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_BLACKOILMODEL_HEADER_INCLUDED
25#define OPM_BLACKOILMODEL_HEADER_INCLUDED
26
27#include <fmt/format.h>
28
29#include <opm/common/ErrorMacros.hpp>
30#include <opm/common/Exceptions.hpp>
31#include <opm/common/OpmLog/OpmLog.hpp>
32
33#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
34#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
35
36#include <opm/simulators/aquifers/AquiferGridUtils.hpp>
37#include <opm/simulators/aquifers/BlackoilAquiferModel.hpp>
38#include <opm/simulators/flow/BlackoilModelNldd.hpp>
39#include <opm/simulators/flow/BlackoilModelParameters.hpp>
40#include <opm/simulators/flow/countGlobalCells.hpp>
42#include <opm/simulators/flow/NonlinearSolver.hpp>
43#include <opm/simulators/flow/RSTConv.hpp>
44#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
45#include <opm/simulators/timestepping/ConvergenceReport.hpp>
46#include <opm/simulators/timestepping/SimulatorReport.hpp>
47#include <opm/simulators/timestepping/SimulatorTimer.hpp>
48
49#include <opm/simulators/wells/BlackoilWellModel.hpp>
50
51#include <opm/simulators/utils/ComponentName.hpp>
52#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
53#include <opm/simulators/utils/ParallelCommunication.hpp>
54#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
55
56#include <dune/common/timer.hh>
57
58#include <fmt/format.h>
59
60#include <algorithm>
61#include <cassert>
62#include <cmath>
63#include <filesystem>
64#include <fstream>
65#include <functional>
66#include <iomanip>
67#include <ios>
68#include <limits>
69#include <memory>
70#include <numeric>
71#include <sstream>
72#include <tuple>
73#include <utility>
74#include <vector>
75
76namespace Opm::Properties {
77
78namespace TTag {
79
80struct FlowProblem { using InheritsFrom = std::tuple<FlowBaseProblemBlackoil, BlackOilModel>; };
81
82}
83
84// default in flow is to formulate the equations in surface volumes
85template<class TypeTag>
87{ static constexpr bool value = true; };
88
89template<class TypeTag>
90struct UseVolumetricResidual<TypeTag, TTag::FlowProblem>
91{ static constexpr bool value = false; };
92
93template<class TypeTag>
94struct AquiferModel<TypeTag, TTag::FlowProblem>
96
97// disable all extensions supported by black oil model. this should not really be
98// necessary but it makes things a bit more explicit
99template<class TypeTag>
100struct EnablePolymer<TypeTag, TTag::FlowProblem>
101{ static constexpr bool value = false; };
102
103template<class TypeTag>
104struct EnableSolvent<TypeTag, TTag::FlowProblem>
105{ static constexpr bool value = false; };
106
107template<class TypeTag>
108struct EnableTemperature<TypeTag, TTag::FlowProblem>
109{ static constexpr bool value = true; };
110
111template<class TypeTag>
112struct EnableEnergy<TypeTag, TTag::FlowProblem>
113{ static constexpr bool value = false; };
114
115template<class TypeTag>
116struct EnableFoam<TypeTag, TTag::FlowProblem>
117{ static constexpr bool value = false; };
118
119template<class TypeTag>
120struct EnableBrine<TypeTag, TTag::FlowProblem>
121{ static constexpr bool value = false; };
122
123template<class TypeTag>
125{ static constexpr bool value = false; };
126
127template<class TypeTag>
128struct EnableMICP<TypeTag, TTag::FlowProblem>
129{ static constexpr bool value = false; };
130
131template<class TypeTag>
132struct EnableDispersion<TypeTag, TTag::FlowProblem>
133{ static constexpr bool value = false; };
134
135template<class TypeTag>
137{ static constexpr bool value = true; };
138
139template<class TypeTag>
140struct WellModel<TypeTag, TTag::FlowProblem>
142
143template<class TypeTag>
144struct LinearSolverSplice<TypeTag, TTag::FlowProblem>
145{ using type = TTag::FlowIstlSolver; };
146
147template<class TypeTag>
148struct EnableDebuggingChecks<TypeTag, TTag::FlowProblem>
149{ static constexpr bool value = false; };
150
151} // namespace Opm::Properties
152
153namespace Opm {
154
161 template <class TypeTag>
162 class BlackoilModel
163 {
164 public:
165 // --------- Types and enums ---------
178 using ModelParameters = BlackoilModelParameters<Scalar>;
179
180 static constexpr int numEq = Indices::numEq;
181 static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
182 static constexpr int contiZfracEqIdx = Indices::contiZfracEqIdx;
183 static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
184 static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
185 static constexpr int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
186 static constexpr int contiFoamEqIdx = Indices::contiFoamEqIdx;
187 static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx;
188 static constexpr int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
189 static constexpr int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
190 static constexpr int contiUreaEqIdx = Indices::contiUreaEqIdx;
191 static constexpr int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
192 static constexpr int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
193 static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
194 static constexpr int zFractionIdx = Indices::zFractionIdx;
195 static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
196 static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
197 static constexpr int temperatureIdx = Indices::temperatureIdx;
198 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
199 static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx;
200 static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
201 static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
202 static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
203 static constexpr int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
204 static constexpr int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
205
206 using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
207 using MatrixBlockType = typename SparseMatrixAdapter::MatrixBlock;
208 using Mat = typename SparseMatrixAdapter::IstlMatrix;
209 using BVector = Dune::BlockVector<VectorBlockType>;
210
212
213 // --------- Public methods ---------
214
225 BlackoilModel(Simulator& simulator,
226 const ModelParameters& param,
227 BlackoilWellModel<TypeTag>& well_model,
228 const bool terminal_output)
229 : simulator_(simulator)
230 , grid_(simulator_.vanguard().grid())
231 , phaseUsage_(phaseUsageFromDeck(eclState()))
232 , param_( param )
233 , well_model_ (well_model)
234 , rst_conv_(simulator_.problem().eclWriter()->collectOnIORank().localIdxToGlobalIdxMapping(),
235 grid_.comm())
236 , terminal_output_ (terminal_output)
237 , current_relaxation_(1.0)
238 , dx_old_(simulator_.model().numGridDof())
239 {
240 // compute global sum of number of cells
241 global_nc_ = detail::countGlobalCells(grid_);
242 convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves.
243 // TODO: remember to fix!
244 if (param_.nonlinear_solver_ == "nldd") {
245 if (terminal_output) {
246 OpmLog::info("Using Non-Linear Domain Decomposition solver (nldd).");
247 }
248 nlddSolver_ = std::make_unique<BlackoilModelNldd<TypeTag>>(*this);
249 } else if (param_.nonlinear_solver_ == "newton") {
250 if (terminal_output) {
251 OpmLog::info("Using Newton nonlinear solver.");
252 }
253 } else {
254 OPM_THROW(std::runtime_error, "Unknown nonlinear solver option: " + param_.nonlinear_solver_);
255 }
256 }
257
258
259 bool isParallel() const
260 { return grid_.comm().size() > 1; }
261
262
263 const EclipseState& eclState() const
264 { return simulator_.vanguard().eclState(); }
265
266
270 {
272 Dune::Timer perfTimer;
273 perfTimer.start();
274 // update the solution variables in the model
275 if ( timer.lastStepFailed() ) {
276 simulator_.model().updateFailed();
277 } else {
278 simulator_.model().advanceTimeLevel();
279 }
280
281 // Set the timestep size, episode index, and non-linear iteration index
282 // for the model explicitly. The model needs to know the report step/episode index
283 // because of timing dependent data despite the fact that flow uses its
284 // own time stepper. (The length of the episode does not matter, though.)
285 simulator_.setTime(timer.simulationTimeElapsed());
286 simulator_.setTimeStepSize(timer.currentStepLength());
287 simulator_.model().newtonMethod().setIterationIndex(0);
288
289 simulator_.problem().beginTimeStep();
290
291 unsigned numDof = simulator_.model().numGridDof();
292 wasSwitched_.resize(numDof);
293 std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
294
295 if (param_.update_equations_scaling_) {
296 OpmLog::error("Equation scaling not supported");
297 //updateEquationsScaling();
298 }
299
300 if (nlddSolver_) {
301 nlddSolver_->prepareStep();
302 }
303
304 report.pre_post_time += perfTimer.stop();
305
306 auto getIdx = [](unsigned phaseIdx) -> int
307 {
308 if (FluidSystem::phaseIsActive(phaseIdx)) {
309 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
310 return Indices::canonicalToActiveComponentIndex(sIdx);
311 }
312
313 return -1;
314 };
315 const auto& schedule = simulator_.vanguard().schedule();
316 rst_conv_.init(simulator_.vanguard().globalNumCells(),
317 schedule[timer.reportStepNum()].rst_config(),
318 {getIdx(FluidSystem::oilPhaseIdx),
319 getIdx(FluidSystem::gasPhaseIdx),
320 getIdx(FluidSystem::waterPhaseIdx),
321 contiPolymerEqIdx,
322 contiBrineEqIdx,
323 contiSolventEqIdx});
324
325 return report;
326 }
327
328
329 void initialLinearization(SimulatorReportSingle& report,
330 const int iteration,
331 const int minIter,
332 const int maxIter,
333 const SimulatorTimerInterface& timer)
334 {
335 // ----------- Set up reports and timer -----------
336 failureReport_ = SimulatorReportSingle();
337 Dune::Timer perfTimer;
338
339 perfTimer.start();
340 report.total_linearizations = 1;
341
342 // ----------- Assemble -----------
343 try {
344 report += assembleReservoir(timer, iteration);
345 report.assemble_time += perfTimer.stop();
346 }
347 catch (...) {
348 report.assemble_time += perfTimer.stop();
349 failureReport_ += report;
350 throw; // continue throwing the stick
351 }
352
353 // ----------- Check if converged -----------
354 std::vector<Scalar> residual_norms;
355 perfTimer.reset();
356 perfTimer.start();
357 // the step is not considered converged until at least minIter iterations is done
358 {
359 auto convrep = getConvergence(timer, iteration, maxIter, residual_norms);
360 report.converged = convrep.converged() && iteration >= minIter;
361 ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
362 convergence_reports_.back().report.push_back(std::move(convrep));
363
364 // Throw if any NaN or too large residual found.
365 if (severity == ConvergenceReport::Severity::NotANumber) {
366 failureReport_ += report;
367 OPM_THROW_PROBLEM(NumericalProblem, "NaN residual found!");
368 } else if (severity == ConvergenceReport::Severity::TooLarge) {
369 failureReport_ += report;
370 OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!");
371 } else if (severity == ConvergenceReport::Severity::ConvergenceMonitorFailure) {
372 failureReport_ += report;
373 OPM_THROW_PROBLEM(ConvergenceMonitorFailure, "Total penalty count exceeded cut-off-limit of " + std::to_string(param_.convergence_monitoring_cutoff_));
374 }
375 }
376 report.update_time += perfTimer.stop();
377 residual_norms_history_.push_back(residual_norms);
378 }
379
380
388 template <class NonlinearSolverType>
390 const SimulatorTimerInterface& timer,
391 NonlinearSolverType& nonlinear_solver)
392 {
393 if (iteration == 0) {
394 // For each iteration we store in a vector the norms of the residual of
395 // the mass balance for each active phase, the well flux and the well equations.
396 residual_norms_history_.clear();
397 total_penaltyCard_.reset();
398 prev_above_tolerance_ = 0;
399 prev_distance_ = std::numeric_limits<double>::infinity();
400 current_relaxation_ = 1.0;
401 dx_old_ = 0.0;
402 convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}});
403 convergence_reports_.back().report.reserve(11);
404 }
405
407 if ((this->param_.nonlinear_solver_ != "nldd") ||
408 (iteration < this->param_.nldd_num_initial_newton_iter_))
409 {
410 result = this->nonlinearIterationNewton(iteration, timer, nonlinear_solver);
411 }
412 else {
413 result = this->nlddSolver_->nonlinearIterationNldd(iteration, timer, nonlinear_solver);
414 }
415
416 rst_conv_.update(simulator_.model().linearizer().residual());
417
418 return result;
419 }
420
421
422 template <class NonlinearSolverType>
423 SimulatorReportSingle nonlinearIterationNewton(const int iteration,
424 const SimulatorTimerInterface& timer,
425 NonlinearSolverType& nonlinear_solver)
426 {
427
428 // ----------- Set up reports and timer -----------
430 Dune::Timer perfTimer;
431
432 this->initialLinearization(report, iteration, nonlinear_solver.minIter(), nonlinear_solver.maxIter(), timer);
433
434 // ----------- If not converged, solve linear system and do Newton update -----------
435 if (!report.converged) {
436 perfTimer.reset();
437 perfTimer.start();
438 report.total_newton_iterations = 1;
439
440 // Compute the nonlinear update.
441 unsigned nc = simulator_.model().numGridDof();
442 BVector x(nc);
443
444 // Solve the linear system.
445 linear_solve_setup_time_ = 0.0;
446 try {
447 // Apply the Schur complement of the well model to
448 // the reservoir linearized equations.
449 // Note that linearize may throw for MSwells.
450 wellModel().linearize(simulator().model().linearizer().jacobian(),
451 simulator().model().linearizer().residual());
452
453 // ---- Solve linear system ----
455
456 report.linear_solve_setup_time += linear_solve_setup_time_;
457 report.linear_solve_time += perfTimer.stop();
458 report.total_linear_iterations += linearIterationsLastSolve();
459 }
460 catch (...) {
461 report.linear_solve_setup_time += linear_solve_setup_time_;
462 report.linear_solve_time += perfTimer.stop();
463 report.total_linear_iterations += linearIterationsLastSolve();
464
465 failureReport_ += report;
466 throw; // re-throw up
467 }
468
469 perfTimer.reset();
470 perfTimer.start();
471
472 // handling well state update before oscillation treatment is a decision based
473 // on observation to avoid some big performance degeneration under some circumstances.
474 // there is no theorectical explanation which way is better for sure.
475 wellModel().postSolve(x);
476
477 if (param_.use_update_stabilization_) {
478 // Stabilize the nonlinear update.
479 bool isOscillate = false;
480 bool isStagnate = false;
481 nonlinear_solver.detectOscillations(residual_norms_history_, residual_norms_history_.size() - 1, isOscillate, isStagnate);
482 if (isOscillate) {
483 current_relaxation_ -= nonlinear_solver.relaxIncrement();
484 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
485 if (terminalOutputEnabled()) {
486 std::string msg = " Oscillating behavior detected: Relaxation set to "
487 + std::to_string(current_relaxation_);
488 OpmLog::info(msg);
489 }
490 }
491 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
492 }
493
494 // ---- Newton update ----
495 // Apply the update, with considering model-dependent limitations and
496 // chopping of the update.
498
499 report.update_time += perfTimer.stop();
500 }
501
502 return report;
503 }
504
505
510 {
512 Dune::Timer perfTimer;
513 perfTimer.start();
514 simulator_.problem().endTimeStep();
515 simulator_.problem().setConvData(rst_conv_.getData());
516 report.pre_post_time += perfTimer.stop();
517 return report;
518 }
519
522 const int iterationIdx)
523 {
524 // -------- Mass balance equations --------
525 simulator_.model().newtonMethod().setIterationIndex(iterationIdx);
526 simulator_.problem().beginIteration();
527 simulator_.model().linearizer().linearizeDomain();
528 simulator_.problem().endIteration();
529 return wellModel().lastReport();
530 }
531
532 // compute the "relative" change of the solution between time steps
533 Scalar relativeChange() const
534 {
535 Scalar resultDelta = 0.0;
536 Scalar resultDenom = 0.0;
537
538 const auto& elemMapper = simulator_.model().elementMapper();
539 const auto& gridView = simulator_.gridView();
540 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
541 unsigned globalElemIdx = elemMapper.index(elem);
542 const auto& priVarsNew = simulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
543
544 Scalar pressureNew;
545 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
546
547 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
548 Scalar oilSaturationNew = 1.0;
549 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
550 FluidSystem::numActivePhases() > 1 &&
551 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
552 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
553 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
554 }
555
556 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
557 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
558 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
559 assert(Indices::compositionSwitchIdx >= 0 );
560 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
561 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
562 }
563
564 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
565 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
566 }
567
568 const auto& priVarsOld = simulator_.model().solution(/*timeIdx=*/1)[globalElemIdx];
569
570 Scalar pressureOld;
571 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
572
573 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
574 Scalar oilSaturationOld = 1.0;
575
576 // NB fix me! adding pressures changes to satutation changes does not make sense
577 Scalar tmp = pressureNew - pressureOld;
578 resultDelta += tmp*tmp;
579 resultDenom += pressureNew*pressureNew;
580
581 if (FluidSystem::numActivePhases() > 1) {
582 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
583 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSwitchIdx];
584 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
585 }
586
587 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
588 {
589 assert(Indices::compositionSwitchIdx >= 0 );
590 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
591 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
592 }
593
594 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
595 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
596 }
597 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
598 Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
599 resultDelta += tmpSat*tmpSat;
600 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
601 assert(std::isfinite(resultDelta));
602 assert(std::isfinite(resultDenom));
603 }
604 }
605 }
606
607 resultDelta = gridView.comm().sum(resultDelta);
608 resultDenom = gridView.comm().sum(resultDenom);
609
610 if (resultDenom > 0.0)
611 return resultDelta/resultDenom;
612 return 0.0;
613 }
614
615
618 {
619 return simulator_.model().newtonMethod().linearSolver().iterations ();
620 }
621
622
623 // Obtain reference to linear solver setup time
624 double& linearSolveSetupTime()
625 {
626 return linear_solve_setup_time_;
627 }
628
629
632 void solveJacobianSystem(BVector& x)
633 {
634 auto& jacobian = simulator_.model().linearizer().jacobian().istlMatrix();
635 auto& residual = simulator_.model().linearizer().residual();
636 auto& linSolver = simulator_.model().newtonMethod().linearSolver();
637
638 const int numSolvers = linSolver.numAvailableSolvers();
639 if ((numSolvers > 1) && (linSolver.getSolveCount() % 100 == 0)) {
640
641 if ( terminal_output_ ) {
642 OpmLog::debug("\nRunning speed test for comparing available linear solvers.");
643 }
644
645 Dune::Timer perfTimer;
646 std::vector<double> times(numSolvers);
647 std::vector<double> setupTimes(numSolvers);
648
649 x = 0.0;
650 std::vector<BVector> x_trial(numSolvers, x);
651 for (int solver = 0; solver < numSolvers; ++solver) {
652 BVector x0(x);
653 linSolver.setActiveSolver(solver);
654 perfTimer.start();
655 linSolver.prepare(jacobian, residual);
656 setupTimes[solver] = perfTimer.stop();
657 perfTimer.reset();
658 linSolver.setResidual(residual);
659 perfTimer.start();
660 linSolver.solve(x_trial[solver]);
661 times[solver] = perfTimer.stop();
662 perfTimer.reset();
663 if (terminal_output_) {
664 OpmLog::debug(fmt::format("Solver time {}: {}", solver, times[solver]));
665 }
666 }
667
668 int fastest_solver = std::min_element(times.begin(), times.end()) - times.begin();
669 // Use timing on rank 0 to determine fastest, must be consistent across ranks.
670 grid_.comm().broadcast(&fastest_solver, 1, 0);
671 linear_solve_setup_time_ = setupTimes[fastest_solver];
672 x = x_trial[fastest_solver];
673 linSolver.setActiveSolver(fastest_solver);
674 } else {
675 // set initial guess
676 x = 0.0;
677
678 Dune::Timer perfTimer;
679 perfTimer.start();
680 linSolver.prepare(jacobian, residual);
681 linear_solve_setup_time_ = perfTimer.stop();
682 linSolver.setResidual(residual);
683 // actually, the error needs to be calculated after setResidual in order to
684 // account for parallelization properly. since the residual of ECFV
685 // discretizations does not need to be synchronized across processes to be
686 // consistent, this is not relevant for OPM-flow...
687 linSolver.solve(x);
688 }
689 }
690
691
693 void updateSolution(const BVector& dx)
694 {
695 OPM_TIMEBLOCK(updateSolution);
696 auto& newtonMethod = simulator_.model().newtonMethod();
697 SolutionVector& solution = simulator_.model().solution(/*timeIdx=*/0);
698
699 newtonMethod.update_(/*nextSolution=*/solution,
700 /*curSolution=*/solution,
701 /*update=*/dx,
702 /*resid=*/dx); // the update routines of the black
703 // oil model do not care about the
704 // residual
705
706 // if the solution is updated, the intensive quantities need to be recalculated
707 {
708 OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities);
709 simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
710 }
711 }
712
715 {
716 return terminal_output_;
717 }
718
719 std::tuple<Scalar,Scalar> convergenceReduction(Parallel::Communication comm,
720 const Scalar pvSumLocal,
721 const Scalar numAquiferPvSumLocal,
722 std::vector< Scalar >& R_sum,
723 std::vector< Scalar >& maxCoeff,
724 std::vector< Scalar >& B_avg)
725 {
726 OPM_TIMEBLOCK(convergenceReduction);
727 // Compute total pore volume (use only owned entries)
728 Scalar pvSum = pvSumLocal;
729 Scalar numAquiferPvSum = numAquiferPvSumLocal;
730
731 if( comm.size() > 1 )
732 {
733 // global reduction
734 std::vector< Scalar > sumBuffer;
735 std::vector< Scalar > maxBuffer;
736 const int numComp = B_avg.size();
737 sumBuffer.reserve( 2*numComp + 2 ); // +2 for (numAquifer)pvSum
738 maxBuffer.reserve( numComp );
739 for( int compIdx = 0; compIdx < numComp; ++compIdx )
740 {
741 sumBuffer.push_back( B_avg[ compIdx ] );
742 sumBuffer.push_back( R_sum[ compIdx ] );
743 maxBuffer.push_back( maxCoeff[ compIdx ] );
744 }
745
746 // Compute total pore volume
747 sumBuffer.push_back( pvSum );
748 sumBuffer.push_back( numAquiferPvSum );
749
750 // compute global sum
751 comm.sum( sumBuffer.data(), sumBuffer.size() );
752
753 // compute global max
754 comm.max( maxBuffer.data(), maxBuffer.size() );
755
756 // restore values to local variables
757 for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
758 {
759 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
760 ++buffIdx;
761
762 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
763 }
764
765 for( int compIdx = 0; compIdx < numComp; ++compIdx )
766 {
767 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
768 }
769
770 // restore global pore volume
771 pvSum = sumBuffer[sumBuffer.size()-2];
772 numAquiferPvSum = sumBuffer.back();
773 }
774
775 // return global pore volume
776 return {pvSum, numAquiferPvSum};
777 }
778
782 std::pair<Scalar,Scalar> localConvergenceData(std::vector<Scalar>& R_sum,
783 std::vector<Scalar>& maxCoeff,
784 std::vector<Scalar>& B_avg,
785 std::vector<int>& maxCoeffCell)
786 {
787 OPM_TIMEBLOCK(localConvergenceData);
788 Scalar pvSumLocal = 0.0;
789 Scalar numAquiferPvSumLocal = 0.0;
790 const auto& model = simulator_.model();
791 const auto& problem = simulator_.problem();
792
793 const auto& residual = simulator_.model().linearizer().residual();
794
795 ElementContext elemCtx(simulator_);
796 const auto& gridView = simulator().gridView();
797 IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
798 OPM_BEGIN_PARALLEL_TRY_CATCH();
799 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
800 elemCtx.updatePrimaryStencil(elem);
801 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
802
803 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
804 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
805 const auto& fs = intQuants.fluidState();
806
807 const auto pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) *
808 model.dofTotalVolume(cell_idx);
809 pvSumLocal += pvValue;
810
811 if (isNumericalAquiferCell(elem))
812 {
813 numAquiferPvSumLocal += pvValue;
814 }
815
816 this->getMaxCoeff(cell_idx, intQuants, fs, residual, pvValue,
817 B_avg, R_sum, maxCoeff, maxCoeffCell);
818 }
819
820 OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::localConvergenceData() failed: ", grid_.comm());
821
822 // compute local average in terms of global number of elements
823 const int bSize = B_avg.size();
824 for ( int i = 0; i<bSize; ++i )
825 {
826 B_avg[ i ] /= Scalar( global_nc_ );
827 }
828
829 return {pvSumLocal, numAquiferPvSumLocal};
830 }
831
832
836 std::pair<std::vector<double>, std::vector<int>>
837 characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt)
838 {
839 OPM_TIMEBLOCK(computeCnvErrorPv);
840
841 // 0: cnv <= tolerance_cnv
842 // 1: tolerance_cnv < cnv <= tolerance_cnv_relaxed
843 // 2: tolerance_cnv_relaxed < cnv
844 constexpr auto numPvGroups = std::vector<double>::size_type{3};
845
846 auto cnvPvSplit = std::pair<std::vector<double>, std::vector<int>> {
847 std::piecewise_construct,
848 std::forward_as_tuple(numPvGroups),
849 std::forward_as_tuple(numPvGroups)
850 };
851
852 auto maxCNV = [&B_avg, dt](const auto& residual, const double pvol)
853 {
854 return (dt / pvol) *
855 std::inner_product(residual.begin(), residual.end(),
856 B_avg.begin(), Scalar{0},
857 [](const Scalar m, const auto& x)
858 {
859 using std::abs;
860 return std::max(m, abs(x));
861 }, std::multiplies<>{});
862 };
863
864 auto& [splitPV, cellCntPV] = cnvPvSplit;
865
866 const auto& model = this->simulator().model();
867 const auto& problem = this->simulator().problem();
868 const auto& residual = model.linearizer().residual();
869 const auto& gridView = this->simulator().gridView();
870
871 const IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
872
873 ElementContext elemCtx(this->simulator());
874
875 OPM_BEGIN_PARALLEL_TRY_CATCH();
876 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
877 // Skip cells of numerical Aquifer
878 if (isNumericalAquiferCell(elem)) {
879 continue;
880 }
881
882 elemCtx.updatePrimaryStencil(elem);
883
884 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
885 const auto pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0)
886 * model.dofTotalVolume(cell_idx);
887
888 const auto maxCnv = maxCNV(residual[cell_idx], pvValue);
889
890 const auto ix = (maxCnv > this->param_.tolerance_cnv_)
891 + (maxCnv > this->param_.tolerance_cnv_relaxed_);
892
893 splitPV[ix] += static_cast<double>(pvValue);
894 ++cellCntPV[ix];
895 }
896
897 OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::characteriseCnvPvSplit() failed: ",
898 this->grid_.comm());
899
900 this->grid_.comm().sum(splitPV .data(), splitPV .size());
901 this->grid_.comm().sum(cellCntPV.data(), cellCntPV.size());
902
903 return cnvPvSplit;
904 }
905
906
907 void updateTUNING(const Tuning& tuning)
908 {
909 this->param_.tolerance_mb_ = tuning.XXXMBE;
910
911 if (terminal_output_) {
912 OpmLog::debug(fmt::format("Setting BlackoilModel mass "
913 "balance limit (XXXMBE) to {:.2e}",
914 tuning.XXXMBE));
915 }
916 }
917
918
919 ConvergenceReport getReservoirConvergence(const double reportTime,
920 const double dt,
921 const int iteration,
922 const int maxIter,
923 std::vector<Scalar>& B_avg,
924 std::vector<Scalar>& residual_norms)
925 {
926 OPM_TIMEBLOCK(getReservoirConvergence);
927 using Vector = std::vector<Scalar>;
928
929 ConvergenceReport report{reportTime};
930
931 const int numComp = numEq;
932
933 Vector R_sum(numComp, Scalar{0});
934 Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest());
935 std::vector<int> maxCoeffCell(numComp, -1);
936
937 const auto [pvSumLocal, numAquiferPvSumLocal] =
938 this->localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell);
939
940 // compute global sum and max of quantities
941 const auto& [pvSum, numAquiferPvSum] =
942 this->convergenceReduction(this->grid_.comm(),
943 pvSumLocal,
944 numAquiferPvSumLocal,
945 R_sum, maxCoeff, B_avg);
946
947 report.setCnvPoreVolSplit(this->characteriseCnvPvSplit(B_avg, dt),
948 pvSum - numAquiferPvSum);
949
950 // For each iteration, we need to determine whether to use the
951 // relaxed tolerances. To disable the usage of relaxed
952 // tolerances, you can set the relaxed tolerances as the strict
953 // tolerances. If min_strict_mb_iter = -1 (default) we use a
954 // relaxed tolerance for the mass balance for the last
955 // iterations. For positive values we use the relaxed tolerance
956 // after the given number of iterations
957 const bool relax_final_iteration_mb =
958 (this->param_.min_strict_mb_iter_ < 0)
959 && (iteration == maxIter);
960
961 const bool use_relaxed_mb = relax_final_iteration_mb ||
962 ((this->param_.min_strict_mb_iter_ >= 0) &&
963 (iteration >= this->param_.min_strict_mb_iter_));
964
965 // If min_strict_cnv_iter = -1 we use a relaxed tolerance for
966 // the cnv for the last iterations. For positive values we use
967 // the relaxed tolerance after the given number of iterations.
968 // We also use relaxed tolerances for cells with total
969 // pore-volume less than relaxed_max_pv_fraction_. Default
970 // value of relaxed_max_pv_fraction_ is 0.03
971 const bool relax_final_iteration_cnv =
972 (this->param_.min_strict_cnv_iter_ < 0)
973 && (iteration == maxIter);
974
975 const bool relax_iter_cnv = (this->param_.min_strict_cnv_iter_ >= 0)
976 && (iteration >= this->param_.min_strict_cnv_iter_);
977
978 // Note trailing parentheses here, just before the final
979 // semicolon. This is an immediately invoked function
980 // expression which calculates a single boolean value.
981 const auto relax_pv_fraction_cnv =
982 [&report, this, eligible = pvSum - numAquiferPvSum]()
983 {
984 const auto& cnvPvSplit = report.cnvPvSplit().first;
985
986 // [1]: tol < cnv <= relaxed
987 // [2]: relaxed < cnv
988 return static_cast<Scalar>(cnvPvSplit[1] + cnvPvSplit[2]) <
989 this->param_.relaxed_max_pv_fraction_ * eligible;
990 }();
991
992 const bool use_relaxed_cnv = relax_final_iteration_cnv
993 || relax_pv_fraction_cnv
994 || relax_iter_cnv;
995
996 if ((relax_final_iteration_mb || relax_final_iteration_cnv) &&
997 this->terminal_output_)
998 {
999 std::string message =
1000 "Number of newton iterations reached its maximum "
1001 "try to continue with relaxed tolerances:";
1002
1003 if (relax_final_iteration_mb) {
1004 message += fmt::format(" MB: {:.1e}", param_.tolerance_mb_relaxed_);
1005 }
1006
1007 if (relax_final_iteration_cnv) {
1008 message += fmt::format(" CNV: {:.1e}", param_.tolerance_cnv_relaxed_);
1009 }
1010
1011 OpmLog::debug(message);
1012 }
1013
1014 const auto tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
1015 const auto tol_mb = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_;
1016 const auto tol_cnv_energy = use_relaxed_cnv ? param_.tolerance_cnv_energy_relaxed_ : param_.tolerance_cnv_energy_;
1017 const auto tol_eb = use_relaxed_mb ? param_.tolerance_energy_balance_relaxed_ : param_.tolerance_energy_balance_;
1018
1019 // Finish computation
1020 std::vector<Scalar> CNV(numComp);
1021 std::vector<Scalar> mass_balance_residual(numComp);
1022 for (int compIdx = 0; compIdx < numComp; ++compIdx)
1023 {
1024 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
1025 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
1026 residual_norms.push_back(CNV[compIdx]);
1027 }
1028
1029 using CR = ConvergenceReport;
1030 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
1031 const Scalar res[2] = {
1032 mass_balance_residual[compIdx], CNV[compIdx],
1033 };
1034
1035 const CR::ReservoirFailure::Type types[2] = {
1036 CR::ReservoirFailure::Type::MassBalance,
1037 CR::ReservoirFailure::Type::Cnv,
1038 };
1039
1040 Scalar tol[2] = { tol_mb, tol_cnv, };
1041 if (has_energy_ && compIdx == contiEnergyEqIdx) {
1042 tol[0] = tol_eb;
1043 tol[1] = tol_cnv_energy;
1044 }
1045
1046 for (int ii : {0, 1}) {
1047 if (std::isnan(res[ii])) {
1048 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
1049 if (this->terminal_output_) {
1050 OpmLog::debug("NaN residual for " + this->compNames_.name(compIdx) + " equation.");
1051 }
1052 }
1053 else if (res[ii] > maxResidualAllowed()) {
1054 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
1055 if (this->terminal_output_) {
1056 OpmLog::debug("Too large residual for " + this->compNames_.name(compIdx) + " equation.");
1057 }
1058 }
1059 else if (res[ii] < 0.0) {
1060 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
1061 if (this->terminal_output_) {
1062 OpmLog::debug("Negative residual for " + this->compNames_.name(compIdx) + " equation.");
1063 }
1064 }
1065 else if (res[ii] > tol[ii]) {
1066 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
1067 }
1068
1069 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]);
1070 }
1071 }
1072
1073 // Output of residuals.
1074 if (this->terminal_output_) {
1075 // Only rank 0 does print to std::cout
1076 if (iteration == 0) {
1077 std::string msg = "Iter";
1078 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
1079 msg += " MB(";
1080 msg += this->compNames_.name(compIdx)[0];
1081 msg += ") ";
1082 }
1083
1084 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
1085 msg += " CNV(";
1086 msg += this->compNames_.name(compIdx)[0];
1087 msg += ") ";
1088 }
1089
1090 OpmLog::debug(msg);
1091 }
1092
1093 std::ostringstream ss;
1094 const std::streamsize oprec = ss.precision(3);
1095 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
1096
1097 ss << std::setw(4) << iteration;
1098 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
1099 ss << std::setw(11) << mass_balance_residual[compIdx];
1100 }
1101
1102 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
1103 ss << std::setw(11) << CNV[compIdx];
1104 }
1105
1106 ss.precision(oprec);
1107 ss.flags(oflags);
1108
1109 OpmLog::debug(ss.str());
1110 }
1111
1112 return report;
1113 }
1114
1115 void checkCardPenalty(ConvergenceReport& report, int iteration)
1116 {
1117
1118 const auto& current_metrics = report.reservoirConvergence();
1119 auto distances = std::vector<double>(current_metrics.size(), 0.0);
1120 int current_above_tolerance = 0;
1121
1122 for (size_t i = 0; i < current_metrics.size(); ++i) {
1123 distances[i] = std::max(std::log10(current_metrics[i].value()/current_metrics[i].tolerance()), 0.0);
1124 // Count number of metrics above tolerance
1125 if (current_metrics[i].value() > current_metrics[i].tolerance()) {
1126 current_above_tolerance++;
1127 }
1128 }
1129
1130 // use L1 norm of the distances vector
1131 double current_distance = std::accumulate(distances.begin(), distances.end(), 0.0);
1132
1133 if (iteration > 0) {
1134 // Add penalty if number of metrics above tolerance has increased
1135 if (current_above_tolerance > prev_above_tolerance_) {
1136 report.addNonConvergedPenalty();
1137 }
1138
1139 if (current_distance > param_.convergence_monitoring_decay_factor_ * prev_distance_) {
1140 report.addDistanceDecayPenalty();
1141 }
1142 }
1143
1144 prev_distance_ = current_distance;
1145 prev_above_tolerance_ = current_above_tolerance;
1146
1147 if (report.wellFailures().size() > 0) {
1148 report.addLargeWellResidualsPenalty();
1149 }
1150
1151 total_penaltyCard_ += report.getPenaltyCard();
1152
1153 if (param_.convergence_monitoring_ && (total_penaltyCard_.total() > param_.convergence_monitoring_cutoff_)) {
1154 report.setReservoirFailed({ConvergenceReport::ReservoirFailure::Type::ConvergenceMonitorFailure,
1155 ConvergenceReport::Severity::ConvergenceMonitorFailure,
1156 -1}); // -1 indicates it's not specific to any component
1157 }
1158 }
1159
1167 const int iteration,
1168 const int maxIter,
1169 std::vector<Scalar>& residual_norms)
1170 {
1171 OPM_TIMEBLOCK(getConvergence);
1172 // Get convergence reports for reservoir and wells.
1173 std::vector<Scalar> B_avg(numEq, 0.0);
1174 auto report = getReservoirConvergence(timer.simulationTimeElapsed(),
1175 timer.currentStepLength(),
1176 iteration, maxIter, B_avg, residual_norms);
1177 {
1178 OPM_TIMEBLOCK(getWellConvergence);
1179 report += wellModel().getWellConvergence(B_avg, /*checkWellGroupControls*/report.converged());
1180 }
1181
1182 checkCardPenalty(report, iteration);
1183
1184 return report;
1185 }
1186
1187
1189 int numPhases() const
1190 {
1191 return phaseUsage_.num_phases;
1192 }
1193
1195 template<class T>
1196 std::vector<std::vector<Scalar> >
1197 computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
1198 {
1199 return computeFluidInPlace(fipnum);
1200 }
1201
1203 std::vector<std::vector<Scalar> >
1204 computeFluidInPlace(const std::vector<int>& /*fipnum*/) const
1205 {
1206 OPM_TIMEBLOCK(computeFluidInPlace);
1207 //assert(true)
1208 //return an empty vector
1209 std::vector<std::vector<Scalar> > regionValues(0, std::vector<Scalar>(0,0.0));
1210 return regionValues;
1211 }
1212
1213 const Simulator& simulator() const
1214 { return simulator_; }
1215
1216 Simulator& simulator()
1217 { return simulator_; }
1218
1221 { return failureReport_; }
1222
1225 {
1226 return nlddSolver_ ? nlddSolver_->localAccumulatedReports()
1228 }
1229
1230 const std::vector<StepReport>& stepReports() const
1231 {
1232 return convergence_reports_;
1233 }
1234
1235 void writePartitions(const std::filesystem::path& odir) const
1236 {
1237 if (this->nlddSolver_ != nullptr) {
1238 this->nlddSolver_->writePartitions(odir);
1239 return;
1240 }
1241
1242 const auto& elementMapper = this->simulator().model().elementMapper();
1243 const auto& cartMapper = this->simulator().vanguard().cartesianIndexMapper();
1244
1245 const auto& grid = this->simulator().vanguard().grid();
1246 const auto& comm = grid.comm();
1247 const auto nDigit = 1 + static_cast<int>(std::floor(std::log10(comm.size())));
1248
1249 std::ofstream pfile { odir / fmt::format("{1:0>{0}}", nDigit, comm.rank()) };
1250
1251 for (const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
1252 pfile << comm.rank() << ' '
1253 << cartMapper.cartesianIndex(elementMapper.index(cell)) << ' '
1254 << comm.rank() << '\n';
1255 }
1256 }
1257
1258 const std::vector<std::vector<int>>& getConvCells() const
1259 { return rst_conv_.getData(); }
1260
1261 protected:
1262 // --------- Data members ---------
1263
1264 Simulator& simulator_;
1265 const Grid& grid_;
1266 const PhaseUsage phaseUsage_;
1267 static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
1268 static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
1269 static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
1270 static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
1271 static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
1272 static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
1273 static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
1274 static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
1275
1276 ModelParameters param_;
1277 SimulatorReportSingle failureReport_;
1278
1279 // Well Model
1280 BlackoilWellModel<TypeTag>& well_model_;
1281
1283
1287 long int global_nc_;
1288
1289 std::vector<std::vector<Scalar>> residual_norms_history_;
1290 Scalar current_relaxation_;
1291 BVector dx_old_;
1292
1293 std::vector<StepReport> convergence_reports_;
1294 ComponentName compNames_{};
1295
1296 std::unique_ptr<BlackoilModelNldd<TypeTag>> nlddSolver_;
1297
1298 public:
1301 wellModel() { return well_model_; }
1302
1304 wellModel() const { return well_model_; }
1305
1306 void beginReportStep()
1307 {
1308 simulator_.problem().beginEpisode();
1309 }
1310
1311 void endReportStep()
1312 {
1313 simulator_.problem().endEpisode();
1314 }
1315
1316 template<class FluidState, class Residual>
1317 void getMaxCoeff(const unsigned cell_idx,
1318 const IntensiveQuantities& intQuants,
1319 const FluidState& fs,
1320 const Residual& modelResid,
1321 const Scalar pvValue,
1322 std::vector<Scalar>& B_avg,
1323 std::vector<Scalar>& R_sum,
1324 std::vector<Scalar>& maxCoeff,
1325 std::vector<int>& maxCoeffCell)
1326 {
1327 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1328 {
1329 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1330 continue;
1331 }
1332
1333 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1334
1335 B_avg[compIdx] += 1.0 / fs.invB(phaseIdx).value();
1336 const auto R2 = modelResid[cell_idx][compIdx];
1337
1338 R_sum[compIdx] += R2;
1339 const Scalar Rval = std::abs(R2) / pvValue;
1340 if (Rval > maxCoeff[compIdx]) {
1341 maxCoeff[compIdx] = Rval;
1342 maxCoeffCell[compIdx] = cell_idx;
1343 }
1344 }
1345
1346 if constexpr (has_solvent_) {
1347 B_avg[contiSolventEqIdx] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
1348 const auto R2 = modelResid[cell_idx][contiSolventEqIdx];
1349 R_sum[contiSolventEqIdx] += R2;
1350 maxCoeff[contiSolventEqIdx] = std::max(maxCoeff[contiSolventEqIdx],
1351 std::abs(R2) / pvValue);
1352 }
1353 if constexpr (has_extbo_) {
1354 B_avg[contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1355 const auto R2 = modelResid[cell_idx][contiZfracEqIdx];
1356 R_sum[ contiZfracEqIdx ] += R2;
1357 maxCoeff[contiZfracEqIdx] = std::max(maxCoeff[contiZfracEqIdx],
1358 std::abs(R2) / pvValue);
1359 }
1360 if constexpr (has_polymer_) {
1361 B_avg[contiPolymerEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1362 const auto R2 = modelResid[cell_idx][contiPolymerEqIdx];
1363 R_sum[contiPolymerEqIdx] += R2;
1364 maxCoeff[contiPolymerEqIdx] = std::max(maxCoeff[contiPolymerEqIdx],
1365 std::abs(R2) / pvValue);
1366 }
1367 if constexpr (has_foam_) {
1368 B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1369 const auto R2 = modelResid[cell_idx][contiFoamEqIdx];
1370 R_sum[contiFoamEqIdx] += R2;
1371 maxCoeff[contiFoamEqIdx] = std::max(maxCoeff[contiFoamEqIdx],
1372 std::abs(R2) / pvValue);
1373 }
1374 if constexpr (has_brine_) {
1375 B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1376 const auto R2 = modelResid[cell_idx][contiBrineEqIdx];
1377 R_sum[contiBrineEqIdx] += R2;
1378 maxCoeff[contiBrineEqIdx] = std::max(maxCoeff[contiBrineEqIdx],
1379 std::abs(R2) / pvValue);
1380 }
1381
1382 if constexpr (has_polymermw_) {
1383 static_assert(has_polymer_);
1384
1385 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1386 // the residual of the polymer molecular equation is scaled down by a 100, since molecular weight
1387 // can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations
1388 // TODO: there should be a more general way to determine the scaling-down coefficient
1389 const auto R2 = modelResid[cell_idx][contiPolymerMWEqIdx] / 100.;
1390 R_sum[contiPolymerMWEqIdx] += R2;
1391 maxCoeff[contiPolymerMWEqIdx] = std::max(maxCoeff[contiPolymerMWEqIdx],
1392 std::abs(R2) / pvValue);
1393 }
1394
1395 if constexpr (has_energy_) {
1396 B_avg[contiEnergyEqIdx] += 1.0 / (4.182e1); // converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water
1397 const auto R2 = modelResid[cell_idx][contiEnergyEqIdx];
1398 R_sum[contiEnergyEqIdx] += R2;
1399 maxCoeff[contiEnergyEqIdx] = std::max(maxCoeff[contiEnergyEqIdx],
1400 std::abs(R2) / pvValue);
1401 }
1402
1403 if constexpr (has_micp_) {
1404 B_avg[contiMicrobialEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1405 const auto R1 = modelResid[cell_idx][contiMicrobialEqIdx];
1406 R_sum[contiMicrobialEqIdx] += R1;
1407 maxCoeff[contiMicrobialEqIdx] = std::max(maxCoeff[contiMicrobialEqIdx],
1408 std::abs(R1) / pvValue);
1409 B_avg[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1410 const auto R2 = modelResid[cell_idx][contiOxygenEqIdx];
1411 R_sum[contiOxygenEqIdx] += R2;
1412 maxCoeff[contiOxygenEqIdx] = std::max(maxCoeff[contiOxygenEqIdx],
1413 std::abs(R2) / pvValue);
1414 B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1415 const auto R3 = modelResid[cell_idx][contiUreaEqIdx];
1416 R_sum[contiUreaEqIdx] += R3;
1417 maxCoeff[contiUreaEqIdx] = std::max(maxCoeff[contiUreaEqIdx],
1418 std::abs(R3) / pvValue);
1419 B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1420 const auto R4 = modelResid[cell_idx][contiBiofilmEqIdx];
1421 R_sum[contiBiofilmEqIdx] += R4;
1422 maxCoeff[contiBiofilmEqIdx] = std::max(maxCoeff[contiBiofilmEqIdx],
1423 std::abs(R4) / pvValue);
1424 B_avg[contiCalciteEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1425 const auto R5 = modelResid[cell_idx][contiCalciteEqIdx];
1426 R_sum[contiCalciteEqIdx] += R5;
1427 maxCoeff[contiCalciteEqIdx] = std::max(maxCoeff[contiCalciteEqIdx],
1428 std::abs(R5) / pvValue);
1429 }
1430 }
1431
1433 const ModelParameters& param() const
1434 {
1435 return param_;
1436 }
1437
1440 {
1441 return compNames_;
1442 }
1443
1444 private:
1445 Scalar dpMaxRel() const { return param_.dp_max_rel_; }
1446 Scalar dsMax() const { return param_.ds_max_; }
1447 Scalar drMaxRel() const { return param_.dr_max_rel_; }
1448 Scalar maxResidualAllowed() const { return param_.max_residual_allowed_; }
1449 double linear_solve_setup_time_;
1450 ConvergenceReport::PenaltyCard total_penaltyCard_;
1451 double prev_distance_ = std::numeric_limits<double>::infinity();
1452 int prev_above_tolerance_ = 0;
1453 public:
1454 std::vector<bool> wasSwitched_;
1455 };
1456
1457} // namespace Opm
1458
1459#endif // OPM_BLACKOILMODEL_HEADER_INCLUDED
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Class for handling the blackoil aquifer model.
Definition BlackoilAquiferModel.hpp:50
A model implementation for three-phase black oil.
Definition BlackoilModelNldd.hpp:75
BlackoilModel(Simulator &simulator, const ModelParameters &param, BlackoilWellModel< TypeTag > &well_model, const bool terminal_output)
Construct the model.
Definition BlackoilModel.hpp:225
std::vector< std::vector< Scalar > > computeFluidInPlace(const T &, const std::vector< int > &fipnum) const
Wrapper required due to not following generic API.
Definition BlackoilModel.hpp:1197
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition BlackoilModel.hpp:617
SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Called once per nonlinear iteration.
Definition BlackoilModel.hpp:389
SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &, const int iterationIdx)
Assemble the residual and Jacobian of the nonlinear system.
Definition BlackoilModel.hpp:521
ConvergenceReport getConvergence(const SimulatorTimerInterface &timer, const int iteration, const int maxIter, std::vector< Scalar > &residual_norms)
Compute convergence based on total mass balance (tol_mb) and maximum residual mass balance (tol_cnv).
Definition BlackoilModel.hpp:1166
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition BlackoilModel.hpp:714
int numPhases() const
The number of active fluid phases in the model.
Definition BlackoilModel.hpp:1189
SimulatorReportSingle afterStep(const SimulatorTimerInterface &)
Called once after each time step.
Definition BlackoilModel.hpp:509
const ComponentName & compNames() const
Returns const reference to component names.
Definition BlackoilModel.hpp:1439
void updateSolution(const BVector &dx)
Apply an update to the primary variables.
Definition BlackoilModel.hpp:693
long int global_nc_
The number of cells of the global grid.
Definition BlackoilModel.hpp:1287
const ModelParameters & param() const
Returns const reference to model parameters.
Definition BlackoilModel.hpp:1433
std::unique_ptr< BlackoilModelNldd< TypeTag > > nlddSolver_
Non-linear DD solver.
Definition BlackoilModel.hpp:1296
std::pair< std::vector< double >, std::vector< int > > characteriseCnvPvSplit(const std::vector< Scalar > &B_avg, const double dt)
Compute pore-volume/cell count split among "converged", "relaxed converged", "unconverged" cells base...
Definition BlackoilModel.hpp:837
const SimulatorReportSingle & failureReport() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModel.hpp:1220
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition BlackoilModel.hpp:1301
std::vector< std::vector< Scalar > > computeFluidInPlace(const std::vector< int > &) const
Should not be called.
Definition BlackoilModel.hpp:1204
SimulatorReportSingle localAccumulatedReports() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModel.hpp:1224
std::pair< Scalar, Scalar > 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 BlackoilModel.hpp:782
void solveJacobianSystem(BVector &x)
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition BlackoilModel.hpp:632
SimulatorReportSingle prepareStep(const SimulatorTimerInterface &timer)
Called once before each time step.
Definition BlackoilModel.hpp:269
RSTConv rst_conv_
Helper class for RPTRST CONV.
Definition BlackoilModel.hpp:1282
bool terminal_output_
Whether we print something to std::cout.
Definition BlackoilModel.hpp:1285
Class for handling the blackoil well model.
Definition WellInterface.hpp:31
Definition ComponentName.hpp:34
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:92
Class computing RPTRST CONV output.
Definition RSTConv.hpp:34
const std::vector< std::vector< int > > & getData() const
Obtain a const-ref to the accumulated data.
Definition RSTConv.hpp:58
void update(const ResidualVector &resid)
Adds the CONV output for given residual vector.
Definition RSTConv.cpp:58
void init(const std::size_t numCells, const RSTConfig &rst_config, const std::array< int, 6 > &compIdx)
Init state at beginning of step.
Definition RSTConv.cpp:35
Interface class for SimulatorTimer objects, to be improved.
Definition SimulatorTimerInterface.hpp:34
virtual int reportStepNum() const
Current report step number. This might differ from currentStepNum in case of sub stepping.
Definition SimulatorTimerInterface.hpp:50
virtual bool lastStepFailed() const =0
Return true if last time step failed.
virtual double currentStepLength() const =0
Current step length.
virtual double simulationTimeElapsed() const =0
Time elapsed since the start of the simulation until the beginning of the current time step [s].
virtual int currentStepNum() const =0
Current step number.
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition phaseUsageFromDeck.cpp:137
Solver parameters for the BlackoilModel.
Definition BlackoilModelParameters.hpp:162
Scalar tolerance_mb_relaxed_
Relaxed mass balance tolerance (can be used when iter >= min_strict_mb_iter_).
Definition BlackoilModelParameters.hpp:182
Scalar tolerance_energy_balance_
Relative energy balance tolerance (total energy balance error).
Definition BlackoilModelParameters.hpp:184
bool update_equations_scaling_
Update scaling factors for mass balance equations.
Definition BlackoilModelParameters.hpp:249
Scalar tolerance_energy_balance_relaxed_
Relaxed energy balance tolerance (can be used when iter >= min_strict_mb_iter_).
Definition BlackoilModelParameters.hpp:186
Scalar convergence_monitoring_decay_factor_
Decay factor used in convergence monitoring.
Definition BlackoilModelParameters.hpp:316
int min_strict_mb_iter_
Minimum number of Newton iterations before we can use relaxed MB convergence criterion.
Definition BlackoilModelParameters.hpp:243
int min_strict_cnv_iter_
Minimum number of Newton iterations before we can use relaxed CNV convergence criterion.
Definition BlackoilModelParameters.hpp:240
bool use_update_stabilization_
Try to detect oscillation or stagnation.
Definition BlackoilModelParameters.hpp:252
Scalar tolerance_cnv_energy_
Local energy convergence tolerance (max of local energy errors).
Definition BlackoilModelParameters.hpp:192
Scalar tolerance_cnv_energy_relaxed_
Relaxed local energy convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViola...
Definition BlackoilModelParameters.hpp:194
Scalar tolerance_mb_
Relative mass balance tolerance (total mass balance error).
Definition BlackoilModelParameters.hpp:180
bool convergence_monitoring_
Whether to enable convergence monitoring.
Definition BlackoilModelParameters.hpp:312
Scalar tolerance_cnv_
Local convergence tolerance (max of local saturation errors).
Definition BlackoilModelParameters.hpp:188
Scalar tolerance_cnv_relaxed_
Relaxed local convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViolatedPV <...
Definition BlackoilModelParameters.hpp:190
int convergence_monitoring_cutoff_
Cut-off limit for convergence monitoring.
Definition BlackoilModelParameters.hpp:314
std::string nonlinear_solver_
Nonlinear solver type: newton or nldd.
Definition BlackoilModelParameters.hpp:294
Definition AquiferGridUtils.hpp:35
Definition FlowBaseProblemProperties.hpp:62
Enable surface volume scaling.
Definition blackoilproperties.hh:54
Enable the ECL-blackoil extension for salt.
Definition blackoilproperties.hh:60
Enable convective mixing?
Definition multiphasebaseproperties.hh:85
Definition FlowBaseProblemProperties.hpp:73
Enable dispersive fluxes?
Definition multiphasebaseproperties.hh:82
Specify whether energy should be considered as a conservation quantity or not.
Definition multiphasebaseproperties.hh:76
Enable the ECL-blackoil extension for foam.
Definition blackoilproperties.hh:57
Enable the ECL-blackoil extension for MICP.
Definition blackoilproperties.hh:72
Enable the ECL-blackoil extension for polymer.
Definition blackoilproperties.hh:48
Enable the ECL-blackoil extension for salt precipitation.
Definition blackoilproperties.hh:63
Enable the ECL-blackoil extension for solvents. ("Second gas")
Definition blackoilproperties.hh:42
Allow the spatial and temporal domains to exhibit non-constant temperature in the black-oil model.
Definition blackoilproperties.hh:78
Definition fvbaseproperties.hh:53
Definition ISTLSolver.hpp:63
Definition BlackoilModel.hpp:80
Specify whether to use volumetric residuals or not.
Definition fvbaseproperties.hh:237
Definition ISTLSolver.hpp:69
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34