My Project
Loading...
Searching...
No Matches
BlackoilModel_impl.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_IMPL_HEADER_INCLUDED
25#define OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED
26
27#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED
28#include <config.h>
29#include <opm/simulators/flow/BlackoilModel.hpp>
30#endif
31
32#include <dune/common/timer.hh>
33
34#include <opm/common/ErrorMacros.hpp>
35#include <opm/common/OpmLog/OpmLog.hpp>
36
37#include <opm/simulators/flow/countGlobalCells.hpp>
38
39#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
40
41#include <algorithm>
42#include <iomanip>
43#include <limits>
44#include <stdexcept>
45#include <sstream>
46
47#include <fmt/format.h>
48
49namespace Opm {
50
51template <class TypeTag>
53BlackoilModel(Simulator& simulator,
54 const ModelParameters& param,
56 const bool terminal_output)
57 : simulator_(simulator)
58 , grid_(simulator_.vanguard().grid())
59 , phaseUsage_(phaseUsageFromDeck(eclState()))
60 , param_( param )
61 , well_model_ (well_model)
62 , terminal_output_ (terminal_output)
63 , current_relaxation_(1.0)
64 , dx_old_(simulator_.model().numGridDof())
65 , conv_monitor_(param_.monitor_params_)
66{
67 // compute global sum of number of cells
68 global_nc_ = detail::countGlobalCells(grid_);
69 convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves.
70 // TODO: remember to fix!
71 if (param_.nonlinear_solver_ == "nldd") {
72 if (terminal_output) {
73 OpmLog::info("Using Non-Linear Domain Decomposition solver (nldd).");
74 }
75 nlddSolver_ = std::make_unique<BlackoilModelNldd<TypeTag>>(*this);
76 } else if (param_.nonlinear_solver_ == "newton") {
77 if (terminal_output) {
78 OpmLog::info("Using Newton nonlinear solver.");
79 }
80 } else {
81 OPM_THROW(std::runtime_error, "Unknown nonlinear solver option: " +
82 param_.nonlinear_solver_);
83 }
84}
85
86template <class TypeTag>
90{
92 Dune::Timer perfTimer;
93 perfTimer.start();
94 // update the solution variables in the model
95 int lastStepFailed = timer.lastStepFailed();
96 if (grid_.comm().size() > 1 && lastStepFailed != grid_.comm().min(lastStepFailed)) {
97 OPM_THROW(std::runtime_error,
98 fmt::format("Misalignment of the parallel simulation run in prepareStep "
99 "- the previous step succeeded on rank {} but failed on the "
100 "other ranks.", grid_.comm().rank()));
101 }
102 if (lastStepFailed) {
103 simulator_.model().updateFailed();
104 }
105 else {
106 simulator_.model().advanceTimeLevel();
107 }
108
109 // Set the timestep size, episode index, and non-linear iteration index
110 // for the model explicitly. The model needs to know the report step/episode index
111 // because of timing dependent data despite the fact that flow uses its
112 // own time stepper. (The length of the episode does not matter, though.)
113 simulator_.setTime(timer.simulationTimeElapsed());
114 simulator_.setTimeStepSize(timer.currentStepLength());
115 simulator_.model().newtonMethod().setIterationIndex(0);
116
117 simulator_.problem().beginTimeStep();
118
119 unsigned numDof = simulator_.model().numGridDof();
120 wasSwitched_.resize(numDof);
121 std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
122
123 if (param_.update_equations_scaling_) {
124 OpmLog::error("Equation scaling not supported");
125 //updateEquationsScaling();
126 }
127
128 if (hasNlddSolver()) {
129 nlddSolver_->prepareStep();
130 }
131
132 report.pre_post_time += perfTimer.stop();
133
134 auto getIdx = [](unsigned phaseIdx) -> int
135 {
136 if (FluidSystem::phaseIsActive(phaseIdx)) {
137 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
138 return Indices::canonicalToActiveComponentIndex(sIdx);
139 }
140
141 return -1;
142 };
143 const auto& schedule = simulator_.vanguard().schedule();
144 auto& rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv();
145 rst_conv.init(simulator_.vanguard().globalNumCells(),
146 schedule[timer.reportStepNum()].rst_config(),
147 {getIdx(FluidSystem::oilPhaseIdx),
148 getIdx(FluidSystem::gasPhaseIdx),
149 getIdx(FluidSystem::waterPhaseIdx),
150 contiPolymerEqIdx,
151 contiBrineEqIdx,
152 contiSolventEqIdx});
153
154 return report;
155}
156
157template <class TypeTag>
158void
161 const int iteration,
162 const int minIter,
163 const int maxIter,
164 const SimulatorTimerInterface& timer)
165{
166 // ----------- Set up reports and timer -----------
167 failureReport_ = SimulatorReportSingle();
168 Dune::Timer perfTimer;
169
170 perfTimer.start();
171 report.total_linearizations = 1;
172
173 // ----------- Assemble -----------
174 try {
175 report += assembleReservoir(timer, iteration);
176 report.assemble_time += perfTimer.stop();
177 }
178 catch (...) {
179 report.assemble_time += perfTimer.stop();
180 failureReport_ += report;
181 throw; // continue throwing the stick
182 }
183
184 // ----------- Check if converged -----------
185 std::vector<Scalar> residual_norms;
186 perfTimer.reset();
187 perfTimer.start();
188 // the step is not considered converged until at least minIter iterations is done
189 {
190 auto convrep = getConvergence(timer, iteration, maxIter, residual_norms);
191 report.converged = convrep.converged() && iteration >= minIter;
192 ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
193 convergence_reports_.back().report.push_back(std::move(convrep));
194
195 // Throw if any NaN or too large residual found.
196 if (severity == ConvergenceReport::Severity::NotANumber) {
197 failureReport_ += report;
198 OPM_THROW_PROBLEM(NumericalProblem, "NaN residual found!");
199 } else if (severity == ConvergenceReport::Severity::TooLarge) {
200 failureReport_ += report;
201 OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!");
202 } else if (severity == ConvergenceReport::Severity::ConvergenceMonitorFailure) {
203 failureReport_ += report;
204 OPM_THROW_PROBLEM(ConvergenceMonitorFailure,
205 fmt::format(
206 "Total penalty count exceeded cut-off-limit of {}",
207 param_.monitor_params_.cutoff_
208 ));
209 }
210 }
211 report.update_time += perfTimer.stop();
212 residual_norms_history_.push_back(residual_norms);
213}
214
215template <class TypeTag>
216template <class NonlinearSolverType>
217SimulatorReportSingle
220 const SimulatorTimerInterface& timer,
222{
223 if (iteration == 0) {
224 // For each iteration we store in a vector the norms of the residual of
225 // the mass balance for each active phase, the well flux and the well equations.
226 residual_norms_history_.clear();
227 conv_monitor_.reset();
228 current_relaxation_ = 1.0;
229 dx_old_ = 0.0;
230 convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}});
231 convergence_reports_.back().report.reserve(11);
232 }
233
235 if ((this->param_.nonlinear_solver_ != "nldd") ||
236 (iteration < this->param_.nldd_num_initial_newton_iter_))
237 {
238 result = this->nonlinearIterationNewton(iteration,
239 timer,
241 }
242 else {
243 result = this->nlddSolver_->nonlinearIterationNldd(iteration,
244 timer,
246 }
247
248 auto& rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv();
249 rst_conv.update(simulator_.model().linearizer().residual());
250
251 return result;
252}
253
254template <class TypeTag>
255template <class NonlinearSolverType>
259 const SimulatorTimerInterface& timer,
261{
262 // ----------- Set up reports and timer -----------
264 Dune::Timer perfTimer;
265
266 this->initialLinearization(report,
267 iteration,
268 nonlinear_solver.minIter(),
269 nonlinear_solver.maxIter(),
270 timer);
271
272 // ----------- If not converged, solve linear system and do Newton update -----------
273 if (!report.converged) {
274 perfTimer.reset();
275 perfTimer.start();
276 report.total_newton_iterations = 1;
277
278 // Compute the nonlinear update.
279 unsigned nc = simulator_.model().numGridDof();
280 BVector x(nc);
281
282 // Solve the linear system.
283 linear_solve_setup_time_ = 0.0;
284 try {
285 // Apply the Schur complement of the well model to
286 // the reservoir linearized equations.
287 // Note that linearize may throw for MSwells.
288 wellModel().linearize(simulator().model().linearizer().jacobian(),
289 simulator().model().linearizer().residual());
290
291 // ---- Solve linear system ----
292 solveJacobianSystem(x);
293
294 report.linear_solve_setup_time += linear_solve_setup_time_;
295 report.linear_solve_time += perfTimer.stop();
296 report.total_linear_iterations += linearIterationsLastSolve();
297 }
298 catch (...) {
299 report.linear_solve_setup_time += linear_solve_setup_time_;
300 report.linear_solve_time += perfTimer.stop();
301 report.total_linear_iterations += linearIterationsLastSolve();
302
303 failureReport_ += report;
304 throw; // re-throw up
305 }
306
307 perfTimer.reset();
308 perfTimer.start();
309
310 // handling well state update before oscillation treatment is a decision based
311 // on observation to avoid some big performance degeneration under some circumstances.
312 // there is no theorectical explanation which way is better for sure.
313 wellModel().postSolve(x);
314
315 if (param_.use_update_stabilization_) {
316 // Stabilize the nonlinear update.
317 bool isOscillate = false;
318 bool isStagnate = false;
319 nonlinear_solver.detectOscillations(residual_norms_history_,
320 residual_norms_history_.size() - 1,
322 isStagnate);
323 if (isOscillate) {
324 current_relaxation_ -= nonlinear_solver.relaxIncrement();
325 current_relaxation_ = std::max(current_relaxation_,
326 nonlinear_solver.relaxMax());
327 if (terminalOutputEnabled()) {
328 std::string msg = " Oscillating behavior detected: Relaxation set to "
329 + std::to_string(current_relaxation_);
330 OpmLog::info(msg);
331 }
332 }
333 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
334 }
335
336 // ---- Newton update ----
337 // Apply the update, with considering model-dependent limitations and
338 // chopping of the update.
339 updateSolution(x);
340
341 report.update_time += perfTimer.stop();
342 }
343
344 return report;
345}
346
347template <class TypeTag>
348SimulatorReportSingle
351{
353 Dune::Timer perfTimer;
354 perfTimer.start();
355 simulator_.problem().endTimeStep();
356 report.pre_post_time += perfTimer.stop();
357 return report;
358}
359
360template <class TypeTag>
364 const int iterationIdx)
365{
366 // -------- Mass balance equations --------
367 simulator_.model().newtonMethod().setIterationIndex(iterationIdx);
368 simulator_.problem().beginIteration();
369 simulator_.model().linearizer().linearizeDomain();
370 simulator_.problem().endIteration();
371 return wellModel().lastReport();
372}
373
374template <class TypeTag>
375typename BlackoilModel<TypeTag>::Scalar
377relativeChange() const
378{
379 Scalar resultDelta = 0.0;
380 Scalar resultDenom = 0.0;
381
382 const auto& elemMapper = simulator_.model().elementMapper();
383 const auto& gridView = simulator_.gridView();
384 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
385 unsigned globalElemIdx = elemMapper.index(elem);
386 const auto& priVarsNew = simulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
387
388 Scalar pressureNew;
389 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
390
391 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
392 Scalar oilSaturationNew = 1.0;
393 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
394 FluidSystem::numActivePhases() > 1 &&
395 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
396 {
397 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
398 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
399 }
400
401 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
402 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
403 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
404 {
405 assert(Indices::compositionSwitchIdx >= 0);
406 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
407 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
408 }
409
410 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
411 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
412 }
413
414 const auto& priVarsOld = simulator_.model().solution(/*timeIdx=*/1)[globalElemIdx];
415
416 Scalar pressureOld;
417 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
418
419 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
420 Scalar oilSaturationOld = 1.0;
421
422 // NB fix me! adding pressures changes to saturation changes does not make sense
423 Scalar tmp = pressureNew - pressureOld;
424 resultDelta += tmp*tmp;
426
427 if (FluidSystem::numActivePhases() > 1) {
428 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
429 saturationsOld[FluidSystem::waterPhaseIdx] =
430 priVarsOld[Indices::waterSwitchIdx];
431 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
432 }
433
434 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
435 {
436 assert(Indices::compositionSwitchIdx >= 0 );
437 saturationsOld[FluidSystem::gasPhaseIdx] =
438 priVarsOld[Indices::compositionSwitchIdx];
439 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
440 }
441
442 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
443 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
444 }
445 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
449 assert(std::isfinite(resultDelta));
450 assert(std::isfinite(resultDenom));
451 }
452 }
453 }
454
455 resultDelta = gridView.comm().sum(resultDelta);
456 resultDenom = gridView.comm().sum(resultDenom);
457
458 return resultDenom > 0.0 ? resultDelta / resultDenom : 0.0;
459}
460
461template <class TypeTag>
462void
464solveJacobianSystem(BVector& x)
465{
466 auto& jacobian = simulator_.model().linearizer().jacobian().istlMatrix();
467 auto& residual = simulator_.model().linearizer().residual();
468 auto& linSolver = simulator_.model().newtonMethod().linearSolver();
469
470 const int numSolvers = linSolver.numAvailableSolvers();
471 if (numSolvers > 1 && (linSolver.getSolveCount() % 100 == 0)) {
472 if (terminal_output_) {
473 OpmLog::debug("\nRunning speed test for comparing available linear solvers.");
474 }
475
476 Dune::Timer perfTimer;
477 std::vector<double> times(numSolvers);
478 std::vector<double> setupTimes(numSolvers);
479
480 x = 0.0;
481 std::vector<BVector> x_trial(numSolvers, x);
482 for (int solver = 0; solver < numSolvers; ++solver) {
483 BVector x0(x);
484 linSolver.setActiveSolver(solver);
485 perfTimer.start();
486 linSolver.prepare(jacobian, residual);
487 setupTimes[solver] = perfTimer.stop();
488 perfTimer.reset();
489 linSolver.setResidual(residual);
490 perfTimer.start();
491 linSolver.solve(x_trial[solver]);
492 times[solver] = perfTimer.stop();
493 perfTimer.reset();
494 if (terminal_output_) {
495 OpmLog::debug(fmt::format("Solver time {}: {}", solver, times[solver]));
496 }
497 }
498
499 int fastest_solver = std::min_element(times.begin(), times.end()) - times.begin();
500 // Use timing on rank 0 to determine fastest, must be consistent across ranks.
501 grid_.comm().broadcast(&fastest_solver, 1, 0);
502 linear_solve_setup_time_ = setupTimes[fastest_solver];
504 linSolver.setActiveSolver(fastest_solver);
505 }
506 else {
507 // set initial guess
508 x = 0.0;
509
510 Dune::Timer perfTimer;
511 perfTimer.start();
512 linSolver.prepare(jacobian, residual);
513 linear_solve_setup_time_ = perfTimer.stop();
514 linSolver.setResidual(residual);
515 // actually, the error needs to be calculated after setResidual in order to
516 // account for parallelization properly. since the residual of ECFV
517 // discretizations does not need to be synchronized across processes to be
518 // consistent, this is not relevant for OPM-flow...
519 linSolver.solve(x);
520 }
521}
522
523template <class TypeTag>
524void
526updateSolution(const BVector& dx)
527{
528 OPM_TIMEBLOCK(updateSolution);
529 auto& newtonMethod = simulator_.model().newtonMethod();
530 SolutionVector& solution = simulator_.model().solution(/*timeIdx=*/0);
531
532 newtonMethod.update_(/*nextSolution=*/solution,
533 /*curSolution=*/solution,
534 /*update=*/dx,
535 /*resid=*/dx); // the update routines of the black
536 // oil model do not care about the
537 // residual
538
539 // if the solution is updated, the intensive quantities need to be recalculated
540 {
541 OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities);
542 simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
543 }
544}
545
546template <class TypeTag>
547std::tuple<typename BlackoilModel<TypeTag>::Scalar,
548 typename BlackoilModel<TypeTag>::Scalar>
550convergenceReduction(Parallel::Communication comm,
551 const Scalar pvSumLocal,
552 const Scalar numAquiferPvSumLocal,
553 std::vector< Scalar >& R_sum,
554 std::vector< Scalar >& maxCoeff,
555 std::vector< Scalar >& B_avg)
556{
557 OPM_TIMEBLOCK(convergenceReduction);
558 // Compute total pore volume (use only owned entries)
559 Scalar pvSum = pvSumLocal;
561
562 if (comm.size() > 1) {
563 // global reduction
564 std::vector< Scalar > sumBuffer;
565 std::vector< Scalar > maxBuffer;
566 const int numComp = B_avg.size();
567 sumBuffer.reserve( 2*numComp + 2 ); // +2 for (numAquifer)pvSum
568 maxBuffer.reserve( numComp );
569 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
570 sumBuffer.push_back(B_avg[compIdx]);
571 sumBuffer.push_back(R_sum[compIdx]);
572 maxBuffer.push_back(maxCoeff[ compIdx]);
573 }
574
575 // Compute total pore volume
576 sumBuffer.push_back( pvSum );
577 sumBuffer.push_back( numAquiferPvSum );
578
579 // compute global sum
580 comm.sum( sumBuffer.data(), sumBuffer.size() );
581
582 // compute global max
583 comm.max( maxBuffer.data(), maxBuffer.size() );
584
585 // restore values to local variables
586 for (int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx) {
588 ++buffIdx;
589
591 }
592
593 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
595 }
596
597 // restore global pore volume
598 pvSum = sumBuffer[sumBuffer.size()-2];
599 numAquiferPvSum = sumBuffer.back();
600 }
601
602 // return global pore volume
603 return {pvSum, numAquiferPvSum};
604}
605
606template <class TypeTag>
607std::pair<typename BlackoilModel<TypeTag>::Scalar,
608 typename BlackoilModel<TypeTag>::Scalar>
610localConvergenceData(std::vector<Scalar>& R_sum,
611 std::vector<Scalar>& maxCoeff,
612 std::vector<Scalar>& B_avg,
613 std::vector<int>& maxCoeffCell)
614{
615 OPM_TIMEBLOCK(localConvergenceData);
616 Scalar pvSumLocal = 0.0;
617 Scalar numAquiferPvSumLocal = 0.0;
618 const auto& model = simulator_.model();
619 const auto& problem = simulator_.problem();
620
621 const auto& residual = simulator_.model().linearizer().residual();
622
623 ElementContext elemCtx(simulator_);
624 const auto& gridView = simulator().gridView();
626 OPM_BEGIN_PARALLEL_TRY_CATCH();
627 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
628 elemCtx.updatePrimaryStencil(elem);
629 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
630
631 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
632 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
633 const auto& fs = intQuants.fluidState();
634
635 const auto pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) *
636 model.dofTotalVolume(cell_idx);
638
641 }
642
643 this->getMaxCoeff(cell_idx, intQuants, fs, residual, pvValue,
645 }
646
647 OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::localConvergenceData() failed: ", grid_.comm());
648
649 // compute local average in terms of global number of elements
650 const int bSize = B_avg.size();
651 for (int i = 0; i < bSize; ++i) {
652 B_avg[i] /= Scalar(global_nc_);
653 }
654
656}
657
658template <class TypeTag>
659std::pair<std::vector<double>, std::vector<int>>
661characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt)
662{
664
665 // 0: cnv <= tolerance_cnv
666 // 1: tolerance_cnv < cnv <= tolerance_cnv_relaxed
667 // 2: tolerance_cnv_relaxed < cnv
668 constexpr auto numPvGroups = std::vector<double>::size_type{3};
669
670 auto cnvPvSplit = std::pair<std::vector<double>, std::vector<int>> {
671 std::piecewise_construct,
672 std::forward_as_tuple(numPvGroups),
673 std::forward_as_tuple(numPvGroups)
674 };
675
676 auto maxCNV = [&B_avg, dt](const auto& residual, const double pvol)
677 {
678 return (dt / pvol) *
679 std::inner_product(residual.begin(), residual.end(),
680 B_avg.begin(), Scalar{0},
681 [](const Scalar m, const auto& x)
682 {
683 using std::abs;
684 return std::max(m, abs(x));
685 }, std::multiplies<>{});
686 };
687
688 auto& [splitPV, cellCntPV] = cnvPvSplit;
689
690 const auto& model = this->simulator().model();
691 const auto& problem = this->simulator().problem();
692 const auto& residual = model.linearizer().residual();
693 const auto& gridView = this->simulator().gridView();
694
695 const IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
696
697 ElementContext elemCtx(this->simulator());
698
699 OPM_BEGIN_PARALLEL_TRY_CATCH();
700 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
701 // Skip cells of numerical Aquifer
703 continue;
704 }
705
706 elemCtx.updatePrimaryStencil(elem);
707
708 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
709 const auto pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0)
710 * model.dofTotalVolume(cell_idx);
711
712 const auto maxCnv = maxCNV(residual[cell_idx], pvValue);
713
714 const auto ix = (maxCnv > this->param_.tolerance_cnv_)
715 + (maxCnv > this->param_.tolerance_cnv_relaxed_);
716
717 splitPV[ix] += static_cast<double>(pvValue);
718 ++cellCntPV[ix];
719 }
720
721 OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::characteriseCnvPvSplit() failed: ",
722 this->grid_.comm());
723
724 this->grid_.comm().sum(splitPV .data(), splitPV .size());
725 this->grid_.comm().sum(cellCntPV.data(), cellCntPV.size());
726
727 return cnvPvSplit;
728}
729
730template <class TypeTag>
731void
734{
735 this->param_.tolerance_mb_ = tuning.XXXMBE;
736
737 if (terminal_output_) {
738 OpmLog::debug(fmt::format("Setting BlackoilModel mass "
739 "balance limit (XXXMBE) to {:.2e}",
740 tuning.XXXMBE));
741 }
742}
743
744template <class TypeTag>
745ConvergenceReport
746BlackoilModel<TypeTag>::
747getReservoirConvergence(const double reportTime,
748 const double dt,
749 const int iteration,
750 const int maxIter,
751 std::vector<Scalar>& B_avg,
752 std::vector<Scalar>& residual_norms)
753{
754 OPM_TIMEBLOCK(getReservoirConvergence);
755 using Vector = std::vector<Scalar>;
756
757 ConvergenceReport report{reportTime};
758
759 const int numComp = numEq;
760
761 Vector R_sum(numComp, Scalar{0});
762 Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest());
763 std::vector<int> maxCoeffCell(numComp, -1);
764
765 const auto [pvSumLocal, numAquiferPvSumLocal] =
766 this->localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell);
767
768 // compute global sum and max of quantities
769 const auto& [pvSum, numAquiferPvSum] =
770 this->convergenceReduction(this->grid_.comm(),
774
775 report.setCnvPoreVolSplit(this->characteriseCnvPvSplit(B_avg, dt),
777
778 // For each iteration, we need to determine whether to use the
779 // relaxed tolerances. To disable the usage of relaxed
780 // tolerances, you can set the relaxed tolerances as the strict
781 // tolerances. If min_strict_mb_iter = -1 (default) we use a
782 // relaxed tolerance for the mass balance for the last
783 // iterations. For positive values we use the relaxed tolerance
784 // after the given number of iterations
785 const bool relax_final_iteration_mb =
786 this->param_.min_strict_mb_iter_ < 0 && iteration == maxIter;
787
789 (this->param_.min_strict_mb_iter_ >= 0 &&
790 iteration >= this->param_.min_strict_mb_iter_);
791
792 // If min_strict_cnv_iter = -1 we use a relaxed tolerance for
793 // the cnv for the last iterations. For positive values we use
794 // the relaxed tolerance after the given number of iterations.
795 // We also use relaxed tolerances for cells with total
796 // pore-volume less than relaxed_max_pv_fraction_. Default
797 // value of relaxed_max_pv_fraction_ is 0.03
798 const bool relax_final_iteration_cnv =
799 this->param_.min_strict_cnv_iter_ < 0 && iteration == maxIter;
800
801 const bool relax_iter_cnv =
802 this->param_.min_strict_cnv_iter_ >= 0 &&
803 iteration >= this->param_.min_strict_cnv_iter_;
804
805 // Note trailing parentheses here, just before the final
806 // semicolon. This is an immediately invoked function
807 // expression which calculates a single boolean value.
808 const auto relax_pv_fraction_cnv =
809 [&report, this, eligible = pvSum - numAquiferPvSum]()
810 {
811 const auto& cnvPvSplit = report.cnvPvSplit().first;
812
813 // [1]: tol < cnv <= relaxed
814 // [2]: relaxed < cnv
815 return static_cast<Scalar>(cnvPvSplit[1] + cnvPvSplit[2]) <
816 this->param_.relaxed_max_pv_fraction_ * eligible;
817 }();
818
822
824 this->terminal_output_)
825 {
826 std::string message =
827 "Number of newton iterations reached its maximum "
828 "try to continue with relaxed tolerances:";
829
831 message += fmt::format(" MB: {:.1e}", param_.tolerance_mb_relaxed_);
832 }
833
835 message += fmt::format(" CNV: {:.1e}", param_.tolerance_cnv_relaxed_);
836 }
837
838 OpmLog::debug(message);
839 }
840
841 const auto tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
842 const auto tol_mb = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_;
843 const auto tol_cnv_energy = use_relaxed_cnv ? param_.tolerance_cnv_energy_relaxed_ : param_.tolerance_cnv_energy_;
844 const auto tol_eb = use_relaxed_mb ? param_.tolerance_energy_balance_relaxed_ : param_.tolerance_energy_balance_;
845
846 // Finish computation
847 std::vector<Scalar> CNV(numComp);
848 std::vector<Scalar> mass_balance_residual(numComp);
849 for (int compIdx = 0; compIdx < numComp; ++compIdx)
850 {
853 residual_norms.push_back(CNV[compIdx]);
854 }
855
856 using CR = ConvergenceReport;
857 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
858 const Scalar res[2] = {
860 };
861
862 const CR::ReservoirFailure::Type types[2] = {
863 CR::ReservoirFailure::Type::MassBalance,
864 CR::ReservoirFailure::Type::Cnv,
865 };
866
867 Scalar tol[2] = { tol_mb, tol_cnv, };
868 if (has_energy_ && compIdx == contiEnergyEqIdx) {
869 tol[0] = tol_eb;
870 tol[1] = tol_cnv_energy;
871 }
872
873 for (int ii : {0, 1}) {
874 if (std::isnan(res[ii])) {
875 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
876 if (this->terminal_output_) {
877 OpmLog::debug("NaN residual for " + this->compNames_.name(compIdx) + " equation.");
878 }
879 }
880 else if (res[ii] > maxResidualAllowed()) {
881 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
882 if (this->terminal_output_) {
883 OpmLog::debug("Too large residual for " + this->compNames_.name(compIdx) + " equation.");
884 }
885 }
886 else if (res[ii] < 0.0) {
887 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
888 if (this->terminal_output_) {
889 OpmLog::debug("Negative residual for " + this->compNames_.name(compIdx) + " equation.");
890 }
891 }
892 else if (res[ii] > tol[ii]) {
893 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
894 }
895
896 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]);
897 }
898 }
899
900 // Output of residuals.
901 if (this->terminal_output_) {
902 // Only rank 0 does print to std::cout
903 if (iteration == 0) {
904 std::string msg = "Iter";
905 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
906 msg += " MB(";
907 msg += this->compNames_.name(compIdx)[0];
908 msg += ") ";
909 }
910
911 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
912 msg += " CNV(";
913 msg += this->compNames_.name(compIdx)[0];
914 msg += ") ";
915 }
916
917 OpmLog::debug(msg);
918 }
919
920 std::ostringstream ss;
921 const std::streamsize oprec = ss.precision(3);
922 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
923
924 ss << std::setw(4) << iteration;
925 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
926 ss << std::setw(11) << mass_balance_residual[compIdx];
927 }
928
929 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
930 ss << std::setw(11) << CNV[compIdx];
931 }
932
933 ss.precision(oprec);
934 ss.flags(oflags);
935
936 OpmLog::debug(ss.str());
937 }
938
939 return report;
940}
941
942template <class TypeTag>
943ConvergenceReport
946 const int iteration,
947 const int maxIter,
948 std::vector<Scalar>& residual_norms)
949{
950 OPM_TIMEBLOCK(getConvergence);
951 // Get convergence reports for reservoir and wells.
952 std::vector<Scalar> B_avg(numEq, 0.0);
953 auto report = getReservoirConvergence(timer.simulationTimeElapsed(),
954 timer.currentStepLength(),
955 iteration, maxIter, B_avg, residual_norms);
956 {
957 OPM_TIMEBLOCK(getWellConvergence);
958 report += wellModel().getWellConvergence(B_avg,
959 /*checkWellGroupControls*/report.converged());
960 }
961
962 conv_monitor_.checkPenaltyCard(report, iteration);
963
964 return report;
965}
966
967template <class TypeTag>
968std::vector<std::vector<typename BlackoilModel<TypeTag>::Scalar> >
970computeFluidInPlace(const std::vector<int>& /*fipnum*/) const
971{
972 OPM_TIMEBLOCK(computeFluidInPlace);
973 // assert(true)
974 // return an empty vector
975 std::vector<std::vector<Scalar> > regionValues(0, std::vector<Scalar>(0,0.0));
976 return regionValues;
977}
978
979template <class TypeTag>
980const SimulatorReport&
983{
984 if (!hasNlddSolver()) {
985 OPM_THROW(std::runtime_error, "Cannot get local reports from a model without NLDD solver");
986 }
987 return nlddSolver_->localAccumulatedReports();
988}
989
990template <class TypeTag>
991const std::vector<SimulatorReport>&
994{
995 if (!nlddSolver_)
996 OPM_THROW(std::runtime_error, "Cannot get domain reports from a model without NLDD solver");
997 return nlddSolver_->domainAccumulatedReports();
998}
999
1000template <class TypeTag>
1001void
1003writeNonlinearIterationsPerCell(const std::filesystem::path& odir) const
1004{
1005 if (hasNlddSolver()) {
1006 nlddSolver_->writeNonlinearIterationsPerCell(odir);
1007 }
1008}
1009
1010template <class TypeTag>
1011void
1013writePartitions(const std::filesystem::path& odir) const
1014{
1015 if (hasNlddSolver()) {
1016 nlddSolver_->writePartitions(odir);
1017 return;
1018 }
1019
1020 const auto& elementMapper = this->simulator().model().elementMapper();
1021 const auto& cartMapper = this->simulator().vanguard().cartesianIndexMapper();
1022
1023 const auto& grid = this->simulator().vanguard().grid();
1024 const auto& comm = grid.comm();
1025 const auto nDigit = 1 + static_cast<int>(std::floor(std::log10(comm.size())));
1026
1027 std::ofstream pfile {odir / fmt::format("{1:0>{0}}", nDigit, comm.rank())};
1028
1029 for (const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
1030 pfile << comm.rank() << ' '
1031 << cartMapper.cartesianIndex(elementMapper.index(cell)) << ' '
1032 << comm.rank() << '\n';
1033 }
1034}
1035
1036template <class TypeTag>
1037template<class FluidState, class Residual>
1038void
1039BlackoilModel<TypeTag>::
1040getMaxCoeff(const unsigned cell_idx,
1041 const IntensiveQuantities& intQuants,
1042 const FluidState& fs,
1043 const Residual& modelResid,
1044 const Scalar pvValue,
1045 std::vector<Scalar>& B_avg,
1046 std::vector<Scalar>& R_sum,
1047 std::vector<Scalar>& maxCoeff,
1048 std::vector<int>& maxCoeffCell)
1049{
1050 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1051 {
1052 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1053 continue;
1054 }
1055
1056 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
1057 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(sIdx);
1058
1059 B_avg[compIdx] += 1.0 / fs.invB(phaseIdx).value();
1060 const auto R2 = modelResid[cell_idx][compIdx];
1061
1062 R_sum[compIdx] += R2;
1063 const Scalar Rval = std::abs(R2) / pvValue;
1064 if (Rval > maxCoeff[compIdx]) {
1067 }
1068 }
1069
1070 if constexpr (has_solvent_) {
1071 B_avg[contiSolventEqIdx] +=
1072 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
1073 const auto R2 = modelResid[cell_idx][contiSolventEqIdx];
1074 R_sum[contiSolventEqIdx] += R2;
1075 maxCoeff[contiSolventEqIdx] = std::max(maxCoeff[contiSolventEqIdx],
1076 std::abs(R2) / pvValue);
1077 }
1078 if constexpr (has_extbo_) {
1079 B_avg[contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1080 const auto R2 = modelResid[cell_idx][contiZfracEqIdx];
1081 R_sum[ contiZfracEqIdx ] += R2;
1082 maxCoeff[contiZfracEqIdx] = std::max(maxCoeff[contiZfracEqIdx],
1083 std::abs(R2) / pvValue);
1084 }
1085 if constexpr (has_polymer_) {
1086 B_avg[contiPolymerEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1087 const auto R2 = modelResid[cell_idx][contiPolymerEqIdx];
1088 R_sum[contiPolymerEqIdx] += R2;
1089 maxCoeff[contiPolymerEqIdx] = std::max(maxCoeff[contiPolymerEqIdx],
1090 std::abs(R2) / pvValue);
1091 }
1092 if constexpr (has_foam_) {
1093 B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1094 const auto R2 = modelResid[cell_idx][contiFoamEqIdx];
1095 R_sum[contiFoamEqIdx] += R2;
1096 maxCoeff[contiFoamEqIdx] = std::max(maxCoeff[contiFoamEqIdx],
1097 std::abs(R2) / pvValue);
1098 }
1099 if constexpr (has_brine_) {
1100 B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1101 const auto R2 = modelResid[cell_idx][contiBrineEqIdx];
1102 R_sum[contiBrineEqIdx] += R2;
1103 maxCoeff[contiBrineEqIdx] = std::max(maxCoeff[contiBrineEqIdx],
1104 std::abs(R2) / pvValue);
1105 }
1106
1107 if constexpr (has_polymermw_) {
1108 static_assert(has_polymer_);
1109
1110 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1111 // the residual of the polymer molecular equation is scaled down by a 100, since molecular weight
1112 // can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations
1113 // TODO: there should be a more general way to determine the scaling-down coefficient
1114 const auto R2 = modelResid[cell_idx][contiPolymerMWEqIdx] / 100.;
1115 R_sum[contiPolymerMWEqIdx] += R2;
1116 maxCoeff[contiPolymerMWEqIdx] = std::max(maxCoeff[contiPolymerMWEqIdx],
1117 std::abs(R2) / pvValue);
1118 }
1119
1120 if constexpr (has_energy_) {
1121 B_avg[contiEnergyEqIdx] += 1.0 / (4.182e1); // converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water
1122 const auto R2 = modelResid[cell_idx][contiEnergyEqIdx];
1123 R_sum[contiEnergyEqIdx] += R2;
1124 maxCoeff[contiEnergyEqIdx] = std::max(maxCoeff[contiEnergyEqIdx],
1125 std::abs(R2) / pvValue);
1126 }
1127
1128 if constexpr (has_micp_) {
1129 B_avg[contiMicrobialEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1130 const auto R1 = modelResid[cell_idx][contiMicrobialEqIdx];
1131 R_sum[contiMicrobialEqIdx] += R1;
1132 maxCoeff[contiMicrobialEqIdx] = std::max(maxCoeff[contiMicrobialEqIdx],
1133 std::abs(R1) / pvValue);
1134 B_avg[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1135 const auto R2 = modelResid[cell_idx][contiOxygenEqIdx];
1136 R_sum[contiOxygenEqIdx] += R2;
1137 maxCoeff[contiOxygenEqIdx] = std::max(maxCoeff[contiOxygenEqIdx],
1138 std::abs(R2) / pvValue);
1139 B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1140 const auto R3 = modelResid[cell_idx][contiUreaEqIdx];
1141 R_sum[contiUreaEqIdx] += R3;
1142 maxCoeff[contiUreaEqIdx] = std::max(maxCoeff[contiUreaEqIdx],
1143 std::abs(R3) / pvValue);
1144 B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1145 const auto R4 = modelResid[cell_idx][contiBiofilmEqIdx];
1146 R_sum[contiBiofilmEqIdx] += R4;
1147 maxCoeff[contiBiofilmEqIdx] = std::max(maxCoeff[contiBiofilmEqIdx],
1148 std::abs(R4) / pvValue);
1149 B_avg[contiCalciteEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1150 const auto R5 = modelResid[cell_idx][contiCalciteEqIdx];
1151 R_sum[contiCalciteEqIdx] += R5;
1152 maxCoeff[contiCalciteEqIdx] = std::max(maxCoeff[contiCalciteEqIdx],
1153 std::abs(R5) / pvValue);
1154 }
1155}
1156
1157} // namespace Opm
1158
1159#endif // OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED
A model implementation for three-phase black oil.
Definition BlackoilModel.hpp:62
BlackoilModel(Simulator &simulator, const ModelParameters &param, BlackoilWellModel< TypeTag > &well_model, const bool terminal_output)
Construct the model.
Definition BlackoilModel_impl.hpp:53
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:245
SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Called once per nonlinear iteration.
Definition BlackoilModel_impl.hpp:219
void writeNonlinearIterationsPerCell(const std::filesystem::path &odir) const
Write the number of nonlinear iterations per cell to a file in ResInsight compatible format.
Definition BlackoilModel_impl.hpp:1003
SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &, const int iterationIdx)
Assemble the residual and Jacobian of the nonlinear system.
Definition BlackoilModel_impl.hpp:363
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_impl.hpp:945
const std::vector< SimulatorReport > & domainAccumulatedReports() const
return the statistics of local solves accumulated for each domain on this rank
Definition BlackoilModel_impl.hpp:993
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_impl.hpp:610
const SimulatorReport & localAccumulatedReports() const
return the statistics of local solves accumulated for this rank
Definition BlackoilModel_impl.hpp:982
SimulatorReportSingle afterStep(const SimulatorTimerInterface &)
Called once after each time step.
Definition BlackoilModel_impl.hpp:350
void updateSolution(const BVector &dx)
Apply an update to the primary variables.
Definition BlackoilModel_impl.hpp:526
long int global_nc_
The number of cells of the global grid.
Definition BlackoilModel.hpp:337
std::unique_ptr< BlackoilModelNldd< TypeTag > > nlddSolver_
Non-linear DD solver.
Definition BlackoilModel.hpp:346
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_impl.hpp:661
void solveJacobianSystem(BVector &x)
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition BlackoilModel_impl.hpp:464
SimulatorReportSingle prepareStep(const SimulatorTimerInterface &timer)
Called once before each time step.
Definition BlackoilModel_impl.hpp:89
Class for handling the blackoil well model.
Definition BlackoilWellModel.hpp:94
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:109
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.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
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:174
std::string nonlinear_solver_
Nonlinear solver type: newton or nldd.
Definition BlackoilModelParameters.hpp:321
Definition AquiferGridUtils.hpp:35
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34
Definition SimulatorReport.hpp:119