00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OPM_BLACKOILREORDERINGTRANSPORTMODEL_HEADER_INCLUDED
00021 #define OPM_BLACKOILREORDERINGTRANSPORTMODEL_HEADER_INCLUDED
00022
00023 #include <opm/autodiff/BlackoilModelBase.hpp>
00024 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
00025 #include <opm/autodiff/BlackoilModelParameters.hpp>
00026 #include <opm/autodiff/DebugTimeReport.hpp>
00027 #include <opm/autodiff/multiPhaseUpwind.hpp>
00028 #include <opm/core/grid.h>
00029 #include <opm/core/transport/reorder/reordersequence.h>
00030 #include <opm/core/simulator/BlackoilState.hpp>
00031
00032 #include <opm/autodiff/BlackoilTransportModel.hpp>
00033
00034 namespace Opm {
00035
00036
00037 namespace detail
00038 {
00039
00040 template <typename Scalar>
00041 struct CreateVariable
00042 {
00043 Scalar operator()(double value, int index)
00044 {
00045 return Scalar::createVariable(value, index);
00046 }
00047 };
00048
00049
00050
00051 template <>
00052 struct CreateVariable<double>
00053 {
00054 double operator()(double value, int)
00055 {
00056 return value;
00057 }
00058 };
00059
00060
00061
00062 template <typename Scalar>
00063 struct CreateConstant
00064 {
00065 Scalar operator()(double value)
00066 {
00067 return Scalar::createConstant(value);
00068 }
00069 };
00070
00071
00072
00073 template <>
00074 struct CreateConstant<double>
00075 {
00076 double operator()(double value)
00077 {
00078 return value;
00079 }
00080 };
00081
00082
00083
00084
00085 struct Connection
00086 {
00087 Connection(const int ind, const double s) : index(ind), sign(s) {}
00088 int index;
00089 double sign;
00090 };
00091
00092
00093
00094
00095 class Connections;
00096
00097
00098
00099 class ConnectivityGraph
00100 {
00101 public:
00102 explicit ConnectivityGraph(const HelperOps& ops)
00103 : grad_(ops.grad)
00104 , div_(ops.div)
00105 {
00106 grad_ia_ = grad_.outerIndexPtr();
00107 grad_ja_ = grad_.innerIndexPtr();
00108 grad_sign_ = grad_.valuePtr();
00109 div_ia_ = div_.outerIndexPtr();
00110 div_ja_ = div_.innerIndexPtr();
00111 div_sign_ = div_.valuePtr();
00112 }
00113
00114 Connections cellConnections(const int cell) const;
00115
00116 std::array<int, 2> connectionCells(const int connection) const
00117 {
00118 const int pos = div_ia_[connection];
00119 assert(div_ia_[connection + 1] == pos + 2);
00120 const double sign1 = div_sign_[pos];
00121 assert(div_sign_[pos + 1] == -sign1);
00122 if (sign1 > 0.0) {
00123 return {{ div_ja_[pos], div_ja_[pos + 1] }};
00124 } else {
00125 return {{ div_ja_[pos + 1], div_ja_[pos] }};
00126 }
00127 }
00128
00129 private:
00130 friend class Connections;
00131 typedef HelperOps::M M;
00132 const M& grad_;
00133 const M& div_;
00134 const int* grad_ia_;
00135 const int* grad_ja_;
00136 const double* grad_sign_;
00137 const int* div_ia_;
00138 const int* div_ja_;
00139 const double* div_sign_;
00140 };
00141
00142
00143 class Connections
00144 {
00145 public:
00146 Connections(const ConnectivityGraph& cg, const int cell) : cg_(cg), cell_(cell) {}
00147 int size() const
00148 {
00149 return cg_.grad_ia_[cell_ + 1] - cg_.grad_ia_[cell_];
00150 }
00151 class Iterator
00152 {
00153 public:
00154 Iterator(const Connections& c, const int index) : c_(c), index_(index) {}
00155 Iterator& operator++()
00156 {
00157 ++index_;
00158 return *this;
00159 }
00160 bool operator!=(const Iterator& other)
00161 {
00162 assert(&c_ == &other.c_);
00163 return index_ != other.index_;
00164 }
00165 Connection operator*()
00166 {
00167 assert(index_ >= 0 && index_ < c_.size());
00168 const int pos = c_.cg_.grad_ia_[c_.cell_] + index_;
00169 return Connection(c_.cg_.grad_ja_[pos], -c_.cg_.grad_sign_[pos]);
00170 }
00171 private:
00172 const Connections& c_;
00173 int index_;
00174 };
00175 Iterator begin() const { return Iterator(*this, 0); }
00176 Iterator end() const { return Iterator(*this, size()); }
00177 private:
00178 friend class Iterator;
00179 const ConnectivityGraph& cg_;
00180 const int cell_;
00181 };
00182
00183
00184 inline Connections ConnectivityGraph::cellConnections(const int cell) const
00185 {
00186 return Connections(*this, cell);
00187 }
00188
00189
00190
00191 }
00192
00193
00194
00195
00196
00197
00198
00200 template<class Grid, class WellModel>
00201 class BlackoilReorderingTransportModel
00202 : public BlackoilModelBase<Grid, WellModel, BlackoilReorderingTransportModel<Grid, WellModel> >
00203 {
00204 public:
00205 typedef BlackoilModelBase<Grid, WellModel, BlackoilReorderingTransportModel<Grid, WellModel> > Base;
00206 friend Base;
00207
00208 typedef typename Base::ReservoirState ReservoirState;
00209 typedef typename Base::WellState WellState;
00210 typedef typename Base::SolutionState SolutionState;
00211 typedef typename Base::V V;
00212
00213
00228 BlackoilReorderingTransportModel(const typename Base::ModelParameters& param,
00229 const Grid& grid,
00230 const BlackoilPropsAdFromDeck& fluid,
00231 const DerivedGeology& geo,
00232 const RockCompressibility* rock_comp_props,
00233 const StandardWells& std_wells,
00234 const NewtonIterationBlackoilInterface& linsolver,
00235 std::shared_ptr<const EclipseState> eclState,
00236 const bool has_disgas,
00237 const bool has_vapoil,
00238 const bool terminal_output)
00239 : Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
00240 eclState, has_disgas, has_vapoil, terminal_output)
00241 , graph_(Base::ops_)
00242 , props_(dynamic_cast<const BlackoilPropsAdFromDeck&>(fluid))
00243 , state0_{ ReservoirState(0, 0, 0), WellState(), V(), V() }
00244 , state_{ ReservoirState(0, 0, 0), WellState(), V(), V() }
00245 , tr_model_(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
00246 eclState, has_disgas, has_vapoil, terminal_output)
00247 {
00248
00249
00250 const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
00251 const V trans_nnc = ops_.nnc_trans;
00252 trans_all_ = V::Zero(transi.size() + trans_nnc.size());
00253 trans_all_ << transi, trans_nnc;
00254 gdz_ = geo_.gravity()[2] * (ops_.grad * geo_.z().matrix());
00255 rhos_ = DataBlock::Zero(ops_.div.rows(), 3);
00256 rhos_.col(Water) = props_.surfaceDensity(Water, Base::cells_);
00257 rhos_.col(Oil) = props_.surfaceDensity(Oil, Base::cells_);
00258 rhos_.col(Gas) = props_.surfaceDensity(Gas, Base::cells_);
00259 }
00260
00261
00262
00263
00264
00265 void prepareStep(const SimulatorTimerInterface& timer,
00266 const ReservoirState& reservoir_state,
00267 const WellState& well_state)
00268 {
00269 tr_model_.prepareStep(timer, reservoir_state, well_state);
00270 Base::prepareStep(timer, reservoir_state, well_state);
00271 Base::param_.solve_welleq_initially_ = false;
00272 state0_.reservoir_state = reservoir_state;
00273 state0_.well_state = well_state;
00274
00275
00276 const std::vector<double>& p = reservoir_state.pressure();
00277 state0_.tr_mult = Base::transMult(ADB::constant(Eigen::Map<const V>(p.data(), p.size()))).value();
00278 state0_.pv_mult = Base::poroMult(ADB::constant(Eigen::Map<const V>(p.data(), p.size()))).value();
00279 const int num_cells = p.size();
00280 cstate0_.resize(num_cells);
00281 for (int cell = 0; cell < num_cells; ++cell) {
00282 computeCellState(cell, state0_, cstate0_[cell]);
00283 }
00284 cstate_ = cstate0_;
00285 }
00286
00287
00288
00289
00290
00291 template <class NonlinearSolverType>
00292 SimulatorReport nonlinearIteration(const int iteration,
00293 const SimulatorTimerInterface& timer,
00294 NonlinearSolverType& nonlinear_solver,
00295 ReservoirState& reservoir_state,
00296 const WellState& well_state)
00297 {
00298
00299 {
00300 DebugTimeReport tr("Extracting fluxes");
00301 extractFluxes(reservoir_state, well_state);
00302 extractState(reservoir_state, well_state);
00303 }
00304
00305
00306 {
00307 DebugTimeReport tr("Topological sort");
00308 computeOrdering();
00309 }
00310
00311
00312 {
00313 DebugTimeReport tr("Solving all components");
00314 for (int ii = 0; ii < 5; ++ii) {
00315 DebugTimeReport tr2("Solving components single sweep.");
00316 solveComponents();
00317 }
00318 }
00319
00320
00321 reservoir_state = state_.reservoir_state;
00322
00323
00324 {
00325 auto rs = reservoir_state;
00326 auto ws = well_state;
00327 tr_model_.nonlinearIteration( 0, timer, nonlinear_solver, rs, ws);
00328 }
00329
00330
00331 SimulatorReport report;
00332 report.converged = true;
00333 return report;
00334 }
00335
00336
00337
00338
00339
00340 void afterStep(const SimulatorTimerInterface& ,
00341 const ReservoirState& ,
00342 const WellState& )
00343 {
00344
00345 }
00346
00347 using Base::numPhases;
00348
00349
00350 protected:
00351
00352
00353
00354 using Vec2 = Dune::FieldVector<double, 2>;
00355 using Mat22 = Dune::FieldMatrix<double, 2, 2>;
00356 using Eval = DenseAd::Evaluation<double, 2>;
00357
00358
00359
00360
00361 struct State
00362 {
00363 ReservoirState reservoir_state;
00364 WellState well_state;
00365 V tr_mult;
00366 V pv_mult;
00367 };
00368
00369
00370
00371
00372
00373 template <typename ScalarT>
00374 struct CellState
00375 {
00376 using Scalar = ScalarT;
00377
00378 Scalar s[3];
00379 Scalar rs;
00380 Scalar rv;
00381 Scalar p[3];
00382 Scalar kr[3];
00383 Scalar pc[3];
00384 Scalar temperature;
00385 Scalar mu[3];
00386 Scalar b[3];
00387 Scalar lambda[3];
00388 Scalar rho[3];
00389 Scalar rssat;
00390 Scalar rvsat;
00391
00392
00393 const Scalar& saturation(int phaseIdx) const
00394 {
00395 return s[phaseIdx];
00396 }
00397
00398 template <typename T>
00399 CellState<T> flatten() const
00400 {
00401 return CellState<T>{
00402 { s[0].value(), s[1].value(), s[2].value() },
00403 rs.value(),
00404 rv.value(),
00405 { p[0].value(), p[1].value(), p[2].value() },
00406 { kr[0].value(), kr[1].value(), kr[2].value() },
00407 { pc[0].value(), pc[1].value(), pc[2].value() },
00408 temperature.value(),
00409 { mu[0].value(), mu[1].value(), mu[2].value() },
00410 { b[0].value(), b[1].value(), b[2].value() },
00411 { lambda[0].value(), lambda[1].value(), lambda[2].value() },
00412 { rho[0].value(), rho[1].value(), rho[2].value() },
00413 rssat.value(),
00414 rvsat.value()
00415 };
00416 }
00417 };
00418
00419
00420
00421
00422
00423 using Base::grid_;
00424 using Base::geo_;
00425 using Base::ops_;
00426
00427 const detail::ConnectivityGraph graph_;
00428
00429 const BlackoilPropsAdFromDeck& props_;
00430
00431 State state0_;
00432 State state_;
00433
00434 std::vector<CellState<double>> cstate0_;
00435 std::vector<CellState<double>> cstate_;
00436
00437 V total_flux_;
00438 V total_wellperf_flux_;
00439 DataBlock comp_wellperf_flux_;
00440 V total_wellflux_cell_;
00441 V oil_wellflux_cell_;
00442 V gas_wellflux_cell_;
00443 std::vector<int> sequence_;
00444 std::vector<int> components_;
00445 V trans_all_;
00446 V gdz_;
00447 DataBlock rhos_;
00448
00449 std::array<double, 2> max_abs_dx_;
00450 std::array<int, 2> max_abs_dx_cell_;
00451
00452
00453 BlackoilTransportModel<Grid, WellModel> tr_model_;
00454
00455
00456
00457
00458
00459
00460
00461
00462 template <typename Scalar>
00463 void computeCellState(const int cell, const State& state, CellState<Scalar>& cstate) const
00464 {
00465 assert(numPhases() == 3);
00466
00467
00468 const auto hcstate = state.reservoir_state.hydroCarbonState()[cell];
00469 const bool is_sg = (hcstate == HydroCarbonState::GasAndOil);
00470 const bool is_rs = (hcstate == HydroCarbonState::OilOnly);
00471 const bool is_rv = (hcstate == HydroCarbonState::GasOnly);
00472 const double swval = state.reservoir_state.saturation()[3*cell + Water];
00473 const double sgval = state.reservoir_state.saturation()[3*cell + Gas];
00474 const double rsval = state.reservoir_state.gasoilratio()[cell];
00475 const double rvval = state.reservoir_state.rv()[cell];
00476 const double poval = state.reservoir_state.pressure()[cell];
00477 const int pvt_region = props_.pvtRegions()[cell];
00478
00479
00480 const auto& waterpvt = props_.waterProps();
00481 const auto& oilpvt = props_.oilProps();
00482 const auto& gaspvt = props_.gasProps();
00483 const auto& satfunc = props_.materialLaws();
00484
00485
00486 detail::CreateVariable<Scalar> variable;
00487 detail::CreateConstant<Scalar> constant;
00488 cstate.s[Water] = variable(swval, 0);
00489 cstate.s[Gas] = is_sg ? variable(sgval, 1) : constant(sgval);
00490 cstate.s[Oil] = 1.0 - cstate.s[Water] - cstate.s[Gas];
00491 cstate.rs = is_rs ? variable(rsval, 1) : constant(rsval);
00492 cstate.rv = is_rv ? variable(rvval, 1) : constant(rvval);
00493
00494
00495 const auto& params = satfunc.materialLawParams(cell);
00496 typedef BlackoilPropsAdFromDeck::MaterialLawManager::MaterialLaw MaterialLaw;
00497 MaterialLaw::relativePermeabilities(cstate.kr, params, cstate);
00498 MaterialLaw::capillaryPressures(cstate.pc, params, cstate);
00499
00500
00501 cstate.p[Oil] = constant(poval);
00502 cstate.p[Water] = cstate.p[Oil] + cstate.pc[Water];
00503 cstate.p[Gas] = cstate.p[Oil] + cstate.pc[Gas];
00504
00505
00506 cstate.temperature = constant(0.0);
00507 cstate.mu[Water] = waterpvt.viscosity(pvt_region, cstate.temperature, cstate.p[Water]);
00508 cstate.mu[Oil] = is_sg
00509 ? oilpvt.saturatedViscosity(pvt_region, cstate.temperature, cstate.p[Oil])
00510 : oilpvt.viscosity(pvt_region, cstate.temperature, cstate.p[Oil], cstate.rs);
00511 cstate.mu[Gas] = is_sg
00512 ? gaspvt.saturatedViscosity(pvt_region, cstate.temperature, cstate.p[Gas])
00513 : gaspvt.viscosity(pvt_region, cstate.temperature, cstate.p[Gas], cstate.rv);
00514 cstate.b[Water] = waterpvt.inverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Water]);
00515 cstate.b[Oil] = is_sg
00516 ? oilpvt.saturatedInverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Oil])
00517 : oilpvt.inverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Oil], cstate.rs);
00518 cstate.b[Gas] = is_sg
00519 ? gaspvt.saturatedInverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Gas])
00520 : gaspvt.inverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Gas], cstate.rv);
00521
00522
00523 for (int phase = 0; phase < 3; ++phase) {
00524 cstate.lambda[phase] = cstate.kr[phase] / cstate.mu[phase];
00525 }
00526
00527
00528 cstate.rho[Water] = rhos_(cell, Water) * cstate.b[Water];
00529 cstate.rho[Oil] = (rhos_(cell, Oil) + cstate.rs*rhos_(cell, Gas)) * cstate.b[Oil];
00530 cstate.rho[Gas] = (rhos_(cell, Gas) + cstate.rv*rhos_(cell, Oil)) * cstate.b[Gas];
00531
00532
00533 cstate.rssat = oilpvt.saturatedGasDissolutionFactor(pvt_region, cstate.temperature, cstate.p[Oil]);
00534 cstate.rvsat = gaspvt.saturatedOilVaporizationFactor(pvt_region, cstate.temperature, cstate.p[Gas]);
00535
00536 }
00537
00538
00539
00540
00541 void extractFluxes(const ReservoirState& reservoir_state,
00542 const WellState& well_state)
00543 {
00544
00545 total_flux_ = Eigen::Map<const V>(reservoir_state.faceflux().data(),
00546 reservoir_state.faceflux().size());
00547 total_wellperf_flux_ = Eigen::Map<const V>(well_state.perfRates().data(),
00548 well_state.perfRates().size());
00549 comp_wellperf_flux_ = Eigen::Map<const DataBlock>(well_state.perfPhaseRates().data(),
00550 well_state.perfRates().size(),
00551 numPhases());
00552 const int num_cells = reservoir_state.pressure().size();
00553 total_wellflux_cell_ = superset(total_wellperf_flux_, Base::wellModel().wellOps().well_cells, num_cells);
00554 assert(Base::numPhases() == 3);
00555 V oilflux = comp_wellperf_flux_.col(1);
00556 V gasflux = comp_wellperf_flux_.col(2);
00557 oil_wellflux_cell_ = superset(oilflux, Base::wellModel().wellOps().well_cells, num_cells);
00558 gas_wellflux_cell_ = superset(gasflux, Base::wellModel().wellOps().well_cells, num_cells);
00559 assert(numPhases() * well_state.perfRates().size() == well_state.perfPhaseRates().size());
00560 }
00561
00562
00563
00564
00565
00566 void extractState(const ReservoirState& reservoir_state,
00567 const WellState& well_state)
00568 {
00569 state_.reservoir_state = reservoir_state;
00570 state_.well_state = well_state;
00571 const std::vector<double>& p = reservoir_state.pressure();
00572 state_.tr_mult = Base::transMult(ADB::constant(Eigen::Map<const V>(p.data(), p.size()))).value();
00573 state_.pv_mult = Base::poroMult(ADB::constant(Eigen::Map<const V>(p.data(), p.size()))).value();
00574 }
00575
00576
00577
00578
00579
00580 void computeOrdering()
00581 {
00582 assert(!geo_.nnc().hasNNC());
00583 static_assert(std::is_same<Grid, UnstructuredGrid>::value,
00584 "compute_sequence() is written in C and therefore requires an UnstructuredGrid, "
00585 "it must be rewritten to use other grid classes such as CpGrid");
00586 using namespace Opm::AutoDiffGrid;
00587 const int num_cells = numCells(grid_);
00588 sequence_.resize(num_cells);
00589 components_.resize(num_cells + 1);
00590 int num_components = -1;
00591
00592
00593 using namespace Opm::AutoDiffGrid;
00594 const int num_faces = numFaces(grid_);
00595 V flux_on_all_faces = superset(total_flux_, ops_.internal_faces, num_faces);
00596 compute_sequence(&grid_, flux_on_all_faces.data(), sequence_.data(), components_.data(), &num_components);
00597 OpmLog::debug(std::string("Number of components: ") + std::to_string(num_components));
00598 components_.resize(num_components + 1);
00599 }
00600
00601
00602
00603
00604 void solveComponents()
00605 {
00606
00607 max_abs_dx_[0] = 0.0;
00608 max_abs_dx_[1] = 0.0;
00609 max_abs_dx_cell_[0] = -1;
00610 max_abs_dx_cell_[1] = -1;
00611
00612
00613 const int num_components = components_.size() - 1;
00614 for (int comp = 0; comp < num_components; ++comp) {
00615 const int comp_size = components_[comp + 1] - components_[comp];
00616 if (comp_size == 1) {
00617 solveSingleCell(sequence_[components_[comp]]);
00618 } else {
00619 solveMultiCell(comp_size, &sequence_[components_[comp]]);
00620 }
00621 }
00622
00623
00624 {
00625 std::ostringstream os;
00626 os << "=== Max abs dx[0]: " << max_abs_dx_[0] << " (cell " << max_abs_dx_cell_[0]
00627 <<") dx[1]: " << max_abs_dx_[1] << " (cell " << max_abs_dx_cell_[1] << ")";
00628 OpmLog::debug(os.str());
00629 }
00630 }
00631
00632
00633
00634
00635
00636 void solveSingleCell(const int cell)
00637 {
00638
00639 Vec2 res;
00640 Mat22 jac;
00641 assembleSingleCell(cell, res, jac);
00642
00643
00644 int iter = 0;
00645 const int max_iter = 100;
00646 double relaxation = 1.0;
00647 while (!getConvergence(cell, res) && iter < max_iter) {
00648 Vec2 dx;
00649 jac.solve(dx, res);
00650 dx *= relaxation;
00651
00652 updateState(cell, -dx);
00653
00654 assembleSingleCell(cell, res, jac);
00655 ++iter;
00656 if (iter > 10) {
00657 relaxation = 0.85;
00658 if (iter > 15) {
00659 relaxation = 0.70;
00660 }
00661 if (iter > 20) {
00662 relaxation = 0.55;
00663 }
00664 if (iter > 25) {
00665 relaxation = 0.40;
00666 }
00667 if (iter > 30) {
00668 relaxation = 0.25;
00669 }
00670
00671
00672
00673
00674
00675 }
00676 }
00677 if (iter == max_iter) {
00678 std::ostringstream os;
00679 os << "Failed to converge in cell " << cell << ", residual = " << res
00680 << ", cell values { s = ( " << cstate_[cell].s[Water] << ", " << cstate_[cell].s[Oil] << ", " << cstate_[cell].s[Gas]
00681 << " ), rs = " << cstate_[cell].rs << ", rv = " << cstate_[cell].rv << " }";
00682 OpmLog::debug(os.str());
00683 }
00684 }
00685
00686
00687
00688
00689
00690 void solveMultiCell(const int comp_size, const int* cell_array)
00691 {
00692
00693 for (int ii = 0; ii < comp_size; ++ii) {
00694 solveSingleCell(cell_array[ii]);
00695 }
00696 }
00697
00698
00699
00700
00701 template <typename Scalar>
00702 Scalar oilAccumulation(const CellState<Scalar>& cs)
00703 {
00704 return cs.b[Oil]*cs.s[Oil] + cs.rv*cs.b[Gas]*cs.s[Gas];
00705 }
00706
00707
00708
00709
00710 template <typename Scalar>
00711 Scalar gasAccumulation(const CellState<Scalar>& cs)
00712 {
00713 return cs.b[Gas]*cs.s[Gas] + cs.rs*cs.b[Oil]*cs.s[Oil];
00714 }
00715
00716
00717
00718
00719 void applyThresholdPressure(const int connection, Eval& dp)
00720 {
00721 const double thres_press = Base::threshold_pressures_by_connection_[connection];
00722 if (std::fabs(dp.value()) < thres_press) {
00723 dp.setValue(0.0);
00724 } else {
00725 dp -= dp.value() > 0.0 ? thres_press : -thres_press;
00726 }
00727 }
00728
00729
00730
00731
00732 void assembleSingleCell(const int cell, Vec2& res, Mat22& jac)
00733 {
00734 assert(numPhases() == 3);
00735
00736 CellState<Eval> st;
00737 computeCellState(cell, state_, st);
00738 cstate_[cell] = st.template flatten<double>();
00739
00740
00741 const double pvm0 = state0_.pv_mult[cell];
00742 const double pvm = state_.pv_mult[cell];
00743 const double ao0 = oilAccumulation(cstate0_[cell]) * pvm0;
00744 const Eval ao = oilAccumulation(st) * pvm;
00745 const double ag0 = gasAccumulation(cstate0_[cell]) * pvm0;
00746 const Eval ag = gasAccumulation(st) * pvm;
00747
00748
00749 Eval div_oilflux = Eval::createConstant(0.0);
00750 Eval div_gasflux = Eval::createConstant(0.0);
00751 for (auto conn : graph_.cellConnections(cell)) {
00752 auto conn_cells = graph_.connectionCells(conn.index);
00753 const int from = conn_cells[0];
00754 const int to = conn_cells[1];
00755 if (from < 0 || to < 0) {
00756 continue;
00757 }
00758 assert((from == cell) == (conn.sign > 0.0));
00759 const int other = from == cell ? to : from;
00760 const double vt = conn.sign * total_flux_[conn.index];
00761 const double gdz = conn.sign * gdz_[conn.index];
00762
00763
00764
00765
00766
00767
00768 Eval dh[3];
00769 Eval dh_sat[3];
00770 const Eval grad_oil_press = cstate_[other].p[Oil] - st.p[Oil];
00771 for (int phase : { Water, Oil, Gas }) {
00772 const Eval gradp = cstate_[other].p[phase] - st.p[phase];
00773 const Eval rhoavg = 0.5 * (st.rho[phase] + cstate_[other].rho[phase]);
00774 dh[phase] = gradp - rhoavg * gdz;
00775 if (Base::use_threshold_pressure_) {
00776 applyThresholdPressure(conn.index, dh[phase]);
00777 }
00778 dh_sat[phase] = grad_oil_press - dh[phase];
00779 }
00780 const double tran = trans_all_[conn.index];
00781 const auto& m1 = st.lambda;
00782 const auto& m2 = cstate_[other].lambda;
00783 const auto upw = connectionMultiPhaseUpwind({{ dh_sat[Water].value(), dh_sat[Oil].value(), dh_sat[Gas].value() }},
00784 {{ m1[Water].value(), m1[Oil].value(), m1[Gas].value() }},
00785 {{ m2[Water], m2[Oil], m2[Gas] }},
00786 tran, vt);
00787
00788
00789
00790 Eval b[3];
00791 Eval mob[3];
00792 Eval tot_mob = Eval::createConstant(0.0);
00793 for (int phase : { Water, Oil, Gas }) {
00794 b[phase] = upw[phase] > 0.0 ? st.b[phase] : cstate_[other].b[phase];
00795 mob[phase] = upw[phase] > 0.0 ? m1[phase] : m2[phase];
00796 tot_mob += mob[phase];
00797 }
00798 Eval rs = upw[Oil] > 0.0 ? st.rs : cstate_[other].rs;
00799 Eval rv = upw[Gas] > 0.0 ? st.rv : cstate_[other].rv;
00800
00801 Eval flux[3];
00802 for (int phase : { Oil, Gas }) {
00803 Eval gflux = Eval::createConstant(0.0);
00804 for (int other_phase : { Water, Oil, Gas }) {
00805 if (phase != other_phase) {
00806 gflux += mob[other_phase] * (dh_sat[phase] - dh_sat[other_phase]);
00807 }
00808 }
00809 flux[phase] = b[phase] * (mob[phase] / tot_mob) * (vt + tran*gflux);
00810 }
00811 div_oilflux += flux[Oil] + rv*flux[Gas];
00812 div_gasflux += flux[Gas] + rs*flux[Oil];
00813 }
00814
00815
00816 if (total_wellflux_cell_[cell] > 0.0) {
00817
00818 assert(oil_wellflux_cell_[cell] >= 0.0);
00819 assert(gas_wellflux_cell_[cell] >= 0.0);
00820 div_oilflux -= oil_wellflux_cell_[cell];
00821 div_gasflux -= gas_wellflux_cell_[cell];
00822 } else if (total_wellflux_cell_[cell] < 0.0) {
00823
00824 Eval totmob = st.lambda[Water] + st.lambda[Oil] + st.lambda[Gas];
00825 Eval oilflux = st.b[Oil] * (st.lambda[Oil]/totmob) * total_wellflux_cell_[cell];
00826 Eval gasflux = st.b[Gas] * (st.lambda[Gas]/totmob) * total_wellflux_cell_[cell];
00827 div_oilflux -= (oilflux + st.rv * gasflux);
00828 div_gasflux -= (gasflux + st.rs * oilflux);
00829 }
00830
00831 const Eval oileq = Base::pvdt_[cell]*(ao - ao0) + div_oilflux;
00832 const Eval gaseq = Base::pvdt_[cell]*(ag - ag0) + div_gasflux;
00833
00834 res[0] = oileq.value();
00835 res[1] = gaseq.value();
00836 jac[0][0] = oileq.derivative(0);
00837 jac[0][1] = oileq.derivative(1);
00838 jac[1][0] = gaseq.derivative(0);
00839 jac[1][1] = gaseq.derivative(1);
00840 }
00841
00842
00843
00844
00845
00846 bool getConvergence(const int cell, const Vec2& res)
00847 {
00848 const double tol = 1e-7;
00849
00850 double sres[] = { res[0] / (cstate_[cell].b[Oil] * Base::pvdt_[cell]),
00851 res[1] / (cstate_[cell].b[Gas] * Base::pvdt_[cell]) };
00852 return std::fabs(sres[0]) < tol && std::fabs(sres[1]) < tol;
00853 }
00854
00855
00856
00857
00858 void updateState(const int cell,
00859 const Vec2& dx)
00860 {
00861 if (std::fabs(dx[0]) > max_abs_dx_[0]) {
00862 max_abs_dx_cell_[0] = cell;
00863 }
00864 if (std::fabs(dx[1]) > max_abs_dx_[1]) {
00865 max_abs_dx_cell_[1] = cell;
00866 }
00867 max_abs_dx_[0] = std::max(max_abs_dx_[0], std::fabs(dx[0]));
00868 max_abs_dx_[1] = std::max(max_abs_dx_[1], std::fabs(dx[1]));
00869
00870
00871 const double dsw = dx[0];
00872 double dsg = 0.0;
00873 auto& hcstate = state_.reservoir_state.hydroCarbonState()[cell];
00874 if (hcstate == HydroCarbonState::GasAndOil) {
00875 dsg = dx[1];
00876 } else if (hcstate == HydroCarbonState::GasOnly) {
00877 dsg = -dsw;
00878 }
00879 const double dso = -(dsw + dsg);
00880
00881
00882 const double maxval = std::max(std::fabs(dsw), std::max(std::fabs(dso), std::fabs(dsg)));
00883 const double sfactor = std::min(1.0, Base::dsMax() / maxval);
00884 double* s = state_.reservoir_state.saturation().data() + 3*cell;
00885 s[Water] += sfactor*dsw;
00886 s[Gas] += sfactor*dsg;
00887 s[Oil] = 1.0 - s[Water] - s[Gas];
00888
00889
00890 for (int phase : { Gas, Oil, Water }) {
00891 if (s[phase] < 0.0) {
00892 for (int other_phase : { Water, Oil, Gas }) {
00893 if (phase != other_phase) {
00894 s[other_phase] /= (1.0 - s[phase]);
00895 }
00896 }
00897 s[phase] = 0.0;
00898 }
00899 }
00900
00901
00902 double& rs = state_.reservoir_state.gasoilratio()[cell];
00903 const double rs_old = rs;
00904 if (hcstate == HydroCarbonState::OilOnly) {
00905
00906 const double drs = dx[1];
00907
00908
00909 rs += drs;
00910 rs = std::max(rs, 0.0);
00911 }
00912
00913
00914 double& rv = state_.reservoir_state.rv()[cell];
00915 const double rv_old = rv;
00916 if (hcstate == HydroCarbonState::GasOnly) {
00917
00918 const double drv = dx[1];
00919
00920
00921 rv += drv;
00922 rv = std::max(rv, 0.0);
00923 }
00924
00925 const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
00926 const bool water_only = s[Water] > (1 - epsilon);
00927 const auto old_hcstate = hcstate;
00928 hcstate = HydroCarbonState::GasAndOil;
00929
00930 {
00931 const double rssat_old = cstate_[cell].rssat;
00932 const double rssat = rssat_old;
00933 const bool is_rs = old_hcstate == HydroCarbonState::OilOnly;
00934 const bool has_gas = (s[Gas] > 0.0 && !is_rs);
00935 const bool gas_vaporized = ( (rs > rssat * (1+epsilon) && is_rs ) && (rs_old > rssat_old * (1-epsilon)) );
00936 if (water_only || has_gas || gas_vaporized) {
00937 rs = rssat;
00938 } else {
00939 hcstate = HydroCarbonState::OilOnly;
00940 }
00941 }
00942
00943
00944 {
00945 const double rvsat_old = cstate_[cell].rvsat;
00946 const double rvsat = rvsat_old;
00947 const bool is_rv = old_hcstate == HydroCarbonState::GasOnly;
00948 const bool has_oil = (s[Oil] > 0.0 && !is_rv);
00949 const bool oil_condensed = ( (rv > rvsat * (1+epsilon) && is_rv) && (rv_old > rvsat_old * (1-epsilon)) );
00950 if (water_only || has_oil || oil_condensed) {
00951 rv = rvsat;
00952 } else {
00953 hcstate = HydroCarbonState::GasOnly;
00954 }
00955 }
00956 }
00957 };
00958
00959
00960
00961
00962
00963
00964
00965
00967 template <class Grid, class WellModel>
00968 struct ModelTraits< BlackoilReorderingTransportModel<Grid, WellModel> >
00969 {
00970 typedef BlackoilState ReservoirState;
00971 typedef WellStateFullyImplicitBlackoil WellState;
00972 typedef BlackoilModelParameters ModelParameters;
00973 typedef DefaultBlackoilSolutionState SolutionState;
00974 };
00975
00976 }
00977
00978
00979
00980
00981 #endif // OPM_BLACKOILREORDERINGTRANSPORTMODEL_HEADER_INCLUDED