00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifndef OPM_BLACKOILMODELBASE_HEADER_INCLUDED
00024 #define OPM_BLACKOILMODELBASE_HEADER_INCLUDED
00025
00026 #include <cassert>
00027
00028 #include <opm/autodiff/AutoDiffBlock.hpp>
00029 #include <opm/autodiff/AutoDiffHelpers.hpp>
00030 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
00031 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
00032 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
00033 #include <opm/autodiff/BlackoilModelEnums.hpp>
00034 #include <opm/autodiff/VFPProperties.hpp>
00035 #include <opm/autodiff/RateConverter.hpp>
00036 #include <opm/autodiff/IterationReport.hpp>
00037 #include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
00038 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
00039 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
00040 #include <opm/core/simulator/SimulatorReport.hpp>
00041
00042 #include <opm/common/data/SimulationDataContainer.hpp>
00043
00044 #include <array>
00045
00046 struct Wells;
00047
00048 namespace Opm {
00049
00050 class ParameterGroup;
00051 class DerivedGeology;
00052 class RockCompressibility;
00053 class NewtonIterationBlackoilInterface;
00054 class VFPProperties;
00055
00059 template <class ConcreteModel>
00060 struct ModelTraits;
00061
00062
00075 template<class Grid, class WellModel, class Implementation>
00076 class BlackoilModelBase
00077 {
00078 public:
00079
00080 typedef AutoDiffBlock<double> ADB;
00081 typedef ADB::V V;
00082 typedef ADB::M M;
00083
00084 struct ReservoirResidualQuant {
00085 ReservoirResidualQuant();
00086 std::vector<ADB> accum;
00087 ADB mflux;
00088 ADB b;
00089 ADB mu;
00090 ADB rho;
00091 ADB kr;
00092 ADB dh;
00093 ADB mob;
00094 };
00095
00096 struct SimulatorData : public Opm::FIPDataEnums {
00097 SimulatorData(int num_phases);
00098
00099 using Opm::FIPDataEnums :: FipId ;
00100 using Opm::FIPDataEnums :: fipValues ;
00101
00102 std::vector<ReservoirResidualQuant> rq;
00103 ADB rsSat;
00104 ADB rvSat;
00105
00106 std::vector<double> soMax;
00107
00108 std::vector<double> Pb;
00109 std::vector<double> Pd;
00110
00111
00112 std::vector<double> krnswdc_ow;
00113 std::vector<double> krnswdc_go;
00114 std::vector<double> pcswmdc_ow;
00115 std::vector<double> pcswmdc_go;
00116
00117 std::array<V, fipValues> fip;
00118 };
00119
00120 typedef Opm::FIPData FIPDataType;
00121
00122 typedef typename ModelTraits<Implementation>::ReservoirState ReservoirState;
00123 typedef typename ModelTraits<Implementation>::WellState WellState;
00124 typedef typename ModelTraits<Implementation>::ModelParameters ModelParameters;
00125 typedef typename ModelTraits<Implementation>::SolutionState SolutionState;
00126
00127
00128 using RateConverterType = RateConverter::
00129 SurfaceToReservoirVoidage<BlackoilPropsAdFromDeck::FluidSystem, std::vector<int> >;
00130
00131
00132
00148 BlackoilModelBase(const ModelParameters& param,
00149 const Grid& grid ,
00150 const BlackoilPropsAdFromDeck& fluid,
00151 const DerivedGeology& geo ,
00152 const RockCompressibility* rock_comp_props,
00153 const WellModel& well_model,
00154 const NewtonIterationBlackoilInterface& linsolver,
00155 std::shared_ptr< const EclipseState > eclState,
00156 const bool has_disgas,
00157 const bool has_vapoil,
00158 const bool terminal_output);
00159
00168 void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
00169
00174 void prepareStep(const SimulatorTimerInterface& timer,
00175 const ReservoirState& reservoir_state,
00176 const WellState& well_state);
00177
00187 template <class NonlinearSolverType>
00188 SimulatorReport nonlinearIteration(const int iteration,
00189 const SimulatorTimerInterface& timer,
00190 NonlinearSolverType& nonlinear_solver,
00191 ReservoirState& reservoir_state,
00192 WellState& well_state);
00193
00199 void afterStep(const SimulatorTimerInterface& timer,
00200 ReservoirState& reservoir_state,
00201 WellState& well_state);
00202
00207 SimulatorReport
00208 assemble(const ReservoirState& reservoir_state,
00209 WellState& well_state,
00210 const bool initial_assembly);
00211
00216 std::vector<double> computeResidualNorms() const;
00217
00219
00220 double relativeChange( const SimulationDataContainer& previous, const SimulationDataContainer& current ) const;
00221
00223 int sizeNonLinear() const;
00224
00226 int linearIterationsLastSolve() const;
00227
00230 V solveJacobianSystem() const;
00231
00236 void updateState(const V& dx,
00237 ReservoirState& reservoir_state,
00238 WellState& well_state);
00239
00241 bool isParallel() const;
00242
00244 bool terminalOutputEnabled() const;
00245
00250 bool getConvergence(const SimulatorTimerInterface& timer, const int iteration);
00251
00253 int numPhases() const;
00254
00258 int numMaterials() const;
00259
00262 const std::string& materialName(int material_index) const;
00263
00265 void updateEquationsScaling();
00266
00268 WellModel& wellModel() { return well_model_; }
00269 const WellModel& wellModel() const { return well_model_; }
00270
00272 const SimulatorData& getSimulatorData(const SimulationDataContainer&) const {
00273 return sd_;
00274 }
00275
00277 FIPDataType getFIPData() const {
00278 return FIPDataType( sd_.fip );
00279 }
00280
00285 std::vector<std::vector<double> >
00286 computeFluidInPlace(const ReservoirState& x,
00287 const std::vector<int>& fipnum);
00288
00292 void computeWellVoidageRates(const ReservoirState& reservoir_state,
00293 const WellState& well_state,
00294 std::vector<double>& well_voidage_rates,
00295 std::vector<double>& voidage_conversion_coeffs);
00296
00297
00298 void applyVREPGroupControl(const ReservoirState& reservoir_state,
00299 WellState& well_state);
00300
00305 const SimulatorReport& failureReport() const
00306 { return failureReport_; }
00307
00308 protected:
00309
00310
00311
00312 typedef Eigen::Array<double,
00313 Eigen::Dynamic,
00314 Eigen::Dynamic,
00315 Eigen::RowMajor> DataBlock;
00316
00317
00318
00319
00320 SimulatorReport failureReport_;
00321 const Grid& grid_;
00322 const BlackoilPropsAdFromDeck& fluid_;
00323 const DerivedGeology& geo_;
00324 const RockCompressibility* rock_comp_props_;
00325 VFPProperties vfp_properties_;
00326 const NewtonIterationBlackoilInterface& linsolver_;
00327
00328 const std::vector<bool> active_;
00329
00330 const std::vector<int> canph_;
00331 const std::vector<int> cells_;
00332 HelperOps ops_;
00333 const bool has_disgas_;
00334 const bool has_vapoil_;
00335
00336 ModelParameters param_;
00337 bool use_threshold_pressure_;
00338 V threshold_pressures_by_connection_;
00339
00340 mutable SimulatorData sd_;
00341 std::vector<PhasePresence> phaseCondition_;
00342
00343
00344 WellModel well_model_;
00345
00346 V isRs_;
00347 V isRv_;
00348 V isSg_;
00349
00350 LinearisedBlackoilResidual residual_;
00351
00353 bool terminal_output_;
00355 int global_nc_;
00356
00357 V pvdt_;
00358 std::vector<std::string> material_name_;
00359 std::vector<std::vector<double>> residual_norms_history_;
00360 double current_relaxation_;
00361 V dx_old_;
00362
00363
00364 RateConverterType rate_converter_;
00365
00366
00367
00370 Implementation& asImpl()
00371 {
00372 return static_cast<Implementation&>(*this);
00373 }
00374
00377 const Implementation& asImpl() const
00378 {
00379 return static_cast<const Implementation&>(*this);
00380 }
00381
00383 const Wells& wells() const { return well_model_.wells(); }
00384
00386 bool wellsActive() const { return well_model_.wellsActive(); }
00387
00389 bool localWellsActive() const { return well_model_.localWellsActive(); }
00390
00391 void
00392 makeConstantState(SolutionState& state) const;
00393
00394 SolutionState
00395 variableState(const ReservoirState& x,
00396 const WellState& xw) const;
00397
00398 std::vector<V>
00399 variableStateInitials(const ReservoirState& x,
00400 const WellState& xw) const;
00401 void
00402 variableReservoirStateInitials(const ReservoirState& x,
00403 std::vector<V>& vars0) const;
00404
00405 std::vector<int>
00406 variableStateIndices() const;
00407
00408 SolutionState
00409 variableStateExtractVars(const ReservoirState& x,
00410 const std::vector<int>& indices,
00411 std::vector<ADB>& vars) const;
00412
00413 void
00414 computeAccum(const SolutionState& state,
00415 const int aix );
00416
00417 void
00418 assembleMassBalanceEq(const SolutionState& state);
00419
00420
00421 SimulatorReport
00422 solveWellEq(const std::vector<ADB>& mob_perfcells,
00423 const std::vector<ADB>& b_perfcells,
00424 const ReservoirState& reservoir_state,
00425 SolutionState& state,
00426 WellState& well_state);
00427
00428 void
00429 addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
00430 const SolutionState& state,
00431 const WellState& xw);
00432
00433 bool getWellConvergence(const int iteration);
00434
00435 bool isVFPActive() const;
00436
00437 std::vector<ADB>
00438 computePressures(const ADB& po,
00439 const ADB& sw,
00440 const ADB& so,
00441 const ADB& sg) const;
00442
00443 V
00444 computeGasPressure(const V& po,
00445 const V& sw,
00446 const V& so,
00447 const V& sg) const;
00448
00449 std::vector<ADB>
00450 computeRelPerm(const SolutionState& state) const;
00451
00452 void
00453 computeMassFlux(const int actph ,
00454 const V& transi,
00455 const ADB& kr ,
00456 const ADB& mu ,
00457 const ADB& rho ,
00458 const ADB& p ,
00459 const SolutionState& state );
00460
00461 void applyThresholdPressures(ADB& dp);
00462
00463 ADB
00464 fluidViscosity(const int phase,
00465 const ADB& p ,
00466 const ADB& temp ,
00467 const ADB& rs ,
00468 const ADB& rv ,
00469 const std::vector<PhasePresence>& cond) const;
00470
00471 ADB
00472 fluidReciprocFVF(const int phase,
00473 const ADB& p ,
00474 const ADB& temp ,
00475 const ADB& rs ,
00476 const ADB& rv ,
00477 const std::vector<PhasePresence>& cond) const;
00478
00479 ADB
00480 fluidDensity(const int phase,
00481 const ADB& b,
00482 const ADB& rs,
00483 const ADB& rv) const;
00484
00485 V
00486 fluidRsSat(const V& p,
00487 const V& so,
00488 const std::vector<int>& cells) const;
00489
00490 ADB
00491 fluidRsSat(const ADB& p,
00492 const ADB& so,
00493 const std::vector<int>& cells) const;
00494
00495 V
00496 fluidRvSat(const V& p,
00497 const V& so,
00498 const std::vector<int>& cells) const;
00499
00500 ADB
00501 fluidRvSat(const ADB& p,
00502 const ADB& so,
00503 const std::vector<int>& cells) const;
00504
00505 ADB
00506 poroMult(const ADB& p) const;
00507
00508 ADB
00509 transMult(const ADB& p) const;
00510
00511 const std::vector<PhasePresence>
00512 phaseCondition() const {return phaseCondition_;}
00513
00514 void
00515 classifyCondition(const ReservoirState& state);
00516
00517
00520 void
00521 updatePrimalVariableFromState(const ReservoirState& state);
00522
00525 void
00526 updatePhaseCondFromPrimalVariable(const ReservoirState& state);
00527
00528
00529
00530 void
00531 computeWellConnectionPressures(const SolutionState& state,
00532 const WellState& well_state);
00533
00554 double
00555 convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
00556 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
00557 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
00558 std::vector<double>& R_sum,
00559 std::vector<double>& maxCoeff,
00560 std::vector<double>& B_avg,
00561 std::vector<double>& maxNormWell,
00562 int nc) const;
00563
00565 void
00566 setupGroupControl(const ReservoirState& reservoir_state,
00567 WellState& well_state);
00568
00569 double dpMaxRel() const { return param_.dp_max_rel_; }
00570 double dbhpMaxRel() const {return param_.dbhp_max_rel_; }
00571 double dsMax() const { return param_.ds_max_; }
00572 double drMaxRel() const { return param_.dr_max_rel_; }
00573 double maxResidualAllowed() const { return param_.max_residual_allowed_; }
00574
00575 };
00576 }
00577
00578 #include "BlackoilModelBase_impl.hpp"
00579
00580 #endif // OPM_BLACKOILMODELBASE_HEADER_INCLUDED