24 #ifndef OPM_BLACKOILMODELEBOS_HEADER_INCLUDED 25 #define OPM_BLACKOILMODELEBOS_HEADER_INCLUDED 27 #include <ebos/eclproblem.hh> 28 #include <ewoms/common/start.hh> 30 #include <opm/autodiff/BlackoilModelParameters.hpp> 31 #include <opm/autodiff/BlackoilWellModel.hpp> 32 #include <opm/autodiff/AutoDiffBlock.hpp> 33 #include <opm/autodiff/AutoDiffHelpers.hpp> 34 #include <opm/autodiff/GridHelpers.hpp> 35 #include <opm/autodiff/WellHelpers.hpp> 36 #include <opm/autodiff/GeoProps.hpp> 37 #include <opm/autodiff/WellDensitySegmented.hpp> 38 #include <opm/autodiff/VFPProperties.hpp> 39 #include <opm/autodiff/VFPProdProperties.hpp> 40 #include <opm/autodiff/VFPInjProperties.hpp> 41 #include <opm/autodiff/BlackoilDetails.hpp> 42 #include <opm/autodiff/BlackoilModelEnums.hpp> 43 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp> 46 #include <opm/core/grid.h> 47 #include <opm/core/simulator/SimulatorReport.hpp> 48 #include <opm/core/linalg/LinearSolverInterface.hpp> 49 #include <opm/core/linalg/ParallelIstlInformation.hpp> 50 #include <opm/core/props/phaseUsageFromDeck.hpp> 51 #include <opm/common/ErrorMacros.hpp> 52 #include <opm/common/Exceptions.hpp> 53 #include <opm/common/OpmLog/OpmLog.hpp> 54 #include <opm/parser/eclipse/Units/Units.hpp> 55 #include <opm/core/well_controls.h> 56 #include <opm/simulators/timestepping/SimulatorTimer.hpp> 57 #include <opm/core/utility/parameters/ParameterGroup.hpp> 58 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp> 59 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp> 61 #include <opm/autodiff/ISTLSolver.hpp> 62 #include <opm/common/data/SimulationDataContainer.hpp> 64 #include <dune/istl/owneroverlapcopy.hh> 65 #include <dune/common/parallel/collectivecommunication.hh> 66 #include <dune/common/timer.hh> 67 #include <dune/common/unused.hh> 81 namespace Properties {
82 NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
83 SET_BOOL_PROP(EclFlowProblem, DisableWells,
true);
84 SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks,
false);
85 SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility,
true);
87 SET_BOOL_PROP(EclFlowProblem, BlackoilConserveSurfaceVolume,
true);
88 SET_BOOL_PROP(EclFlowProblem, UseVolumetricResidual,
false);
93 SET_BOOL_PROP(EclFlowProblem, EnableSwatinit,
false);
103 template <
class TypeTag>
108 typedef BlackoilState ReservoirState;
112 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
113 typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
114 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
115 typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
116 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables ;
117 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
118 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
119 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
120 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
122 typedef double Scalar;
123 static const int numEq = Indices::numEq;
124 static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
125 static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
126 static const int solventSaturationIdx = Indices::solventSaturationIdx;
127 static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
129 typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
130 typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
131 typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
132 typedef Dune::BlockVector<VectorBlockType> BVector;
139 SurfaceToReservoirVoidage<BlackoilPropsAdFromDeck::FluidSystem, std::vector<int> >;
160 const bool terminal_output
162 : ebosSimulator_(ebosSimulator)
163 , grid_(ebosSimulator_.gridManager().grid())
164 , istlSolver_( dynamic_cast< const
ISTLSolverType* > (&linsolver) )
165 , phaseUsage_(phaseUsageFromDeck(eclState()))
167 eclState().getTableManager().getVFPInjTables(),
168 eclState().getTableManager().getVFPProdTables())
169 , active_(detail::activePhases(phaseUsage_))
170 , has_disgas_(FluidSystem::enableDissolvedGas())
171 , has_vapoil_(FluidSystem::enableVaporizedOil())
172 , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
173 , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
175 , well_model_ (well_model)
177 , rate_converter_(rate_converter)
178 , current_relaxation_(1.0)
179 , dx_old_(AutoDiffGrid::numCells(grid_))
180 , isBeginReportStep_(false)
185 wellsActive = grid_.comm().max(wellsActive);
186 wellModel().setWellsActive( wellsActive );
189 wellModel().setVFPProperties(&vfp_properties_);
192 OPM_THROW(std::logic_error,
"solver down cast to ISTLSolver failed");
196 bool isParallel()
const 197 {
return grid_.comm().size() > 1; }
199 const EclipseState& eclState()
const 200 {
return ebosSimulator_.gridManager().eclState(); }
207 const ReservoirState& ,
210 if (
wellModel().wellCollection()->havingVREPGroups() ) {
211 updateRateConverter();
219 ebosSimulator_.model().solution( 0 ) = ebosSimulator_.model().solution( 1 );
220 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
223 ebosSimulator_.model().solution( 1 ) = ebosSimulator_.model().solution( 0 );
226 unsigned numDof = ebosSimulator_.model().numGridDof();
227 wasSwitched_.resize(numDof);
228 std::fill(wasSwitched_.begin(), wasSwitched_.end(),
false);
241 template <
class NonlinearSolverType>
244 NonlinearSolverType& nonlinear_solver,
248 SimulatorReport report;
249 failureReport_ = SimulatorReport();
250 Dune::Timer perfTimer;
253 if (iteration == 0) {
256 residual_norms_history_.clear();
257 current_relaxation_ = 1.0;
261 report.total_linearizations = 1;
264 report +=
assemble(timer, iteration, well_state);
265 report.assemble_time += perfTimer.stop();
268 report.assemble_time += perfTimer.stop();
269 failureReport_ += report;
274 std::vector<double> residual_norms;
278 report.converged =
getConvergence(timer, iteration,residual_norms) && iteration > nonlinear_solver.minIter();
281 if (
wellModel().wellCollection()->groupControlActive()) {
282 report.converged = report.converged &&
wellModel().wellCollection()->groupTargetConverged(well_state.wellRates());
285 report.update_time += perfTimer.stop();
286 residual_norms_history_.push_back(residual_norms);
287 if (!report.converged) {
290 report.total_newton_iterations = 1;
296 const int nc = AutoDiffGrid::numCells(grid_);
297 const int nw = numWells();
302 report.linear_solve_time += perfTimer.stop();
306 report.linear_solve_time += perfTimer.stop();
309 failureReport_ += report;
322 wellModel().recoverWellSolutionAndUpdateWellState(x, well_state);
327 bool isOscillate =
false;
328 bool isStagnate =
false;
329 nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
331 current_relaxation_ -= nonlinear_solver.relaxIncrement();
332 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
334 std::string msg =
" Oscillating behavior detected: Relaxation set to " 335 + std::to_string(current_relaxation_);
339 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
346 report.update_time += perfTimer.stop();
352 void printIf(
int c,
double x,
double y,
double eps, std::string type) {
353 if (std::abs(x-y) > eps) {
354 std::cout << type <<
" " <<c <<
": "<<x <<
" " << y << std::endl;
365 const ReservoirState& reservoir_state,
368 DUNE_UNUSED_PARAMETER(timer);
369 DUNE_UNUSED_PARAMETER(reservoir_state);
370 DUNE_UNUSED_PARAMETER(well_state);
378 const int iterationIdx,
383 SimulatorReport report;
386 if (
wellModel().wellCollection()->havingVREPGroups() ) {
387 updateRateConverter();
391 assembleMassBalanceEq(timer, iterationIdx);
400 report =
wellModel().assemble(ebosSimulator_, iterationIdx, dt, well_state);
402 catch (
const Dune::FMatrixError& e )
404 OPM_THROW(Opm::NumericalProblem,
"Error encounted when solving well equations");
411 template <
class Dummy>
412 double relativeChange(
const Dummy&,
const Dummy&)
const 414 Scalar resultDelta = 0.0;
415 Scalar resultDenom = 0.0;
417 const auto& elemMapper = ebosSimulator_.model().elementMapper();
418 const auto& gridView = ebosSimulator_.gridView();
419 auto elemIt = gridView.template begin<0>();
420 const auto& elemEndIt = gridView.template end<0>();
421 for (; elemIt != elemEndIt; ++elemIt) {
422 const auto& elem = *elemIt;
423 if (elem.partitionType() != Dune::InteriorEntity)
426 unsigned globalElemIdx = elemMapper.index(elem);
427 const auto& priVarsNew = ebosSimulator_.model().solution(0)[globalElemIdx];
430 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
432 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
433 Scalar oilSaturationNew = 1.0;
434 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
435 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx];
436 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
439 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
440 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
441 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
444 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
445 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
448 const auto& priVarsOld = ebosSimulator_.model().solution(1)[globalElemIdx];
451 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
453 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
454 Scalar oilSaturationOld = 1.0;
455 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
456 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
457 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
460 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
461 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
462 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
465 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
466 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
469 Scalar tmp = pressureNew - pressureOld;
470 resultDelta += tmp*tmp;
471 resultDenom += pressureNew*pressureNew;
473 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
474 Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
475 resultDelta += tmp*tmp;
476 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
480 resultDelta = gridView.comm().sum(resultDelta);
481 resultDenom = gridView.comm().sum(resultDenom);
483 if (resultDenom > 0.0)
484 return resultDelta/resultDenom;
492 const int nc = Opm::AutoDiffGrid::numCells(grid_);
493 const int nw = numWells();
494 return numComponents() * (nc + nw);
507 const auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
508 auto& ebosResid = ebosSimulator_.model().linearizer().residual();
529 Operator opA(ebosJac, well_model_, istlSolver().parallelInformation() );
530 assert( opA.comm() );
531 istlSolver().
solve( opA, x, ebosResid, *(opA.comm()) );
536 Operator opA(ebosJac, well_model_);
537 istlSolver().
solve( opA, x, ebosResid );
550 template<
class M,
class X,
class Y,
class WellModel,
bool overlapping >
553 typedef Dune::AssembledLinearOperator<M,X,Y> BaseType;
557 typedef X domain_type;
558 typedef Y range_type;
559 typedef typename X::field_type field_type;
562 typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
564 typedef Dune::CollectiveCommunication< Grid > communication_type;
570 Dune::SolverCategory::overlapping :
571 Dune::SolverCategory::sequential
576 : A_( A ), wellMod_( wellMod ), comm_()
579 if( parallelInformation.type() ==
typeid(ParallelISTLInformation) )
581 const ParallelISTLInformation& info =
582 boost::any_cast<
const ParallelISTLInformation&>( parallelInformation);
583 comm_.reset(
new communication_type( info.communicator() ) );
588 virtual void apply(
const X& x, Y& y )
const 592 wellMod_.apply(x, y );
601 virtual void applyscaleadd (field_type alpha,
const X& x, Y& y)
const 605 wellMod_.applyScaleAdd( alpha, x, y );
613 virtual const matrix_type& getmat()
const {
return A_; }
615 communication_type* comm()
617 return comm_.operator->();
621 const matrix_type& A_ ;
622 const WellModel& wellMod_;
623 std::unique_ptr< communication_type > comm_;
634 const auto& ebosProblem = ebosSimulator_.problem();
636 unsigned numSwitched = 0;
637 ElementContext elemCtx( ebosSimulator_ );
638 const auto& gridView = ebosSimulator_.gridView();
639 const auto& elemEndIt = gridView.template end<0>();
640 SolutionVector& solution = ebosSimulator_.model().solution( 0 );
642 for (
auto elemIt = gridView.template begin</*codim=*/0>();
646 const auto& elem = *elemIt;
648 elemCtx.updatePrimaryStencil(elem);
649 elemCtx.updatePrimaryIntensiveQuantities(0);
650 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
652 PrimaryVariables& priVars = solution[ cell_idx ];
654 const double& dp = dx[cell_idx][Indices::pressureSwitchIdx];
655 double& p = priVars[Indices::pressureSwitchIdx];
656 const double& dp_rel_max = dpMaxRel();
657 const int sign_dp = dp > 0 ? 1: -1;
658 p -= sign_dp * std::min(std::abs(dp), std::abs(p)*dp_rel_max);
659 p = std::max(p, 0.0);
662 const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0;
663 const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0;
671 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
674 else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
678 assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv);
684 const double dss = has_solvent_ ? dx[cell_idx][Indices::solventSaturationIdx] : 0.0;
687 const double dc = has_polymer_ ? dx[cell_idx][Indices::polymerConcentrationIdx] : 0.0;
690 dso = - (dsw + dsg + dss);
695 maxVal = std::max(std::abs(dsw),maxVal);
696 maxVal = std::max(std::abs(dsg),maxVal);
697 maxVal = std::max(std::abs(dso),maxVal);
698 maxVal = std::max(std::abs(dss),maxVal);
700 double satScaleFactor = 1.0;
701 if (maxVal > dsMax()) {
702 satScaleFactor = dsMax()/maxVal;
705 if (active_[Water]) {
706 double& sw = priVars[Indices::waterSaturationIdx];
707 sw -= satScaleFactor * dsw;
711 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
712 double& sg = priVars[Indices::compositionSwitchIdx];
713 sg -= satScaleFactor * dsg;
718 double& ss = priVars[Indices::solventSaturationIdx];
719 ss -= satScaleFactor * dss;
720 ss = std::min(std::max(ss, 0.0),1.0);
723 double& c = priVars[Indices::polymerConcentrationIdx];
724 c -= satScaleFactor * dc;
725 c = std::max(c, 0.0);
729 if (active_[Gas] && active_[Oil] ) {
730 unsigned pvtRegionIdx = ebosSimulator_.problem().pvtRegionIndex(cell_idx);
731 const double drmaxrel = drMaxRel();
733 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
735 FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx, 300.0, p);
737 double& rs = priVars[Indices::compositionSwitchIdx];
738 rs -= ((drs<0)?-1:1)*std::min(std::abs(drs), RsSat*drmaxrel);
739 rs = std::max(rs, 0.0);
744 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
746 FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx, 300.0, p);
748 double& rv = priVars[Indices::compositionSwitchIdx];
749 rv -= ((drv<0)?-1:1)*std::min(std::abs(drv), RvSat*drmaxrel);
750 rv = std::max(rv, 0.0);
756 if (wasSwitched_[cell_idx])
757 wasSwitched_[cell_idx] = priVars.adaptPrimaryVariables(ebosProblem, cell_idx, 1e-5);
759 wasSwitched_[cell_idx] = priVars.adaptPrimaryVariables(ebosProblem, cell_idx);
761 if (wasSwitched_[cell_idx])
766 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
776 template <
class CollectiveCommunication>
777 double convergenceReduction(
const CollectiveCommunication& comm,
778 const double pvSumLocal,
779 std::vector< Scalar >& R_sum,
780 std::vector< Scalar >& maxCoeff,
781 std::vector< Scalar >& B_avg)
784 double pvSum = pvSumLocal;
786 if( comm.size() > 1 )
789 std::vector< Scalar > sumBuffer;
790 std::vector< Scalar > maxBuffer;
791 const int numComp = B_avg.size();
792 sumBuffer.reserve( 2*numComp + 1 );
793 maxBuffer.reserve( numComp );
794 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
796 sumBuffer.push_back( B_avg[ compIdx ] );
797 sumBuffer.push_back( R_sum[ compIdx ] );
798 maxBuffer.push_back( maxCoeff[ compIdx ] );
802 sumBuffer.push_back( pvSum );
805 comm.sum( sumBuffer.data(), sumBuffer.size() );
808 comm.max( maxBuffer.data(), maxBuffer.size() );
811 for(
int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
813 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
816 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
819 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
821 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
825 pvSum = sumBuffer.back();
839 typedef std::vector< Scalar > Vector;
846 const int numComp = numComponents();
848 Vector R_sum(numComp, 0.0 );
849 Vector B_avg(numComp, 0.0 );
850 Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
852 const auto& ebosModel = ebosSimulator_.model();
853 const auto& ebosProblem = ebosSimulator_.problem();
855 const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
857 ElementContext elemCtx(ebosSimulator_);
858 const auto& gridView = ebosSimulator().gridView();
859 const auto& elemEndIt = gridView.template end<0, Dune::Interior_Partition>();
861 double pvSumLocal = 0.0;
862 for (
auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
866 const auto& elem = *elemIt;
867 elemCtx.updatePrimaryStencil(elem);
868 elemCtx.updatePrimaryIntensiveQuantities(0);
869 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
870 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
871 const auto& fs = intQuants.fluidState();
873 const double pvValue = ebosProblem.porosity(cell_idx) * ebosModel.dofTotalVolume( cell_idx );
874 pvSumLocal += pvValue;
876 for (
int phaseIdx = 0; phaseIdx < np; ++phaseIdx )
878 const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx);
879 const int ebosCompIdx = flowPhaseToEbosCompIdx(phaseIdx);
881 B_avg[ phaseIdx ] += 1.0 / fs.invB(ebosPhaseIdx).value();
882 const auto R2 = ebosResid[cell_idx][ebosCompIdx];
884 R_sum[ phaseIdx ] += R2;
885 maxCoeff[ phaseIdx ] = std::max( maxCoeff[ phaseIdx ], std::abs( R2 ) / pvValue );
888 if ( has_solvent_ ) {
889 B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
890 const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
891 R_sum[ contiSolventEqIdx ] += R2;
892 maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue );
895 B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
896 const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
897 R_sum[ contiPolymerEqIdx ] += R2;
898 maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
904 const int bSize = B_avg.size();
905 for (
int i = 0; i<bSize; ++i )
915 const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal,
916 R_sum, maxCoeff, B_avg);
919 Vector mass_balance_residual(numComp);
921 bool converged_MB =
true;
922 bool converged_CNV =
true;
924 for (
int compIdx = 0; compIdx < numComp; ++compIdx )
926 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
927 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
928 converged_MB = converged_MB && (mass_balance_residual[compIdx] < tol_mb);
929 converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv);
931 residual_norms.push_back(CNV[compIdx]);
934 const bool converged_Well =
wellModel().getWellConvergence(ebosSimulator_, B_avg);
936 bool converged = converged_MB && converged_Well;
941 converged = converged && converged_CNV;
946 if (iteration == 0) {
947 std::string msg =
"Iter";
949 std::vector< std::string > key( numComp );
950 for (
int phaseIdx = 0; phaseIdx <
numPhases(); ++phaseIdx) {
951 const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
952 key[ phaseIdx ] = std::toupper( phaseName.front() );
955 key[ solventSaturationIdx ] =
"S";
959 key[ polymerConcentrationIdx ] =
"P";
962 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
963 msg +=
" MB(" + key[ compIdx ] +
") ";
965 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
966 msg +=
" CNV(" + key[ compIdx ] +
") ";
970 std::ostringstream ss;
971 const std::streamsize oprec = ss.precision(3);
972 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
973 ss << std::setw(4) << iteration;
974 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
975 ss << std::setw(11) << mass_balance_residual[compIdx];
977 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
978 ss << std::setw(11) << CNV[compIdx];
982 OpmLog::note(ss.str());
985 for (
int phaseIdx = 0; phaseIdx <
numPhases(); ++phaseIdx) {
986 const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
988 if (std::isnan(mass_balance_residual[phaseIdx])
989 || std::isnan(CNV[phaseIdx])) {
990 OPM_THROW(Opm::NumericalProblem,
"NaN residual for phase " << phaseName);
992 if (mass_balance_residual[phaseIdx] > maxResidualAllowed()
993 || CNV[phaseIdx] > maxResidualAllowed()) {
994 OPM_THROW(Opm::NumericalProblem,
"Too large residual for phase " << phaseName);
1005 return phaseUsage_.num_phases;
1008 int numComponents()
const 1013 int numComp = FluidSystem::numComponents;
1025 std::vector<std::vector<double> >
1031 std::vector<std::vector<double> >
1034 const auto& comm = grid_.comm();
1035 const auto& gridView = ebosSimulator().gridView();
1036 const int nc = gridView.size(0);
1037 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
1038 int ntFip = *std::max_element(fipnum.begin(), fipnum.end());
1039 ntFip = comm.max(ntFip);
1041 std::vector<double> tpv(ntFip, 0.0);
1042 std::vector<double> hcpv(ntFip, 0.0);
1044 std::vector<std::vector<double> > regionValues(ntFip, std::vector<double>(FIPDataType::fipValues,0.0));
1046 for (
int i = 0; i<FIPDataType::fipValues; i++) {
1047 fip_.fip[i].resize(nc,0.0);
1050 ElementContext elemCtx(ebosSimulator_);
1051 const auto& elemEndIt = gridView.template end<0, Dune::Interior_Partition>();
1052 for (
auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
1053 elemIt != elemEndIt;
1056 elemCtx.updatePrimaryStencil(*elemIt);
1057 elemCtx.updatePrimaryIntensiveQuantities(0);
1059 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
1060 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
1061 const auto& fs = intQuants.fluidState();
1063 const int regionIdx = fipnum[cellIdx] - 1;
1064 if (regionIdx < 0) {
1077 ebosSimulator_.model().dofTotalVolume(cellIdx)
1078 * intQuants.porosity().value();
1080 for (
int phase = 0; phase < maxnp; ++phase) {
1081 const double b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value();
1082 const double s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value();
1084 fip_.fip[phase][cellIdx] = b * s * pv;
1086 if (active_[ phase ]) {
1087 regionValues[regionIdx][phase] += fip_.fip[phase][cellIdx];
1091 if (active_[ Oil ] && active_[ Gas ]) {
1093 fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][cellIdx] = fs.Rs().value() * fip_.fip[FIPDataType::FIP_LIQUID][cellIdx];
1094 fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][cellIdx] = fs.Rv().value() * fip_.fip[FIPDataType::FIP_VAPOUR][cellIdx];
1096 regionValues[regionIdx][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][cellIdx];
1097 regionValues[regionIdx][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][cellIdx];
1100 const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
1102 tpv[regionIdx] += pv;
1103 hcpv[regionIdx] += pv * hydrocarbon;
1108 comm.sum(tpv.data(), tpv.size());
1109 comm.sum(hcpv.data(), hcpv.size());
1111 for (
auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
1112 elemIt != elemEndIt;
1115 const auto& elem = *elemIt;
1117 elemCtx.updatePrimaryStencil(elem);
1118 elemCtx.updatePrimaryIntensiveQuantities(0);
1120 unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
1121 const int regionIdx = fipnum[cellIdx] - 1;
1122 if (regionIdx < 0) {
1127 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
1128 const auto& fs = intQuants.fluidState();
1138 ebosSimulator_.model().dofTotalVolume(cellIdx)
1139 * intQuants.porosity().value();
1141 fip_.fip[FIPDataType::FIP_PV][cellIdx] = pv;
1142 const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
1146 if (hcpv[regionIdx] > 1e-10) {
1147 fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[regionIdx];
1149 fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() / tpv[regionIdx];
1152 regionValues[regionIdx][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][cellIdx];
1153 regionValues[regionIdx][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx];
1157 for(
int regionIdx=0; regionIdx < ntFip; ++regionIdx) {
1158 comm.sum(regionValues[regionIdx].data(), regionValues[regionIdx].size());
1161 return regionValues;
1164 SimulationDataContainer getSimulatorData (
const SimulationDataContainer& localState)
const 1166 typedef std::vector<double> VectorType;
1168 const auto& ebosModel = ebosSimulator().model();
1169 const auto& phaseUsage = phaseUsage_;
1172 const int numCells = ebosModel.numGridDof();
1175 SimulationDataContainer simData( numCells, 0, num_phases );
1178 const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
1179 const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
1180 const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
1182 const int aqua_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Aqua ];
1183 const int liquid_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Liquid ];
1184 const int vapour_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Vapour ];
1188 VectorType& pressureOil = simData.pressure();
1189 VectorType& temperature = simData.temperature();
1190 VectorType& saturation = simData.saturation();
1194 simData.registerCellData(
"1OVERBW", 1 );
1195 simData.registerCellData(
"WAT_DEN", 1 );
1196 simData.registerCellData(
"WAT_VISC", 1 );
1197 simData.registerCellData(
"WATKR", 1 );
1200 VectorType& bWater = aqua_active ? simData.getCellData(
"1OVERBW" ) : zero;
1201 VectorType& rhoWater = aqua_active ? simData.getCellData(
"WAT_DEN" ) : zero;
1202 VectorType& muWater = aqua_active ? simData.getCellData(
"WAT_VISC" ) : zero;
1203 VectorType& krWater = aqua_active ? simData.getCellData(
"WATKR" ) : zero;
1206 if( liquid_active ) {
1207 simData.registerCellData(
"1OVERBO", 1 );
1208 simData.registerCellData(
"OIL_DEN", 1 );
1209 simData.registerCellData(
"OIL_VISC", 1 );
1210 simData.registerCellData(
"OILKR", 1 );
1213 VectorType& bOil = liquid_active ? simData.getCellData(
"1OVERBO" ) : zero;
1214 VectorType& rhoOil = liquid_active ? simData.getCellData(
"OIL_DEN" ) : zero;
1215 VectorType& muOil = liquid_active ? simData.getCellData(
"OIL_VISC" ) : zero;
1216 VectorType& krOil = liquid_active ? simData.getCellData(
"OILKR" ) : zero;
1219 if( vapour_active ) {
1220 simData.registerCellData(
"1OVERBG", 1 );
1221 simData.registerCellData(
"GAS_DEN", 1 );
1222 simData.registerCellData(
"GAS_VISC", 1 );
1223 simData.registerCellData(
"GASKR", 1 );
1226 VectorType& bGas = vapour_active ? simData.getCellData(
"1OVERBG" ) : zero;
1227 VectorType& rhoGas = vapour_active ? simData.getCellData(
"GAS_DEN" ) : zero;
1228 VectorType& muGas = vapour_active ? simData.getCellData(
"GAS_VISC" ) : zero;
1229 VectorType& krGas = vapour_active ? simData.getCellData(
"GASKR" ) : zero;
1231 simData.registerCellData( BlackoilState::GASOILRATIO, 1 );
1232 simData.registerCellData( BlackoilState::RV, 1 );
1233 simData.registerCellData(
"RSSAT", 1 );
1234 simData.registerCellData(
"RVSAT", 1 );
1236 VectorType& Rs = simData.getCellData( BlackoilState::GASOILRATIO );
1237 VectorType& Rv = simData.getCellData( BlackoilState::RV );
1238 VectorType& RsSat = simData.getCellData(
"RSSAT" );
1239 VectorType& RvSat = simData.getCellData(
"RVSAT" );
1241 simData.registerCellData(
"PBUB", 1 );
1242 simData.registerCellData(
"PDEW", 1 );
1244 VectorType& Pb = simData.getCellData(
"PBUB" );
1245 VectorType& Pd = simData.getCellData(
"PDEW" );
1247 simData.registerCellData(
"SOMAX", 1 );
1248 VectorType& somax = simData.getCellData(
"SOMAX" );
1252 simData.registerCellData(
"PCSWMDC_GO", 1 );
1253 simData.registerCellData(
"KRNSWMDC_GO", 1 );
1255 simData.registerCellData(
"PCSWMDC_OW", 1 );
1256 simData.registerCellData(
"KRNSWMDC_OW", 1 );
1258 VectorType& pcSwMdc_go = simData.getCellData(
"PCSWMDC_GO" );
1259 VectorType& krnSwMdc_go = simData.getCellData(
"KRNSWMDC_GO" );
1261 VectorType& pcSwMdc_ow = simData.getCellData(
"PCSWMDC_OW" );
1262 VectorType& krnSwMdc_ow = simData.getCellData(
"KRNSWMDC_OW" );
1265 simData.registerCellData(
"SSOL", 1 );
1267 VectorType& ssol = has_solvent_ ? simData.getCellData(
"SSOL" ) : zero;
1270 simData.registerCellData(
"POLYMER", 1 );
1272 VectorType& cpolymer = has_polymer_ ? simData.getCellData(
"POLYMER" ) : zero;
1274 std::vector<int> failed_cells_pb;
1275 std::vector<int> failed_cells_pd;
1276 const auto& gridView = ebosSimulator().gridView();
1277 auto elemIt = gridView.template begin< 0, Dune::Interior_Partition>();
1278 const auto& elemEndIt = gridView.template end< 0, Dune::Interior_Partition>();
1279 ElementContext elemCtx(ebosSimulator());
1281 for (; elemIt != elemEndIt; ++elemIt) {
1282 const auto& elem = *elemIt;
1284 elemCtx.updatePrimaryStencil(elem);
1285 elemCtx.updatePrimaryIntensiveQuantities(0);
1287 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
1288 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
1290 const auto& fs = intQuants.fluidState();
1292 const int satIdx = cellIdx * num_phases;
1294 pressureOil[cellIdx] = fs.pressure(FluidSystem::oilPhaseIdx).value();
1296 temperature[cellIdx] = fs.temperature(FluidSystem::oilPhaseIdx).value();
1298 somax[cellIdx] = ebosSimulator().model().maxOilSaturation(cellIdx);
1300 const auto& matLawManager = ebosSimulator().problem().materialLawManager();
1301 if (matLawManager->enableHysteresis()) {
1302 matLawManager->oilWaterHysteresisParams(
1303 pcSwMdc_ow[cellIdx],
1304 krnSwMdc_ow[cellIdx],
1306 matLawManager->gasOilHysteresisParams(
1307 pcSwMdc_go[cellIdx],
1308 krnSwMdc_go[cellIdx],
1313 saturation[ satIdx + aqua_pos ] = fs.saturation(FluidSystem::waterPhaseIdx).value();
1314 bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value();
1315 rhoWater[cellIdx] = fs.density(FluidSystem::waterPhaseIdx).value();
1316 muWater[cellIdx] = fs.viscosity(FluidSystem::waterPhaseIdx).value();
1317 krWater[cellIdx] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
1319 if (vapour_active) {
1320 saturation[ satIdx + vapour_pos ] = fs.saturation(FluidSystem::gasPhaseIdx).value();
1321 bGas[cellIdx] = fs.invB(FluidSystem::gasPhaseIdx).value();
1322 rhoGas[cellIdx] = fs.density(FluidSystem::gasPhaseIdx).value();
1323 muGas[cellIdx] = fs.viscosity(FluidSystem::gasPhaseIdx).value();
1324 krGas[cellIdx] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
1325 Rs[cellIdx] = fs.Rs().value();
1326 Rv[cellIdx] = fs.Rv().value();
1327 RsSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
1328 FluidSystem::oilPhaseIdx,
1329 intQuants.pvtRegionIndex(),
1331 RvSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
1332 FluidSystem::gasPhaseIdx,
1333 intQuants.pvtRegionIndex(),
1336 Pb[cellIdx] = FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()).value();
1338 catch (
const NumericalProblem& e) {
1339 const auto globalIdx = ebosSimulator_.gridManager().grid().globalCell()[cellIdx];
1340 failed_cells_pb.push_back(globalIdx);
1343 Pd[cellIdx] = FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()).value();
1345 catch (
const NumericalProblem& e) {
1346 const auto globalIdx = ebosSimulator_.gridManager().grid().globalCell()[cellIdx];
1347 failed_cells_pd.push_back(globalIdx);
1352 saturation[ satIdx + liquid_pos ] = fs.saturation(FluidSystem::oilPhaseIdx).value();
1353 bOil[cellIdx] = fs.invB(FluidSystem::oilPhaseIdx).value();
1354 rhoOil[cellIdx] = fs.density(FluidSystem::oilPhaseIdx).value();
1355 muOil[cellIdx] = fs.viscosity(FluidSystem::oilPhaseIdx).value();
1356 krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
1361 ssol[cellIdx] = intQuants.solventSaturation().value();
1366 cpolymer[cellIdx] = intQuants.polymerConcentration().value();
1374 if (ebosSimulator_.episodeIndex() < 0 && vapour_active && liquid_active ) {
1376 Rs[cellIdx] = localState.getCellData( BlackoilState::GASOILRATIO )[cellIdx];
1377 Rv[cellIdx] = localState.getCellData( BlackoilState::RV)[cellIdx];
1380 auto fs_updated = fs;
1381 auto rs_eval = fs_updated.Rs();
1382 rs_eval.setValue( Rs[cellIdx] );
1383 fs_updated.setRs(rs_eval);
1384 auto rv_eval = fs_updated.Rv();
1385 rv_eval.setValue( Rv[cellIdx] );
1386 fs_updated.setRv(rv_eval);
1389 rhoOil[cellIdx] = FluidSystem::density(fs_updated,
1390 FluidSystem::oilPhaseIdx,
1391 intQuants.pvtRegionIndex()).value();
1392 rhoGas[cellIdx] = FluidSystem::density(fs_updated,
1393 FluidSystem::gasPhaseIdx,
1394 intQuants.pvtRegionIndex()).value();
1396 bOil[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated,
1397 FluidSystem::oilPhaseIdx,
1398 intQuants.pvtRegionIndex()).value();
1399 bGas[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated,
1400 FluidSystem::gasPhaseIdx,
1401 intQuants.pvtRegionIndex()).value();
1403 muOil[cellIdx] = FluidSystem::viscosity(fs_updated,
1404 FluidSystem::oilPhaseIdx,
1405 intQuants.pvtRegionIndex()).value();
1406 muGas[cellIdx] = FluidSystem::viscosity(fs_updated,
1407 FluidSystem::gasPhaseIdx,
1408 intQuants.pvtRegionIndex()).value();
1413 const size_t max_num_cells_faillog = 20;
1415 int pb_size = failed_cells_pb.size(), pd_size = failed_cells_pd.size();
1416 std::vector<int> displ_pb, displ_pd, recv_len_pb, recv_len_pd;
1417 const auto& comm = grid_.comm();
1419 if ( comm.rank() == 0 )
1421 displ_pb.resize(comm.size()+1, 0);
1422 displ_pd.resize(comm.size()+1, 0);
1423 recv_len_pb.resize(comm.size());
1424 recv_len_pd.resize(comm.size());
1427 comm.gather(&pb_size, recv_len_pb.data(), 1, 0);
1428 comm.gather(&pd_size, recv_len_pd.data(), 1, 0);
1429 std::partial_sum(recv_len_pb.begin(), recv_len_pb.end(), displ_pb.begin()+1);
1430 std::partial_sum(recv_len_pd.begin(), recv_len_pd.end(), displ_pd.begin()+1);
1431 std::vector<int> global_failed_cells_pb, global_failed_cells_pd;
1433 if ( comm.rank() == 0 )
1435 global_failed_cells_pb.resize(displ_pb.back());
1436 global_failed_cells_pd.resize(displ_pd.back());
1439 comm.gatherv(failed_cells_pb.data(),
static_cast<int>(failed_cells_pb.size()),
1440 global_failed_cells_pb.data(), recv_len_pb.data(),
1441 displ_pb.data(), 0);
1442 comm.gatherv(failed_cells_pd.data(),
static_cast<int>(failed_cells_pd.size()),
1443 global_failed_cells_pd.data(), recv_len_pd.data(),
1444 displ_pd.data(), 0);
1445 std::sort(global_failed_cells_pb.begin(), global_failed_cells_pb.end());
1446 std::sort(global_failed_cells_pd.begin(), global_failed_cells_pd.end());
1448 if (global_failed_cells_pb.size() > 0) {
1449 std::stringstream errlog;
1450 errlog <<
"Finding the bubble point pressure failed for " << global_failed_cells_pb.size() <<
" cells [";
1451 errlog << global_failed_cells_pb[0];
1452 const size_t max_elems = std::min(max_num_cells_faillog, failed_cells_pb.size());
1453 for (
size_t i = 1; i < max_elems; ++i) {
1454 errlog <<
", " << global_failed_cells_pb[i];
1456 if (global_failed_cells_pb.size() > max_num_cells_faillog) {
1460 OpmLog::warning(
"Bubble point numerical problem", errlog.str());
1462 if (global_failed_cells_pd.size() > 0) {
1463 std::stringstream errlog;
1464 errlog <<
"Finding the dew point pressure failed for " << global_failed_cells_pd.size() <<
" cells [";
1465 errlog << global_failed_cells_pd[0];
1466 const size_t max_elems = std::min(max_num_cells_faillog, global_failed_cells_pd.size());
1467 for (
size_t i = 1; i < max_elems; ++i) {
1468 errlog <<
", " << global_failed_cells_pd[i];
1470 if (global_failed_cells_pd.size() > max_num_cells_faillog) {
1474 OpmLog::warning(
"Dew point numerical problem", errlog.str());
1480 const FIPDataType& getFIPData()
const {
1484 const Simulator& ebosSimulator()
const 1485 {
return ebosSimulator_; }
1489 {
return failureReport_; }
1492 const ISTLSolverType& istlSolver()
const 1494 assert( istlSolver_ );
1495 return *istlSolver_;
1500 Simulator& ebosSimulator_;
1502 const ISTLSolverType* istlSolver_;
1503 const PhaseUsage phaseUsage_;
1504 VFPProperties vfp_properties_;
1506 const std::vector<bool> active_;
1508 const std::vector<int> cells_;
1509 const bool has_disgas_;
1510 const bool has_vapoil_;
1511 const bool has_solvent_;
1512 const bool has_polymer_;
1514 ModelParameters param_;
1515 SimulatorReport failureReport_;
1518 BlackoilWellModel<TypeTag>& well_model_;
1528 std::vector<std::vector<double>> residual_norms_history_;
1529 double current_relaxation_;
1539 wellModel()
const {
return well_model_; }
1541 int numWells()
const {
return well_model_.numWells(); }
1550 int flowPhaseToEbosCompIdx(
const int phaseIdx )
const 1552 const auto& pu = phaseUsage_;
1553 if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
1554 return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1555 if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
1556 return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1557 if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
1558 return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1564 int flowPhaseToEbosPhaseIdx(
const int phaseIdx )
const 1566 const auto& pu = phaseUsage_;
1567 if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
1568 return FluidSystem::waterPhaseIdx;
1569 if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
1570 return FluidSystem::oilPhaseIdx;
1571 if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
1572 return FluidSystem::gasPhaseIdx;
1574 assert(phaseIdx < 3);
1581 void updateRateConverter()
1583 rate_converter_.
defineState<ElementContext>(ebosSimulator_);
1588 void beginReportStep()
1590 isBeginReportStep_ =
true;
1593 void endReportStep()
1595 ebosSimulator_.problem().endEpisode();
1599 void assembleMassBalanceEq(
const SimulatorTimerInterface& timer,
1600 const int iterationIdx)
1602 ebosSimulator_.startNextEpisode( timer.currentStepLength() );
1603 ebosSimulator_.setEpisodeIndex( timer.reportStepNum() );
1604 ebosSimulator_.setTimeStepIndex( timer.reportStepNum() );
1605 ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
1607 static int prevEpisodeIdx = 10000;
1610 if (isBeginReportStep_) {
1611 isBeginReportStep_ =
false;
1612 ebosSimulator_.problem().beginEpisode();
1618 if (ebosSimulator_.model().newtonMethod().numIterations() == 0
1619 && prevEpisodeIdx < timer.reportStepNum())
1621 ebosSimulator_.problem().endTimeStep();
1624 ebosSimulator_.setTimeStepSize( timer.currentStepLength() );
1625 if (ebosSimulator_.model().newtonMethod().numIterations() == 0)
1627 ebosSimulator_.problem().beginTimeStep();
1630 ebosSimulator_.problem().beginIteration();
1631 ebosSimulator_.model().linearizer().linearize();
1632 ebosSimulator_.problem().endIteration();
1634 prevEpisodeIdx = ebosSimulator_.episodeIndex();
1636 if (param_.update_equations_scaling_) {
1637 std::cout <<
"equation scaling not suported yet" << std::endl;
1642 double dpMaxRel()
const {
return param_.dp_max_rel_; }
1643 double dsMax()
const {
return param_.ds_max_; }
1644 double drMaxRel()
const {
return param_.dr_max_rel_; }
1645 double maxResidualAllowed()
const {
return param_.max_residual_allowed_; }
1648 bool isBeginReportStep_;
1649 std::vector<bool> wasSwitched_;
1653 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &, const WellState &)
Called once before each time step.
Definition: BlackoilModelEbos.hpp:206
bool getConvergence(const SimulatorTimerInterface &timer, const int iteration, std::vector< double > &residual_norms)
Compute convergence based on total mass balance (tol_mb) and maximum residual mass balance (tol_cnv)...
Definition: BlackoilModelEbos.hpp:837
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition: BlackoilModelEbos.hpp:498
Definition: BlackoilModelEnums.hpp:66
SimulatorReport assemble(const SimulatorTimerInterface &timer, const int iterationIdx, WellState &well_state)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilModelEbos.hpp:377
int max_strict_iter_
Maximum number of Newton iterations before we give up on the CNV convergence criterion.
Definition: BlackoilModelParameters.hpp:73
int sizeNonLinear() const
The size (number of unknowns) of the nonlinear system of equations.
Definition: BlackoilModelEbos.hpp:490
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:66
std::vector< std::vector< double > > computeFluidInPlace(const T &, const std::vector< int > &fipnum) const
Wrapper required due to not following generic API.
Definition: BlackoilModelEbos.hpp:1026
Adapter to turn a matrix into a linear operator.
Definition: BlackoilModelEbos.hpp:551
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition: BlackoilModelEbos.hpp:1536
Definition: BlackoilModelEbos.hpp:80
void solve(Operator &opA, Vector &x, Vector &istlb, ScalarProd &sp, Precond &precond, Dune::InverseOperatorResult &result) const
Solve the system using the given preconditioner and scalar product.
Definition: ISTLSolver.hpp:487
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelEbos.hpp:1003
double tolerance_cnv_
Local convergence tolerance (max of local saturation errors).
Definition: BlackoilModelParameters.hpp:48
void solveJacobianSystem(BVector &x) const
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilModelEbos.hpp:505
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolver.hpp:323
A model implementation for three-phase black oil.
Definition: BlackoilModelEbos.hpp:104
AutoDiffMatrix is a wrapper class that optimizes matrix operations.
Definition: AutoDiffMatrix.hpp:43
double tolerance_mb_
Relative mass balance tolerance (total mass balance error).
Definition: BlackoilModelParameters.hpp:46
void afterStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, WellState &well_state)
Called once after each time step.
Definition: BlackoilModelEbos.hpp:364
SimulatorReport nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver, ReservoirState &, WellState &well_state)
Called once per nonlinear iteration.
Definition: BlackoilModelEbos.hpp:242
WellModelMatrixAdapter(const M &A, const WellModel &wellMod, const boost::any ¶llelInformation=boost::any())
constructor: just store a reference to a matrix
Definition: BlackoilModelEbos.hpp:575
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
Facility for converting component rates at surface conditions to phase (voidage) rates at reservoir c...
virtual bool lastStepFailed() const =0
Return true if last time step failed.
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:31
virtual double currentStepLength() const =0
Current step length.
BlackoilModelEbos(Simulator &ebosSimulator, const ModelParameters ¶m, BlackoilWellModel< TypeTag > &well_model, RateConverterType &rate_converter, const NewtonIterationBlackoilInterface &linsolver, const bool terminal_output)
Construct the model.
Definition: BlackoilModelEbos.hpp:155
Definition: GridHelpers.cpp:27
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition: BlackoilModelEbos.hpp:771
bool localWellsActive() const
return true if wells are available on this process
Definition: BlackoilModelEbos.hpp:1544
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
const SimulatorReport & failureReport() const
return the statistics if the nonlinearIteration() method failed
Definition: BlackoilModelEbos.hpp:1488
The solver category.
Definition: BlackoilModelEbos.hpp:569
void updateState(const BVector &dx)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilModelEbos.hpp:630
bool terminal_output_
Whether we print something to std::cout.
Definition: BlackoilModelEbos.hpp:1521
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
long int global_nc_
The number of cells of the global grid.
Definition: BlackoilModelEbos.hpp:1523
void defineState(const BlackoilState &state, const boost::any &info=boost::any())
Compute average hydrocarbon pressure and maximum dissolution and evaporation at average hydrocarbon p...
Definition: RateConverter.hpp:440
bool use_update_stabilization_
Try to detect oscillation or stagnation.
Definition: BlackoilModelParameters.hpp:82
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolver.hpp:374