00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #ifndef OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
00025 #define OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
00026
00027 #include <ebos/eclproblem.hh>
00028 #include <ewoms/common/start.hh>
00029
00030 #include <opm/autodiff/BlackoilModelParameters.hpp>
00031 #include <opm/autodiff/BlackoilWellModel.hpp>
00032 #include <opm/autodiff/AutoDiffBlock.hpp>
00033 #include <opm/autodiff/AutoDiffHelpers.hpp>
00034 #include <opm/autodiff/GridHelpers.hpp>
00035 #include <opm/autodiff/WellHelpers.hpp>
00036 #include <opm/autodiff/GeoProps.hpp>
00037 #include <opm/autodiff/WellDensitySegmented.hpp>
00038 #include <opm/autodiff/VFPProperties.hpp>
00039 #include <opm/autodiff/VFPProdProperties.hpp>
00040 #include <opm/autodiff/VFPInjProperties.hpp>
00041 #include <opm/autodiff/BlackoilDetails.hpp>
00042 #include <opm/autodiff/BlackoilModelEnums.hpp>
00043 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
00044 #include <opm/autodiff/RateConverter.hpp>
00045
00046 #include <opm/core/grid.h>
00047 #include <opm/core/simulator/SimulatorReport.hpp>
00048 #include <opm/core/linalg/LinearSolverInterface.hpp>
00049 #include <opm/core/linalg/ParallelIstlInformation.hpp>
00050 #include <opm/core/props/phaseUsageFromDeck.hpp>
00051 #include <opm/common/ErrorMacros.hpp>
00052 #include <opm/common/Exceptions.hpp>
00053 #include <opm/common/OpmLog/OpmLog.hpp>
00054 #include <opm/parser/eclipse/Units/Units.hpp>
00055 #include <opm/core/well_controls.h>
00056 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
00057 #include <opm/core/utility/parameters/ParameterGroup.hpp>
00058 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00059 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
00060
00061 #include <opm/autodiff/ISTLSolver.hpp>
00062 #include <opm/common/data/SimulationDataContainer.hpp>
00063
00064 #include <dune/istl/owneroverlapcopy.hh>
00065 #include <dune/common/parallel/collectivecommunication.hh>
00066 #include <dune/common/timer.hh>
00067 #include <dune/common/unused.hh>
00068
00069 #include <cassert>
00070 #include <cmath>
00071 #include <iostream>
00072 #include <iomanip>
00073 #include <limits>
00074 #include <vector>
00075 #include <algorithm>
00076
00077
00078
00079
00080 namespace Ewoms {
00081 namespace Properties {
00082 NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
00083 SET_BOOL_PROP(EclFlowProblem, DisableWells, true);
00084 SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false);
00085 SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility, true);
00086
00087 SET_BOOL_PROP(EclFlowProblem, BlackoilConserveSurfaceVolume, true);
00088 SET_BOOL_PROP(EclFlowProblem, UseVolumetricResidual, false);
00089
00090
00091
00092
00093 SET_BOOL_PROP(EclFlowProblem, EnableSwatinit, false);
00094 }}
00095
00096 namespace Opm {
00103 template <class TypeTag>
00104 class BlackoilModelEbos
00105 {
00106 public:
00107
00108 typedef BlackoilState ReservoirState;
00109 typedef WellStateFullyImplicitBlackoil WellState;
00110 typedef BlackoilModelParameters ModelParameters;
00111
00112 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
00113 typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
00114 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
00115 typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
00116 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables ;
00117 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
00118 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
00119 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
00120 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
00121
00122 typedef double Scalar;
00123 static const int numEq = Indices::numEq;
00124 static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
00125 static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
00126 static const int solventSaturationIdx = Indices::solventSaturationIdx;
00127 static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
00128
00129 typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
00130 typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
00131 typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
00132 typedef Dune::BlockVector<VectorBlockType> BVector;
00133
00134 typedef ISTLSolver< MatrixBlockType, VectorBlockType, Indices::pressureSwitchIdx > ISTLSolverType;
00135
00136
00137
00138 using RateConverterType = RateConverter::
00139 SurfaceToReservoirVoidage<BlackoilPropsAdFromDeck::FluidSystem, std::vector<int> >;
00140
00141 typedef Opm::FIPData FIPDataType;
00142
00143
00144
00155 BlackoilModelEbos(Simulator& ebosSimulator,
00156 const ModelParameters& param,
00157 BlackoilWellModel<TypeTag>& well_model,
00158 RateConverterType& rate_converter,
00159 const NewtonIterationBlackoilInterface& linsolver,
00160 const bool terminal_output
00161 )
00162 : ebosSimulator_(ebosSimulator)
00163 , grid_(ebosSimulator_.gridManager().grid())
00164 , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
00165 , phaseUsage_(phaseUsageFromDeck(eclState()))
00166 , vfp_properties_(
00167 eclState().getTableManager().getVFPInjTables(),
00168 eclState().getTableManager().getVFPProdTables())
00169 , active_(detail::activePhases(phaseUsage_))
00170 , has_disgas_(FluidSystem::enableDissolvedGas())
00171 , has_vapoil_(FluidSystem::enableVaporizedOil())
00172 , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
00173 , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
00174 , param_( param )
00175 , well_model_ (well_model)
00176 , terminal_output_ (terminal_output)
00177 , rate_converter_(rate_converter)
00178 , current_relaxation_(1.0)
00179 , dx_old_(AutoDiffGrid::numCells(grid_))
00180 , isBeginReportStep_(false)
00181 {
00182
00183
00184 int wellsActive = localWellsActive() ? 1 : 0;
00185 wellsActive = grid_.comm().max(wellsActive);
00186 wellModel().setWellsActive( wellsActive );
00187
00188 global_nc_ = detail::countGlobalCells(grid_);
00189 wellModel().setVFPProperties(&vfp_properties_);
00190 if (!istlSolver_)
00191 {
00192 OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed");
00193 }
00194 }
00195
00196 bool isParallel() const
00197 { return grid_.comm().size() > 1; }
00198
00199 const EclipseState& eclState() const
00200 { return ebosSimulator_.gridManager().eclState(); }
00201
00206 void prepareStep(const SimulatorTimerInterface& timer,
00207 const ReservoirState& ,
00208 const WellState& )
00209 {
00210 if ( wellModel().wellCollection()->havingVREPGroups() ) {
00211 updateRateConverter();
00212 }
00213
00214
00215
00216
00217
00218 if ( timer.lastStepFailed() ) {
00219 ebosSimulator_.model().solution( 0 ) = ebosSimulator_.model().solution( 1 );
00220 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
00221 } else {
00222
00223 ebosSimulator_.model().solution( 1 ) = ebosSimulator_.model().solution( 0 );
00224 }
00225
00226 unsigned numDof = ebosSimulator_.model().numGridDof();
00227 wasSwitched_.resize(numDof);
00228 std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
00229 }
00230
00231
00241 template <class NonlinearSolverType>
00242 SimulatorReport nonlinearIteration(const int iteration,
00243 const SimulatorTimerInterface& timer,
00244 NonlinearSolverType& nonlinear_solver,
00245 ReservoirState& ,
00246 WellState& well_state)
00247 {
00248 SimulatorReport report;
00249 failureReport_ = SimulatorReport();
00250 Dune::Timer perfTimer;
00251
00252 perfTimer.start();
00253 if (iteration == 0) {
00254
00255
00256 residual_norms_history_.clear();
00257 current_relaxation_ = 1.0;
00258 dx_old_ = 0.0;
00259 }
00260
00261 report.total_linearizations = 1;
00262
00263 try {
00264 report += assemble(timer, iteration, well_state);
00265 report.assemble_time += perfTimer.stop();
00266 }
00267 catch (...) {
00268 report.assemble_time += perfTimer.stop();
00269 failureReport_ += report;
00270
00271 throw;
00272 }
00273
00274 std::vector<double> residual_norms;
00275 perfTimer.reset();
00276 perfTimer.start();
00277
00278 report.converged = getConvergence(timer, iteration,residual_norms) && iteration > nonlinear_solver.minIter();
00279
00280
00281 if (wellModel().wellCollection()->groupControlActive()) {
00282 report.converged = report.converged && wellModel().wellCollection()->groupTargetConverged(well_state.wellRates());
00283 }
00284
00285 report.update_time += perfTimer.stop();
00286 residual_norms_history_.push_back(residual_norms);
00287 if (!report.converged) {
00288 perfTimer.reset();
00289 perfTimer.start();
00290 report.total_newton_iterations = 1;
00291
00292
00293
00294
00295
00296 const int nc = AutoDiffGrid::numCells(grid_);
00297 const int nw = numWells();
00298 BVector x(nc);
00299
00300 try {
00301 solveJacobianSystem(x);
00302 report.linear_solve_time += perfTimer.stop();
00303 report.total_linear_iterations += linearIterationsLastSolve();
00304 }
00305 catch (...) {
00306 report.linear_solve_time += perfTimer.stop();
00307 report.total_linear_iterations += linearIterationsLastSolve();
00308
00309 failureReport_ += report;
00310 throw;
00311 }
00312
00313 perfTimer.reset();
00314 perfTimer.start();
00315
00316
00317
00318
00319
00320 if( nw > 0 )
00321 {
00322 wellModel().recoverWellSolutionAndUpdateWellState(x, well_state);
00323 }
00324
00325 if (param_.use_update_stabilization_) {
00326
00327 bool isOscillate = false;
00328 bool isStagnate = false;
00329 nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
00330 if (isOscillate) {
00331 current_relaxation_ -= nonlinear_solver.relaxIncrement();
00332 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
00333 if (terminalOutputEnabled()) {
00334 std::string msg = " Oscillating behavior detected: Relaxation set to "
00335 + std::to_string(current_relaxation_);
00336 OpmLog::info(msg);
00337 }
00338 }
00339 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
00340 }
00341
00342
00343
00344 updateState(x);
00345
00346 report.update_time += perfTimer.stop();
00347 }
00348
00349 return report;
00350 }
00351
00352 void printIf(int c, double x, double y, double eps, std::string type) {
00353 if (std::abs(x-y) > eps) {
00354 std::cout << type << " " <<c << ": "<<x << " " << y << std::endl;
00355 }
00356 }
00357
00358
00364 void afterStep(const SimulatorTimerInterface& timer,
00365 const ReservoirState& reservoir_state,
00366 WellState& well_state)
00367 {
00368 DUNE_UNUSED_PARAMETER(timer);
00369 DUNE_UNUSED_PARAMETER(reservoir_state);
00370 DUNE_UNUSED_PARAMETER(well_state);
00371 }
00372
00377 SimulatorReport assemble(const SimulatorTimerInterface& timer,
00378 const int iterationIdx,
00379 WellState& well_state)
00380 {
00381 using namespace Opm::AutoDiffGrid;
00382
00383 SimulatorReport report;
00384
00385
00386 if ( wellModel().wellCollection()->havingVREPGroups() ) {
00387 updateRateConverter();
00388 }
00389
00390
00391 assembleMassBalanceEq(timer, iterationIdx);
00392
00393
00394 double dt = timer.currentStepLength();
00395
00396 try
00397 {
00398
00399
00400 report = wellModel().assemble(ebosSimulator_, iterationIdx, dt, well_state);
00401 }
00402 catch ( const Dune::FMatrixError& e )
00403 {
00404 OPM_THROW(Opm::NumericalProblem,"Error encounted when solving well equations");
00405 }
00406
00407 return report;
00408 }
00409
00410
00411 template <class Dummy>
00412 double relativeChange(const Dummy&, const Dummy&) const
00413 {
00414 Scalar resultDelta = 0.0;
00415 Scalar resultDenom = 0.0;
00416
00417 const auto& elemMapper = ebosSimulator_.model().elementMapper();
00418 const auto& gridView = ebosSimulator_.gridView();
00419 auto elemIt = gridView.template begin<0>();
00420 const auto& elemEndIt = gridView.template end<0>();
00421 for (; elemIt != elemEndIt; ++elemIt) {
00422 const auto& elem = *elemIt;
00423 if (elem.partitionType() != Dune::InteriorEntity)
00424 continue;
00425
00426 unsigned globalElemIdx = elemMapper.index(elem);
00427 const auto& priVarsNew = ebosSimulator_.model().solution(0)[globalElemIdx];
00428
00429 Scalar pressureNew;
00430 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
00431
00432 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
00433 Scalar oilSaturationNew = 1.0;
00434 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
00435 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx];
00436 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
00437 }
00438
00439 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
00440 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
00441 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
00442 }
00443
00444 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
00445 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
00446 }
00447
00448 const auto& priVarsOld = ebosSimulator_.model().solution(1)[globalElemIdx];
00449
00450 Scalar pressureOld;
00451 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
00452
00453 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
00454 Scalar oilSaturationOld = 1.0;
00455 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
00456 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
00457 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
00458 }
00459
00460 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
00461 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
00462 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
00463 }
00464
00465 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
00466 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
00467 }
00468
00469 Scalar tmp = pressureNew - pressureOld;
00470 resultDelta += tmp*tmp;
00471 resultDenom += pressureNew*pressureNew;
00472
00473 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
00474 Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
00475 resultDelta += tmp*tmp;
00476 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
00477 }
00478 }
00479
00480 resultDelta = gridView.comm().sum(resultDelta);
00481 resultDenom = gridView.comm().sum(resultDenom);
00482
00483 if (resultDenom > 0.0)
00484 return resultDelta/resultDenom;
00485 return 0.0;
00486 }
00487
00488
00490 int sizeNonLinear() const
00491 {
00492 const int nc = Opm::AutoDiffGrid::numCells(grid_);
00493 const int nw = numWells();
00494 return numComponents() * (nc + nw);
00495 }
00496
00498 int linearIterationsLastSolve() const
00499 {
00500 return istlSolver().iterations();
00501 }
00502
00505 void solveJacobianSystem(BVector& x) const
00506 {
00507 const auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
00508 auto& ebosResid = ebosSimulator_.model().linearizer().residual();
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520 wellModel().apply(ebosResid);
00521
00522
00523 x = 0.0;
00524
00525
00526 if( isParallel() )
00527 {
00528 typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
00529 Operator opA(ebosJac, well_model_, istlSolver().parallelInformation() );
00530 assert( opA.comm() );
00531 istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
00532 }
00533 else
00534 {
00535 typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
00536 Operator opA(ebosJac, well_model_);
00537 istlSolver().solve( opA, x, ebosResid );
00538 }
00539 }
00540
00541
00542
00543
00544
00550 template<class M, class X, class Y, class WellModel, bool overlapping >
00551 class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
00552 {
00553 typedef Dune::AssembledLinearOperator<M,X,Y> BaseType;
00554
00555 public:
00556 typedef M matrix_type;
00557 typedef X domain_type;
00558 typedef Y range_type;
00559 typedef typename X::field_type field_type;
00560
00561 #if HAVE_MPI
00562 typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
00563 #else
00564 typedef Dune::CollectiveCommunication< Grid > communication_type;
00565 #endif
00566
00567 enum {
00569 category = overlapping ?
00570 Dune::SolverCategory::overlapping :
00571 Dune::SolverCategory::sequential
00572 };
00573
00575 WellModelMatrixAdapter (const M& A, const WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
00576 : A_( A ), wellMod_( wellMod ), comm_()
00577 {
00578 #if HAVE_MPI
00579 if( parallelInformation.type() == typeid(ParallelISTLInformation) )
00580 {
00581 const ParallelISTLInformation& info =
00582 boost::any_cast<const ParallelISTLInformation&>( parallelInformation);
00583 comm_.reset( new communication_type( info.communicator() ) );
00584 }
00585 #endif
00586 }
00587
00588 virtual void apply( const X& x, Y& y ) const
00589 {
00590 A_.mv( x, y );
00591
00592 wellMod_.apply(x, y );
00593
00594 #if HAVE_MPI
00595 if( comm_ )
00596 comm_->project( y );
00597 #endif
00598 }
00599
00600
00601 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
00602 {
00603 A_.usmv(alpha,x,y);
00604
00605 wellMod_.applyScaleAdd( alpha, x, y );
00606
00607 #if HAVE_MPI
00608 if( comm_ )
00609 comm_->project( y );
00610 #endif
00611 }
00612
00613 virtual const matrix_type& getmat() const { return A_; }
00614
00615 communication_type* comm()
00616 {
00617 return comm_.operator->();
00618 }
00619
00620 protected:
00621 const matrix_type& A_ ;
00622 const WellModel& wellMod_;
00623 std::unique_ptr< communication_type > comm_;
00624 };
00625
00630 void updateState(const BVector& dx)
00631 {
00632 using namespace Opm::AutoDiffGrid;
00633
00634 const auto& ebosProblem = ebosSimulator_.problem();
00635
00636 unsigned numSwitched = 0;
00637 ElementContext elemCtx( ebosSimulator_ );
00638 const auto& gridView = ebosSimulator_.gridView();
00639 const auto& elemEndIt = gridView.template end<0>();
00640 SolutionVector& solution = ebosSimulator_.model().solution( 0 );
00641
00642 for (auto elemIt = gridView.template begin</*codim=*/0>();
00643 elemIt != elemEndIt;
00644 ++elemIt)
00645 {
00646 const auto& elem = *elemIt;
00647
00648 elemCtx.updatePrimaryStencil(elem);
00649 elemCtx.updatePrimaryIntensiveQuantities(0);
00650 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
00651
00652 PrimaryVariables& priVars = solution[ cell_idx ];
00653
00654 const double& dp = dx[cell_idx][Indices::pressureSwitchIdx];
00655 double& p = priVars[Indices::pressureSwitchIdx];
00656 const double& dp_rel_max = dpMaxRel();
00657 const int sign_dp = dp > 0 ? 1: -1;
00658 p -= sign_dp * std::min(std::abs(dp), std::abs(p)*dp_rel_max);
00659 p = std::max(p, 0.0);
00660
00661
00662 const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0;
00663 const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0;
00664
00665 double dso = 0.0;
00666 double dsg = 0.0;
00667 double drs = 0.0;
00668 double drv = 0.0;
00669
00670
00671 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
00672 dsg = dxvar;
00673 }
00674 else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
00675 drs = dxvar;
00676 }
00677 else {
00678 assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv);
00679 drv = dxvar;
00680 dsg = 0.0;
00681 }
00682
00683
00684 const double dss = has_solvent_ ? dx[cell_idx][Indices::solventSaturationIdx] : 0.0;
00685
00686
00687 const double dc = has_polymer_ ? dx[cell_idx][Indices::polymerConcentrationIdx] : 0.0;
00688
00689
00690 dso = - (dsw + dsg + dss);
00691
00692
00693
00694 double maxVal = 0.0;
00695 maxVal = std::max(std::abs(dsw),maxVal);
00696 maxVal = std::max(std::abs(dsg),maxVal);
00697 maxVal = std::max(std::abs(dso),maxVal);
00698 maxVal = std::max(std::abs(dss),maxVal);
00699
00700 double satScaleFactor = 1.0;
00701 if (maxVal > dsMax()) {
00702 satScaleFactor = dsMax()/maxVal;
00703 }
00704
00705 if (active_[Water]) {
00706 double& sw = priVars[Indices::waterSaturationIdx];
00707 sw -= satScaleFactor * dsw;
00708 }
00709
00710 if (active_[Gas]) {
00711 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
00712 double& sg = priVars[Indices::compositionSwitchIdx];
00713 sg -= satScaleFactor * dsg;
00714 }
00715 }
00716
00717 if (has_solvent_) {
00718 double& ss = priVars[Indices::solventSaturationIdx];
00719 ss -= satScaleFactor * dss;
00720 ss = std::min(std::max(ss, 0.0),1.0);
00721 }
00722 if (has_polymer_) {
00723 double& c = priVars[Indices::polymerConcentrationIdx];
00724 c -= satScaleFactor * dc;
00725 c = std::max(c, 0.0);
00726 }
00727
00728
00729 if (active_[Gas] && active_[Oil] ) {
00730 unsigned pvtRegionIdx = ebosSimulator_.problem().pvtRegionIndex(cell_idx);
00731 const double drmaxrel = drMaxRel();
00732 if (has_disgas_) {
00733 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
00734 Scalar RsSat =
00735 FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx, 300.0, p);
00736
00737 double& rs = priVars[Indices::compositionSwitchIdx];
00738 rs -= ((drs<0)?-1:1)*std::min(std::abs(drs), RsSat*drmaxrel);
00739 rs = std::max(rs, 0.0);
00740 }
00741
00742 }
00743 if (has_vapoil_) {
00744 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
00745 Scalar RvSat =
00746 FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx, 300.0, p);
00747
00748 double& rv = priVars[Indices::compositionSwitchIdx];
00749 rv -= ((drv<0)?-1:1)*std::min(std::abs(drv), RvSat*drmaxrel);
00750 rv = std::max(rv, 0.0);
00751 }
00752 }
00753 }
00754
00755
00756 if (wasSwitched_[cell_idx])
00757 wasSwitched_[cell_idx] = priVars.adaptPrimaryVariables(ebosProblem, cell_idx, 1e-5);
00758 else
00759 wasSwitched_[cell_idx] = priVars.adaptPrimaryVariables(ebosProblem, cell_idx);
00760
00761 if (wasSwitched_[cell_idx])
00762 ++numSwitched;
00763 }
00764
00765
00766 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
00767
00768 }
00769
00771 bool terminalOutputEnabled() const
00772 {
00773 return terminal_output_;
00774 }
00775
00776 template <class CollectiveCommunication>
00777 double convergenceReduction(const CollectiveCommunication& comm,
00778 const double pvSumLocal,
00779 std::vector< Scalar >& R_sum,
00780 std::vector< Scalar >& maxCoeff,
00781 std::vector< Scalar >& B_avg)
00782 {
00783
00784 double pvSum = pvSumLocal;
00785
00786 if( comm.size() > 1 )
00787 {
00788
00789 std::vector< Scalar > sumBuffer;
00790 std::vector< Scalar > maxBuffer;
00791 const int numComp = B_avg.size();
00792 sumBuffer.reserve( 2*numComp + 1 );
00793 maxBuffer.reserve( numComp );
00794 for( int compIdx = 0; compIdx < numComp; ++compIdx )
00795 {
00796 sumBuffer.push_back( B_avg[ compIdx ] );
00797 sumBuffer.push_back( R_sum[ compIdx ] );
00798 maxBuffer.push_back( maxCoeff[ compIdx ] );
00799 }
00800
00801
00802 sumBuffer.push_back( pvSum );
00803
00804
00805 comm.sum( sumBuffer.data(), sumBuffer.size() );
00806
00807
00808 comm.max( maxBuffer.data(), maxBuffer.size() );
00809
00810
00811 for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
00812 {
00813 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
00814 ++buffIdx;
00815
00816 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
00817 }
00818
00819 for( int compIdx = 0; compIdx < numComp; ++compIdx )
00820 {
00821 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
00822 }
00823
00824
00825 pvSum = sumBuffer.back();
00826 }
00827
00828
00829 return pvSum;
00830 }
00831
00837 bool getConvergence(const SimulatorTimerInterface& timer, const int iteration, std::vector<double>& residual_norms)
00838 {
00839 typedef std::vector< Scalar > Vector;
00840
00841 const double dt = timer.currentStepLength();
00842 const double tol_mb = param_.tolerance_mb_;
00843 const double tol_cnv = param_.tolerance_cnv_;
00844
00845 const int np = numPhases();
00846 const int numComp = numComponents();
00847
00848 Vector R_sum(numComp, 0.0 );
00849 Vector B_avg(numComp, 0.0 );
00850 Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
00851
00852 const auto& ebosModel = ebosSimulator_.model();
00853 const auto& ebosProblem = ebosSimulator_.problem();
00854
00855 const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
00856
00857 ElementContext elemCtx(ebosSimulator_);
00858 const auto& gridView = ebosSimulator().gridView();
00859 const auto& elemEndIt = gridView.template end<0, Dune::Interior_Partition>();
00860
00861 double pvSumLocal = 0.0;
00862 for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
00863 elemIt != elemEndIt;
00864 ++elemIt)
00865 {
00866 const auto& elem = *elemIt;
00867 elemCtx.updatePrimaryStencil(elem);
00868 elemCtx.updatePrimaryIntensiveQuantities(0);
00869 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
00870 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
00871 const auto& fs = intQuants.fluidState();
00872
00873 const double pvValue = ebosProblem.porosity(cell_idx) * ebosModel.dofTotalVolume( cell_idx );
00874 pvSumLocal += pvValue;
00875
00876 for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx )
00877 {
00878 const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx);
00879 const int ebosCompIdx = flowPhaseToEbosCompIdx(phaseIdx);
00880
00881 B_avg[ phaseIdx ] += 1.0 / fs.invB(ebosPhaseIdx).value();
00882 const auto R2 = ebosResid[cell_idx][ebosCompIdx];
00883
00884 R_sum[ phaseIdx ] += R2;
00885 maxCoeff[ phaseIdx ] = std::max( maxCoeff[ phaseIdx ], std::abs( R2 ) / pvValue );
00886 }
00887
00888 if ( has_solvent_ ) {
00889 B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
00890 const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
00891 R_sum[ contiSolventEqIdx ] += R2;
00892 maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue );
00893 }
00894 if (has_polymer_ ) {
00895 B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
00896 const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
00897 R_sum[ contiPolymerEqIdx ] += R2;
00898 maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
00899 }
00900
00901 }
00902
00903
00904 const int bSize = B_avg.size();
00905 for ( int i = 0; i<bSize; ++i )
00906 {
00907 B_avg[ i ] /= Scalar( global_nc_ );
00908 }
00909
00910
00911
00912
00913
00914
00915 const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal,
00916 R_sum, maxCoeff, B_avg);
00917
00918 Vector CNV(numComp);
00919 Vector mass_balance_residual(numComp);
00920
00921 bool converged_MB = true;
00922 bool converged_CNV = true;
00923
00924 for ( int compIdx = 0; compIdx < numComp; ++compIdx )
00925 {
00926 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
00927 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
00928 converged_MB = converged_MB && (mass_balance_residual[compIdx] < tol_mb);
00929 converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv);
00930
00931 residual_norms.push_back(CNV[compIdx]);
00932 }
00933
00934 const bool converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg);
00935
00936 bool converged = converged_MB && converged_Well;
00937
00938
00939
00940 if (iteration < param_.max_strict_iter_)
00941 converged = converged && converged_CNV;
00942
00943 if ( terminal_output_ )
00944 {
00945
00946 if (iteration == 0) {
00947 std::string msg = "Iter";
00948
00949 std::vector< std::string > key( numComp );
00950 for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) {
00951 const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
00952 key[ phaseIdx ] = std::toupper( phaseName.front() );
00953 }
00954 if (has_solvent_) {
00955 key[ solventSaturationIdx ] = "S";
00956 }
00957
00958 if (has_polymer_) {
00959 key[ polymerConcentrationIdx ] = "P";
00960 }
00961
00962 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
00963 msg += " MB(" + key[ compIdx ] + ") ";
00964 }
00965 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
00966 msg += " CNV(" + key[ compIdx ] + ") ";
00967 }
00968 OpmLog::note(msg);
00969 }
00970 std::ostringstream ss;
00971 const std::streamsize oprec = ss.precision(3);
00972 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
00973 ss << std::setw(4) << iteration;
00974 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
00975 ss << std::setw(11) << mass_balance_residual[compIdx];
00976 }
00977 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
00978 ss << std::setw(11) << CNV[compIdx];
00979 }
00980 ss.precision(oprec);
00981 ss.flags(oflags);
00982 OpmLog::note(ss.str());
00983 }
00984
00985 for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) {
00986 const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
00987
00988 if (std::isnan(mass_balance_residual[phaseIdx])
00989 || std::isnan(CNV[phaseIdx])) {
00990 OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName);
00991 }
00992 if (mass_balance_residual[phaseIdx] > maxResidualAllowed()
00993 || CNV[phaseIdx] > maxResidualAllowed()) {
00994 OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName);
00995 }
00996 }
00997
00998 return converged;
00999 }
01000
01001
01003 int numPhases() const
01004 {
01005 return phaseUsage_.num_phases;
01006 }
01007
01008 int numComponents() const
01009 {
01010 if (numPhases() == 2) {
01011 return 2;
01012 }
01013 int numComp = FluidSystem::numComponents;
01014 if (has_solvent_)
01015 numComp ++;
01016
01017 if (has_polymer_)
01018 numComp ++;
01019
01020 return numComp;
01021 }
01022
01024 template<class T>
01025 std::vector<std::vector<double> >
01026 computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
01027 {
01028 return computeFluidInPlace(fipnum);
01029 }
01030
01031 std::vector<std::vector<double> >
01032 computeFluidInPlace(const std::vector<int>& fipnum) const
01033 {
01034 const auto& comm = grid_.comm();
01035 const auto& gridView = ebosSimulator().gridView();
01036 const int nc = gridView.size(0);
01037 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
01038 int ntFip = *std::max_element(fipnum.begin(), fipnum.end());
01039 ntFip = comm.max(ntFip);
01040
01041 std::vector<double> tpv(ntFip, 0.0);
01042 std::vector<double> hcpv(ntFip, 0.0);
01043
01044 std::vector<std::vector<double> > regionValues(ntFip, std::vector<double>(FIPDataType::fipValues,0.0));
01045
01046 for (int i = 0; i<FIPDataType::fipValues; i++) {
01047 fip_.fip[i].resize(nc,0.0);
01048 }
01049
01050 ElementContext elemCtx(ebosSimulator_);
01051 const auto& elemEndIt = gridView.template end<0, Dune::Interior_Partition>();
01052 for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
01053 elemIt != elemEndIt;
01054 ++elemIt)
01055 {
01056 elemCtx.updatePrimaryStencil(*elemIt);
01057 elemCtx.updatePrimaryIntensiveQuantities(0);
01058
01059 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
01060 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
01061 const auto& fs = intQuants.fluidState();
01062
01063 const int regionIdx = fipnum[cellIdx] - 1;
01064 if (regionIdx < 0) {
01065
01066 continue;
01067 }
01068
01069
01070
01071
01072
01073
01074
01075
01076 const double pv =
01077 ebosSimulator_.model().dofTotalVolume(cellIdx)
01078 * intQuants.porosity().value();
01079
01080 for (int phase = 0; phase < maxnp; ++phase) {
01081 const double b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value();
01082 const double s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value();
01083
01084 fip_.fip[phase][cellIdx] = b * s * pv;
01085
01086 if (active_[ phase ]) {
01087 regionValues[regionIdx][phase] += fip_.fip[phase][cellIdx];
01088 }
01089 }
01090
01091 if (active_[ Oil ] && active_[ Gas ]) {
01092
01093 fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][cellIdx] = fs.Rs().value() * fip_.fip[FIPDataType::FIP_LIQUID][cellIdx];
01094 fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][cellIdx] = fs.Rv().value() * fip_.fip[FIPDataType::FIP_VAPOUR][cellIdx];
01095
01096 regionValues[regionIdx][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][cellIdx];
01097 regionValues[regionIdx][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][cellIdx];
01098 }
01099
01100 const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
01101
01102 tpv[regionIdx] += pv;
01103 hcpv[regionIdx] += pv * hydrocarbon;
01104 }
01105
01106
01107
01108 comm.sum(tpv.data(), tpv.size());
01109 comm.sum(hcpv.data(), hcpv.size());
01110
01111 for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
01112 elemIt != elemEndIt;
01113 ++elemIt)
01114 {
01115 const auto& elem = *elemIt;
01116
01117 elemCtx.updatePrimaryStencil(elem);
01118 elemCtx.updatePrimaryIntensiveQuantities(0);
01119
01120 unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
01121 const int regionIdx = fipnum[cellIdx] - 1;
01122 if (regionIdx < 0) {
01123
01124 continue;
01125 }
01126
01127 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
01128 const auto& fs = intQuants.fluidState();
01129
01130
01131
01132
01133
01134
01135
01136
01137 const double pv =
01138 ebosSimulator_.model().dofTotalVolume(cellIdx)
01139 * intQuants.porosity().value();
01140
01141 fip_.fip[FIPDataType::FIP_PV][cellIdx] = pv;
01142 const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
01143
01144
01145
01146 if (hcpv[regionIdx] > 1e-10) {
01147 fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[regionIdx];
01148 } else {
01149 fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() / tpv[regionIdx];
01150 }
01151
01152 regionValues[regionIdx][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][cellIdx];
01153 regionValues[regionIdx][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx];
01154 }
01155
01156
01157 for(int regionIdx=0; regionIdx < ntFip; ++regionIdx) {
01158 comm.sum(regionValues[regionIdx].data(), regionValues[regionIdx].size());
01159 }
01160
01161 return regionValues;
01162 }
01163
01164 SimulationDataContainer getSimulatorData ( const SimulationDataContainer& localState) const
01165 {
01166 typedef std::vector<double> VectorType;
01167
01168 const auto& ebosModel = ebosSimulator().model();
01169 const auto& phaseUsage = phaseUsage_;
01170
01171
01172 const int numCells = ebosModel.numGridDof();
01173 const int num_phases = numPhases();
01174
01175 SimulationDataContainer simData( numCells, 0, num_phases );
01176
01177
01178 const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
01179 const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
01180 const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
01181
01182 const int aqua_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Aqua ];
01183 const int liquid_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Liquid ];
01184 const int vapour_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Vapour ];
01185
01186 VectorType zero;
01187
01188 VectorType& pressureOil = simData.pressure();
01189 VectorType& temperature = simData.temperature();
01190 VectorType& saturation = simData.saturation();
01191
01192
01193 if( aqua_active ) {
01194 simData.registerCellData( "1OVERBW", 1 );
01195 simData.registerCellData( "WAT_DEN", 1 );
01196 simData.registerCellData( "WAT_VISC", 1 );
01197 simData.registerCellData( "WATKR", 1 );
01198 }
01199
01200 VectorType& bWater = aqua_active ? simData.getCellData( "1OVERBW" ) : zero;
01201 VectorType& rhoWater = aqua_active ? simData.getCellData( "WAT_DEN" ) : zero;
01202 VectorType& muWater = aqua_active ? simData.getCellData( "WAT_VISC" ) : zero;
01203 VectorType& krWater = aqua_active ? simData.getCellData( "WATKR" ) : zero;
01204
01205
01206 if( liquid_active ) {
01207 simData.registerCellData( "1OVERBO", 1 );
01208 simData.registerCellData( "OIL_DEN", 1 );
01209 simData.registerCellData( "OIL_VISC", 1 );
01210 simData.registerCellData( "OILKR", 1 );
01211 }
01212
01213 VectorType& bOil = liquid_active ? simData.getCellData( "1OVERBO" ) : zero;
01214 VectorType& rhoOil = liquid_active ? simData.getCellData( "OIL_DEN" ) : zero;
01215 VectorType& muOil = liquid_active ? simData.getCellData( "OIL_VISC" ) : zero;
01216 VectorType& krOil = liquid_active ? simData.getCellData( "OILKR" ) : zero;
01217
01218
01219 if( vapour_active ) {
01220 simData.registerCellData( "1OVERBG", 1 );
01221 simData.registerCellData( "GAS_DEN", 1 );
01222 simData.registerCellData( "GAS_VISC", 1 );
01223 simData.registerCellData( "GASKR", 1 );
01224 }
01225
01226 VectorType& bGas = vapour_active ? simData.getCellData( "1OVERBG" ) : zero;
01227 VectorType& rhoGas = vapour_active ? simData.getCellData( "GAS_DEN" ) : zero;
01228 VectorType& muGas = vapour_active ? simData.getCellData( "GAS_VISC" ) : zero;
01229 VectorType& krGas = vapour_active ? simData.getCellData( "GASKR" ) : zero;
01230
01231 simData.registerCellData( BlackoilState::GASOILRATIO, 1 );
01232 simData.registerCellData( BlackoilState::RV, 1 );
01233 simData.registerCellData( "RSSAT", 1 );
01234 simData.registerCellData( "RVSAT", 1 );
01235
01236 VectorType& Rs = simData.getCellData( BlackoilState::GASOILRATIO );
01237 VectorType& Rv = simData.getCellData( BlackoilState::RV );
01238 VectorType& RsSat = simData.getCellData( "RSSAT" );
01239 VectorType& RvSat = simData.getCellData( "RVSAT" );
01240
01241 simData.registerCellData( "PBUB", 1 );
01242 simData.registerCellData( "PDEW", 1 );
01243
01244 VectorType& Pb = simData.getCellData( "PBUB" );
01245 VectorType& Pd = simData.getCellData( "PDEW" );
01246
01247 simData.registerCellData( "SOMAX", 1 );
01248 VectorType& somax = simData.getCellData( "SOMAX" );
01249
01250
01251
01252 simData.registerCellData( "PCSWMDC_GO", 1 );
01253 simData.registerCellData( "KRNSWMDC_GO", 1 );
01254
01255 simData.registerCellData( "PCSWMDC_OW", 1 );
01256 simData.registerCellData( "KRNSWMDC_OW", 1 );
01257
01258 VectorType& pcSwMdc_go = simData.getCellData( "PCSWMDC_GO" );
01259 VectorType& krnSwMdc_go = simData.getCellData( "KRNSWMDC_GO" );
01260
01261 VectorType& pcSwMdc_ow = simData.getCellData( "PCSWMDC_OW" );
01262 VectorType& krnSwMdc_ow = simData.getCellData( "KRNSWMDC_OW" );
01263
01264 if (has_solvent_) {
01265 simData.registerCellData( "SSOL", 1 );
01266 }
01267 VectorType& ssol = has_solvent_ ? simData.getCellData( "SSOL" ) : zero;
01268
01269 if (has_polymer_) {
01270 simData.registerCellData( "POLYMER", 1 );
01271 }
01272 VectorType& cpolymer = has_polymer_ ? simData.getCellData( "POLYMER" ) : zero;
01273
01274 std::vector<int> failed_cells_pb;
01275 std::vector<int> failed_cells_pd;
01276 const auto& gridView = ebosSimulator().gridView();
01277 auto elemIt = gridView.template begin< 0, Dune::Interior_Partition>();
01278 const auto& elemEndIt = gridView.template end< 0, Dune::Interior_Partition>();
01279 ElementContext elemCtx(ebosSimulator());
01280
01281 for (; elemIt != elemEndIt; ++elemIt) {
01282 const auto& elem = *elemIt;
01283
01284 elemCtx.updatePrimaryStencil(elem);
01285 elemCtx.updatePrimaryIntensiveQuantities(0);
01286
01287 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
01288 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
01289
01290 const auto& fs = intQuants.fluidState();
01291
01292 const int satIdx = cellIdx * num_phases;
01293
01294 pressureOil[cellIdx] = fs.pressure(FluidSystem::oilPhaseIdx).value();
01295
01296 temperature[cellIdx] = fs.temperature(FluidSystem::oilPhaseIdx).value();
01297
01298 somax[cellIdx] = ebosSimulator().model().maxOilSaturation(cellIdx);
01299
01300 const auto& matLawManager = ebosSimulator().problem().materialLawManager();
01301 if (matLawManager->enableHysteresis()) {
01302 matLawManager->oilWaterHysteresisParams(
01303 pcSwMdc_ow[cellIdx],
01304 krnSwMdc_ow[cellIdx],
01305 cellIdx);
01306 matLawManager->gasOilHysteresisParams(
01307 pcSwMdc_go[cellIdx],
01308 krnSwMdc_go[cellIdx],
01309 cellIdx);
01310 }
01311
01312 if (aqua_active) {
01313 saturation[ satIdx + aqua_pos ] = fs.saturation(FluidSystem::waterPhaseIdx).value();
01314 bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value();
01315 rhoWater[cellIdx] = fs.density(FluidSystem::waterPhaseIdx).value();
01316 muWater[cellIdx] = fs.viscosity(FluidSystem::waterPhaseIdx).value();
01317 krWater[cellIdx] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
01318 }
01319 if (vapour_active) {
01320 saturation[ satIdx + vapour_pos ] = fs.saturation(FluidSystem::gasPhaseIdx).value();
01321 bGas[cellIdx] = fs.invB(FluidSystem::gasPhaseIdx).value();
01322 rhoGas[cellIdx] = fs.density(FluidSystem::gasPhaseIdx).value();
01323 muGas[cellIdx] = fs.viscosity(FluidSystem::gasPhaseIdx).value();
01324 krGas[cellIdx] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
01325 Rs[cellIdx] = fs.Rs().value();
01326 Rv[cellIdx] = fs.Rv().value();
01327 RsSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
01328 FluidSystem::oilPhaseIdx,
01329 intQuants.pvtRegionIndex(),
01330 1.0).value();
01331 RvSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
01332 FluidSystem::gasPhaseIdx,
01333 intQuants.pvtRegionIndex(),
01334 1.0).value();
01335 try {
01336 Pb[cellIdx] = FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()).value();
01337 }
01338 catch (const NumericalProblem& e) {
01339 const auto globalIdx = ebosSimulator_.gridManager().grid().globalCell()[cellIdx];
01340 failed_cells_pb.push_back(globalIdx);
01341 }
01342 try {
01343 Pd[cellIdx] = FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()).value();
01344 }
01345 catch (const NumericalProblem& e) {
01346 const auto globalIdx = ebosSimulator_.gridManager().grid().globalCell()[cellIdx];
01347 failed_cells_pd.push_back(globalIdx);
01348 }
01349 }
01350 if( liquid_active )
01351 {
01352 saturation[ satIdx + liquid_pos ] = fs.saturation(FluidSystem::oilPhaseIdx).value();
01353 bOil[cellIdx] = fs.invB(FluidSystem::oilPhaseIdx).value();
01354 rhoOil[cellIdx] = fs.density(FluidSystem::oilPhaseIdx).value();
01355 muOil[cellIdx] = fs.viscosity(FluidSystem::oilPhaseIdx).value();
01356 krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
01357 }
01358
01359 if (has_solvent_)
01360 {
01361 ssol[cellIdx] = intQuants.solventSaturation().value();
01362 }
01363
01364 if (has_polymer_)
01365 {
01366 cpolymer[cellIdx] = intQuants.polymerConcentration().value();
01367 }
01368
01369
01370
01371
01372
01373
01374 if (ebosSimulator_.episodeIndex() < 0 && vapour_active && liquid_active ) {
01375
01376 Rs[cellIdx] = localState.getCellData( BlackoilState::GASOILRATIO )[cellIdx];
01377 Rv[cellIdx] = localState.getCellData( BlackoilState::RV)[cellIdx];
01378
01379
01380 auto fs_updated = fs;
01381 auto rs_eval = fs_updated.Rs();
01382 rs_eval.setValue( Rs[cellIdx] );
01383 fs_updated.setRs(rs_eval);
01384 auto rv_eval = fs_updated.Rv();
01385 rv_eval.setValue( Rv[cellIdx] );
01386 fs_updated.setRv(rv_eval);
01387
01388
01389 rhoOil[cellIdx] = FluidSystem::density(fs_updated,
01390 FluidSystem::oilPhaseIdx,
01391 intQuants.pvtRegionIndex()).value();
01392 rhoGas[cellIdx] = FluidSystem::density(fs_updated,
01393 FluidSystem::gasPhaseIdx,
01394 intQuants.pvtRegionIndex()).value();
01395
01396 bOil[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated,
01397 FluidSystem::oilPhaseIdx,
01398 intQuants.pvtRegionIndex()).value();
01399 bGas[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated,
01400 FluidSystem::gasPhaseIdx,
01401 intQuants.pvtRegionIndex()).value();
01402
01403 muOil[cellIdx] = FluidSystem::viscosity(fs_updated,
01404 FluidSystem::oilPhaseIdx,
01405 intQuants.pvtRegionIndex()).value();
01406 muGas[cellIdx] = FluidSystem::viscosity(fs_updated,
01407 FluidSystem::gasPhaseIdx,
01408 intQuants.pvtRegionIndex()).value();
01409
01410 }
01411 }
01412
01413 const size_t max_num_cells_faillog = 20;
01414
01415 int pb_size = failed_cells_pb.size(), pd_size = failed_cells_pd.size();
01416 std::vector<int> displ_pb, displ_pd, recv_len_pb, recv_len_pd;
01417 const auto& comm = grid_.comm();
01418
01419 if ( comm.rank() == 0 )
01420 {
01421 displ_pb.resize(comm.size()+1, 0);
01422 displ_pd.resize(comm.size()+1, 0);
01423 recv_len_pb.resize(comm.size());
01424 recv_len_pd.resize(comm.size());
01425 }
01426
01427 comm.gather(&pb_size, recv_len_pb.data(), 1, 0);
01428 comm.gather(&pd_size, recv_len_pd.data(), 1, 0);
01429 std::partial_sum(recv_len_pb.begin(), recv_len_pb.end(), displ_pb.begin()+1);
01430 std::partial_sum(recv_len_pd.begin(), recv_len_pd.end(), displ_pd.begin()+1);
01431 std::vector<int> global_failed_cells_pb, global_failed_cells_pd;
01432
01433 if ( comm.rank() == 0 )
01434 {
01435 global_failed_cells_pb.resize(displ_pb.back());
01436 global_failed_cells_pd.resize(displ_pd.back());
01437 }
01438
01439 comm.gatherv(failed_cells_pb.data(), static_cast<int>(failed_cells_pb.size()),
01440 global_failed_cells_pb.data(), recv_len_pb.data(),
01441 displ_pb.data(), 0);
01442 comm.gatherv(failed_cells_pd.data(), static_cast<int>(failed_cells_pd.size()),
01443 global_failed_cells_pd.data(), recv_len_pd.data(),
01444 displ_pd.data(), 0);
01445 std::sort(global_failed_cells_pb.begin(), global_failed_cells_pb.end());
01446 std::sort(global_failed_cells_pd.begin(), global_failed_cells_pd.end());
01447
01448 if (global_failed_cells_pb.size() > 0) {
01449 std::stringstream errlog;
01450 errlog << "Finding the bubble point pressure failed for " << global_failed_cells_pb.size() << " cells [";
01451 errlog << global_failed_cells_pb[0];
01452 const size_t max_elems = std::min(max_num_cells_faillog, failed_cells_pb.size());
01453 for (size_t i = 1; i < max_elems; ++i) {
01454 errlog << ", " << global_failed_cells_pb[i];
01455 }
01456 if (global_failed_cells_pb.size() > max_num_cells_faillog) {
01457 errlog << ", ...";
01458 }
01459 errlog << "]";
01460 OpmLog::warning("Bubble point numerical problem", errlog.str());
01461 }
01462 if (global_failed_cells_pd.size() > 0) {
01463 std::stringstream errlog;
01464 errlog << "Finding the dew point pressure failed for " << global_failed_cells_pd.size() << " cells [";
01465 errlog << global_failed_cells_pd[0];
01466 const size_t max_elems = std::min(max_num_cells_faillog, global_failed_cells_pd.size());
01467 for (size_t i = 1; i < max_elems; ++i) {
01468 errlog << ", " << global_failed_cells_pd[i];
01469 }
01470 if (global_failed_cells_pd.size() > max_num_cells_faillog) {
01471 errlog << ", ...";
01472 }
01473 errlog << "]";
01474 OpmLog::warning("Dew point numerical problem", errlog.str());
01475 }
01476
01477 return simData;
01478 }
01479
01480 const FIPDataType& getFIPData() const {
01481 return fip_;
01482 }
01483
01484 const Simulator& ebosSimulator() const
01485 { return ebosSimulator_; }
01486
01488 const SimulatorReport& failureReport() const
01489 { return failureReport_; }
01490
01491 protected:
01492 const ISTLSolverType& istlSolver() const
01493 {
01494 assert( istlSolver_ );
01495 return *istlSolver_;
01496 }
01497
01498
01499
01500 Simulator& ebosSimulator_;
01501 const Grid& grid_;
01502 const ISTLSolverType* istlSolver_;
01503 const PhaseUsage phaseUsage_;
01504 VFPProperties vfp_properties_;
01505
01506 const std::vector<bool> active_;
01507
01508 const std::vector<int> cells_;
01509 const bool has_disgas_;
01510 const bool has_vapoil_;
01511 const bool has_solvent_;
01512 const bool has_polymer_;
01513
01514 ModelParameters param_;
01515 SimulatorReport failureReport_;
01516
01517
01518 BlackoilWellModel<TypeTag>& well_model_;
01519
01521 bool terminal_output_;
01523 long int global_nc_;
01524
01525
01526 RateConverterType& rate_converter_;
01527
01528 std::vector<std::vector<double>> residual_norms_history_;
01529 double current_relaxation_;
01530 BVector dx_old_;
01531 mutable FIPDataType fip_;
01532
01533 public:
01535 BlackoilWellModel<TypeTag>&
01536 wellModel() { return well_model_; }
01537
01538 const BlackoilWellModel<TypeTag>&
01539 wellModel() const { return well_model_; }
01540
01541 int numWells() const { return well_model_.numWells(); }
01542
01544 bool localWellsActive() const { return well_model_.localWellsActive(); }
01545
01546
01547
01548
01549 public:
01550 int flowPhaseToEbosCompIdx( const int phaseIdx ) const
01551 {
01552 const auto& pu = phaseUsage_;
01553 if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
01554 return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
01555 if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
01556 return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
01557 if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
01558 return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
01559
01560
01561 return phaseIdx;
01562 }
01563
01564 int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
01565 {
01566 const auto& pu = phaseUsage_;
01567 if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
01568 return FluidSystem::waterPhaseIdx;
01569 if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
01570 return FluidSystem::oilPhaseIdx;
01571 if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
01572 return FluidSystem::gasPhaseIdx;
01573
01574 assert(phaseIdx < 3);
01575
01576 return phaseIdx;
01577 }
01578
01579 private:
01580
01581 void updateRateConverter()
01582 {
01583 rate_converter_.defineState<ElementContext>(ebosSimulator_);
01584 }
01585
01586
01587 public:
01588 void beginReportStep()
01589 {
01590 isBeginReportStep_ = true;
01591 }
01592
01593 void endReportStep()
01594 {
01595 ebosSimulator_.problem().endEpisode();
01596 }
01597
01598 private:
01599 void assembleMassBalanceEq(const SimulatorTimerInterface& timer,
01600 const int iterationIdx)
01601 {
01602 ebosSimulator_.startNextEpisode( timer.currentStepLength() );
01603 ebosSimulator_.setEpisodeIndex( timer.reportStepNum() );
01604 ebosSimulator_.setTimeStepIndex( timer.reportStepNum() );
01605 ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
01606
01607 static int prevEpisodeIdx = 10000;
01608
01609
01610 if (isBeginReportStep_) {
01611 isBeginReportStep_ = false;
01612 ebosSimulator_.problem().beginEpisode();
01613 }
01614
01615
01616
01617
01618 if (ebosSimulator_.model().newtonMethod().numIterations() == 0
01619 && prevEpisodeIdx < timer.reportStepNum())
01620 {
01621 ebosSimulator_.problem().endTimeStep();
01622 }
01623
01624 ebosSimulator_.setTimeStepSize( timer.currentStepLength() );
01625 if (ebosSimulator_.model().newtonMethod().numIterations() == 0)
01626 {
01627 ebosSimulator_.problem().beginTimeStep();
01628 }
01629
01630 ebosSimulator_.problem().beginIteration();
01631 ebosSimulator_.model().linearizer().linearize();
01632 ebosSimulator_.problem().endIteration();
01633
01634 prevEpisodeIdx = ebosSimulator_.episodeIndex();
01635
01636 if (param_.update_equations_scaling_) {
01637 std::cout << "equation scaling not suported yet" << std::endl;
01638
01639 }
01640 }
01641
01642 double dpMaxRel() const { return param_.dp_max_rel_; }
01643 double dsMax() const { return param_.ds_max_; }
01644 double drMaxRel() const { return param_.dr_max_rel_; }
01645 double maxResidualAllowed() const { return param_.max_residual_allowed_; }
01646
01647 public:
01648 bool isBeginReportStep_;
01649 std::vector<bool> wasSwitched_;
01650 };
01651 }
01652
01653 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED