00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED
00022 #define OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED
00023
00024 #include <opm/autodiff/BlackoilModelBase.hpp>
00025 #include <opm/core/simulator/BlackoilState.hpp>
00026 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
00027 #include <opm/autodiff/BlackoilModelParameters.hpp>
00028 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
00029 #include <opm/autodiff/multiPhaseUpwind.hpp>
00030
00031 namespace Opm {
00032
00034 template<class Grid, class WellModel>
00035 class BlackoilTransportModel : public BlackoilModelBase<Grid, WellModel, BlackoilTransportModel<Grid, WellModel> >
00036 {
00037 public:
00038 typedef BlackoilModelBase<Grid, WellModel, BlackoilTransportModel<Grid, WellModel> > Base;
00039 friend Base;
00040
00041 typedef typename Base::ReservoirState ReservoirState;
00042 typedef typename Base::WellState WellState;
00043 typedef typename Base::SolutionState SolutionState;
00044 typedef typename Base::V V;
00045
00046
00061 BlackoilTransportModel(const typename Base::ModelParameters& param,
00062 const Grid& grid,
00063 const BlackoilPropsAdFromDeck& fluid,
00064 const DerivedGeology& geo,
00065 const RockCompressibility* rock_comp_props,
00066 const StandardWells& std_wells,
00067 const NewtonIterationBlackoilInterface& linsolver,
00068 std::shared_ptr<const EclipseState> eclState,
00069 const bool has_disgas,
00070 const bool has_vapoil,
00071 const bool terminal_output)
00072 : Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
00073 eclState, has_disgas, has_vapoil, terminal_output)
00074 {
00075 }
00076
00077
00078
00079
00080
00081 void prepareStep(const SimulatorTimerInterface& timer,
00082 const ReservoirState& reservoir_state,
00083 const WellState& well_state)
00084 {
00085 Base::prepareStep(timer, reservoir_state, well_state);
00086 Base::param_.solve_welleq_initially_ = false;
00087 SolutionState state0 = variableState(reservoir_state, well_state);
00088 asImpl().makeConstantState(state0);
00089 asImpl().computeAccum(state0, 0);
00090 }
00091
00092
00093
00094
00095
00096 SimulatorReport
00097 assemble(const ReservoirState& reservoir_state,
00098 WellState& well_state,
00099 const bool initial_assembly)
00100 {
00101 using namespace Opm::AutoDiffGrid;
00102
00103 SimulatorReport report;
00104
00105
00106
00107
00108 if (isVFPActive()) {
00109 SolutionState state = asImpl().variableState(reservoir_state, well_state);
00110 SolutionState state0 = state;
00111 asImpl().makeConstantState(state0);
00112 asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
00113 }
00114
00115
00116
00117 asImpl().wellModel().updateWellControls(well_state);
00118
00119 if (initial_assembly) {
00120
00121 const_cast<V&>(total_flux_)
00122 = Eigen::Map<const V>(reservoir_state.faceflux().data(), reservoir_state.faceflux().size());
00123 const_cast<V&>(total_wellperf_flux_)
00124 = Eigen::Map<const V>(well_state.perfRates().data(), well_state.perfRates().size());
00125 const_cast<DataBlock&>(comp_wellperf_flux_)
00126 = Eigen::Map<const DataBlock>(well_state.perfPhaseRates().data(), well_state.perfRates().size(), numPhases());
00127 assert(numPhases() * well_state.perfRates().size() == well_state.perfPhaseRates().size());
00128 asImpl().updatePhaseCondFromPrimalVariable(reservoir_state);
00129 }
00130
00131
00132 SolutionState state = asImpl().variableState(reservoir_state, well_state);
00133
00134 if (initial_assembly) {
00135 SolutionState state0 = state;
00136 asImpl().makeConstantState(state0);
00137 asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
00138 }
00139
00140
00141
00142 asImpl().assembleMassBalanceEq(state);
00143
00144
00145 if ( ! wellsActive() ) {
00146 return report;
00147 }
00148
00149 std::vector<ADB> mob_perfcells;
00150 std::vector<ADB> b_perfcells;
00151 asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
00152 if (param_.solve_welleq_initially_ && initial_assembly) {
00153
00154 report += asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
00155 }
00156 V aliveWells;
00157 std::vector<ADB> cq_s;
00158
00159
00160
00161 asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
00162
00163 asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
00164 asImpl().wellModel().addWellFluxEq(cq_s, state, residual_);
00165 asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
00166 asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
00167
00168 return report;
00169 }
00170
00171
00172
00173
00176 V solveJacobianSystem() const
00177 {
00178 const int n_transport = residual_.material_balance_eq[1].size();
00179 const int n_full = residual_.sizeNonLinear();
00180 const auto& mb = residual_.material_balance_eq;
00181 LinearisedBlackoilResidual transport_res = {
00182 {
00183
00184 ADB::function(mb[1].value(), { mb[1].derivative()[1], mb[1].derivative()[2] }),
00185 ADB::function(mb[2].value(), { mb[2].derivative()[1], mb[2].derivative()[2] })
00186 },
00187 ADB::null(),
00188 ADB::null(),
00189 residual_.matbalscale,
00190 residual_.singlePrecision
00191 };
00192 assert(transport_res.sizeNonLinear() == 2*n_transport);
00193 V dx_transport = linsolver_.computeNewtonIncrement(transport_res);
00194 assert(dx_transport.size() == 2*n_transport);
00195 V dx_full = V::Zero(n_full);
00196 for (int i = 0; i < 2*n_transport; ++i) {
00197 dx_full(n_transport + i) = dx_transport(i);
00198 }
00199 return dx_full;
00200 }
00201
00202
00203
00204
00205
00206 using Base::numPhases;
00207 using Base::numMaterials;
00208
00209 protected:
00210 using Base::asImpl;
00211 using Base::materialName;
00212 using Base::convergenceReduction;
00213 using Base::maxResidualAllowed;
00214
00215 using Base::linsolver_;
00216 using Base::residual_;
00217 using Base::sd_;
00218 using Base::geo_;
00219 using Base::ops_;
00220 using Base::grid_;
00221 using Base::use_threshold_pressure_;
00222 using Base::canph_;
00223 using Base::active_;
00224 using Base::pvdt_;
00225 using Base::fluid_;
00226 using Base::param_;
00227 using Base::terminal_output_;
00228
00229 using Base::isVFPActive;
00230 using Base::phaseCondition;
00231 using Base::vfp_properties_;
00232 using Base::wellsActive;
00233
00234 V total_flux_;
00235 V total_wellperf_flux_;
00236 DataBlock comp_wellperf_flux_;
00237 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind_flags_;
00238
00239
00240
00241
00242
00243 SolutionState
00244 variableState(const ReservoirState& x,
00245 const WellState& xw) const
00246 {
00247
00248 std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
00249 std::vector<ADB> vars = ADB::variables(vars0);
00250 const std::vector<int> indices = asImpl().variableStateIndices();
00251 vars[indices[Pressure]] = ADB::constant(vars[indices[Pressure]].value());
00252 vars[indices[Qs]] = ADB::constant(vars[indices[Qs]].value());
00253 vars[indices[Bhp]] = ADB::constant(vars[indices[Bhp]].value());
00254 return asImpl().variableStateExtractVars(x, indices, vars);
00255 }
00256
00257
00258
00259
00260
00261
00262
00263 void assembleMassBalanceEq(const SolutionState& state)
00264 {
00265
00266
00267
00268
00269
00270
00271 asImpl().computeAccum(state, 1);
00272
00273
00274
00275 const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
00276 const V trans_nnc = ops_.nnc_trans;
00277 V trans_all(transi.size() + trans_nnc.size());
00278 trans_all << transi, trans_nnc;
00279 const ADB tr_mult = asImpl().transMult(state.pressure);
00280 const V gdz = geo_.gravity()[2] * (ops_.grad * geo_.z().matrix());
00281
00282
00283 const std::vector<PhasePresence>& cond = asImpl().phaseCondition();
00284 const std::vector<ADB> kr = asImpl().computeRelPerm(state);
00285 #pragma omp parallel for schedule(static)
00286 for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
00287
00288 const int canonical_phase_idx = canph_[ phase_idx ];
00289 const ADB& phase_pressure = state.canonical_phase_pressures[canonical_phase_idx];
00290 sd_.rq[phase_idx].mu = asImpl().fluidViscosity(canonical_phase_idx, phase_pressure, state.temperature, state.rs, state.rv, cond);
00291 sd_.rq[phase_idx].kr = kr[canonical_phase_idx];
00292
00293
00294 sd_.rq[ phase_idx ].mob = tr_mult * sd_.rq[phase_idx].kr / sd_.rq[phase_idx].mu;
00295
00296
00297 sd_.rq[phase_idx].rho = asImpl().fluidDensity(canonical_phase_idx, sd_.rq[phase_idx].b, state.rs, state.rv);
00298 const ADB rhoavg = ops_.caver * sd_.rq[phase_idx].rho;
00299 sd_.rq[ phase_idx ].dh = ops_.grad * phase_pressure - rhoavg * gdz;
00300
00301 if (use_threshold_pressure_) {
00302 asImpl().applyThresholdPressures(sd_.rq[ phase_idx ].dh);
00303 }
00304 }
00305
00306
00307 const ADB gradp = ops_.grad * state.pressure;
00308 std::vector<ADB> dh_sat(numPhases(), ADB::null());
00309 for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
00310 dh_sat[phase_idx] = gradp - sd_.rq[phase_idx].dh;
00311 }
00312
00313
00314 upwind_flags_ = multiPhaseUpwind(dh_sat, trans_all);
00315
00316
00317
00318 std::vector<ADB> mob(numPhases(), ADB::null());
00319 std::vector<ADB> b(numPhases(), ADB::null());
00320 ADB rs = ADB::null();
00321 ADB rv = ADB::null();
00322 ADB tot_mob = ADB::constant(V::Zero(gdz.size()));
00323 for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
00324 UpwindSelector<double> upwind(grid_, ops_, upwind_flags_.col(phase_idx));
00325 mob[phase_idx] = upwind.select(sd_.rq[phase_idx].mob);
00326 tot_mob += mob[phase_idx];
00327 b[phase_idx] = upwind.select(sd_.rq[phase_idx].b);
00328 if (canph_[phase_idx] == Oil) {
00329 rs = upwind.select(state.rs);
00330 }
00331 if (canph_[phase_idx] == Gas) {
00332 rv = upwind.select(state.rv);
00333 }
00334 }
00335
00336
00337 for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
00338 ADB gflux = ADB::constant(V::Zero(gdz.size()));
00339 for (int other_phase = 0; other_phase < numPhases(); ++other_phase) {
00340 if (phase_idx != other_phase) {
00341 gflux += mob[other_phase] * (dh_sat[phase_idx] - dh_sat[other_phase]);
00342 }
00343 }
00344 sd_.rq[phase_idx].mflux = b[phase_idx] * (mob[phase_idx] / tot_mob) * (total_flux_ + trans_all * gflux);
00345 }
00346
00347 #pragma omp parallel for schedule(static)
00348 for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
00349
00350
00351
00352
00353
00354 residual_.material_balance_eq[ phase_idx ] =
00355 pvdt_ * (sd_.rq[phase_idx].accum[1] - sd_.rq[phase_idx].accum[0])
00356 + ops_.div*sd_.rq[phase_idx].mflux;
00357 }
00358
00359
00360
00361
00362
00363
00364 if (active_[ Oil ] && active_[ Gas ]) {
00365 const int po = fluid_.phaseUsage().phase_pos[ Oil ];
00366 const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
00367 residual_.material_balance_eq[ pg ] += ops_.div * (rs * sd_.rq[po].mflux);
00368 residual_.material_balance_eq[ po ] += ops_.div * (rv * sd_.rq[pg].mflux);
00369 }
00370
00371 if (param_.update_equations_scaling_) {
00372 asImpl().updateEquationsScaling();
00373 }
00374
00375 }
00376
00377
00378
00379
00380
00381 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> multiPhaseUpwind(const std::vector<ADB>& head_diff,
00382 const V& transmissibility)
00383 {
00384 assert(numPhases() == 3);
00385 const int num_connections = head_diff[0].size();
00386 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind(num_connections, numPhases());
00387 for (int conn = 0; conn < num_connections; ++conn) {
00388 const double q = total_flux_[conn];
00389 const double t = transmissibility[conn];
00390 const int a = ops_.connection_cells(conn, 0);
00391 const int b = ops_.connection_cells(conn, 1);
00392 auto up = connectionMultiPhaseUpwind(
00393 {{ head_diff[0].value()[conn], head_diff[1].value()[conn], head_diff[2].value()[conn] }},
00394 {{ sd_.rq[0].mob.value()[a], sd_.rq[1].mob.value()[a], sd_.rq[2].mob.value()[a]}},
00395 {{ sd_.rq[0].mob.value()[b], sd_.rq[1].mob.value()[b], sd_.rq[2].mob.value()[b]}},
00396 t,
00397 q);
00398 for (int ii = 0; ii < numPhases(); ++ii) {
00399 upwind(conn, ii) = up[ii];
00400 }
00401 }
00402 return upwind;
00403 }
00404
00405
00406
00407
00408
00409 void computeWellFlux(const SolutionState& state,
00410 const std::vector<ADB>& mob_perfcells,
00411 const std::vector<ADB>& b_perfcells,
00412 V& ,
00413 std::vector<ADB>& cq_s) const
00414 {
00415
00416
00417 if( ! asImpl().localWellsActive() ) return ;
00418
00419 const int np = asImpl().wells().number_of_phases;
00420 const int nw = asImpl().wells().number_of_wells;
00421 const int nperf = asImpl().wells().well_connpos[nw];
00422 const Opm::PhaseUsage& pu = asImpl().fluid_.phaseUsage();
00423
00424
00425 ADB totmob_perfcells = ADB::constant(V::Zero(nperf));
00426 for (int phase = 0; phase < numPhases(); ++phase) {
00427 totmob_perfcells += mob_perfcells[phase];
00428 }
00429
00430
00431 std::vector<ADB> frac_flow(np, ADB::null());
00432 for (int phase = 0; phase < np; ++phase) {
00433 frac_flow[phase] = mob_perfcells[phase] / totmob_perfcells;
00434 }
00435
00436
00437 V is_inj = V::Zero(nperf);
00438 V is_prod = V::Zero(nperf);
00439 for (int c = 0; c < nperf; ++c){
00440 if (total_wellperf_flux_[c] > 0.0) {
00441 is_inj[c] = 1;
00442 } else {
00443 is_prod[c] = 1;
00444 }
00445 }
00446
00447
00448 std::vector<ADB> cq_s_prod(3, ADB::null());
00449 for (int phase = 0; phase < np; ++phase) {
00450
00451 cq_s_prod[phase] = b_perfcells[phase] * frac_flow[phase] * total_wellperf_flux_;
00452 }
00453 if (asImpl().has_disgas_ || asImpl().has_vapoil_) {
00454 const int oilpos = pu.phase_pos[Oil];
00455 const int gaspos = pu.phase_pos[Gas];
00456 const ADB cq_s_prod_oil = cq_s_prod[oilpos];
00457 const ADB cq_s_prod_gas = cq_s_prod[gaspos];
00458 cq_s_prod[gaspos] += subset(state.rs, Base::well_model_.wellOps().well_cells) * cq_s_prod_oil;
00459 cq_s_prod[oilpos] += subset(state.rv, Base::well_model_.wellOps().well_cells) * cq_s_prod_gas;
00460 }
00461
00462
00463 cq_s.resize(np, ADB::null());
00464 for (int phase = 0; phase < np; ++phase) {
00465 const int pos = pu.phase_pos[phase];
00466
00467 const V cq_s_inj = comp_wellperf_flux_.col(pos);
00468 cq_s[phase] = is_prod * cq_s_prod[phase] + is_inj * cq_s_inj;
00469 }
00470 }
00471
00472
00473
00474
00475
00476 bool getConvergence(const SimulatorTimerInterface& timer, const int iteration)
00477 {
00478 const double dt = timer.currentStepLength();
00479 const double tol_mb = param_.tolerance_mb_;
00480 const double tol_cnv = param_.tolerance_cnv_;
00481
00482 const int nc = Opm::AutoDiffGrid::numCells(grid_);
00483 const int np = asImpl().numPhases();
00484 const int nm = asImpl().numMaterials();
00485 assert(int(sd_.rq.size()) == nm);
00486
00487 const V& pv = geo_.poreVolume();
00488
00489 std::vector<double> R_sum(nm);
00490 std::vector<double> B_avg(nm);
00491 std::vector<double> maxCoeff(nm);
00492 std::vector<double> maxNormWell(np);
00493 Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
00494 Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
00495 Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
00496
00497 for ( int idx = 0; idx < nm; ++idx )
00498 {
00499 const ADB& tempB = sd_.rq[idx].b;
00500 B.col(idx) = 1./tempB.value();
00501 R.col(idx) = residual_.material_balance_eq[idx].value();
00502 tempV.col(idx) = R.col(idx).abs()/pv;
00503 }
00504
00505 const double pvSum = convergenceReduction(B, tempV, R,
00506 R_sum, maxCoeff, B_avg, maxNormWell,
00507 nc);
00508
00509 std::vector<double> CNV(nm);
00510 std::vector<double> mass_balance_residual(nm);
00511 std::vector<double> well_flux_residual(np);
00512
00513 bool converged_MB = true;
00514 bool converged_CNV = true;
00515
00516 for ( int idx = 1; idx < nm; ++idx ) {
00517 CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
00518 mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
00519 converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
00520 converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
00521 assert(nm >= np);
00522 }
00523
00524 const bool converged = converged_MB && converged_CNV;
00525
00526 for (int idx = 0; idx < nm; ++idx) {
00527 if (std::isnan(mass_balance_residual[idx])
00528 || std::isnan(CNV[idx])
00529 || (idx < np && std::isnan(well_flux_residual[idx]))) {
00530 OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx));
00531 }
00532 if (mass_balance_residual[idx] > maxResidualAllowed()
00533 || CNV[idx] > maxResidualAllowed()
00534 || (idx < np && well_flux_residual[idx] > maxResidualAllowed())) {
00535 OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx));
00536 }
00537 }
00538
00539 if ( terminal_output_ ) {
00540
00541 std::ostringstream os;
00542 if (iteration == 0) {
00543 os << "\nIter";
00544 for (int idx = 1; idx < nm; ++idx) {
00545 os << " MB(" << materialName(idx).substr(0, 3) << ") ";
00546 }
00547 for (int idx = 1; idx < nm; ++idx) {
00548 os << " CNV(" << materialName(idx).substr(0, 1) << ") ";
00549 }
00550 os << '\n';
00551 }
00552 os.precision(3);
00553 os.setf(std::ios::scientific);
00554 os << std::setw(4) << iteration;
00555 for (int idx = 1; idx < nm; ++idx) {
00556 os << std::setw(11) << mass_balance_residual[idx];
00557 }
00558 for (int idx = 1; idx < nm; ++idx) {
00559 os << std::setw(11) << CNV[idx];
00560 }
00561 OpmLog::info(os.str());
00562 }
00563 return converged;
00564 }
00565 };
00566
00567
00569 template <class Grid, class WellModel>
00570 struct ModelTraits< BlackoilTransportModel<Grid, WellModel> >
00571 {
00572 typedef BlackoilState ReservoirState;
00573 typedef WellStateFullyImplicitBlackoil WellState;
00574 typedef BlackoilModelParameters ModelParameters;
00575 typedef DefaultBlackoilSolutionState SolutionState;
00576 };
00577
00578 }
00579
00580
00581
00582
00583 #endif // OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED