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_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
00025 #define OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
00026
00027 #include <opm/autodiff/BlackoilModelBase.hpp>
00028 #include <opm/autodiff/BlackoilDetails.hpp>
00029 #include <opm/autodiff/BlackoilLegacyDetails.hpp>
00030
00031 #include <opm/autodiff/AutoDiffBlock.hpp>
00032 #include <opm/autodiff/AutoDiffHelpers.hpp>
00033 #include <opm/autodiff/GridHelpers.hpp>
00034 #include <opm/autodiff/WellHelpers.hpp>
00035 #include <opm/autodiff/BlackoilPropsAdFromDeck.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
00042 #include <opm/core/grid.h>
00043 #include <opm/core/linalg/LinearSolverInterface.hpp>
00044 #include <opm/core/linalg/ParallelIstlInformation.hpp>
00045 #include <opm/core/props/rock/RockCompressibility.hpp>
00046 #include <opm/common/ErrorMacros.hpp>
00047 #include <opm/common/Exceptions.hpp>
00048 #include <opm/common/OpmLog/OpmLog.hpp>
00049 #include <opm/parser/eclipse/Units/Units.hpp>
00050 #include <opm/core/well_controls.h>
00051 #include <opm/core/utility/parameters/ParameterGroup.hpp>
00052 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00053 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
00054
00055 #include <dune/common/timer.hh>
00056
00057 #include <cassert>
00058 #include <cmath>
00059 #include <iostream>
00060 #include <iomanip>
00061 #include <limits>
00062 #include <vector>
00063 #include <algorithm>
00064
00065
00066
00067 #define OPM_AD_DUMP(foo) \
00068 do { \
00069 std::cout << "==========================================\n" \
00070 << #foo ":\n" \
00071 << collapseJacs(foo) << std::endl; \
00072 } while (0)
00073
00074 #define OPM_AD_DUMPVAL(foo) \
00075 do { \
00076 std::cout << "==========================================\n" \
00077 << #foo ":\n" \
00078 << foo.value() << std::endl; \
00079 } while (0)
00080
00081 #define OPM_AD_DISKVAL(foo) \
00082 do { \
00083 std::ofstream os(#foo); \
00084 os.precision(16); \
00085 os << foo.value() << std::endl; \
00086 } while (0)
00087
00088
00089 namespace Opm {
00090
00091 typedef AutoDiffBlock<double> ADB;
00092 typedef ADB::V V;
00093 typedef ADB::M M;
00094 typedef Eigen::Array<double,
00095 Eigen::Dynamic,
00096 Eigen::Dynamic,
00097 Eigen::RowMajor> DataBlock;
00098
00099 template <class Grid, class WellModel, class Implementation>
00100 BlackoilModelBase<Grid, WellModel, Implementation>::
00101 BlackoilModelBase(const ModelParameters& param,
00102 const Grid& grid ,
00103 const BlackoilPropsAdFromDeck& fluid,
00104 const DerivedGeology& geo ,
00105 const RockCompressibility* rock_comp_props,
00106 const WellModel& well_model,
00107 const NewtonIterationBlackoilInterface& linsolver,
00108 std::shared_ptr< const Opm::EclipseState > eclState,
00109 const bool has_disgas,
00110 const bool has_vapoil,
00111 const bool terminal_output)
00112 : grid_ (grid)
00113 , fluid_ (fluid)
00114 , geo_ (geo)
00115 , rock_comp_props_(rock_comp_props)
00116 , vfp_properties_(
00117 eclState->getTableManager().getVFPInjTables(),
00118 eclState->getTableManager().getVFPProdTables())
00119 , linsolver_ (linsolver)
00120 , active_(detail::activePhases(fluid.phaseUsage()))
00121 , canph_ (detail::active2Canonical(fluid.phaseUsage()))
00122 , cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
00123 , ops_ (grid, geo.nnc())
00124 , has_disgas_(has_disgas)
00125 , has_vapoil_(has_vapoil)
00126 , param_( param )
00127 , use_threshold_pressure_(false)
00128 , sd_ (fluid.numPhases())
00129 , phaseCondition_(AutoDiffGrid::numCells(grid))
00130 , well_model_ (well_model)
00131 , isRs_(V::Zero(AutoDiffGrid::numCells(grid)))
00132 , isRv_(V::Zero(AutoDiffGrid::numCells(grid)))
00133 , isSg_(V::Zero(AutoDiffGrid::numCells(grid)))
00134 , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
00135 ADB::null(),
00136 ADB::null(),
00137 { 1.1169, 1.0031, 0.0031 },
00138 false } )
00139 , terminal_output_ (terminal_output)
00140 , material_name_(0)
00141 , current_relaxation_(1.0)
00142
00143
00144
00145
00146 , rate_converter_(fluid_.phaseUsage(), std::vector<int>(AutoDiffGrid::numCells(grid_),0))
00147 {
00148 if (active_[Water]) {
00149 material_name_.push_back("Water");
00150 }
00151 if (active_[Oil]) {
00152 material_name_.push_back("Oil");
00153 }
00154 if (active_[Gas]) {
00155 material_name_.push_back("Gas");
00156 }
00157
00158 assert(numMaterials() == std::accumulate(active_.begin(), active_.end(), 0));
00159
00160 const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
00161 const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_);
00162
00163 well_model_.init(&fluid_, &active_, &phaseCondition_, &vfp_properties_, gravity, depth);
00164
00165
00166
00167 #if HAVE_MPI
00168 const Wells* wells_arg = asImpl().well_model_.wellsPointer();
00169
00170 if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
00171 {
00172 const ParallelISTLInformation& info =
00173 boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
00174 if ( terminal_output_ ) {
00175
00176 terminal_output_ = (info.communicator().rank()==0);
00177 }
00178 int local_number_of_wells = localWellsActive() ? wells().number_of_wells : 0;
00179 int global_number_of_wells = info.communicator().sum(local_number_of_wells);
00180 const bool wells_active = ( wells_arg && global_number_of_wells > 0 );
00181 wellModel().setWellsActive(wells_active);
00182
00183 std::vector<int> v( Opm::AutoDiffGrid::numCells(grid_), 1);
00184 global_nc_ = 0;
00185 info.computeReduction(v, Opm::Reduction::makeGlobalSumFunctor<int>(), global_nc_);
00186 }else
00187 #endif
00188 {
00189 wellModel().setWellsActive( localWellsActive() );
00190 global_nc_ = Opm::AutoDiffGrid::numCells(grid_);
00191 }
00192 }
00193
00194
00195 template <class Grid, class WellModel, class Implementation>
00196 bool
00197 BlackoilModelBase<Grid, WellModel, Implementation>::
00198 isParallel() const
00199 {
00200 #if HAVE_MPI
00201 if ( linsolver_.parallelInformation().type() !=
00202 typeid(ParallelISTLInformation) )
00203 {
00204 return false;
00205 }
00206 else
00207 {
00208 const auto& comm =boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation()).communicator();
00209 return comm.size() > 1;
00210 }
00211 #else
00212 return false;
00213 #endif
00214 }
00215
00216
00217 template <class Grid, class WellModel, class Implementation>
00218 void
00219 BlackoilModelBase<Grid, WellModel, Implementation>::
00220 prepareStep(const SimulatorTimerInterface& timer,
00221 const ReservoirState& reservoir_state,
00222 const WellState& )
00223 {
00224 const double dt = timer.currentStepLength();
00225
00226 pvdt_ = geo_.poreVolume() / dt;
00227 if (active_[Gas]) {
00228 updatePrimalVariableFromState(reservoir_state);
00229 }
00230 }
00231
00232
00233
00234
00235
00236 template <class Grid, class WellModel, class Implementation>
00237 template <class NonlinearSolverType>
00238 SimulatorReport
00239 BlackoilModelBase<Grid, WellModel, Implementation>::
00240 nonlinearIteration(const int iteration,
00241 const SimulatorTimerInterface& timer,
00242 NonlinearSolverType& nonlinear_solver,
00243 ReservoirState& reservoir_state,
00244 WellState& well_state)
00245 {
00246 SimulatorReport report;
00247 Dune::Timer perfTimer;
00248
00249 perfTimer.start();
00250 const double dt = timer.currentStepLength();
00251
00252 if (iteration == 0) {
00253
00254
00255 residual_norms_history_.clear();
00256 current_relaxation_ = 1.0;
00257 dx_old_ = V::Zero(sizeNonLinear());
00258 }
00259 try {
00260 report += asImpl().assemble(reservoir_state, well_state, iteration == 0);
00261 report.assemble_time += perfTimer.stop();
00262 }
00263 catch (...) {
00264 report.assemble_time += perfTimer.stop();
00265 throw;
00266 }
00267
00268 report.total_linearizations = 1;
00269 perfTimer.reset();
00270 perfTimer.start();
00271 report.converged = asImpl().getConvergence(timer, iteration);
00272 residual_norms_history_.push_back(asImpl().computeResidualNorms());
00273 report.update_time += perfTimer.stop();
00274
00275 const bool must_solve = (iteration < nonlinear_solver.minIter()) || (!report.converged);
00276 if (must_solve) {
00277 perfTimer.reset();
00278 perfTimer.start();
00279 report.total_newton_iterations = 1;
00280
00281
00282 residual_.singlePrecision = ( dt < param_.maxSinglePrecisionTimeStep_ );
00283
00284
00285 V dx;
00286 try {
00287 dx = asImpl().solveJacobianSystem();
00288 report.linear_solve_time += perfTimer.stop();
00289 report.total_linear_iterations += linearIterationsLastSolve();
00290 }
00291 catch (...) {
00292 report.linear_solve_time += perfTimer.stop();
00293 report.total_linear_iterations += linearIterationsLastSolve();
00294 throw;
00295 }
00296
00297 perfTimer.reset();
00298 perfTimer.start();
00299
00300 if (param_.use_update_stabilization_) {
00301
00302 bool isOscillate = false;
00303 bool isStagnate = false;
00304 nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
00305 if (isOscillate) {
00306 current_relaxation_ -= nonlinear_solver.relaxIncrement();
00307 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
00308 if (terminalOutputEnabled()) {
00309 std::string msg = " Oscillating behavior detected: Relaxation set to "
00310 + std::to_string(current_relaxation_);
00311 OpmLog::info(msg);
00312 }
00313 }
00314 nonlinear_solver.stabilizeNonlinearUpdate(dx, dx_old_, current_relaxation_);
00315 }
00316
00317
00318
00319 asImpl().updateState(dx, reservoir_state, well_state);
00320 report.update_time += perfTimer.stop();
00321 }
00322
00323 return report;
00324 }
00325
00326
00327
00328
00329
00330 template <class Grid, class WellModel, class Implementation>
00331 void
00332 BlackoilModelBase<Grid, WellModel, Implementation>::
00333 afterStep(const SimulatorTimerInterface& ,
00334 ReservoirState& ,
00335 WellState& )
00336 {
00337
00338 }
00339
00340
00341
00342
00343
00344 template <class Grid, class WellModel, class Implementation>
00345 int
00346 BlackoilModelBase<Grid, WellModel, Implementation>::
00347 sizeNonLinear() const
00348 {
00349 return residual_.sizeNonLinear();
00350 }
00351
00352
00353
00354
00355
00356 template <class Grid, class WellModel, class Implementation>
00357 int
00358 BlackoilModelBase<Grid, WellModel, Implementation>::
00359 linearIterationsLastSolve() const
00360 {
00361 return linsolver_.iterations();
00362 }
00363
00364
00365
00366
00367
00368 template <class Grid, class WellModel, class Implementation>
00369 bool
00370 BlackoilModelBase<Grid, WellModel, Implementation>::
00371 terminalOutputEnabled() const
00372 {
00373 return terminal_output_;
00374 }
00375
00376
00377
00378
00379
00380 template <class Grid, class WellModel, class Implementation>
00381 int
00382 BlackoilModelBase<Grid, WellModel, Implementation>::
00383 numPhases() const
00384 {
00385 return fluid_.numPhases();
00386 }
00387
00388
00389
00390
00391 template <class Grid, class WellModel, class Implementation>
00392 int
00393 BlackoilModelBase<Grid, WellModel, Implementation>::
00394 numMaterials() const
00395 {
00396 return material_name_.size();
00397 }
00398
00399
00400
00401
00402
00403 template <class Grid, class WellModel, class Implementation>
00404 const std::string&
00405 BlackoilModelBase<Grid, WellModel, Implementation>::
00406 materialName(int material_index) const
00407 {
00408 assert(material_index < numMaterials());
00409 return material_name_[material_index];
00410 }
00411
00412
00413
00414
00415
00416 template <class Grid, class WellModel, class Implementation>
00417 void
00418 BlackoilModelBase<Grid, WellModel, Implementation>::
00419 setThresholdPressures(const std::vector<double>& threshold_pressures)
00420 {
00421 const int num_faces = AutoDiffGrid::numFaces(grid_);
00422 const int num_nnc = geo_.nnc().numNNC();
00423 const int num_connections = num_faces + num_nnc;
00424 if (int(threshold_pressures.size()) != num_connections) {
00425 OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures input ( " << threshold_pressures.size()
00426 << " ), must be equal to number of faces + nncs ( " << num_faces << " + " << num_nnc << " ).");
00427 }
00428 use_threshold_pressure_ = true;
00429
00430 const int num_ifaces = ops_.internal_faces.size();
00431 threshold_pressures_by_connection_.resize(num_ifaces + num_nnc);
00432 for (int ii = 0; ii < num_ifaces; ++ii) {
00433 threshold_pressures_by_connection_[ii] = threshold_pressures[ops_.internal_faces[ii]];
00434 }
00435
00436
00437 for (int ii = 0; ii < num_nnc; ++ii) {
00438 threshold_pressures_by_connection_[ii + num_ifaces] = threshold_pressures[ii + num_faces];
00439 }
00440
00441 }
00442
00443
00444
00445
00446 template <class Grid, class WellModel, class Implementation>
00447 BlackoilModelBase<Grid, WellModel, Implementation>::
00448 ReservoirResidualQuant::ReservoirResidualQuant()
00449 : accum(2, ADB::null())
00450 , mflux( ADB::null())
00451 , b ( ADB::null())
00452 , mu ( ADB::null())
00453 , rho ( ADB::null())
00454 , kr ( ADB::null())
00455 , dh ( ADB::null())
00456 , mob ( ADB::null())
00457 {
00458 }
00459
00460 template <class Grid, class WellModel, class Implementation>
00461 BlackoilModelBase<Grid, WellModel, Implementation>::
00462 SimulatorData::SimulatorData(int num_phases)
00463 : rq(num_phases)
00464 , rsSat(ADB::null())
00465 , rvSat(ADB::null())
00466 , soMax()
00467 , Pb()
00468 , Pd()
00469 , krnswdc_ow()
00470 , krnswdc_go()
00471 , pcswmdc_ow()
00472 , pcswmdc_go()
00473 , fip()
00474 {
00475 }
00476
00477
00478
00479
00480
00481 template <class Grid, class WellModel, class Implementation>
00482 void
00483 BlackoilModelBase<Grid, WellModel, Implementation>::
00484 makeConstantState(SolutionState& state) const
00485 {
00486
00487
00488
00489
00490 state.pressure = ADB::constant(state.pressure.value());
00491 state.temperature = ADB::constant(state.temperature.value());
00492 state.rs = ADB::constant(state.rs.value());
00493 state.rv = ADB::constant(state.rv.value());
00494 const int num_phases = state.saturation.size();
00495 for (int phaseIdx = 0; phaseIdx < num_phases; ++ phaseIdx) {
00496 state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value());
00497 }
00498 state.qs = ADB::constant(state.qs.value());
00499 state.bhp = ADB::constant(state.bhp.value());
00500 assert(state.canonical_phase_pressures.size() == static_cast<std::size_t>(Opm::BlackoilPhases::MaxNumPhases));
00501 for (int canphase = 0; canphase < Opm::BlackoilPhases::MaxNumPhases; ++canphase) {
00502 ADB& pp = state.canonical_phase_pressures[canphase];
00503 pp = ADB::constant(pp.value());
00504 }
00505 }
00506
00507
00508
00509
00510
00511 template <class Grid, class WellModel, class Implementation>
00512 typename BlackoilModelBase<Grid, WellModel, Implementation>::SolutionState
00513 BlackoilModelBase<Grid, WellModel, Implementation>::
00514 variableState(const ReservoirState& x,
00515 const WellState& xw) const
00516 {
00517 std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
00518 std::vector<ADB> vars = ADB::variables(vars0);
00519 return asImpl().variableStateExtractVars(x, asImpl().variableStateIndices(), vars);
00520 }
00521
00522
00523
00524
00525
00526 template <class Grid, class WellModel, class Implementation>
00527 std::vector<V>
00528 BlackoilModelBase<Grid, WellModel, Implementation>::
00529 variableStateInitials(const ReservoirState& x,
00530 const WellState& xw) const
00531 {
00532 assert(active_[ Oil ]);
00533
00534 const int np = x.numPhases();
00535
00536 std::vector<V> vars0;
00537
00538
00539 vars0.reserve(np + 1);
00540 variableReservoirStateInitials(x, vars0);
00541 asImpl().wellModel().variableWellStateInitials(xw, vars0);
00542 return vars0;
00543 }
00544
00545
00546
00547
00548
00549 template <class Grid, class WellModel, class Implementation>
00550 void
00551 BlackoilModelBase<Grid, WellModel, Implementation>::
00552 variableReservoirStateInitials(const ReservoirState& x, std::vector<V>& vars0) const
00553 {
00554 using namespace Opm::AutoDiffGrid;
00555 const int nc = numCells(grid_);
00556 const int np = x.numPhases();
00557
00558 assert (not x.pressure().empty());
00559 const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
00560 vars0.push_back(p);
00561
00562
00563 assert (not x.saturation().empty());
00564 const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
00565 const Opm::PhaseUsage pu = fluid_.phaseUsage();
00566
00567 assert (active_[ Oil]);
00568 if (active_[ Water ]) {
00569 const V sw = s.col(pu.phase_pos[ Water ]);
00570 vars0.push_back(sw);
00571 }
00572
00573 if (active_[ Gas ]) {
00574
00575 V xvar(nc);
00576 const V sg = s.col(pu.phase_pos[ Gas ]);
00577 const V rs = Eigen::Map<const V>(& x.gasoilratio()[0], x.gasoilratio().size());
00578 const V rv = Eigen::Map<const V>(& x.rv()[0], x.rv().size());
00579 xvar = isRs_*rs + isRv_*rv + isSg_*sg;
00580 vars0.push_back(xvar);
00581 }
00582 }
00583
00584
00585
00586
00587
00588 template <class Grid, class WellModel, class Implementation>
00589 std::vector<int>
00590 BlackoilModelBase<Grid, WellModel, Implementation>::
00591 variableStateIndices() const
00592 {
00593 assert(active_[Oil]);
00594 std::vector<int> indices(5, -1);
00595 int next = 0;
00596 indices[Pressure] = next++;
00597 if (active_[Water]) {
00598 indices[Sw] = next++;
00599 }
00600 if (active_[Gas]) {
00601 indices[Xvar] = next++;
00602 }
00603 asImpl().wellModel().variableStateWellIndices(indices, next);
00604 assert(next == fluid_.numPhases() + 2);
00605 return indices;
00606 }
00607
00608
00609
00610
00611
00612 template <class Grid, class WellModel, class Implementation>
00613 typename BlackoilModelBase<Grid, WellModel, Implementation>::SolutionState
00614 BlackoilModelBase<Grid, WellModel, Implementation>::
00615 variableStateExtractVars(const ReservoirState& x,
00616 const std::vector<int>& indices,
00617 std::vector<ADB>& vars) const
00618 {
00619
00620 const int nc = Opm::AutoDiffGrid::numCells(grid_);
00621 const Opm::PhaseUsage pu = fluid_.phaseUsage();
00622
00623 SolutionState state(fluid_.numPhases());
00624
00625
00626 state.pressure = std::move(vars[indices[Pressure]]);
00627
00628
00629 const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
00630 state.temperature = ADB::constant(temp);
00631
00632
00633 {
00634 ADB so = ADB::constant(V::Ones(nc, 1));
00635
00636 if (active_[ Water ]) {
00637 state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]);
00638 const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
00639 so -= sw;
00640 }
00641
00642 if (active_[ Gas ]) {
00643
00644
00645 const ADB& xvar = vars[indices[Xvar]];
00646 ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
00647 sg = isSg_*xvar + isRv_*so;
00648 so -= sg;
00649
00650
00651 {
00652 const ADB& sw = (active_[ Water ]
00653 ? state.saturation[ pu.phase_pos[ Water ] ]
00654 : ADB::null());
00655 state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
00656 }
00657
00658 if (active_[ Oil ]) {
00659
00660 sd_.rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
00661 if (has_disgas_) {
00662 state.rs = (1-isRs_)*sd_.rsSat + isRs_*xvar;
00663 } else {
00664 state.rs = sd_.rsSat;
00665 }
00666 sd_.rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_);
00667 if (has_vapoil_) {
00668 state.rv = (1-isRv_)*sd_.rvSat + isRv_*xvar;
00669 } else {
00670 state.rv = sd_.rvSat;
00671 }
00672 sd_.soMax = fluid_.satOilMax();
00673 fluid_.getGasOilHystParams(sd_.krnswdc_go, sd_.pcswmdc_go, cells_);
00674 fluid_.getOilWaterHystParams(sd_.krnswdc_ow, sd_.pcswmdc_ow, cells_);
00675
00676 sd_.Pb = fluid_.bubblePointPressure(cells_,
00677 state.temperature.value(),
00678 state.rs.value());
00679 sd_.Pd = fluid_.dewPointPressure(cells_,
00680 state.temperature.value(),
00681 state.rv.value());
00682 }
00683 }
00684 else {
00685
00686 const ADB& sw = (active_[ Water ]
00687 ? state.saturation[ pu.phase_pos[ Water ] ]
00688 : ADB::null());
00689 const ADB& sg = ADB::null();
00690 state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
00691 }
00692
00693 if (active_[ Oil ]) {
00694
00695 state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
00696 }
00697 }
00698
00699 asImpl().wellModel().variableStateExtractWellsVars(indices, vars, state);
00700 return state;
00701 }
00702
00703
00704
00705
00706
00707 template <class Grid, class WellModel, class Implementation>
00708 void
00709 BlackoilModelBase<Grid, WellModel, Implementation>::
00710 computeAccum(const SolutionState& state,
00711 const int aix )
00712 {
00713 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00714
00715 const ADB& press = state.pressure;
00716 const ADB& temp = state.temperature;
00717 const std::vector<ADB>& sat = state.saturation;
00718 const ADB& rs = state.rs;
00719 const ADB& rv = state.rv;
00720
00721 const std::vector<PhasePresence> cond = phaseCondition();
00722
00723 const ADB pv_mult = poroMult(press);
00724
00725 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
00726 for (int phase = 0; phase < maxnp; ++phase) {
00727 if (active_[ phase ]) {
00728 const int pos = pu.phase_pos[ phase ];
00729 sd_.rq[pos].b = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
00730 sd_.rq[pos].accum[aix] = pv_mult * sd_.rq[pos].b * sat[pos];
00731
00732
00733 }
00734 }
00735
00736 if (active_[ Oil ] && active_[ Gas ]) {
00737
00738 const int po = pu.phase_pos[ Oil ];
00739 const int pg = pu.phase_pos[ Gas ];
00740
00741
00742
00743 const ADB accum_gas_copy =sd_.rq[pg].accum[aix];
00744
00745 sd_.rq[pg].accum[aix] += state.rs * sd_.rq[po].accum[aix];
00746 sd_.rq[po].accum[aix] += state.rv * accum_gas_copy;
00747
00748 }
00749 }
00750
00751
00752
00753
00754
00755 template <class Grid, class WellModel, class Implementation>
00756 SimulatorReport
00757 BlackoilModelBase<Grid, WellModel, Implementation>::
00758 assemble(const ReservoirState& reservoir_state,
00759 WellState& well_state,
00760 const bool initial_assembly)
00761 {
00762 using namespace Opm::AutoDiffGrid;
00763
00764 SimulatorReport report;
00765
00766
00767
00768
00769 if (isVFPActive()) {
00770 SolutionState state = asImpl().variableState(reservoir_state, well_state);
00771 SolutionState state0 = state;
00772 asImpl().makeConstantState(state0);
00773 asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
00774 }
00775
00776
00777 if (asImpl().wellModel().wellCollection()->groupControlActive() && initial_assembly) {
00778 setupGroupControl(reservoir_state, well_state);
00779 }
00780
00781
00782
00783 asImpl().wellModel().updateWellControls(well_state);
00784
00785 if (asImpl().wellModel().wellCollection()->groupControlActive()) {
00786
00787 applyVREPGroupControl(reservoir_state, well_state);
00788
00789 asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
00790 }
00791
00792
00793 SolutionState state = asImpl().variableState(reservoir_state, well_state);
00794
00795 if (initial_assembly) {
00796
00797 SolutionState state0 = state;
00798 asImpl().makeConstantState(state0);
00799
00800
00801 asImpl().computeAccum(state0, 0);
00802 asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
00803 }
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815 asImpl().assembleMassBalanceEq(state);
00816
00817
00818 if ( ! wellsActive() ) {
00819 return report;
00820 }
00821
00822 std::vector<ADB> mob_perfcells;
00823 std::vector<ADB> b_perfcells;
00824 asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
00825 if (param_.solve_welleq_initially_ && initial_assembly) {
00826
00827 report += asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
00828 }
00829 V aliveWells;
00830 std::vector<ADB> cq_s;
00831 asImpl().wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
00832 asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
00833 asImpl().wellModel().addWellFluxEq(cq_s, state, residual_);
00834 asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
00835 asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
00836
00837 return report;
00838 }
00839
00840
00841
00842
00843 template <class Grid, class WellModel, class Implementation>
00844 void
00845 BlackoilModelBase<Grid, WellModel, Implementation>::
00846 assembleMassBalanceEq(const SolutionState& state)
00847 {
00848
00849
00850
00851
00852
00853
00854 asImpl().computeAccum(state, 1);
00855
00856
00857
00858 const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
00859 const V trans_nnc = ops_.nnc_trans;
00860 V trans_all(transi.size() + trans_nnc.size());
00861 trans_all << transi, trans_nnc;
00862
00863
00864 {
00865 const std::vector<ADB> kr = asImpl().computeRelPerm(state);
00866 for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
00867 sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
00868 }
00869 }
00870 #pragma omp parallel for schedule(static)
00871 for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
00872 const std::vector<PhasePresence>& cond = phaseCondition();
00873 sd_.rq[phaseIdx].mu = asImpl().fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
00874 sd_.rq[phaseIdx].rho = asImpl().fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
00875 asImpl().computeMassFlux(phaseIdx, trans_all, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
00876
00877 residual_.material_balance_eq[ phaseIdx ] =
00878 pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
00879 + ops_.div*sd_.rq[phaseIdx].mflux;
00880 }
00881
00882
00883
00884
00885
00886
00887 if (active_[ Oil ] && active_[ Gas ]) {
00888 const int po = fluid_.phaseUsage().phase_pos[ Oil ];
00889 const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
00890
00891 const UpwindSelector<double> upwindOil(grid_, ops_,
00892 sd_.rq[po].dh.value());
00893 const ADB rs_face = upwindOil.select(state.rs);
00894
00895 const UpwindSelector<double> upwindGas(grid_, ops_,
00896 sd_.rq[pg].dh.value());
00897 const ADB rv_face = upwindGas.select(state.rv);
00898
00899 residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
00900 residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
00901
00902
00903
00904 }
00905
00906
00907 if (param_.update_equations_scaling_) {
00908 asImpl().updateEquationsScaling();
00909 }
00910
00911 }
00912
00913
00914
00915
00916
00917 template <class Grid, class WellModel, class Implementation>
00918 void
00919 BlackoilModelBase<Grid, WellModel, Implementation>::
00920 updateEquationsScaling() {
00921 ADB::V B;
00922 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00923 for ( int idx=0; idx<MaxNumPhases; ++idx )
00924 {
00925 if (active_[idx]) {
00926 const int pos = pu.phase_pos[idx];
00927 const ADB& temp_b = sd_.rq[pos].b;
00928 B = 1. / temp_b.value();
00929 #if HAVE_MPI
00930 if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
00931 {
00932 const ParallelISTLInformation& real_info =
00933 boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
00934 double B_global_sum = 0;
00935 real_info.computeReduction(B, Reduction::makeGlobalSumFunctor<double>(), B_global_sum);
00936 residual_.matbalscale[idx] = B_global_sum / global_nc_;
00937 }
00938 else
00939 #endif
00940 {
00941 residual_.matbalscale[idx] = B.mean();
00942 }
00943 }
00944 }
00945 }
00946
00947
00948
00949
00950
00951 template <class Grid, class WellModel, class Implementation>
00952 void
00953 BlackoilModelBase<Grid, WellModel, Implementation>::
00954 addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
00955 const SolutionState&,
00956 const WellState&)
00957 {
00958 if ( !asImpl().localWellsActive() )
00959 {
00960
00961
00962 return;
00963 }
00964
00965
00966 const int nc = Opm::AutoDiffGrid::numCells(grid_);
00967 const int np = asImpl().numPhases();
00968 const V& efficiency_factors = wellModel().wellPerfEfficiencyFactors();
00969 for (int phase = 0; phase < np; ++phase) {
00970 residual_.material_balance_eq[phase] -= superset(efficiency_factors * cq_s[phase],
00971 wellModel().wellOps().well_cells, nc);
00972 }
00973 }
00974
00975
00976
00977
00978
00979 template <class Grid, class WellModel, class Implementation>
00980 bool
00981 BlackoilModelBase<Grid, WellModel, Implementation>::
00982 isVFPActive() const
00983 {
00984 if( ! localWellsActive() ) {
00985 return false;
00986 }
00987
00988 if ( vfp_properties_.getProd()->empty() && vfp_properties_.getInj()->empty() ) {
00989 return false;
00990 }
00991
00992 const int nw = wells().number_of_wells;
00993
00994 for (int w = 0; w < nw; ++w) {
00995 const WellControls* wc = wells().ctrls[w];
00996
00997 const int nwc = well_controls_get_num(wc);
00998
00999
01000 for (int c=0; c < nwc; ++c) {
01001 const WellControlType ctrl_type = well_controls_iget_type(wc, c);
01002
01003 if (ctrl_type == THP) {
01004 return true;
01005 }
01006 }
01007 }
01008
01009 return false;
01010 }
01011
01012
01013
01014
01015 template <class Grid, class WellModel, class Implementation>
01016 SimulatorReport
01017 BlackoilModelBase<Grid, WellModel, Implementation>::
01018 solveWellEq(const std::vector<ADB>& mob_perfcells,
01019 const std::vector<ADB>& b_perfcells,
01020 const ReservoirState& reservoir_state,
01021 SolutionState& state,
01022 WellState& well_state)
01023 {
01024 V aliveWells;
01025 const int np = wells().number_of_phases;
01026 std::vector<ADB> cq_s(np, ADB::null());
01027 std::vector<int> indices = asImpl().wellModel().variableWellStateIndices();
01028 SolutionState state0 = state;
01029 WellState well_state0 = well_state;
01030 asImpl().makeConstantState(state0);
01031
01032 std::vector<ADB> mob_perfcells_const(np, ADB::null());
01033 std::vector<ADB> b_perfcells_const(np, ADB::null());
01034
01035 if (asImpl().localWellsActive() ){
01036
01037
01038 for (int phase = 0; phase < np; ++phase) {
01039 mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value());
01040 b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value());
01041 }
01042 }
01043
01044 int it = 0;
01045 bool converged;
01046 do {
01047
01048 std::vector<V> vars0;
01049 vars0.reserve(2);
01050 asImpl().wellModel().variableWellStateInitials(well_state, vars0);
01051 std::vector<ADB> vars = ADB::variables(vars0);
01052
01053 SolutionState wellSolutionState = state0;
01054 asImpl().wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState);
01055 asImpl().wellModel().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
01056 asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
01057 asImpl().wellModel().addWellFluxEq(cq_s, wellSolutionState, residual_);
01058 asImpl().wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_);
01059 converged = getWellConvergence(it);
01060
01061 if (converged) {
01062 break;
01063 }
01064
01065 ++it;
01066 if( localWellsActive() )
01067 {
01068 std::vector<ADB> eqs;
01069 eqs.reserve(2);
01070 eqs.push_back(residual_.well_flux_eq);
01071 eqs.push_back(residual_.well_eq);
01072 ADB total_residual = vertcatCollapseJacs(eqs);
01073 const std::vector<M>& Jn = total_residual.derivative();
01074 typedef Eigen::SparseMatrix<double> Sp;
01075 Sp Jn0;
01076 Jn[0].toSparse(Jn0);
01077 const Eigen::SparseLU< Sp > solver(Jn0);
01078 ADB::V total_residual_v = total_residual.value();
01079 const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
01080 assert(dx.size() == total_residual_v.size());
01081 asImpl().wellModel().updateWellState(dx.array(), dbhpMaxRel(), well_state);
01082 }
01083
01084
01085
01086 asImpl().wellModel().updateWellControls(well_state);
01087
01088 if (asImpl().wellModel().wellCollection()->groupControlActive()) {
01089
01090 applyVREPGroupControl(reservoir_state, well_state);
01091 asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
01092 }
01093 } while (it < 15);
01094
01095 if (converged) {
01096 if (terminalOutputEnabled())
01097 {
01098 OpmLog::note("well converged iter: " + std::to_string(it));
01099 }
01100
01101 const int nw = wells().number_of_wells;
01102 {
01103
01104
01105 ADB::V new_bhp = Eigen::Map<ADB::V>(well_state.bhp().data(), nw);
01106
01107
01108 std::vector<ADB::M> old_derivs = state.bhp.derivative();
01109 state.bhp = ADB::function(std::move(new_bhp), std::move(old_derivs));
01110 }
01111 {
01112
01113
01114
01115 const DataBlock wrates = Eigen::Map<const DataBlock>(well_state.wellRates().data(), nw, np).transpose();
01116 ADB::V new_qs = Eigen::Map<const V>(wrates.data(), nw*np);
01117 std::vector<ADB::M> old_derivs = state.qs.derivative();
01118 state.qs = ADB::function(std::move(new_qs), std::move(old_derivs));
01119 }
01120 asImpl().computeWellConnectionPressures(state, well_state);
01121 }
01122
01123 if (!converged) {
01124 well_state = well_state0;
01125 }
01126
01127 SimulatorReport report;
01128 report.total_well_iterations = it;
01129 report.converged = converged;
01130 return report;
01131 }
01132
01133
01134
01135
01136
01137 template <class Grid, class WellModel, class Implementation>
01138 V
01139 BlackoilModelBase<Grid, WellModel, Implementation>::
01140 solveJacobianSystem() const
01141 {
01142 return linsolver_.computeNewtonIncrement(residual_);
01143 }
01144
01145 template <class Grid, class WellModel, class Implementation>
01146 void
01147 BlackoilModelBase<Grid, WellModel, Implementation>::
01148 updateState(const V& dx,
01149 ReservoirState& reservoir_state,
01150 WellState& well_state)
01151 {
01152 using namespace Opm::AutoDiffGrid;
01153 const int np = fluid_.numPhases();
01154 const int nc = numCells(grid_);
01155 const V null;
01156 assert(null.size() == 0);
01157 const V zero = V::Zero(nc);
01158 const V ones = V::Constant(nc,1.0);
01159
01160
01161 const V dp = subset(dx, Span(nc));
01162 int varstart = nc;
01163 const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
01164 varstart += dsw.size();
01165
01166 const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
01167 varstart += dxvar.size();
01168
01169
01170 const V dwells = subset(dx, Span(asImpl().wellModel().numWellVars(), 1, varstart));
01171 varstart += dwells.size();
01172
01173 assert(varstart == dx.size());
01174
01175
01176 const double dpmaxrel = dpMaxRel();
01177 const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
01178 const V absdpmax = dpmaxrel*p_old.abs();
01179 const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
01180 const V p = (p_old - dp_limited).max(zero);
01181 std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
01182
01183
01184 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
01185 const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
01186 const double dsmax = dsMax();
01187
01188 V so;
01189 V sw;
01190 V sg;
01191
01192 {
01193 V maxVal = zero;
01194 V dso = zero;
01195 if (active_[Water]){
01196 maxVal = dsw.abs().max(maxVal);
01197 dso = dso - dsw;
01198 }
01199
01200 V dsg;
01201 if (active_[Gas]){
01202 dsg = isSg_ * dxvar - isRv_ * dsw;
01203 maxVal = dsg.abs().max(maxVal);
01204 dso = dso - dsg;
01205 }
01206
01207 maxVal = dso.abs().max(maxVal);
01208
01209 V step = dsmax/maxVal;
01210 step = step.min(1.);
01211
01212 if (active_[Water]) {
01213 const int pos = pu.phase_pos[ Water ];
01214 const V sw_old = s_old.col(pos);
01215 sw = sw_old - step * dsw;
01216 }
01217
01218 if (active_[Gas]) {
01219 const int pos = pu.phase_pos[ Gas ];
01220 const V sg_old = s_old.col(pos);
01221 sg = sg_old - step * dsg;
01222 }
01223
01224 assert(active_[Oil]);
01225 const int pos = pu.phase_pos[ Oil ];
01226 const V so_old = s_old.col(pos);
01227 so = so_old - step * dso;
01228 }
01229
01230 if (active_[Gas]) {
01231 auto ixg = sg < 0;
01232 for (int c = 0; c < nc; ++c) {
01233 if (ixg[c]) {
01234 if (active_[Water]) {
01235 sw[c] = sw[c] / (1-sg[c]);
01236 }
01237 so[c] = so[c] / (1-sg[c]);
01238 sg[c] = 0;
01239 }
01240 }
01241 }
01242
01243 if (active_[Oil]) {
01244 auto ixo = so < 0;
01245 for (int c = 0; c < nc; ++c) {
01246 if (ixo[c]) {
01247 if (active_[Water]) {
01248 sw[c] = sw[c] / (1-so[c]);
01249 }
01250 if (active_[Gas]) {
01251 sg[c] = sg[c] / (1-so[c]);
01252 }
01253 so[c] = 0;
01254 }
01255 }
01256 }
01257
01258 if (active_[Water]) {
01259 auto ixw = sw < 0;
01260 for (int c = 0; c < nc; ++c) {
01261 if (ixw[c]) {
01262 so[c] = so[c] / (1-sw[c]);
01263 if (active_[Gas]) {
01264 sg[c] = sg[c] / (1-sw[c]);
01265 }
01266 sw[c] = 0;
01267 }
01268 }
01269 }
01270
01271
01272 const double drmaxrel = drMaxRel();
01273 V rs;
01274 if (has_disgas_) {
01275 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
01276 const V drs = isRs_ * dxvar;
01277 const V drs_limited = sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1.0));
01278 rs = rs_old - drs_limited;
01279 rs = rs.max(zero);
01280 }
01281 V rv;
01282 if (has_vapoil_) {
01283 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
01284 const V drv = isRv_ * dxvar;
01285 const V drv_limited = sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-3));
01286 rv = rv_old - drv_limited;
01287 rv = rv.max(zero);
01288 }
01289
01290
01291 const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
01292 auto watOnly = sw > (1 - epsilon);
01293
01294
01295 std::vector<HydroCarbonState>& hydroCarbonState = reservoir_state.hydroCarbonState();
01296 std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
01297
01298 if (has_disgas_) {
01299 const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
01300 const V rsSat = fluidRsSat(p, so, cells_);
01301 sd_.rsSat = ADB::constant(rsSat);
01302
01303 auto hasGas = (sg > 0 && isRs_ == 0);
01304
01305
01306 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
01307 auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
01308 auto useSg = watOnly || hasGas || gasVaporized;
01309 for (int c = 0; c < nc; ++c) {
01310 if (useSg[c]) {
01311 rs[c] = rsSat[c];
01312 if (watOnly[c]) {
01313 so[c] = 0;
01314 sg[c] = 0;
01315 rs[c] = 0;
01316 }
01317 } else {
01318 hydroCarbonState[c] = HydroCarbonState::OilOnly;
01319 }
01320 }
01321
01322 }
01323
01324
01325 if (has_vapoil_) {
01326
01327
01328 const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
01329 const V gaspress = computeGasPressure(p, sw, so, sg);
01330 const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
01331 const V rvSat = fluidRvSat(gaspress, so, cells_);
01332 sd_.rvSat = ADB::constant(rvSat);
01333
01334
01335 auto hasOil = (so > 0 && isRv_ == 0);
01336
01337
01338 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
01339 auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
01340 auto useSg = watOnly || hasOil || oilCondensed;
01341 for (int c = 0; c < nc; ++c) {
01342 if (useSg[c]) {
01343 rv[c] = rvSat[c];
01344 if (watOnly[c]) {
01345 so[c] = 0;
01346 sg[c] = 0;
01347 rv[c] = 0;
01348 }
01349 } else {
01350 hydroCarbonState[c] = HydroCarbonState::GasOnly;
01351 }
01352 }
01353
01354 }
01355
01356
01357 if (active_[Water]) {
01358 for (int c = 0; c < nc; ++c) {
01359 reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
01360 }
01361 }
01362
01363 if (active_[Gas]) {
01364 for (int c = 0; c < nc; ++c) {
01365 reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
01366 }
01367 }
01368
01369 if (active_[ Oil ]) {
01370 for (int c = 0; c < nc; ++c) {
01371 reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
01372 }
01373 }
01374
01375 if (has_disgas_) {
01376 std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
01377 }
01378
01379 if (has_vapoil_) {
01380 std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
01381 }
01382
01383 asImpl().wellModel().updateWellState(dwells, dbhpMaxRel(), well_state);
01384
01385
01386 updatePhaseCondFromPrimalVariable(reservoir_state);
01387 }
01388
01389
01390
01391
01392
01393 template <class Grid, class WellModel, class Implementation>
01394 std::vector<ADB>
01395 BlackoilModelBase<Grid, WellModel, Implementation>::
01396 computeRelPerm(const SolutionState& state) const
01397 {
01398 using namespace Opm::AutoDiffGrid;
01399 const int nc = numCells(grid_);
01400
01401 const ADB zero = ADB::constant(V::Zero(nc));
01402
01403 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
01404 const ADB& sw = (active_[ Water ]
01405 ? state.saturation[ pu.phase_pos[ Water ] ]
01406 : zero);
01407
01408 const ADB& so = (active_[ Oil ]
01409 ? state.saturation[ pu.phase_pos[ Oil ] ]
01410 : zero);
01411
01412 const ADB& sg = (active_[ Gas ]
01413 ? state.saturation[ pu.phase_pos[ Gas ] ]
01414 : zero);
01415
01416 return fluid_.relperm(sw, so, sg, cells_);
01417 }
01418
01419
01420
01421
01422
01423 template <class Grid, class WellModel, class Implementation>
01424 std::vector<ADB>
01425 BlackoilModelBase<Grid, WellModel, Implementation>::
01426 computePressures(const ADB& po,
01427 const ADB& sw,
01428 const ADB& so,
01429 const ADB& sg) const
01430 {
01431
01432 std::vector<ADB> pressure = fluid_.capPress(sw, so, sg, cells_);
01433 for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
01434
01435 if (phaseIdx == BlackoilPhases::Liquid)
01436 continue;
01437 if (active_[phaseIdx]) {
01438 pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
01439 }
01440 }
01441
01442
01443
01444
01445
01446
01447 for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
01448 if (active_[phaseIdx]) {
01449 if (phaseIdx == BlackoilPhases::Aqua) {
01450 pressure[phaseIdx] = po - pressure[phaseIdx];
01451 } else {
01452 pressure[phaseIdx] += po;
01453 }
01454 }
01455 }
01456
01457 return pressure;
01458 }
01459
01460
01461
01462
01463
01464 template <class Grid, class WellModel, class Implementation>
01465 V
01466 BlackoilModelBase<Grid, WellModel, Implementation>::
01467 computeGasPressure(const V& po,
01468 const V& sw,
01469 const V& so,
01470 const V& sg) const
01471 {
01472 assert (active_[Gas]);
01473 std::vector<ADB> cp = fluid_.capPress(ADB::constant(sw),
01474 ADB::constant(so),
01475 ADB::constant(sg),
01476 cells_);
01477 return cp[Gas].value() + po;
01478 }
01479
01480
01481
01482 template <class Grid, class WellModel, class Implementation>
01483 void
01484 BlackoilModelBase<Grid, WellModel, Implementation>::
01485 computeMassFlux(const int actph ,
01486 const V& transi,
01487 const ADB& kr ,
01488 const ADB& mu ,
01489 const ADB& rho ,
01490 const ADB& phasePressure,
01491 const SolutionState& state)
01492 {
01493
01494 const ADB tr_mult = transMult(state.pressure);
01495 sd_.rq[ actph ].mob = tr_mult * kr / mu;
01496
01497
01498 const ADB rhoavg = ops_.caver * rho;
01499 sd_.rq[ actph ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
01500 if (use_threshold_pressure_) {
01501 applyThresholdPressures(sd_.rq[ actph ].dh);
01502 }
01503
01504
01505 const ADB& b = sd_.rq[ actph ].b;
01506 const ADB& mob = sd_.rq[ actph ].mob;
01507 const ADB& dh = sd_.rq[ actph ].dh;
01508 UpwindSelector<double> upwind(grid_, ops_, dh.value());
01509 sd_.rq[ actph ].mflux = upwind.select(b * mob) * (transi * dh);
01510 }
01511
01512
01513
01514
01515
01516 template <class Grid, class WellModel, class Implementation>
01517 void
01518 BlackoilModelBase<Grid, WellModel, Implementation>::
01519 applyThresholdPressures(ADB& dp)
01520 {
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532 const V high_potential = (dp.value().abs() >= threshold_pressures_by_connection_).template cast<double>();
01533
01534
01535 const M keep_high_potential(high_potential.matrix().asDiagonal());
01536
01537
01538 const V sign_dp = sign(dp.value());
01539 const V threshold_modification = sign_dp * threshold_pressures_by_connection_;
01540
01541
01542 dp = keep_high_potential * (dp - threshold_modification);
01543 }
01544
01545
01546
01547
01548
01549 template <class Grid, class WellModel, class Implementation>
01550 std::vector<double>
01551 BlackoilModelBase<Grid, WellModel, Implementation>::
01552 computeResidualNorms() const
01553 {
01554 std::vector<double> residualNorms;
01555
01556 std::vector<ADB>::const_iterator massBalanceIt = residual_.material_balance_eq.begin();
01557 const std::vector<ADB>::const_iterator endMassBalanceIt = residual_.material_balance_eq.end();
01558
01559 for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) {
01560 const double massBalanceResid = detail::infinityNorm( (*massBalanceIt),
01561 linsolver_.parallelInformation() );
01562 if (!std::isfinite(massBalanceResid)) {
01563 OPM_THROW(Opm::NumericalProblem,
01564 "Encountered a non-finite residual");
01565 }
01566 residualNorms.push_back(massBalanceResid);
01567 }
01568
01569
01570 const double wellFluxResid = detail::infinityNormWell( residual_.well_flux_eq,
01571 linsolver_.parallelInformation() );
01572 if (!std::isfinite(wellFluxResid)) {
01573 OPM_THROW(Opm::NumericalProblem,
01574 "Encountered a non-finite residual");
01575 }
01576 residualNorms.push_back(wellFluxResid);
01577
01578 const double wellResid = detail::infinityNormWell( residual_.well_eq,
01579 linsolver_.parallelInformation() );
01580 if (!std::isfinite(wellResid)) {
01581 OPM_THROW(Opm::NumericalProblem,
01582 "Encountered a non-finite residual");
01583 }
01584 residualNorms.push_back(wellResid);
01585
01586 return residualNorms;
01587 }
01588
01589
01590 template <class Grid, class WellModel, class Implementation>
01591 double
01592 BlackoilModelBase<Grid, WellModel, Implementation>::
01593 relativeChange(const SimulationDataContainer& previous,
01594 const SimulationDataContainer& current ) const
01595 {
01596 std::vector< double > p0 ( previous.pressure() );
01597 std::vector< double > sat0( previous.saturation() );
01598
01599 const std::size_t pSize = p0.size();
01600 const std::size_t satSize = sat0.size();
01601
01602
01603 for( std::size_t i=0; i<pSize; ++i ) {
01604 p0[ i ] -= current.pressure()[ i ];
01605 }
01606
01607 for( std::size_t i=0; i<satSize; ++i ) {
01608 sat0[ i ] -= current.saturation()[ i ];
01609 }
01610
01611
01612 const double stateOld = detail::euclidianNormSquared( p0.begin(), p0.end(), 1, linsolver_.parallelInformation() ) +
01613 detail::euclidianNormSquared( sat0.begin(), sat0.end(),
01614 current.numPhases(),
01615 linsolver_.parallelInformation() );
01616
01617
01618 const double stateNew = detail::euclidianNormSquared( current.pressure().begin(), current.pressure().end(), 1, linsolver_.parallelInformation() ) +
01619 detail::euclidianNormSquared( current.saturation().begin(), current.saturation().end(),
01620 current.numPhases(),
01621 linsolver_.parallelInformation() );
01622
01623 if( stateNew > 0.0 ) {
01624 return stateOld / stateNew ;
01625 }
01626 else {
01627 return 0.0;
01628 }
01629 }
01630
01631 template <class Grid, class WellModel, class Implementation>
01632 double
01633 BlackoilModelBase<Grid, WellModel, Implementation>::
01634 convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
01635 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
01636 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
01637 std::vector<double>& R_sum,
01638 std::vector<double>& maxCoeff,
01639 std::vector<double>& B_avg,
01640 std::vector<double>& maxNormWell,
01641 int nc) const
01642 {
01643 const int np = asImpl().numPhases();
01644 const int nm = asImpl().numMaterials();
01645 const int nw = residual_.well_flux_eq.size() / np;
01646 assert(nw * np == int(residual_.well_flux_eq.size()));
01647
01648
01649 #if HAVE_MPI
01650 if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
01651 {
01652 const ParallelISTLInformation& info =
01653 boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
01654
01655
01656 std::vector<int> v(nc, 1);
01657 auto nc_and_pv = std::tuple<int, double>(0, 0.0);
01658 auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
01659 Opm::Reduction::makeGlobalSumFunctor<double>());
01660 auto nc_and_pv_containers = std::make_tuple(v, geo_.poreVolume());
01661 info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
01662
01663 for ( int idx = 0; idx < nm; ++idx )
01664 {
01665 auto values = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
01666 auto containers = std::make_tuple(B.col(idx),
01667 tempV.col(idx),
01668 R.col(idx));
01669 auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
01670 Opm::Reduction::makeGlobalMaxFunctor<double>(),
01671 Opm::Reduction::makeGlobalSumFunctor<double>());
01672 info.computeReduction(containers, operators, values);
01673 B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
01674 maxCoeff[idx] = std::get<1>(values);
01675 R_sum[idx] = std::get<2>(values);
01676 assert(nm >= np);
01677 if (idx < np) {
01678 maxNormWell[idx] = 0.0;
01679 for ( int w = 0; w < nw; ++w ) {
01680 maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
01681 }
01682 }
01683 }
01684 info.communicator().max(maxNormWell.data(), np);
01685
01686 return std::get<1>(nc_and_pv);
01687 }
01688 else
01689 #endif
01690 {
01691 B_avg.resize(nm);
01692 maxCoeff.resize(nm);
01693 R_sum.resize(nm);
01694 maxNormWell.resize(np);
01695 for ( int idx = 0; idx < nm; ++idx )
01696 {
01697 B_avg[idx] = B.col(idx).sum()/nc;
01698 maxCoeff[idx] = tempV.col(idx).maxCoeff();
01699 R_sum[idx] = R.col(idx).sum();
01700
01701 assert(nm >= np);
01702 if (idx < np) {
01703 maxNormWell[idx] = 0.0;
01704 for ( int w = 0; w < nw; ++w ) {
01705 maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
01706 }
01707 }
01708 }
01709
01710 return geo_.poreVolume().sum();
01711 }
01712 }
01713
01714
01715
01716
01717
01718 template <class Grid, class WellModel, class Implementation>
01719 bool
01720 BlackoilModelBase<Grid, WellModel, Implementation>::
01721 getConvergence(const SimulatorTimerInterface& timer, const int iteration)
01722 {
01723 const double dt = timer.currentStepLength();
01724 const double tol_mb = param_.tolerance_mb_;
01725 const double tol_cnv = param_.tolerance_cnv_;
01726 const double tol_wells = param_.tolerance_wells_;
01727 const double tol_well_control = param_.tolerance_well_control_;
01728
01729 const int nc = Opm::AutoDiffGrid::numCells(grid_);
01730 const int np = asImpl().numPhases();
01731 const int nm = asImpl().numMaterials();
01732 assert(int(sd_.rq.size()) == nm);
01733
01734 const V& pv = geo_.poreVolume();
01735
01736 std::vector<double> R_sum(nm);
01737 std::vector<double> B_avg(nm);
01738 std::vector<double> maxCoeff(nm);
01739 std::vector<double> maxNormWell(np);
01740 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
01741 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
01742 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
01743
01744 for ( int idx = 0; idx < nm; ++idx )
01745 {
01746 const ADB& tempB = sd_.rq[idx].b;
01747 B.col(idx) = 1./tempB.value();
01748 R.col(idx) = residual_.material_balance_eq[idx].value();
01749 tempV.col(idx) = R.col(idx).abs()/pv;
01750 }
01751
01752 const double pvSum = convergenceReduction(B, tempV, R,
01753 R_sum, maxCoeff, B_avg, maxNormWell,
01754 nc);
01755
01756 std::vector<double> CNV(nm);
01757 std::vector<double> mass_balance_residual(nm);
01758 std::vector<double> well_flux_residual(np);
01759
01760 bool converged_MB = true;
01761 bool converged_CNV = true;
01762 bool converged_Well = true;
01763
01764 for ( int idx = 0; idx < nm; ++idx )
01765 {
01766 CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
01767 mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
01768 converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
01769 converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
01770
01771
01772 assert(nm >= np);
01773 if (idx < np) {
01774 well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
01775 converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
01776 }
01777 }
01778
01779 const double residualWell = detail::infinityNormWell(residual_.well_eq,
01780 linsolver_.parallelInformation());
01781 converged_Well = converged_Well && (residualWell < tol_well_control);
01782
01783 const bool converged = converged_MB && converged_CNV && converged_Well;
01784
01785
01786 const double maxWellResidualAllowed = 1000.0 * maxResidualAllowed();
01787
01788 if ( terminal_output_ )
01789 {
01790
01791 if (iteration == 0) {
01792 std::string msg = "Iter";
01793 for (int idx = 0; idx < nm; ++idx) {
01794 msg += " MB(" + materialName(idx).substr(0, 3) + ") ";
01795 }
01796 for (int idx = 0; idx < nm; ++idx) {
01797 msg += " CNV(" + materialName(idx).substr(0, 1) + ") ";
01798 }
01799 for (int idx = 0; idx < np; ++idx) {
01800 msg += " W-FLUX(" + materialName(idx).substr(0, 1) + ")";
01801 }
01802 msg += " WELL-CONT";
01803
01804 OpmLog::note(msg);
01805 }
01806 std::ostringstream ss;
01807 const std::streamsize oprec = ss.precision(3);
01808 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
01809 ss << std::setw(4) << iteration;
01810 for (int idx = 0; idx < nm; ++idx) {
01811 ss << std::setw(11) << mass_balance_residual[idx];
01812 }
01813 for (int idx = 0; idx < nm; ++idx) {
01814 ss << std::setw(11) << CNV[idx];
01815 }
01816 for (int idx = 0; idx < np; ++idx) {
01817 ss << std::setw(11) << well_flux_residual[idx];
01818 }
01819 ss << std::setw(11) << residualWell;
01820
01821 ss.precision(oprec);
01822 ss.flags(oflags);
01823 OpmLog::note(ss.str());
01824 }
01825
01826 for (int idx = 0; idx < nm; ++idx) {
01827 if (std::isnan(mass_balance_residual[idx])
01828 || std::isnan(CNV[idx])
01829 || (idx < np && std::isnan(well_flux_residual[idx]))) {
01830 const auto msg = std::string("NaN residual for phase ") + materialName(idx);
01831 if (terminal_output_) {
01832 OpmLog::bug(msg);
01833 }
01834 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
01835 }
01836 if (mass_balance_residual[idx] > maxResidualAllowed()
01837 || CNV[idx] > maxResidualAllowed()
01838 || (idx < np && well_flux_residual[idx] > maxResidualAllowed())) {
01839 const auto msg = std::string("Too large residual for phase ") + materialName(idx);
01840 if (terminal_output_) {
01841 OpmLog::problem(msg);
01842 }
01843 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
01844 }
01845 }
01846 if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) {
01847 const auto msg = std::string("NaN or too large residual for well control equation");
01848 if (terminal_output_) {
01849 OpmLog::problem(msg);
01850 }
01851 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
01852 }
01853
01854 return converged;
01855 }
01856
01857
01858
01859
01860
01861 template <class Grid, class WellModel, class Implementation>
01862 bool
01863 BlackoilModelBase<Grid, WellModel, Implementation>::
01864 getWellConvergence(const int iteration)
01865 {
01866 const double tol_wells = param_.tolerance_wells_;
01867 const double tol_well_control = param_.tolerance_well_control_;
01868
01869 const int nc = Opm::AutoDiffGrid::numCells(grid_);
01870 const int np = asImpl().numPhases();
01871 const int nm = asImpl().numMaterials();
01872
01873 const V& pv = geo_.poreVolume();
01874 std::vector<double> R_sum(nm);
01875 std::vector<double> B_avg(nm);
01876 std::vector<double> maxCoeff(nm);
01877 std::vector<double> maxNormWell(np);
01878 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
01879 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
01880 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
01881 for ( int idx = 0; idx < nm; ++idx )
01882 {
01883 const ADB& tempB = sd_.rq[idx].b;
01884 B.col(idx) = 1./tempB.value();
01885 R.col(idx) = residual_.material_balance_eq[idx].value();
01886 tempV.col(idx) = R.col(idx).abs()/pv;
01887 }
01888
01889 convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc);
01890
01891 std::vector<double> well_flux_residual(np);
01892 bool converged_Well = true;
01893
01894 for ( int idx = 0; idx < np; ++idx )
01895 {
01896 well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
01897 converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
01898 }
01899
01900 const double residualWell = detail::infinityNormWell(residual_.well_eq,
01901 linsolver_.parallelInformation());
01902 converged_Well = converged_Well && (residualWell < tol_well_control);
01903
01904 const bool converged = converged_Well;
01905
01906
01907 for (int idx = 0; idx < np; ++idx) {
01908 if (std::isnan(well_flux_residual[idx])) {
01909 const auto msg = std::string("NaN residual for phase ") + materialName(idx);
01910 if (terminal_output_) {
01911 OpmLog::bug(msg);
01912 }
01913 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
01914 }
01915 if (well_flux_residual[idx] > maxResidualAllowed()) {
01916 const auto msg = std::string("Too large residual for phase ") + materialName(idx);
01917 if (terminal_output_) {
01918 OpmLog::problem(msg);
01919 }
01920 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
01921 }
01922 }
01923
01924 if ( terminal_output_ )
01925 {
01926
01927 if (iteration == 0) {
01928 std::string msg;
01929 msg = "Iter";
01930 for (int idx = 0; idx < np; ++idx) {
01931 msg += " W-FLUX(" + materialName(idx).substr(0, 1) + ")";
01932 }
01933 msg += " WELL-CONT";
01934 OpmLog::note(msg);
01935 }
01936 std::ostringstream ss;
01937 const std::streamsize oprec = ss.precision(3);
01938 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
01939 ss << std::setw(4) << iteration;
01940 for (int idx = 0; idx < np; ++idx) {
01941 ss << std::setw(11) << well_flux_residual[idx];
01942 }
01943 ss << std::setw(11) << residualWell;
01944 ss.precision(oprec);
01945 ss.flags(oflags);
01946 OpmLog::note(ss.str());
01947 }
01948 return converged;
01949 }
01950
01951
01952
01953
01954
01955 template <class Grid, class WellModel, class Implementation>
01956 ADB
01957 BlackoilModelBase<Grid, WellModel, Implementation>::
01958 fluidViscosity(const int phase,
01959 const ADB& p ,
01960 const ADB& temp ,
01961 const ADB& rs ,
01962 const ADB& rv ,
01963 const std::vector<PhasePresence>& cond) const
01964 {
01965 switch (phase) {
01966 case Water:
01967 return fluid_.muWat(p, temp, cells_);
01968 case Oil:
01969 return fluid_.muOil(p, temp, rs, cond, cells_);
01970 case Gas:
01971 return fluid_.muGas(p, temp, rv, cond, cells_);
01972 default:
01973 OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
01974 }
01975 }
01976
01977
01978
01979
01980
01981 template <class Grid, class WellModel, class Implementation>
01982 ADB
01983 BlackoilModelBase<Grid, WellModel, Implementation>::
01984 fluidReciprocFVF(const int phase,
01985 const ADB& p ,
01986 const ADB& temp ,
01987 const ADB& rs ,
01988 const ADB& rv ,
01989 const std::vector<PhasePresence>& cond) const
01990 {
01991 switch (phase) {
01992 case Water:
01993 return fluid_.bWat(p, temp, cells_);
01994 case Oil:
01995 return fluid_.bOil(p, temp, rs, cond, cells_);
01996 case Gas:
01997 return fluid_.bGas(p, temp, rv, cond, cells_);
01998 default:
01999 OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
02000 }
02001 }
02002
02003
02004
02005
02006
02007 template <class Grid, class WellModel, class Implementation>
02008 ADB
02009 BlackoilModelBase<Grid, WellModel, Implementation>::
02010 fluidDensity(const int phase,
02011 const ADB& b,
02012 const ADB& rs,
02013 const ADB& rv) const
02014 {
02015 const V& rhos = fluid_.surfaceDensity(phase, cells_);
02016 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
02017 ADB rho = rhos * b;
02018 if (phase == Oil && active_[Gas]) {
02019 rho += fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * rs * b;
02020 }
02021 if (phase == Gas && active_[Oil]) {
02022 rho += fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * rv * b;
02023 }
02024 return rho;
02025 }
02026
02027
02028
02029
02030
02031 template <class Grid, class WellModel, class Implementation>
02032 V
02033 BlackoilModelBase<Grid, WellModel, Implementation>::
02034 fluidRsSat(const V& p,
02035 const V& satOil,
02036 const std::vector<int>& cells) const
02037 {
02038 return fluid_.rsSat(ADB::constant(p), ADB::constant(satOil), cells).value();
02039 }
02040
02041
02042
02043
02044
02045 template <class Grid, class WellModel, class Implementation>
02046 ADB
02047 BlackoilModelBase<Grid, WellModel, Implementation>::
02048 fluidRsSat(const ADB& p,
02049 const ADB& satOil,
02050 const std::vector<int>& cells) const
02051 {
02052 return fluid_.rsSat(p, satOil, cells);
02053 }
02054
02055
02056
02057
02058
02059 template <class Grid, class WellModel, class Implementation>
02060 V
02061 BlackoilModelBase<Grid, WellModel, Implementation>::
02062 fluidRvSat(const V& p,
02063 const V& satOil,
02064 const std::vector<int>& cells) const
02065 {
02066 return fluid_.rvSat(ADB::constant(p), ADB::constant(satOil), cells).value();
02067 }
02068
02069
02070
02071
02072
02073 template <class Grid, class WellModel, class Implementation>
02074 ADB
02075 BlackoilModelBase<Grid, WellModel, Implementation>::
02076 fluidRvSat(const ADB& p,
02077 const ADB& satOil,
02078 const std::vector<int>& cells) const
02079 {
02080 return fluid_.rvSat(p, satOil, cells);
02081 }
02082
02083
02084
02085
02086
02087 template <class Grid, class WellModel, class Implementation>
02088 ADB
02089 BlackoilModelBase<Grid, WellModel, Implementation>::
02090 poroMult(const ADB& p) const
02091 {
02092 const int n = p.size();
02093 if (rock_comp_props_ && rock_comp_props_->isActive()) {
02094 V pm(n);
02095 V dpm(n);
02096 #pragma omp parallel for schedule(static)
02097 for (int i = 0; i < n; ++i) {
02098 pm[i] = rock_comp_props_->poroMult(p.value()[i]);
02099 dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]);
02100 }
02101 ADB::M dpm_diag(dpm.matrix().asDiagonal());
02102 const int num_blocks = p.numBlocks();
02103 std::vector<ADB::M> jacs(num_blocks);
02104 #pragma omp parallel for schedule(dynamic)
02105 for (int block = 0; block < num_blocks; ++block) {
02106 fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]);
02107 }
02108 return ADB::function(std::move(pm), std::move(jacs));
02109 } else {
02110 return ADB::constant(V::Constant(n, 1.0));
02111 }
02112 }
02113
02114
02115
02116
02117
02118 template <class Grid, class WellModel, class Implementation>
02119 ADB
02120 BlackoilModelBase<Grid, WellModel, Implementation>::
02121 transMult(const ADB& p) const
02122 {
02123 const int n = p.size();
02124 if (rock_comp_props_ && rock_comp_props_->isActive()) {
02125 V tm(n);
02126 V dtm(n);
02127 #pragma omp parallel for schedule(static)
02128 for (int i = 0; i < n; ++i) {
02129 tm[i] = rock_comp_props_->transMult(p.value()[i]);
02130 dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]);
02131 }
02132 ADB::M dtm_diag(dtm.matrix().asDiagonal());
02133 const int num_blocks = p.numBlocks();
02134 std::vector<ADB::M> jacs(num_blocks);
02135 #pragma omp parallel for schedule(dynamic)
02136 for (int block = 0; block < num_blocks; ++block) {
02137 fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]);
02138 }
02139 return ADB::function(std::move(tm), std::move(jacs));
02140 } else {
02141 return ADB::constant(V::Constant(n, 1.0));
02142 }
02143 }
02144
02145
02146
02147
02148
02149 template <class Grid, class WellModel, class Implementation>
02150 void
02151 BlackoilModelBase<Grid, WellModel, Implementation>::
02152 classifyCondition(const ReservoirState& state)
02153 {
02154 using namespace Opm::AutoDiffGrid;
02155 const int nc = numCells(grid_);
02156 const int np = state.numPhases();
02157
02158 const PhaseUsage& pu = fluid_.phaseUsage();
02159 const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
02160 if (active_[ Gas ]) {
02161
02162 const V so = s.col(pu.phase_pos[ Oil ]);
02163 const V sg = s.col(pu.phase_pos[ Gas ]);
02164
02165 for (V::Index c = 0, e = sg.size(); c != e; ++c) {
02166 if (so[c] > 0) { phaseCondition_[c].setFreeOil (); }
02167 if (sg[c] > 0) { phaseCondition_[c].setFreeGas (); }
02168 if (active_[ Water ]) { phaseCondition_[c].setFreeWater(); }
02169 }
02170 }
02171 else {
02172
02173 assert (active_[ Water ]);
02174
02175 const V so = s.col(pu.phase_pos[ Oil ]);
02176
02177
02178 for (V::Index c = 0, e = so.size(); c != e; ++c) {
02179 phaseCondition_[c].setFreeWater();
02180
02181 if (so[c] > 0) { phaseCondition_[c].setFreeOil(); }
02182 }
02183 }
02184 }
02185
02186
02187
02188
02189
02190 template <class Grid, class WellModel, class Implementation>
02191 void
02192 BlackoilModelBase<Grid, WellModel, Implementation>::
02193 updatePrimalVariableFromState(const ReservoirState& state)
02194 {
02195 updatePhaseCondFromPrimalVariable(state);
02196 }
02197
02198
02199
02200
02201
02203 template <class Grid, class WellModel, class Implementation>
02204 void
02205 BlackoilModelBase<Grid, WellModel, Implementation>::
02206 updatePhaseCondFromPrimalVariable(const ReservoirState& state)
02207 {
02208 const int nc = Opm::AutoDiffGrid::numCells(grid_);
02209 isRs_ = V::Zero(nc);
02210 isRv_ = V::Zero(nc);
02211 isSg_ = V::Zero(nc);
02212
02213 if (! (active_[Gas] && active_[Oil])) {
02214
02215 phaseCondition_.assign(nc, PhasePresence());
02216 return;
02217 }
02218 for (int c = 0; c < nc; ++c) {
02219 phaseCondition_[c] = PhasePresence();
02220 phaseCondition_[c].setFreeWater();
02221 switch (state.hydroCarbonState()[c]) {
02222 case HydroCarbonState::GasAndOil:
02223 phaseCondition_[c].setFreeOil();
02224 phaseCondition_[c].setFreeGas();
02225 isSg_[c] = 1;
02226 break;
02227 case HydroCarbonState::OilOnly:
02228 phaseCondition_[c].setFreeOil();
02229 isRs_[c] = 1;
02230 break;
02231 case HydroCarbonState::GasOnly:
02232 phaseCondition_[c].setFreeGas();
02233 isRv_[c] = 1;
02234 break;
02235 default:
02236 OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << c << ": " << state.hydroCarbonState()[c]);
02237 }
02238 }
02239 }
02240
02241
02242
02243
02244
02245 template <class Grid, class WellModel, class Implementation>
02246 void
02247 BlackoilModelBase<Grid, WellModel, Implementation>::
02248 computeWellConnectionPressures(const SolutionState& state,
02249 const WellState& well_state)
02250 {
02251 asImpl().wellModel().computeWellConnectionPressures(state, well_state);
02252 }
02253
02254
02255
02256
02257
02258 template <class Grid, class WellModel, class Implementation>
02259 std::vector<std::vector<double> >
02260 BlackoilModelBase<Grid, WellModel, Implementation>::
02261 computeFluidInPlace(const ReservoirState& x,
02262 const std::vector<int>& fipnum)
02263 {
02264 using namespace Opm::AutoDiffGrid;
02265 const int nc = numCells(grid_);
02266 std::vector<ADB> saturation(3, ADB::null());
02267 const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, x.numPhases());
02268 const ADB pressure = ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1));
02269 const ADB temperature = ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1));
02270 saturation[Water] = active_[Water] ? ADB::constant(s.col(Water)) : ADB::null();
02271 saturation[Oil] = active_[Oil] ? ADB::constant(s.col(Oil)) : ADB::constant(V::Zero(nc));
02272 saturation[Gas] = active_[Gas] ? ADB::constant(s.col(Gas)) : ADB::constant(V::Zero(nc));
02273 const ADB rs = ADB::constant(Eigen::Map<const V>(& x.gasoilratio()[0], nc, 1));
02274 const ADB rv = ADB::constant(Eigen::Map<const V>(& x.rv()[0], nc, 1));
02275 const auto canonical_phase_pressures = computePressures(pressure, saturation[Water], saturation[Oil], saturation[Gas]);
02276 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
02277 const std::vector<PhasePresence> cond = phaseCondition();
02278
02279 const ADB pv_mult = poroMult(pressure);
02280 const V& pv = geo_.poreVolume();
02281 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
02282 for (int phase = 0; phase < maxnp; ++phase) {
02283 if (active_[ phase ]) {
02284 const int pos = pu.phase_pos[ phase ];
02285 const auto& b = asImpl().fluidReciprocFVF(phase, canonical_phase_pressures[phase], temperature, rs, rv, cond);
02286 sd_.fip[phase] = ((pv_mult * b * saturation[pos] * pv).value());
02287 }
02288 }
02289
02290 if (active_[ Oil ] && active_[ Gas ]) {
02291
02292 sd_.fip[SimulatorData::FIP_DISSOLVED_GAS] = rs.value() * sd_.fip[SimulatorData::FIP_LIQUID];
02293 sd_.fip[SimulatorData::FIP_VAPORIZED_OIL] = rv.value() * sd_.fip[SimulatorData::FIP_VAPOUR];
02294 }
02295
02296
02297 int dims = *std::max_element(fipnum.begin(), fipnum.end());
02298 std::vector<std::vector<double> > values(dims);
02299 for (int i=0; i < dims; ++i) {
02300 values[i].resize(7, 0.0);
02301 }
02302
02303 const V hydrocarbon = saturation[Oil].value() + saturation[Gas].value();
02304 V hcpv;
02305 V pres;
02306
02307 if ( !isParallel() )
02308 {
02309
02310 for (int phase = 0; phase < maxnp; ++phase) {
02311 if (active_[ phase ]) {
02312 for (int c = 0; c < nc; ++c) {
02313 const int region = fipnum[c] - 1;
02314 if (region != -1) {
02315 values[region][phase] += sd_.fip[phase][c];
02316 }
02317 }
02318 }
02319 }
02320
02321
02322 if (active_[ Oil ] && active_[ Gas ]) {
02323 for (int c = 0; c < nc; ++c) {
02324 const int region = fipnum[c] - 1;
02325 if (region != -1) {
02326 values[region][SimulatorData::FIP_DISSOLVED_GAS] += sd_.fip[SimulatorData::FIP_DISSOLVED_GAS][c];
02327 values[region][SimulatorData::FIP_VAPORIZED_OIL] += sd_.fip[SimulatorData::FIP_VAPORIZED_OIL][c];
02328 }
02329 }
02330 }
02331
02332
02333 hcpv = V::Zero(dims);
02334 pres = V::Zero(dims);
02335
02336 for (int c = 0; c < nc; ++c) {
02337 const int region = fipnum[c] - 1;
02338 if (region != -1) {
02339 hcpv[region] += pv[c] * hydrocarbon[c];
02340 pres[region] += pv[c] * pressure.value()[c];
02341 }
02342 }
02343
02344 sd_.fip[SimulatorData::FIP_PV] = V::Zero(nc);
02345 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE] = V::Zero(nc);
02346
02347 for (int c = 0; c < nc; ++c) {
02348 const int region = fipnum[c] - 1;
02349 if (region != -1) {
02350 sd_.fip[SimulatorData::FIP_PV][c] = pv[c];
02351
02352
02353
02354 if (hcpv[region] != 0) {
02355 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * pressure.value()[c] * hydrocarbon[c] / hcpv[region];
02356 } else {
02357 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
02358 }
02359
02360 values[region][SimulatorData::FIP_PV] += sd_.fip[SimulatorData::FIP_PV][c];
02361 values[region][SimulatorData::FIP_WEIGHTED_PRESSURE] += sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c];
02362 }
02363 }
02364 }
02365 else
02366 {
02367 #if HAVE_MPI
02368
02369 const auto & pinfo =
02370 boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
02371 const auto& mask = pinfo.getOwnerMask();
02372 auto comm = pinfo.communicator();
02373
02374 dims = comm.max(dims);
02375 values.resize(dims);
02376 for (int i=0; i < dims; ++i) {
02377 values[i].resize(7);
02378 std::fill(values[i].begin(), values[i].end(), 0.0);
02379 }
02380
02381
02382 for (int phase = 0; phase < maxnp; ++phase) {
02383 for (int c = 0; c < nc; ++c) {
02384 const int region = fipnum[c] - 1;
02385 if (region != -1 && mask[c]) {
02386 values[region][phase] += sd_.fip[phase][c];
02387 }
02388 }
02389 }
02390
02391
02392 if (active_[ Oil ] && active_[ Gas ]) {
02393 for (int c = 0; c < nc; ++c) {
02394 const int region = fipnum[c] - 1;
02395 if (region != -1 && mask[c]) {
02396 values[region][SimulatorData::FIP_DISSOLVED_GAS] += sd_.fip[SimulatorData::FIP_DISSOLVED_GAS][c];
02397 values[region][SimulatorData::FIP_VAPORIZED_OIL] += sd_.fip[SimulatorData::FIP_VAPORIZED_OIL][c];
02398 }
02399 }
02400 }
02401
02402 hcpv = V::Zero(dims);
02403 pres = V::Zero(dims);
02404
02405 for (int c = 0; c < nc; ++c) {
02406 const int region = fipnum[c] - 1;
02407 if (region != -1 && mask[c]) {
02408 hcpv[region] += pv[c] * hydrocarbon[c];
02409 pres[region] += pv[c] * pressure.value()[c];
02410 }
02411 }
02412
02413 comm.sum(hcpv.data(), hcpv.size());
02414 comm.sum(pres.data(), pres.size());
02415
02416 sd_.fip[SimulatorData::FIP_PV] = V::Zero(nc);
02417 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE] = V::Zero(nc);
02418
02419 for (int c = 0; c < nc; ++c) {
02420 const int region = fipnum[c] - 1;
02421 if (region != -1 && mask[c]) {
02422 sd_.fip[SimulatorData::FIP_PV][c] = pv[c];
02423
02424 if (hcpv[region] != 0) {
02425 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * pressure.value()[c] * hydrocarbon[c] / hcpv[region];
02426 } else {
02427 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
02428 }
02429
02430 values[region][SimulatorData::FIP_PV] += sd_.fip[SimulatorData::FIP_PV][c];
02431 values[region][SimulatorData::FIP_WEIGHTED_PRESSURE] += sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c];
02432 }
02433 }
02434
02435
02436
02437
02438 for(int reg=0; reg < dims; ++reg)
02439 {
02440 comm.sum(values[reg].data(), values[reg].size());
02441 }
02442 #else
02443
02444 OPM_THROW(std::logic_error, "HAVE_MPI should be defined if we are running in parallel");
02445 #endif
02446 }
02447
02448 return values;
02449 }
02450
02451
02452
02453
02454
02455 template <class Grid, class WellModel, class Implementation>
02456 void
02457 BlackoilModelBase<Grid, WellModel, Implementation>::
02458 computeWellVoidageRates(const ReservoirState& reservoir_state,
02459 const WellState& well_state,
02460 std::vector<double>& well_voidage_rates,
02461 std::vector<double>& voidage_conversion_coeffs)
02462 {
02463
02464
02465
02466
02467
02468 const int nw = well_state.numWells();
02469 const int np = numPhases();
02470
02471 const Wells* wells = asImpl().wellModel().wellsPointer();
02472
02473
02474 well_voidage_rates.resize(nw, 0);
02475
02476 voidage_conversion_coeffs.resize(nw * np, 1.0);
02477
02478 int global_number_wells = nw;
02479
02480 #if HAVE_MPI
02481 if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
02482 {
02483 const auto& info =
02484 boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
02485 global_number_wells = info.communicator().sum(global_number_wells);
02486 if ( global_number_wells )
02487 {
02488 rate_converter_.defineState(reservoir_state, boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation()));
02489 }
02490 }
02491 else
02492 #endif
02493 {
02494 if ( global_number_wells )
02495 {
02496 rate_converter_.defineState(reservoir_state);
02497 }
02498 }
02499
02500 std::vector<double> well_rates(np, 0.0);
02501 std::vector<double> convert_coeff(np, 1.0);
02502
02503
02504 if ( !well_voidage_rates.empty() ) {
02505 for (int w = 0; w < nw; ++w) {
02506 const bool is_producer = wells->type[w] == PRODUCER;
02507
02508
02509 if (is_producer) {
02510 std::transform(well_state.wellRates().begin() + np * w,
02511 well_state.wellRates().begin() + np * (w + 1),
02512 well_rates.begin(), std::negate<double>());
02513
02514
02515 const int fipreg = 0;
02516 const int well_cell_top = wells->well_cells[wells->well_connpos[w]];
02517 const int pvtreg = fluid_.cellPvtRegionIndex()[well_cell_top];
02518
02519 rate_converter_.calcCoeff(fipreg, pvtreg, convert_coeff);
02520 well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(),
02521 convert_coeff.begin(), 0.0);
02522 } else {
02523
02524
02525 std::copy(well_state.wellRates().begin() + np * w,
02526 well_state.wellRates().begin() + np * (w + 1),
02527 well_rates.begin());
02528
02529 const int fipreg = 0;
02530 const int well_cell_top = wells->well_cells[wells->well_connpos[w]];
02531 const int pvtreg = fluid_.cellPvtRegionIndex()[well_cell_top];
02532 rate_converter_.calcCoeff(fipreg, pvtreg, convert_coeff);
02533 std::copy(convert_coeff.begin(), convert_coeff.end(),
02534 voidage_conversion_coeffs.begin() + np * w);
02535 }
02536 }
02537 }
02538 }
02539
02540
02541
02542
02543
02544 template <class Grid, class WellModel, class Implementation>
02545 void
02546 BlackoilModelBase<Grid, WellModel, Implementation>::
02547 applyVREPGroupControl(const ReservoirState& reservoir_state,
02548 WellState& well_state)
02549 {
02550 if (asImpl().wellModel().wellCollection()->havingVREPGroups()) {
02551 std::vector<double> well_voidage_rates;
02552 std::vector<double> voidage_conversion_coeffs;
02553 computeWellVoidageRates(reservoir_state, well_state, well_voidage_rates, voidage_conversion_coeffs);
02554 asImpl().wellModel().wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs);
02555
02556
02557 for (const WellNode* well_node : asImpl().wellModel().wellCollection()->getLeafNodes()) {
02558 if (well_node->isInjector() && !well_node->individualControl()) {
02559 const int well_index = well_node->selfIndex();
02560 well_state.currentControls()[well_index] = well_node->groupControlIndex();
02561 }
02562 }
02563 }
02564 }
02565
02566
02567
02568
02569
02570 template <class Grid, class WellModel, class Implementation>
02571 void
02572 BlackoilModelBase<Grid, WellModel, Implementation>::
02573 setupGroupControl(const ReservoirState& reservoir_state,
02574 WellState& well_state)
02575 {
02576 if (asImpl().wellModel().wellCollection()->requireWellPotentials()) {
02577 SolutionState state = asImpl().variableState(reservoir_state, well_state);
02578 asImpl().makeConstantState(state);
02579 asImpl().wellModel().computeWellConnectionPressures(state, well_state);
02580
02581 const int np = numPhases();
02582 std::vector<ADB> b(np, ADB::null());
02583 std::vector<ADB> mob(np, ADB::null());
02584
02585 const ADB& press = state.pressure;
02586 const ADB& temp = state.temperature;
02587 const ADB& rs = state.rs;
02588 const ADB& rv = state.rv;
02589
02590 const std::vector<PhasePresence> cond = phaseCondition();
02591
02592 const ADB pv_mult = poroMult(press);
02593 const ADB tr_mult = transMult(press);
02594 const std::vector<ADB> kr = asImpl().computeRelPerm(state);
02595
02596 std::vector<ADB> mob_perfcells(np, ADB::null());
02597 std::vector<ADB> b_perfcells(np, ADB::null());
02598
02599 for (int phase = 0; phase < np; ++phase) {
02600 if (active_[phase]) {
02601 const std::vector<int>& well_cells = asImpl().wellModel().wellOps().well_cells;
02602 const ADB mu = asImpl().fluidViscosity(canph_[phase], state.canonical_phase_pressures[canph_[phase]],
02603 temp, rs, rv, cond);
02604 mob[phase] = tr_mult * kr[canph_[phase]] / mu;
02605 mob_perfcells[phase] = subset(mob[phase], well_cells);
02606
02607 b[phase] = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
02608 b_perfcells[phase] = subset(b[phase], well_cells);
02609 }
02610 }
02611
02612
02613 std::vector<double> well_potentials;
02614 asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, well_state, state, well_potentials);
02615 asImpl().wellModel().wellCollection()->setGuideRatesWithPotentials(asImpl().wellModel().wellsPointer(),
02616 fluid_.phaseUsage(), well_potentials);
02617 }
02618
02619 applyVREPGroupControl(reservoir_state, well_state);
02620
02621 if (asImpl().wellModel().wellCollection()->groupControlApplied()) {
02622 asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
02623 } else {
02624 asImpl().wellModel().wellCollection()->applyGroupControls();
02625
02626
02627
02628 const int nw = wells().number_of_wells;
02629 for (int w = 0; w < nw; ++w) {
02630 const WellNode& well_node = asImpl().wellModel().wellCollection()->findWellNode(wells().name[w]);
02631 if (!well_node.individualControl()) {
02632 well_state.currentControls()[w] = well_node.groupControlIndex();
02633 }
02634 }
02635 }
02636 }
02637
02638 }
02639
02640 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED