From 8a166ebbd61f09b85b7a5f5e3f49011dc7c17628 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 12 Jun 2015 06:49:25 +0200 Subject: [PATCH] Add option for solving the wellEq seperatly The well equations is solved seperatly in order to provide more accurate initinal values to the reservoir model. --- opm/autodiff/BlackoilModelBase.hpp | 10 ++ opm/autodiff/BlackoilModelBase_impl.hpp | 204 ++++++++++++++++++++--- opm/autodiff/BlackoilModelParameters.cpp | 4 +- opm/autodiff/BlackoilModelParameters.hpp | 3 + 4 files changed, 200 insertions(+), 21 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 95543310d..559edb616 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -335,9 +335,17 @@ namespace Opm { const WellState& xw, const V& aliveWells); + void + solveWellEq(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state, + WellState& well_state); + void addWellEq(const SolutionState& state, WellState& xw, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, V& aliveWells, std::vector& cq_s); @@ -354,6 +362,8 @@ namespace Opm { void updateWellState(const V& dx, WellState& well_state); + bool getWellConvergence(const int iteration); + std::vector computePressures(const ADB& po, const ADB& sw, diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index cdaa09039..087fcf080 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -478,10 +478,12 @@ namespace detail { std::vector BlackoilModelBase::variableWellsStateIndices() const { - std::vector indices(2, -1); + // Black oil model standard is 5 equation. + // For the pure well solve, only the well equations are picked. + std::vector indices(5, -1); int next = 0; - indices[0] = next++; - indices[1] = next++; + indices[Qs] = next++; + indices[Bhp] = next++; assert(next == 2); return indices; } @@ -749,13 +751,29 @@ namespace detail { asImpl().assembleMassBalanceEq(state); // -------- Well equations ---------- + V aliveWells; const int np = wells().number_of_phases; std::vector cq_s(np, ADB::null()); - asImpl().addWellEq(state, well_state, aliveWells,cq_s); - addWellContribution2MassBalanceEq(cq_s); - addWellControlEq(state, well_state, aliveWells); + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); + + std::vector mob_perfcells(np, ADB::null()); + std::vector b_perfcells(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + mob_perfcells[phase] = subset(rq_[phase].mob,well_cells); + b_perfcells[phase] = subset(rq_[phase].b,well_cells); + } + if (param_.solve_wellEq_initially_ && initial_assembly) { + // solve the well equations as a pre-processing step + solveWellEq(mob_perfcells,b_perfcells,state,well_state); + } + + asImpl().addWellEq(state, well_state, mob_perfcells, b_perfcells, aliveWells,cq_s); + addWellContribution2MassBalanceEq(cq_s); + addWellControlEq(state, well_state, aliveWells); } @@ -831,8 +849,11 @@ namespace detail { template void BlackoilModelBase::addWellEq(const SolutionState& state, WellState& xw, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, V& aliveWells, - std::vector& cq_s) + std::vector& cq_s + ) { if( ! wellsActive() ) return ; @@ -846,17 +867,10 @@ namespace detail { // pressure diffs computed already (once per step, not changing per iteration) const V& cdp = well_perforation_pressure_diffs_; - // Extract needed quantities for the perforation cells const ADB& p_perfcells = subset(state.pressure, well_cells); const ADB& rv_perfcells = subset(state.rv,well_cells); const ADB& rs_perfcells = subset(state.rs,well_cells); - std::vector mob_perfcells(np, ADB::null()); - std::vector b_perfcells(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - mob_perfcells[phase] = subset(rq_[phase].mob,well_cells); - b_perfcells[phase] = subset(rq_[phase].b,well_cells); - } // Perforation pressure const ADB perfpressure = (wops_.w2p * state.bhp) + cdp; @@ -881,7 +895,6 @@ namespace detail { } // HANDLE FLOW INTO WELLBORE - // compute phase volumetric rates at standard conditions std::vector cq_ps(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { @@ -898,7 +911,6 @@ namespace detail { } // HANDLE FLOW OUT FROM WELLBORE - // Using total mobilities ADB total_mob = mob_perfcells[0]; for (int phase = 1; phase < np; ++phase) { @@ -923,7 +935,6 @@ namespace detail { wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps; wbqt += wbq[phase]; } - // compute wellbore mixture at standard conditions. Selector notDeadWells_selector(wbqt.value(), Selector::Zero); std::vector cmix_s(np, ADB::null()); @@ -1156,6 +1167,88 @@ namespace detail { } } + template + void BlackoilModelBase::solveWellEq(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state, + WellState& well_state) + { + V aliveWells; + const int np = wells().number_of_phases; + std::vector cq_s(np, ADB::null()); + + int it = 0; + std::vector vars0; + //bhp and Q for the wells + vars0.reserve(2); + variableWellStateInitials(well_state,vars0); + std::vector vars = ADB::variables(vars0); + std::vector indices = variableWellsStateIndices(); + SolutionState state0 = state; + asImpl().makeConstantState(state0); + SolutionState wellSolutionState = state0; + variableStateExtractWellsVars(indices,vars,wellSolutionState); + std::vector mob_perfcells_const(np, ADB::null()); + std::vector b_perfcells_const(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value()); + b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value()); + } + asImpl().addWellEq(wellSolutionState, well_state, mob_perfcells_const, b_perfcells_const, aliveWells,cq_s); + addWellControlEq(wellSolutionState, well_state, aliveWells); + bool converged = getWellConvergence(it); + while ( (!converged && (it < 15))) { + + std::vector eqs; + eqs.reserve(2); + eqs.push_back(residual_.well_flux_eq); + eqs.push_back(residual_.well_eq); + ADB total_residual = vertcatCollapseJacs(eqs); + const std::vector& Jn = total_residual.derivative(); + const Eigen::SparseLU< M > solver(Jn[0]); + const Eigen::VectorXd& dx = solver.solve(total_residual.value().matrix()); + const int numeq = well_state.numWells()*(well_state.numPhases()+1); + V dx_V = V(numeq); + std::copy_n(dx.data(),numeq, dx_V.data()); + updateWellState(dx_V,well_state); + updateWellControls(well_state); + //bhp and Q for the wells + vars0.clear(); + variableWellStateInitials(well_state,vars0); + vars = ADB::variables(vars0); + wellSolutionState = state0; + variableStateExtractWellsVars(indices,vars,wellSolutionState); + asImpl().addWellEq(wellSolutionState, well_state, mob_perfcells_const, b_perfcells_const, aliveWells,cq_s); + addWellControlEq(wellSolutionState, well_state, aliveWells); + it++; + converged = getWellConvergence(it); + } + if (converged) { + std::cout << "well converged iter: " << it << std::endl; + const int nw = wells().number_of_wells; + { + // We will set the bhp primary variable to the new ones, + // but we do not change the derivatives here. + ADB::V new_bhp = Eigen::Map(well_state.bhp().data(), nw); + // Avoiding the copy below would require a value setter method + // in AutoDiffBlock. + std::vector old_derivs = state.bhp.derivative(); + state.bhp = ADB::function(std::move(new_bhp), std::move(old_derivs)); + } + { + // Need to reshuffle well rates, from phase running fastest + // to wells running fastest. + // The transpose() below switches the ordering. + const DataBlock wrates = Eigen::Map(well_state.wellRates().data(), nw, np).transpose(); + ADB::V new_qs = Eigen::Map(wrates.data(), nw*np); + std::vector old_derivs = state.qs.derivative(); + state.qs = ADB::function(std::move(new_qs), std::move(old_derivs)); + } + computeWellConnectionPressures(state, well_state); + } + + } + @@ -1539,7 +1632,6 @@ namespace detail { const V dbhp = subset(dx, Span(nw, 1, varstart)); varstart += dbhp.size(); assert(varstart == dx.size()); - const double dpmaxrel = dpMaxRel(); @@ -1559,7 +1651,6 @@ namespace detail { const V bhp = bhp_old - dbhp_limited; std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); } - } @@ -1896,7 +1987,7 @@ namespace detail { mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum; converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb); converged_CNV = converged_CNV && (CNV[idx] < tol_cnv); - well_flux_residual[idx] = B_avg[idx] * dt * maxNormWell[idx]; + well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx]; converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); } @@ -1946,6 +2037,79 @@ namespace detail { return converged; } + template + bool + BlackoilModelBase::getWellConvergence(const int iteration) + { + const double tol_wells = param_.tolerance_wells_; + + const int nc = Opm::AutoDiffGrid::numCells(grid_); + const int nw = wellsActive() ? wells().number_of_wells : 0; + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + + const V pv = geo_.poreVolume(); + std::array R_sum = {{0., 0., 0.}}; + std::array B_avg = {{0., 0., 0.}}; + std::array maxCoeff = {{0., 0., 0.}}; + std::array well_flux_residual = {{0., 0., 0.}}; + std::size_t cols = MaxNumPhases; // needed to pass the correct type to Eigen + Eigen::Array B(nc, cols); + Eigen::Array R(nc, cols); + Eigen::Array tempV(nc, cols); + std::vector maxNormWell(MaxNumPhases); + for ( int idx=0; idx maxResidualAllowed() || + std::isnan(well_flux_residual[Oil]) || well_flux_residual[Oil] > maxResidualAllowed() || + std::isnan(well_flux_residual[Gas]) || well_flux_residual[Gas] > maxResidualAllowed() ) + { + OPM_THROW(Opm::NumericalProblem,"One of the well residuals is NaN or to large!"); + } + + if ( terminal_output_ ) + { + // Only rank 0 does print to std::cout + if (iteration == 0) { + std::cout << "\nIter W-FLUX(W) W-FLUX(O) W-FLUX(G)\n"; + } + const std::streamsize oprec = std::cout.precision(3); + const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific); + std::cout << std::setw(4) << iteration + << std::setw(11) << well_flux_residual[Water] + << std::setw(11) << well_flux_residual[Oil] + << std::setw(11) << well_flux_residual[Gas] + << std::endl; + std::cout.precision(oprec); + std::cout.flags(oflags); + } + return converged; + } + diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index 7a18c2fed..854284f18 100644 --- a/opm/autodiff/BlackoilModelParameters.cpp +++ b/opm/autodiff/BlackoilModelParameters.cpp @@ -46,6 +46,7 @@ namespace Opm tolerance_mb_ = param.getDefault("tolerance_mb", tolerance_mb_); tolerance_cnv_ = param.getDefault("tolerance_cnv", tolerance_cnv_); tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ ); + solve_wellEq_initially_ = param.getDefault("solve_wellEq_initially",solve_wellEq_initially_); } @@ -60,7 +61,8 @@ namespace Opm max_residual_allowed_ = 1e7; tolerance_mb_ = 1.0e-5; tolerance_cnv_ = 1.0e-2; - tolerance_wells_ = 5.0e-1; + tolerance_wells_ = 1.0e-3; + solve_wellEq_initially_ = false; } diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index 054345cbe..c50ce4cce 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -43,6 +43,9 @@ namespace Opm /// Well convergence tolerance. double tolerance_wells_; + /// Solve well equation initially + bool solve_wellEq_initially_; + /// Construct from user parameters or defaults. explicit BlackoilModelParameters( const parameter::ParameterGroup& param );