My Project
Loading...
Searching...
No Matches
BlackoilModelEbos.hpp
1/*
2 Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3 Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
4 Copyright 2014, 2015 Statoil ASA.
5 Copyright 2015 NTNU
6 Copyright 2015, 2016, 2017 IRIS AS
7
8 This file is part of the Open Porous Media project (OPM).
9
10 OPM is free software: you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation, either version 3 of the License, or
13 (at your option) any later version.
14
15 OPM is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with OPM. If not, see <http://www.gnu.org/licenses/>.
22*/
23
24#ifndef OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
25#define OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
26
27#include <ebos/eclproblem.hh>
28#include <opm/models/utils/start.hh>
29
30#include <opm/common/ErrorMacros.hpp>
31#include <opm/common/Exceptions.hpp>
32#include <opm/common/OpmLog/OpmLog.hpp>
33
34#include <opm/core/props/phaseUsageFromDeck.hpp>
35
36#include <opm/grid/UnstructuredGrid.h>
37
38#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
39#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
40
41#include <opm/simulators/aquifers/BlackoilAquiferModel.hpp>
42#include <opm/simulators/flow/countGlobalCells.hpp>
43#include <opm/simulators/flow/NonlinearSolverEbos.hpp>
44#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
45#include <opm/simulators/linalg/ISTLSolverEbos.hpp>
46#include <opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp>
47#include <opm/simulators/timestepping/ConvergenceReport.hpp>
48#include <opm/simulators/timestepping/SimulatorReport.hpp>
49#include <opm/simulators/timestepping/SimulatorTimer.hpp>
50
51#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
52#include <opm/simulators/wells/BlackoilWellModel.hpp>
53#include <opm/simulators/wells/WellConnectionAuxiliaryModule.hpp>
54
55#include <dune/istl/owneroverlapcopy.hh>
56#include <dune/common/parallel/communication.hh>
57#include <dune/common/timer.hh>
58#include <dune/common/unused.hh>
59
60#include <algorithm>
61#include <cassert>
62#include <cmath>
63#include <iomanip>
64#include <iostream>
65#include <limits>
66#include <utility>
67#include <vector>
68
69namespace Opm::Properties {
70
71namespace TTag {
76}
77template<class TypeTag>
78struct OutputDir<TypeTag, TTag::EclFlowProblem> {
79 static constexpr auto value = "";
80};
81template<class TypeTag>
82struct EnableDebuggingChecks<TypeTag, TTag::EclFlowProblem> {
83 static constexpr bool value = false;
84};
85// default in flow is to formulate the equations in surface volumes
86template<class TypeTag>
87struct BlackoilConserveSurfaceVolume<TypeTag, TTag::EclFlowProblem> {
88 static constexpr bool value = true;
89};
90template<class TypeTag>
91struct UseVolumetricResidual<TypeTag, TTag::EclFlowProblem> {
92 static constexpr bool value = false;
93};
94
95template<class TypeTag>
96struct EclAquiferModel<TypeTag, TTag::EclFlowProblem> {
98};
99
100// disable all extensions supported by black oil model. this should not really be
101// necessary but it makes things a bit more explicit
102template<class TypeTag>
103struct EnablePolymer<TypeTag, TTag::EclFlowProblem> {
104 static constexpr bool value = false;
105};
106template<class TypeTag>
107struct EnableSolvent<TypeTag, TTag::EclFlowProblem> {
108 static constexpr bool value = false;
109};
110template<class TypeTag>
111struct EnableTemperature<TypeTag, TTag::EclFlowProblem> {
112 static constexpr bool value = true;
113};
114template<class TypeTag>
115struct EnableEnergy<TypeTag, TTag::EclFlowProblem> {
116 static constexpr bool value = false;
117};
118template<class TypeTag>
119struct EnableFoam<TypeTag, TTag::EclFlowProblem> {
120 static constexpr bool value = false;
121};
122template<class TypeTag>
123struct EnableBrine<TypeTag, TTag::EclFlowProblem> {
124 static constexpr bool value = false;
125};
126template<class TypeTag>
127struct EnableSaltPrecipitation<TypeTag, TTag::EclFlowProblem> {
128 static constexpr bool value = false;
129};
130template<class TypeTag>
131struct EnableMICP<TypeTag, TTag::EclFlowProblem> {
132 static constexpr bool value = false;
133};
134
135template<class TypeTag>
136struct EclWellModel<TypeTag, TTag::EclFlowProblem> {
138};
139template<class TypeTag>
140struct LinearSolverSplice<TypeTag, TTag::EclFlowProblem> {
142};
143
144} // namespace Opm::Properties
145
146namespace Opm {
153 template <class TypeTag>
155 {
156 public:
157 // --------- Types and enums ---------
159
170
171 typedef double Scalar;
172 static const int numEq = Indices::numEq;
173 static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
174 static const int contiZfracEqIdx = Indices::contiZfracEqIdx;
175 static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
176 static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
177 static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
178 static const int contiFoamEqIdx = Indices::contiFoamEqIdx;
179 static const int contiBrineEqIdx = Indices::contiBrineEqIdx;
180 static const int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
181 static const int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
182 static const int contiUreaEqIdx = Indices::contiUreaEqIdx;
183 static const int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
184 static const int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
185 static const int solventSaturationIdx = Indices::solventSaturationIdx;
186 static const int zFractionIdx = Indices::zFractionIdx;
187 static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
188 static const int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
189 static const int temperatureIdx = Indices::temperatureIdx;
190 static const int foamConcentrationIdx = Indices::foamConcentrationIdx;
191 static const int saltConcentrationIdx = Indices::saltConcentrationIdx;
192 static const int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
193 static const int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
194 static const int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
195 static const int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
196 static const int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
197
198 typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
199 typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
200 typedef typename SparseMatrixAdapter::IstlMatrix Mat;
201 typedef Dune::BlockVector<VectorBlockType> BVector;
202
204
206 {
207 public:
209 : names_(numEq)
210 {
211 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
212 if (!FluidSystem::phaseIsActive(phaseIdx)) {
213 continue;
214 }
215
216 const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
217 names_[Indices::canonicalToActiveComponentIndex(canonicalCompIdx)]
218 = FluidSystem::componentName(canonicalCompIdx);
219 }
220
221 if constexpr (has_solvent_) {
222 names_[solventSaturationIdx] = "Solvent";
223 }
224
225 if constexpr (has_extbo_) {
226 names_[zFractionIdx] = "ZFraction";
227 }
228
229 if constexpr (has_polymer_) {
230 names_[polymerConcentrationIdx] = "Polymer";
231 }
232
233 if constexpr (has_polymermw_) {
234 assert(has_polymer_);
235 names_[polymerMoleWeightIdx] = "MolecularWeightP";
236 }
237
238 if constexpr (has_energy_) {
239 names_[temperatureIdx] = "Energy";
240 }
241
242 if constexpr (has_foam_) {
243 names_[foamConcentrationIdx] = "Foam";
244 }
245
246 if constexpr (has_brine_) {
247 names_[saltConcentrationIdx] = "Brine";
248 }
249
250 if constexpr (has_micp_) {
251 names_[microbialConcentrationIdx] = "Microbes";
252 names_[oxygenConcentrationIdx] = "Oxygen";
253 names_[ureaConcentrationIdx] = "Urea";
254 names_[biofilmConcentrationIdx] = "Biofilm";
255 names_[calciteConcentrationIdx] = "Calcite";
256 }
257 }
258
259 const std::string& name(const int compIdx) const
260 {
261 return this->names_[compIdx];
262 }
263
264 private:
265 std::vector<std::string> names_{};
266 };
267
268 //typedef typename SolutionVector :: value_type PrimaryVariables ;
269
270 // --------- Public methods ---------
271
282 BlackoilModelEbos(Simulator& ebosSimulator,
283 const ModelParameters& param,
285 const bool terminal_output)
286 : ebosSimulator_(ebosSimulator)
287 , grid_(ebosSimulator_.vanguard().grid())
288 , phaseUsage_(phaseUsageFromDeck(eclState()))
289 , param_( param )
290 , well_model_ (well_model)
292 , current_relaxation_(1.0)
293 , dx_old_(ebosSimulator_.model().numGridDof())
294 {
295 // compute global sum of number of cells
296 global_nc_ = detail::countGlobalCells(grid_);
297 convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves.
298 }
299
300 bool isParallel() const
301 { return grid_.comm().size() > 1; }
302
303 const EclipseState& eclState() const
304 { return ebosSimulator_.vanguard().eclState(); }
305
309 {
311 Dune::Timer perfTimer;
312 perfTimer.start();
313 // update the solution variables in ebos
314 if ( timer.lastStepFailed() ) {
315 ebosSimulator_.model().updateFailed();
316 } else {
317 ebosSimulator_.model().advanceTimeLevel();
318 }
319
320 // Set the timestep size, episode index, and non-linear iteration index
321 // for ebos explicitly. ebos needs to know the report step/episode index
322 // because of timing dependent data despite the fact that flow uses its
323 // own time stepper. (The length of the episode does not matter, though.)
324 ebosSimulator_.setTime(timer.simulationTimeElapsed());
325 ebosSimulator_.setTimeStepSize(timer.currentStepLength());
326 ebosSimulator_.model().newtonMethod().setIterationIndex(0);
327
328 ebosSimulator_.problem().beginTimeStep();
329
330 unsigned numDof = ebosSimulator_.model().numGridDof();
331 wasSwitched_.resize(numDof);
332 std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
333
334 if (param_.update_equations_scaling_) {
335 std::cout << "equation scaling not supported yet" << std::endl;
336 //updateEquationsScaling();
337 }
338 report.pre_post_time += perfTimer.stop();
339
340 return report;
341 }
342
343
353 template <class NonlinearSolverType>
357 {
359 failureReport_ = SimulatorReportSingle();
360 Dune::Timer perfTimer;
361
362 perfTimer.start();
363 if (iteration == 0) {
364 // For each iteration we store in a vector the norms of the residual of
365 // the mass balance for each active phase, the well flux and the well equations.
366 residual_norms_history_.clear();
367 current_relaxation_ = 1.0;
368 dx_old_ = 0.0;
369 convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}});
370 convergence_reports_.back().report.reserve(11);
371 }
372
373 report.total_linearizations = 1;
374
375 try {
377 report.assemble_time += perfTimer.stop();
378 }
379 catch (...) {
380 report.assemble_time += perfTimer.stop();
381 failureReport_ += report;
382 // todo (?): make the report an attribute of the class
383 throw; // continue throwing the stick
384 }
385
386 std::vector<double> residual_norms;
387 perfTimer.reset();
388 perfTimer.start();
389 // the step is not considered converged until at least minIter iterations is done
390 {
392 report.converged = convrep.converged() && iteration > nonlinear_solver.minIter();;
393 ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
394 convergence_reports_.back().report.push_back(std::move(convrep));
395
396 // Throw if any NaN or too large residual found.
397 if (severity == ConvergenceReport::Severity::NotANumber) {
398 OPM_THROW(NumericalProblem, "NaN residual found!");
399 } else if (severity == ConvergenceReport::Severity::TooLarge) {
400 OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!");
401 }
402 }
403 report.update_time += perfTimer.stop();
404 residual_norms_history_.push_back(residual_norms);
405 if (!report.converged) {
406 perfTimer.reset();
407 perfTimer.start();
408 report.total_newton_iterations = 1;
409
410 // enable single precision for solvers when dt is smaller then 20 days
411 //residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
412
413 // Compute the nonlinear update.
414 unsigned nc = ebosSimulator_.model().numGridDof();
415 BVector x(nc);
416
417 // Solve the linear system.
418 linear_solve_setup_time_ = 0.0;
419 try {
420 // apply the Schur compliment of the well model to the reservoir linearized
421 // equations
422 // Note that linearize may throw for MSwells.
423 wellModel().linearize(ebosSimulator().model().linearizer().jacobian(),
424 ebosSimulator().model().linearizer().residual());
425
427 report.linear_solve_setup_time += linear_solve_setup_time_;
428 report.linear_solve_time += perfTimer.stop();
429 report.total_linear_iterations += linearIterationsLastSolve();
430 }
431 catch (...) {
432 report.linear_solve_setup_time += linear_solve_setup_time_;
433 report.linear_solve_time += perfTimer.stop();
434 report.total_linear_iterations += linearIterationsLastSolve();
435
436 failureReport_ += report;
437 throw; // re-throw up
438 }
439
440 perfTimer.reset();
441 perfTimer.start();
442
443 // handling well state update before oscillation treatment is a decision based
444 // on observation to avoid some big performance degeneration under some circumstances.
445 // there is no theorectical explanation which way is better for sure.
446 wellModel().postSolve(x);
447
448 if (param_.use_update_stabilization_) {
449 // Stabilize the nonlinear update.
450 bool isOscillate = false;
451 bool isStagnate = false;
452 nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
453 if (isOscillate) {
454 current_relaxation_ -= nonlinear_solver.relaxIncrement();
455 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
456 if (terminalOutputEnabled()) {
457 std::string msg = " Oscillating behavior detected: Relaxation set to "
458 + std::to_string(current_relaxation_);
459 OpmLog::info(msg);
460 }
461 }
462 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
463 }
464
465 // Apply the update, with considering model-dependent limitations and
466 // chopping of the update.
468
469 report.update_time += perfTimer.stop();
470 }
471
472 return report;
473 }
474
475 void printIf(int c, double x, double y, double eps, std::string type) {
476 if (std::abs(x-y) > eps) {
477 std::cout << type << " " <<c << ": "<<x << " " << y << std::endl;
478 }
479 }
480
481
486 {
488 Dune::Timer perfTimer;
489 perfTimer.start();
490 ebosSimulator_.problem().endTimeStep();
491 report.pre_post_time += perfTimer.stop();
492 return report;
493 }
494
500 const int iterationIdx)
501 {
502 // -------- Mass balance equations --------
503 ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
504 ebosSimulator_.problem().beginIteration();
505 ebosSimulator_.model().linearizer().linearizeDomain();
506 ebosSimulator_.problem().endIteration();
507
508 return wellModel().lastReport();
509 }
510
511 // compute the "relative" change of the solution between time steps
512 double relativeChange() const
513 {
514 Scalar resultDelta = 0.0;
515 Scalar resultDenom = 0.0;
516
517 const auto& elemMapper = ebosSimulator_.model().elementMapper();
518 const auto& gridView = ebosSimulator_.gridView();
519 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
520 unsigned globalElemIdx = elemMapper.index(elem);
521 const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
522
523 Scalar pressureNew;
524 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
525
526 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
527 Scalar oilSaturationNew = 1.0;
528 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
529 FluidSystem::numActivePhases() > 1 &&
530 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
531 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
532 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
533 }
534
535 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
536 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
537 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
538 assert(Indices::compositionSwitchIdx >= 0 );
539 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
540 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
541 }
542
543 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
544 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
545 }
546
547 const auto& priVarsOld = ebosSimulator_.model().solution(/*timeIdx=*/1)[globalElemIdx];
548
549 Scalar pressureOld;
550 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
551
552 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
553 Scalar oilSaturationOld = 1.0;
554
555 // NB fix me! adding pressures changes to satutation changes does not make sense
556 Scalar tmp = pressureNew - pressureOld;
557 resultDelta += tmp*tmp;
558 resultDenom += pressureNew*pressureNew;
559
560 if (FluidSystem::numActivePhases() > 1) {
561 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
562 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSwitchIdx];
563 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
564 }
565
566 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
567 {
568 assert(Indices::compositionSwitchIdx >= 0 );
569 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
570 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
571 }
572
573 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
574 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
575 }
576 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
577 Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
578 resultDelta += tmpSat*tmpSat;
579 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
580 assert(std::isfinite(resultDelta));
581 assert(std::isfinite(resultDenom));
582 }
583 }
584 }
585
586 resultDelta = gridView.comm().sum(resultDelta);
587 resultDenom = gridView.comm().sum(resultDenom);
588
589 if (resultDenom > 0.0)
590 return resultDelta/resultDenom;
591 return 0.0;
592 }
593
594
597 {
598 return ebosSimulator_.model().newtonMethod().linearSolver().iterations ();
599 }
600
603 void solveJacobianSystem(BVector& x)
604 {
605
606 auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
607 auto& ebosResid = ebosSimulator_.model().linearizer().residual();
608
609 // set initial guess
610 x = 0.0;
611
612 auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
613 Dune::Timer perfTimer;
614 perfTimer.start();
615 ebosSolver.prepare(ebosJac, ebosResid);
616 linear_solve_setup_time_ = perfTimer.stop();
617 ebosSolver.setResidual(ebosResid);
618 // actually, the error needs to be calculated after setResidual in order to
619 // account for parallelization properly. since the residual of ECFV
620 // discretizations does not need to be synchronized across processes to be
621 // consistent, this is not relevant for OPM-flow...
622 ebosSolver.setMatrix(ebosJac);
623 ebosSolver.solve(x);
624 }
625
626
627
629 void updateSolution(const BVector& dx)
630 {
632 auto& ebosNewtonMethod = ebosSimulator_.model().newtonMethod();
633 SolutionVector& solution = ebosSimulator_.model().solution(/*timeIdx=*/0);
634
635 ebosNewtonMethod.update_(/*nextSolution=*/solution,
636 /*curSolution=*/solution,
637 /*update=*/dx,
638 /*resid=*/dx); // the update routines of the black
639 // oil model do not care about the
640 // residual
641
642 // if the solution is updated, the intensive quantities need to be recalculated
643 {
645 ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
646 ebosSimulator_.problem().eclWriter()->mutableEclOutputModule().invalidateLocalData();
647 }
648 }
649
652 {
653 return terminal_output_;
654 }
655
656 template <class CollectiveCommunication>
657 std::tuple<double,double> convergenceReduction(const CollectiveCommunication& comm,
658 const double pvSumLocal,
659 const double numAquiferPvSumLocal,
660 std::vector< Scalar >& R_sum,
661 std::vector< Scalar >& maxCoeff,
662 std::vector< Scalar >& B_avg)
663 {
664 OPM_TIMEBLOCK(convergenceReduction);
665 // Compute total pore volume (use only owned entries)
666 double pvSum = pvSumLocal;
668
669 if( comm.size() > 1 )
670 {
671 // global reduction
672 std::vector< Scalar > sumBuffer;
673 std::vector< Scalar > maxBuffer;
674 const int numComp = B_avg.size();
675 sumBuffer.reserve( 2*numComp + 2 ); // +2 for (numAquifer)pvSum
676 maxBuffer.reserve( numComp );
677 for( int compIdx = 0; compIdx < numComp; ++compIdx )
678 {
679 sumBuffer.push_back( B_avg[ compIdx ] );
680 sumBuffer.push_back( R_sum[ compIdx ] );
681 maxBuffer.push_back( maxCoeff[ compIdx ] );
682 }
683
684 // Compute total pore volume
685 sumBuffer.push_back( pvSum );
686 sumBuffer.push_back( numAquiferPvSum );
687
688 // compute global sum
689 comm.sum( sumBuffer.data(), sumBuffer.size() );
690
691 // compute global max
692 comm.max( maxBuffer.data(), maxBuffer.size() );
693
694 // restore values to local variables
695 for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
696 {
697 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
698 ++buffIdx;
699
700 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
701 }
702
703 for( int compIdx = 0; compIdx < numComp; ++compIdx )
704 {
705 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
706 }
707
708 // restore global pore volume
709 pvSum = sumBuffer[sumBuffer.size()-2];
710 numAquiferPvSum = sumBuffer.back();
711 }
712
713 // return global pore volume
714 return {pvSum, numAquiferPvSum};
715 }
716
720 std::tuple<double,double> localConvergenceData(std::vector<Scalar>& R_sum,
721 std::vector<Scalar>& maxCoeff,
722 std::vector<Scalar>& B_avg)
723 {
725 double pvSumLocal = 0.0;
726 double numAquiferPvSumLocal = 0.0;
727 const auto& ebosModel = ebosSimulator_.model();
728 const auto& ebosProblem = ebosSimulator_.problem();
729
730 const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
731
732 ElementContext elemCtx(ebosSimulator_);
733 const auto& gridView = ebosSimulator().gridView();
734 OPM_BEGIN_PARALLEL_TRY_CATCH();
735 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
736 elemCtx.updatePrimaryStencil(elem);
737 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
738 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
739 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
740 const auto& fs = intQuants.fluidState();
741
742 const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx );
744
745 if (isNumericalAquiferCell(gridView.grid(), elem))
746 {
748 }
749
750 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
751 {
752 if (!FluidSystem::phaseIsActive(phaseIdx)) {
753 continue;
754 }
755
756 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
757
758 B_avg[ compIdx ] += 1.0 / fs.invB(phaseIdx).value();
759 const auto R2 = ebosResid[cell_idx][compIdx];
760
761 R_sum[ compIdx ] += R2;
762 maxCoeff[ compIdx ] = std::max( maxCoeff[ compIdx ], std::abs( R2 ) / pvValue );
763 }
764
765 if constexpr (has_solvent_) {
766 B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
767 const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
768 R_sum[ contiSolventEqIdx ] += R2;
769 maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue );
770 }
771 if constexpr (has_extbo_) {
772 B_avg[ contiZfracEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
773 const auto R2 = ebosResid[cell_idx][contiZfracEqIdx];
774 R_sum[ contiZfracEqIdx ] += R2;
775 maxCoeff[ contiZfracEqIdx ] = std::max( maxCoeff[ contiZfracEqIdx ], std::abs( R2 ) / pvValue );
776 }
777 if constexpr (has_polymer_) {
778 B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
779 const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
780 R_sum[ contiPolymerEqIdx ] += R2;
781 maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
782 }
783 if constexpr (has_foam_) {
784 B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
785 const auto R2 = ebosResid[cell_idx][contiFoamEqIdx];
786 R_sum[ contiFoamEqIdx ] += R2;
787 maxCoeff[ contiFoamEqIdx ] = std::max( maxCoeff[ contiFoamEqIdx ], std::abs( R2 ) / pvValue );
788 }
789 if constexpr (has_brine_) {
790 B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
791 const auto R2 = ebosResid[cell_idx][contiBrineEqIdx];
792 R_sum[ contiBrineEqIdx ] += R2;
793 maxCoeff[ contiBrineEqIdx ] = std::max( maxCoeff[ contiBrineEqIdx ], std::abs( R2 ) / pvValue );
794 }
795
796 if constexpr (has_polymermw_) {
797 static_assert(has_polymer_);
798
799 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
800 // the residual of the polymer molecular equation is scaled down by a 100, since molecular weight
801 // can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations
802 // TODO: there should be a more general way to determine the scaling-down coefficient
803 const auto R2 = ebosResid[cell_idx][contiPolymerMWEqIdx] / 100.;
804 R_sum[contiPolymerMWEqIdx] += R2;
805 maxCoeff[contiPolymerMWEqIdx] = std::max( maxCoeff[contiPolymerMWEqIdx], std::abs( R2 ) / pvValue );
806 }
807
808 if constexpr (has_energy_) {
809 B_avg[ contiEnergyEqIdx ] += 1.0 / (4.182e1); // converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water
810 const auto R2 = ebosResid[cell_idx][contiEnergyEqIdx];
811 R_sum[ contiEnergyEqIdx ] += R2;
812 maxCoeff[ contiEnergyEqIdx ] = std::max( maxCoeff[ contiEnergyEqIdx ], std::abs( R2 ) / pvValue );
813 }
814
815 if constexpr (has_micp_) {
816 B_avg[ contiMicrobialEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
817 const auto R1 = ebosResid[cell_idx][contiMicrobialEqIdx];
818 R_sum[ contiMicrobialEqIdx ] += R1;
819 maxCoeff[ contiMicrobialEqIdx ] = std::max( maxCoeff[ contiMicrobialEqIdx ], std::abs( R1 ) / pvValue );
820 B_avg[ contiOxygenEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
821 const auto R2 = ebosResid[cell_idx][contiOxygenEqIdx];
822 R_sum[ contiOxygenEqIdx ] += R2;
823 maxCoeff[ contiOxygenEqIdx ] = std::max( maxCoeff[ contiOxygenEqIdx ], std::abs( R2 ) / pvValue );
824 B_avg[ contiUreaEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
825 const auto R3 = ebosResid[cell_idx][contiUreaEqIdx];
826 R_sum[ contiUreaEqIdx ] += R3;
827 maxCoeff[ contiUreaEqIdx ] = std::max( maxCoeff[ contiUreaEqIdx ], std::abs( R3 ) / pvValue );
828 B_avg[ contiBiofilmEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
829 const auto R4 = ebosResid[cell_idx][contiBiofilmEqIdx];
830 R_sum[ contiBiofilmEqIdx ] += R4;
831 maxCoeff[ contiBiofilmEqIdx ] = std::max( maxCoeff[ contiBiofilmEqIdx ], std::abs( R4 ) / pvValue );
832 B_avg[ contiCalciteEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
833 const auto R5 = ebosResid[cell_idx][contiCalciteEqIdx];
834 R_sum[ contiCalciteEqIdx ] += R5;
835 maxCoeff[ contiCalciteEqIdx ] = std::max( maxCoeff[ contiCalciteEqIdx ], std::abs( R5 ) / pvValue );
836 }
837 }
838
839 OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm());
840
841 // compute local average in terms of global number of elements
842 const int bSize = B_avg.size();
843 for ( int i = 0; i<bSize; ++i )
844 {
845 B_avg[ i ] /= Scalar( global_nc_ );
846 }
847
849 }
850
853 double computeCnvErrorPv(const std::vector<Scalar>& B_avg, double dt)
854 {
856 double errorPV{};
857 const auto& ebosModel = ebosSimulator_.model();
858 const auto& ebosProblem = ebosSimulator_.problem();
859 const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
860 const auto& gridView = ebosSimulator().gridView();
861 ElementContext elemCtx(ebosSimulator_);
862
863 OPM_BEGIN_PARALLEL_TRY_CATCH();
864
865 for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder))
866 {
867 // Skip cells of numerical Aquifer
868 if (isNumericalAquiferCell(gridView.grid(), elem))
869 {
870 continue;
871 }
872 elemCtx.updatePrimaryStencil(elem);
873 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
874 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
875 const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx );
876 const auto& cellResidual = ebosResid[cell_idx];
877 bool cnvViolated = false;
878
879 for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx)
880 {
881 using std::abs;
883 cnvViolated = cnvViolated || (abs(CNV) > param_.tolerance_cnv_);
884 }
885
886 if (cnvViolated)
887 {
888 errorPV += pvValue;
889 }
890 }
891
892 OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::ComputeCnvError() failed: ", grid_.comm());
893
894 return grid_.comm().sum(errorPV);
895 }
896
897 ConvergenceReport getReservoirConvergence(const double reportTime,
898 const double dt,
899 const int iteration,
900 std::vector<Scalar>& B_avg,
901 std::vector<Scalar>& residual_norms)
902 {
903 OPM_TIMEBLOCK(getReservoirConvergence);
904 typedef std::vector< Scalar > Vector;
905
906 const int numComp = numEq;
907 Vector R_sum(numComp, 0.0 );
908 Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
910
911 // compute global sum and max of quantities
912 const auto [ pvSum, numAquiferPvSum ] =
913 convergenceReduction(grid_.comm(), pvSumLocal,
916
919
920 const double tol_mb = param_.tolerance_mb_;
921 // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0.
922 // For each iteration, we need to determine whether to use the relaxed CNV tolerance.
923 // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0.
924 const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.min_strict_cnv_iter_;
925 const double tol_cnv = use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
926
927 // Finish computation
928 std::vector<Scalar> CNV(numComp);
929 std::vector<Scalar> mass_balance_residual(numComp);
930 for ( int compIdx = 0; compIdx < numComp; ++compIdx )
931 {
934 residual_norms.push_back(CNV[compIdx]);
935 }
936
937 // Create convergence report.
938 ConvergenceReport report{reportTime};
939 using CR = ConvergenceReport;
940 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
941 double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
942 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
943 CR::ReservoirFailure::Type::Cnv };
944 double tol[2] = { tol_mb, tol_cnv };
945 for (int ii : {0, 1}) {
946 if (std::isnan(res[ii])) {
947 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
948 if ( terminal_output_ ) {
949 OpmLog::debug("NaN residual for " + this->compNames_.name(compIdx) + " equation.");
950 }
951 } else if (res[ii] > maxResidualAllowed()) {
952 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
953 if ( terminal_output_ ) {
954 OpmLog::debug("Too large residual for " + this->compNames_.name(compIdx) + " equation.");
955 }
956 } else if (res[ii] < 0.0) {
957 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
958 if ( terminal_output_ ) {
959 OpmLog::debug("Negative residual for " + this->compNames_.name(compIdx) + " equation.");
960 }
961 } else if (res[ii] > tol[ii]) {
962 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
963 }
964 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
965 }
966 }
967
968 // Output of residuals.
969 if ( terminal_output_ )
970 {
971 // Only rank 0 does print to std::cout
972 if (iteration == 0) {
973 std::string msg = "Iter";
974 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
975 msg += " MB(";
976 msg += this->compNames_.name(compIdx)[0];
977 msg += ") ";
978 }
979 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
980 msg += " CNV(";
981 msg += this->compNames_.name(compIdx)[0];
982 msg += ") ";
983 }
984 OpmLog::debug(msg);
985 }
986 std::ostringstream ss;
987 const std::streamsize oprec = ss.precision(3);
988 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
989 ss << std::setw(4) << iteration;
990 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
991 ss << std::setw(11) << mass_balance_residual[compIdx];
992 }
993 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
994 ss << std::setw(11) << CNV[compIdx];
995 }
996 ss.precision(oprec);
997 ss.flags(oflags);
998 OpmLog::debug(ss.str());
999 }
1000
1001 return report;
1002 }
1003
1010 const int iteration,
1011 std::vector<double>& residual_norms)
1012 {
1014 // Get convergence reports for reservoir and wells.
1015 std::vector<Scalar> B_avg(numEq, 0.0);
1016 auto report = getReservoirConvergence(timer.simulationTimeElapsed(),
1017 timer.currentStepLength(),
1019 {
1020 OPM_TIMEBLOCK(getWellConvergence);
1021 report += wellModel().getWellConvergence(B_avg, /*checkWellGroupControls*/report.converged());
1022 }
1023 return report;
1024 }
1025
1026
1028 int numPhases() const
1029 {
1030 return phaseUsage_.num_phases;
1031 }
1032
1034 template<class T>
1035 std::vector<std::vector<double> >
1036 computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
1037 {
1039 }
1040
1042 std::vector<std::vector<double> >
1043 computeFluidInPlace(const std::vector<int>& /*fipnum*/) const
1044 {
1046 //assert(true)
1047 //return an empty vector
1048 std::vector<std::vector<double> > regionValues(0, std::vector<double>(0,0.0));
1049 return regionValues;
1050 }
1051
1052 const Simulator& ebosSimulator() const
1053 { return ebosSimulator_; }
1054
1055 Simulator& ebosSimulator()
1056 { return ebosSimulator_; }
1057
1060 { return failureReport_; }
1061
1062 const std::vector<StepReport>& stepReports() const
1063 {
1064 return convergence_reports_;
1065 }
1066
1067 std::vector<StepReport> getStepReportsDestructively() const
1068 {
1069 return std::move(this->convergence_reports_);
1070 }
1071
1072 protected:
1073 // --------- Data members ---------
1074
1075 Simulator& ebosSimulator_;
1076 const Grid& grid_;
1077 const PhaseUsage phaseUsage_;
1078 static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
1079 static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
1080 static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
1081 static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
1082 static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
1083 static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
1084 static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
1085 static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
1086
1087 ModelParameters param_;
1088 SimulatorReportSingle failureReport_;
1089
1090 // Well Model
1091 BlackoilWellModel<TypeTag>& well_model_;
1092
1096 long int global_nc_;
1097
1098 std::vector<std::vector<double>> residual_norms_history_;
1099 double current_relaxation_;
1100 BVector dx_old_;
1101
1102 std::vector<StepReport> convergence_reports_;
1103 ComponentName compNames_{};
1104
1105 public:
1108 wellModel() { return well_model_; }
1109
1111 wellModel() const { return well_model_; }
1112
1113 void beginReportStep()
1114 {
1115 ebosSimulator_.problem().beginEpisode();
1116 }
1117
1118 void endReportStep()
1119 {
1120 ebosSimulator_.problem().endEpisode();
1121 }
1122
1123 private:
1124 template<class T>
1125 bool isNumericalAquiferCell(const Dune::CpGrid& grid, const T& elem)
1126 {
1127 const auto& aquiferCells = grid.sortedNumAquiferCells();
1128 if (aquiferCells.empty())
1129 {
1130 return false;
1131 }
1132 auto candidate = std::lower_bound(aquiferCells.begin(), aquiferCells.end(),
1133 elem.index());
1134 return candidate != aquiferCells.end() && *candidate == elem.index();
1135 }
1136
1137 template<class G, class T>
1138 typename std::enable_if<!std::is_same<G,Dune::CpGrid>::value, bool>::type
1139 isNumericalAquiferCell(const G&, const T&)
1140 {
1141 return false;
1142 }
1143
1144 double dpMaxRel() const { return param_.dp_max_rel_; }
1145 double dsMax() const { return param_.ds_max_; }
1146 double drMaxRel() const { return param_.dr_max_rel_; }
1147 double maxResidualAllowed() const { return param_.max_residual_allowed_; }
1148 double linear_solve_setup_time_;
1149
1150 public:
1151 std::vector<bool> wasSwitched_;
1152 };
1153} // namespace Opm
1154
1155#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
Definition AquiferInterface.hpp:35
Definition BlackoilModelEbos.hpp:206
A model implementation for three-phase black oil.
Definition BlackoilModelEbos.hpp:155
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition BlackoilModelEbos.hpp:1108
std::vector< std::vector< double > > computeFluidInPlace(const std::vector< int > &) const
Should not be called.
Definition BlackoilModelEbos.hpp:1043
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition BlackoilModelEbos.hpp:596
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition BlackoilModelEbos.hpp:651
const SimulatorReportSingle & failureReport() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModelEbos.hpp:1059
ConvergenceReport getConvergence(const SimulatorTimerInterface &timer, const int iteration, std::vector< double > &residual_norms)
Compute convergence based on total mass balance (tol_mb) and maximum residual mass balance (tol_cnv).
Definition BlackoilModelEbos.hpp:1009
BlackoilModelEbos(Simulator &ebosSimulator, const ModelParameters &param, BlackoilWellModel< TypeTag > &well_model, const bool terminal_output)
Construct the model.
Definition BlackoilModelEbos.hpp:282
bool terminal_output_
Whether we print something to std::cout.
Definition BlackoilModelEbos.hpp:1094
SimulatorReportSingle afterStep(const SimulatorTimerInterface &)
Called once after each time step.
Definition BlackoilModelEbos.hpp:485
std::tuple< double, double > localConvergenceData(std::vector< Scalar > &R_sum, std::vector< Scalar > &maxCoeff, std::vector< Scalar > &B_avg)
Get reservoir quantities on this process needed for convergence calculations.
Definition BlackoilModelEbos.hpp:720
SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Called once per nonlinear iteration.
Definition BlackoilModelEbos.hpp:354
SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &, const int iterationIdx)
Assemble the residual and Jacobian of the nonlinear system.
Definition BlackoilModelEbos.hpp:499
SimulatorReportSingle prepareStep(const SimulatorTimerInterface &timer)
Called once before each time step.
Definition BlackoilModelEbos.hpp:308
long int global_nc_
The number of cells of the global grid.
Definition BlackoilModelEbos.hpp:1096
double computeCnvErrorPv(const std::vector< Scalar > &B_avg, double dt)
Compute the total pore volume of cells violating CNV that are not part of a numerical aquifer.
Definition BlackoilModelEbos.hpp:853
int numPhases() const
The number of active fluid phases in the model.
Definition BlackoilModelEbos.hpp:1028
std::vector< std::vector< double > > computeFluidInPlace(const T &, const std::vector< int > &fipnum) const
Wrapper required due to not following generic API.
Definition BlackoilModelEbos.hpp:1036
void solveJacobianSystem(BVector &x)
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition BlackoilModelEbos.hpp:603
void updateSolution(const BVector &dx)
Apply an update to the primary variables.
Definition BlackoilModelEbos.hpp:629
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
Interface class for SimulatorTimer objects, to be improved.
Definition SimulatorTimerInterface.hpp:34
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition phaseUsageFromDeck.cpp:145
Definition ISTLSolverEbos.hpp:68
Definition BlackoilModelEbos.hpp:72
Definition ISTLSolverEbos.hpp:62
Definition BlackoilModelParametersEbos.hpp:31
Definition NonlinearSolverEbos.hpp:41
Definition AdaptiveTimeSteppingEbos.hpp:49
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34