From 07339be80e77c92ea11b776bf837da227c0df92f Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 29 Sep 2016 14:44:21 +0200 Subject: [PATCH 1/4] BUGFIX. Make frankenstein compile and run - remove outout of initial timestep - pass phaseUsage to report(...) in WellStateFullyImplicitBlackoil.hpp --- opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp | 2 +- opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp | 2 +- opm/autodiff/WellStateFullyImplicitBlackoil.hpp | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 86a9cf76e..9ed1c4360 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -239,7 +239,7 @@ public: // write the inital state at the report stage if (timer.initialStep()) { - output_writer_.writeTimeStep( timer, state, well_state, solver->model() ); + //output_writer_.writeTimeStep( timer, state, well_state, solver->model() ); } if( terminal_output_ ) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp index acb29bf0f..3816544af 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp @@ -167,7 +167,7 @@ namespace Opm substep, timer.simulationTimeElapsed(), simToSolution( state, phaseUsage_ ), - wellState.report(), + wellState.report(phaseUsage_), simProps); } } diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index ab228346b..67a5951c4 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -276,8 +276,8 @@ namespace Opm std::vector& wellSolutions() { return well_solutions_; } const std::vector& wellSolutions() const { return well_solutions_; } - data::Wells report() const override { - data::Wells res = WellState::report(); + data::Wells report(const PhaseUsage& pu) const override { + data::Wells res = WellState::report(pu); const int nw = this->numWells(); // If there are now wells numPhases throws a floating point From 84c33e0e9697062954c925f9cc101e3838ad6ffe Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 29 Sep 2016 15:29:55 +0200 Subject: [PATCH 2/4] Add trivial well equation when there is no flow Fixes case when all perforations in a well are crossflowing while crossflow is not allowed --- opm/autodiff/StandardWellsDense.hpp | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index aff5512ca..873f1e5ef 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -226,6 +226,7 @@ namespace Opm { const double volume = 0.002831684659200; // 0.1 cu ft; for (int w = 0; w < nw; ++w) { + double total_flux = 0.0; for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { @@ -243,6 +244,7 @@ namespace Opm { // subtract sum of phase fluxes in the well equations. resWell_[w][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].value; + total_flux += cq_s[p1].value; // assemble the jacobians for (int p2 = 0; p2 < np; ++p2) { @@ -260,15 +262,23 @@ namespace Opm { // Store the perforation pressure for later usage. well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf]; } - - // add vol * dF/dt + Q to the well equations; - for (int p1 = 0; p1 < np; ++p1) { - EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; - resWell_loc += getQs(w, p1); - for (int p2 = 0; p2 < np; ++p2) { - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivatives[p2+3]; + if (std::abs(total_flux) == 0) { // add a trivial equation for the case when there is no flow + for (int p1 = 0; p1 < np; ++p1) { + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p1)] += 1; + resWell_[w][flowPhaseToEbosCompIdx(p1)] += 0; + } + } + else { + + // add vol * dF/dt + Q to the well equations; + for (int p1 = 0; p1 < np; ++p1) { + EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; + resWell_loc += getQs(w, p1); + for (int p2 = 0; p2 < np; ++p2) { + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivatives[p2+3]; + } + resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value; } - resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value; } } From 69811e9585f3a8c8423c4dea0a14f7033bf69433 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 29 Sep 2016 15:30:59 +0200 Subject: [PATCH 3/4] Use correct pvtregions in rsSat and rvSat calculations --- opm/autodiff/BlackoilModelEbos.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 9c7e5ab43..a46e059a2 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -523,6 +523,8 @@ namespace Opm { // phase translation sg <-> rs const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx]; + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); switch (hydroCarbonState) { case HydroCarbonState::GasAndOil: { double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]]; @@ -534,7 +536,7 @@ namespace Opm { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs sg = 0; so = 1.0 - sw - sg; - double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); double& rs = reservoir_state.gasoilratio()[cell_idx]; rs = rsSat*(1-epsilon); } else if (so <= 0.0 && has_vapoil_) { @@ -543,7 +545,7 @@ namespace Opm { sg = 1.0 - sw - so; double& rv = reservoir_state.rv()[cell_idx]; // use gas pressure? - double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); rv = rvSat*(1-epsilon); } break; @@ -551,8 +553,7 @@ namespace Opm { case HydroCarbonState::OilOnly: { double& rs = reservoir_state.gasoilratio()[cell_idx]; double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; - // TODO:: not hardcode pvtRegion = 0 - double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); if (rs > ( rsSat * (1+epsilon) ) ) { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; sg = epsilon; @@ -563,7 +564,7 @@ namespace Opm { } case HydroCarbonState::GasOnly: { double& rv = reservoir_state.rv()[cell_idx]; - double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); if (rv > rvSat * (1+epsilon) ) { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; so = epsilon; From 8e2657ce4d4a40c84fa44dbdf9c66ff5b1cf655d Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 30 Sep 2016 08:29:19 +0200 Subject: [PATCH 4/4] Some tuning in the Appleyard --- opm/autodiff/BlackoilModelEbos.hpp | 40 ++++++++++++++++++++++++----- opm/autodiff/StandardWellsDense.hpp | 2 +- 2 files changed, 34 insertions(+), 8 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index a46e059a2..5be2eeddf 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -456,7 +456,10 @@ namespace Opm { for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { double dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)]; - reservoir_state.pressure()[cell_idx] -= dp; + //reservoir_state.pressure()[cell_idx] -= dp; + double& p = reservoir_state.pressure()[cell_idx]; + p -= dp; + p = std::max(p, 1e5); // Saturation updates. const double dsw = active_[Water] ? dx[cell_idx][flowPhaseToEbosCompIdx(1)] : 0.0; @@ -512,14 +515,22 @@ namespace Opm { if (has_disgas_) { double& rs = reservoir_state.gasoilratio()[cell_idx]; rs -= drs; + rs = std::max(rs, 0.0); + } if (has_vapoil_) { double& rv = reservoir_state.rv()[cell_idx]; rv -= drv; + rv = std::max(rv, 0.0); } // Sg is used as primal variable for water only cells. const double epsilon = 1e-4; //std::sqrt(std::numeric_limits::epsilon()); + double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]]; + double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; + double& rs = reservoir_state.gasoilratio()[cell_idx]; + double& rv = reservoir_state.rv()[cell_idx]; + // phase translation sg <-> rs const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx]; @@ -527,11 +538,10 @@ namespace Opm { const auto& fs = intQuants.fluidState(); switch (hydroCarbonState) { case HydroCarbonState::GasAndOil: { - double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]]; - double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; if (sw > (1.0 - epsilon)) // water only i.e. do nothing break; + if (sg <= 0.0 && has_disgas_) { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs sg = 0; @@ -551,8 +561,18 @@ namespace Opm { break; } case HydroCarbonState::OilOnly: { - double& rs = reservoir_state.gasoilratio()[cell_idx]; - double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; + if (sw > (1.0 - epsilon)) { + // water only change to Sg + rs = 0; + rv = 0; + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + //std::cout << "watonly rv -> sg" << cell_idx << std::endl; + break; + } + + + + double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); if (rs > ( rsSat * (1+epsilon) ) ) { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; @@ -563,13 +583,19 @@ namespace Opm { break; } case HydroCarbonState::GasOnly: { - double& rv = reservoir_state.rv()[cell_idx]; + if (sw > (1.0 - epsilon)) { + // water only change to Sg + rs = 0; + rv = 0; + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + //std::cout << "watonly rv -> sg" << cell_idx << std::endl; + break; + } double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); if (rv > rvSat * (1+epsilon) ) { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; so = epsilon; rv = rvSat; - double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; sg -= epsilon; } break; diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 873f1e5ef..f71eab732 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -1007,7 +1007,7 @@ namespace Opm { { const int sign1 = dwells[w][flowPhaseToEbosCompIdx(0)] > 0 ? 1: -1; const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(0)]),std::abs(xvar_well_old[w])*dBHPLimit); - well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited; + well_state.wellSolutions()[w] = std::max(xvar_well_old[w] - dx1_limited,1e5); well_state.bhp()[w] = well_state.wellSolutions()[w]; if (well_controls_iget_type(wc, current) == SURFACE_RATE) {