diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index 3a6d24ec5..e3cc97f23 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -218,7 +218,10 @@ try rock_comp->isActive() ? rock_comp.get() : 0, wells, *fis_solver, - grav); + grav, + deck->hasKeyword("DISGAS"), + deck->hasKeyword("VAPOIL") ); + SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); ++simtimer; diff --git a/examples/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp index 842e981dd..44420f8f1 100644 --- a/examples/test_implicit_ad.cpp +++ b/examples/test_implicit_ad.cpp @@ -104,7 +104,7 @@ try Opm::NewtonIterationBlackoilSimple fis_solver(param); - Opm::FullyImplicitBlackoilSolver solver(param, *g, props, geo, 0, *wells, fis_solver); + Opm::FullyImplicitBlackoilSolver solver(param, *g, props, geo, 0, *wells, fis_solver, /*hasDisgas*/ true, /*hasVapoil=*/false); Opm::BlackoilState state; initStateBasic(*g, props0, param, 0.0, state); diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp index 17b382daa..ae5c9e186 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -94,6 +94,8 @@ namespace Opm phase_usage_ = phaseUsageFromDeck(deck); + + // Surface densities. Accounting for different orders in eclipse and our code. Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY"); int numRegions = densityKeyword->size(); diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp index 6677b0bdb..bced83fa9 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.hpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp @@ -360,6 +360,8 @@ namespace Opm std::unique_ptr satprops_; PhaseUsage phase_usage_; + bool has_vapoil_; + bool has_disgas_; // The PVT region which is to be used for each cell std::vector cellPvtRegionIdx_; diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index a85f719dc..381204b64 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -70,7 +70,9 @@ namespace Opm { const DerivedGeology& geo , const RockCompressibility* rock_comp_props, const Wells& wells, - const NewtonIterationBlackoilInterface& linsolver); + const NewtonIterationBlackoilInterface& linsolver, + const bool has_disgas, + const bool has_vapoil ); /// Take a single forward step, modifiying /// state.pressure() @@ -127,7 +129,7 @@ namespace Opm { // the Newton relaxation type enum RelaxType { DAMPEN, SOR }; - + enum PrimalVariables { Sg = 0, RS = 1, RV = 2 }; // Member data const Grid& grid_; @@ -144,6 +146,8 @@ namespace Opm { HelperOps ops_; const WellOps wops_; const M grav_; + const bool has_disgas_; + const bool has_vapoil_; double dp_max_rel_; double ds_max_; double drs_max_rel_; @@ -159,6 +163,8 @@ namespace Opm { LinearisedBlackoilResidual residual_; + std::vector primalVariable_; + // Private methods. SolutionState constantState(const BlackoilState& x, @@ -279,6 +285,12 @@ namespace Opm { void classifyCondition(const BlackoilState& state); + + /// update the primal variable for Sg, Rv or Rs. The Gas phase must + /// be active to call this method. + void + updatePrimalVariableFromState(const BlackoilState& state); + /// Compute convergence based on total mass balance (tol_mb) and maximum /// residual mass balance (tol_cnv). bool getConvergence(const double dt); diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index c4a041dc1..4ac2cb7e2 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -219,7 +219,9 @@ namespace { const DerivedGeology& geo , const RockCompressibility* rock_comp_props, const Wells& wells, - const NewtonIterationBlackoilInterface& linsolver) + const NewtonIterationBlackoilInterface& linsolver, + const bool has_disgas, + const bool has_vapoil) : grid_ (grid) , fluid_ (fluid) , geo_ (geo) @@ -232,6 +234,8 @@ namespace { , ops_ (grid) , wops_ (wells) , grav_ (gravityOperator(grid_, ops_, geo_)) + , has_disgas_(has_disgas) + , has_vapoil_(has_vapoil) , dp_max_rel_ (1.0e9) , ds_max_ (0.2) , drs_max_rel_ (1.0e9) @@ -272,7 +276,8 @@ namespace { { const V pvdt = geo_.poreVolume() / dt; - classifyCondition(x); + if (active_[Gas]) { updatePrimalVariableFromState(x); } + { const SolutionState state = constantState(x, xw); computeAccum(state, 0); @@ -463,26 +468,13 @@ namespace { V isRs = V::Zero(nc,1); V isRv = V::Zero(nc,1); V isSg = V::Zero(nc,1); - bool disgas = false; - bool vapoil = false; if (active_[ Gas ]){ - // this is a temporary hack to find if vapoil or disgas - // is a active component. Should be given directly from - // DISGAS and VAPOIL keywords in the deck. - for (int c = 0; c < nc; c++){ - if(x.rv()[c] > 0) - vapoil = true; - if(x.gasoilratio ()[c] > 0) - disgas = true; - } - for (int c = 0; c < nc ; c++ ) { - const PhasePresence cond = phaseCondition()[c]; - if ( (!cond.hasFreeGas()) && disgas ) { + if ( primalVariable_[c] == PrimalVariables::RS ) { isRs[c] = 1; } - else if ( (!cond.hasFreeOil()) && vapoil ) { + else if ( primalVariable_[c] == PrimalVariables::RV ) { isRv[c] = 1; } else { @@ -546,12 +538,12 @@ namespace { state.saturation[ pu.phase_pos[ Gas ] ] = sg; so = so - sg; - if (disgas) { + if (has_disgas_) { state.rs = (1-isRs) * rsSat + isRs*xvar; } else { state.rs = rsSat; } - if (vapoil) { + if (has_vapoil_) { state.rv = (1-isRv) * rvSat + isRv*xvar; } else { state.rv = rvSat; @@ -1232,31 +1224,17 @@ namespace { V isRs = V::Zero(nc,1); V isRv = V::Zero(nc,1); V isSg = V::Zero(nc,1); - - bool disgas = false; - bool vapoil = false; - - // this is a temporary hack to find if vapoil or disgas - // is a active component. Should be given directly from - // DISGAS and VAPOIL keywords in the deck. - for (int c = 0; c0) - vapoil = true; - if(state.gasoilratio()[c]>0) - disgas = true; - } - - const std::vector conditions = phaseCondition(); - for (int c = 0; c < nc; c++ ) { - const PhasePresence cond = conditions[c]; - if ( (!cond.hasFreeGas()) && disgas ) { - isRs[c] = 1; - } - else if ( (!cond.hasFreeOil()) && vapoil ) { - isRv[c] = 1; - } - else { - isSg[c] = 1; + if (active_[Gas]) { + for (int c = 0; c < nc ; c++ ) { + if ( primalVariable_[c] == PrimalVariables::RS ) { + isRs[c] = 1; + } + else if ( primalVariable_[c] == PrimalVariables::RV ) { + isRv[c] = 1; + } + else { + isSg[c] = 1; + } } } @@ -1330,14 +1308,14 @@ namespace { const double drsmaxrel = drsMaxRel(); const double drvmax = 1e9;//% same as in Mrst V rs; - if (disgas) { + if (has_disgas_) { const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); const V drs = isRs * dxvar; const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel); rs = rs_old - drs_limited; } V rv; - if (vapoil) { + if (has_vapoil_) { const V rv_old = Eigen::Map(&state.rv()[0], nc); const V drv = isRv * dxvar; const V drv_limited = sign(drv) * drv.abs().min(drvmax); @@ -1356,23 +1334,23 @@ namespace { // reset the phase conditions std::vector cond(nc); - if (disgas) { - // The obvioious case - auto ix0 = (sg > 0 && isRs == 0); + std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg); + + if (has_disgas_) { + // The obvious case + auto hasGas = (sg > 0 && isRs == 0); // keep oil saturated if previous sg is sufficient large: const int pos = pu.phase_pos[ Gas ]; - auto ix1 = (sg < 0 && s_old.col(pos) > epsilon); + auto hadGas = (sg < 0 && s_old.col(pos) > epsilon); // Set oil saturated if previous rs is sufficiently large const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); - auto ix2 = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) ); + auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) ); - auto gasPresent = watOnly || ix0 || ix1 || ix2; + auto useSg = watOnly || hasGas || hadGas || gasVaporized; for (int c = 0; c < nc; ++c) { - if (gasPresent[c]) { - rs[c] = rsSat[c]; - cond[c].setFreeGas(); - } + if (useSg[c]) { rs[c] = rsSat[c];} + else { primalVariable_[c] = PrimalVariables::RS; } } } @@ -1381,26 +1359,24 @@ namespace { const V rvSat0 = fluidRvSat(p_old, cells_); const V rvSat = fluidRvSat(p, cells_); - if (vapoil) { + if (has_vapoil_) { // The obvious case - auto ix0 = (so > 0 && isRv == 0); + auto hasOil = (so > 0 && isRv == 0); - // keep oil saturated if previous sg is sufficient large: + // keep oil saturated if previous so is sufficient large: const int pos = pu.phase_pos[ Oil ]; - auto ix1 = (so < 0 && s_old.col(pos) > epsilon ); - // Set oil saturated if previous rs is sufficiently large + auto hadOil = (so < 0 && s_old.col(pos) > epsilon ); + // Set oil saturated if previous rv is sufficiently large const V rv_old = Eigen::Map(&state.rv()[0], nc); - auto ix2 = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) ); - auto oilPresent = watOnly || ix0 || ix1 || ix2; + auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) ); + auto useSg = watOnly || hasOil || hadOil || oilCondensed; for (int c = 0; c < nc; ++c) { - if (oilPresent[c]) { - rv[c] = rvSat[c]; - cond[c].setFreeOil(); - } + if (useSg[c]) { rv[c] = rvSat[c]; } + else {primalVariable_[c] = PrimalVariables::RV; } + } } - std::copy(&cond[0], &cond[0] + nc, phaseCondition_.begin()); auto ixg = sg < 0; for (int c = 0; c < nc; ++c) { @@ -1449,10 +1425,10 @@ namespace { } // Rs and Rv updates - if (disgas) + if (has_disgas_) std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin()); - if (vapoil) + if (has_vapoil_) std::copy(&rv[0], &rv[0] + nc, state.rv().begin()); @@ -2104,4 +2080,45 @@ namespace { } + template + void + FullyImplicitBlackoilSolver::updatePrimalVariableFromState(const BlackoilState& state) + { + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); + const int np = state.numPhases(); + + const PhaseUsage& pu = fluid_.phaseUsage(); + const DataBlock s = Eigen::Map(& state.saturation()[0], nc, np); + + assert (active_[ Gas ]); + + // reset the primary variables if RV and RS is not set Sg is used as primary variable. + primalVariable_.resize(nc); + std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg); + + if (has_disgas_) { + // Oil/Gas or Water/Oil/Gas system + const V sg = s.col(pu.phase_pos[ Gas ]); + const V so = s.col(pu.phase_pos[ Oil ]); + + for (V::Index c = 0, e = sg.size(); c != e; ++c) { + if ( sg[c] <= 0 && so[c] > 0 ) {primalVariable_[c] = PrimalVariables::RS; } + } + } + + if (has_vapoil_) { + // Oil/Gas or Water/Oil/Gas system + const V sg = s.col(pu.phase_pos[ Gas ]); + const V so = s.col(pu.phase_pos[ Oil ]); + + for (V::Index c = 0, e = so.size(); c != e; ++c) { + if (so[c] <= 0 && sg[c] > 0) {primalVariable_[c] = PrimalVariables::RV; } + } + } + + } + + + } // namespace Opm diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp index 67f7cda81..efad51ff6 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp @@ -74,7 +74,9 @@ namespace Opm const RockCompressibility* rock_comp_props, WellsManager& wells_manager, NewtonIterationBlackoilInterface& linsolver, - const double* gravity); + const double* gravity, + const bool disgas, + const bool vapoil ); /// Run the simulation. /// This will run succesive timesteps until timer.done() is true. It will diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 81fa58fce..fcad56038 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -66,7 +66,9 @@ namespace Opm const RockCompressibility* rock_comp_props, WellsManager& wells_manager, NewtonIterationBlackoilInterface& linsolver, - const double* gravity); + const double* gravity, + bool has_disgas, + bool has_vapoil ); SimulatorReport run(SimulatorTimer& timer, BlackoilState& state, @@ -107,10 +109,12 @@ namespace Opm const RockCompressibility* rock_comp_props, WellsManager& wells_manager, NewtonIterationBlackoilInterface& linsolver, - const double* gravity) + const double* gravity, + const bool has_disgas, + const bool has_vapoil ) { - pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity)); + pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity, has_disgas, has_vapoil)); } @@ -192,7 +196,9 @@ namespace Opm const RockCompressibility* rock_comp_props, WellsManager& wells_manager, NewtonIterationBlackoilInterface& linsolver, - const double* gravity) + const double* gravity, + const bool has_disgas, + const bool has_vapoil) : grid_(grid), props_(props), rock_comp_props_(rock_comp_props), @@ -200,7 +206,7 @@ namespace Opm wells_(wells_manager.c_wells()), gravity_(gravity), geo_(grid_, props_, gravity_), - solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver) + solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver, has_disgas, has_vapoil) /* param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_maxiter", 10),