diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index d4de95e30..348afc765 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -294,15 +294,32 @@ namespace Opm { std::vector variableStateInitials(const ReservoirState& x, const WellState& xw) const; + void + variableReservoirStateInitials(const ReservoirState& x, + std::vector& vars0) const; + void + variableWellStateInitials(const WellState& xw, + std::vector& vars0) const; std::vector variableStateIndices() const; + std::vector + variableWellStateIndices() const; + + void + addWellContributionToMassBalanceEq(const std::vector& cq_s); + SolutionState variableStateExtractVars(const ReservoirState& x, const std::vector& indices, std::vector& vars) const; + void + variableStateExtractWellsVars(const std::vector& indices, + std::vector& vars, + SolutionState& state) const; + void computeAccum(const SolutionState& state, const int aix ); @@ -318,10 +335,19 @@ 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, - V& aliveWells); + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + V& aliveWells, + std::vector& cq_s); void extraAddWellEq(const SolutionState& state, @@ -333,6 +359,11 @@ namespace Opm { void updateWellControls(WellState& xw) const; + void updateWellState(const V& dwells, + 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 96b756d72..18381a797 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -187,6 +187,8 @@ namespace detail { + + template void BlackoilModelBase:: @@ -203,6 +205,7 @@ namespace detail { + template void BlackoilModelBase:: @@ -216,6 +219,7 @@ namespace detail { + template int BlackoilModelBase:: @@ -227,6 +231,7 @@ namespace detail { + template int BlackoilModelBase:: @@ -238,6 +243,7 @@ namespace detail { + template bool BlackoilModelBase:: @@ -249,6 +255,7 @@ namespace detail { + template int BlackoilModelBase:: @@ -260,6 +267,7 @@ namespace detail { + template void BlackoilModelBase:: @@ -281,6 +289,7 @@ namespace detail { + template BlackoilModelBase::ReservoirResidualQuant::ReservoirResidualQuant() : accum(2, ADB::null()) @@ -381,13 +390,28 @@ namespace detail { { assert(active_[ Oil ]); - using namespace Opm::AutoDiffGrid; - const int nc = numCells(grid_); const int np = x.numPhases(); std::vector vars0; // p, Sw and Rs, Rv or Sg is used as primary depending on solution conditions + // and bhp and Q for the wells vars0.reserve(np + 1); + variableReservoirStateInitials(x, vars0); + variableWellStateInitials(xw, vars0); + return vars0; + } + + + + + + template + void + BlackoilModelBase::variableReservoirStateInitials(const ReservoirState& x, std::vector& vars0) const + { + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); + const int np = x.numPhases(); // Initial pressure. assert (not x.pressure().empty()); const V p = Eigen::Map(& x.pressure()[0], nc, 1); @@ -413,14 +437,23 @@ namespace detail { xvar = isRs_*rs + isRv_*rv + isSg_*sg; vars0.push_back(xvar); } + } + + + + template + void + BlackoilModelBase::variableWellStateInitials(const WellState& xw, std::vector& vars0) const + { // Initial well rates. if ( wellsActive() ) { // Need to reshuffle well rates, from phase running fastest // to wells running fastest. const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; // The transpose() below switches the ordering. const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); @@ -438,8 +471,6 @@ namespace detail { vars0.push_back(V()); vars0.push_back(V()); } - - return vars0; } @@ -469,6 +500,23 @@ namespace detail { + template + std::vector + BlackoilModelBase::variableWellStateIndices() const + { + // 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[Qs] = next++; + indices[Bhp] = next++; + assert(next == 2); + return indices; + } + + + + template typename BlackoilModelBase::SolutionState @@ -533,14 +581,26 @@ namespace detail { state.saturation[pu.phase_pos[ Oil ]] = std::move(so); } } + // wells + variableStateExtractWellsVars(indices, vars, state); + return state; + } + + + + + template + void + BlackoilModelBase::variableStateExtractWellsVars(const std::vector& indices, + std::vector& vars, + SolutionState& state) const + { // Qs. state.qs = std::move(vars[indices[Qs]]); // Bhp. state.bhp = std::move(vars[indices[Bhp]]); - - return state; } @@ -726,9 +786,33 @@ namespace detail { asImpl().assembleMassBalanceEq(state); // -------- Well equations ---------- + + if ( ! wellsActive() ) { + return; + } + V aliveWells; - asImpl().addWellEq(state, well_state, aliveWells); - addWellControlEq(state, well_state, aliveWells); + const int np = wells().number_of_phases; + std::vector cq_s(np, ADB::null()); + + 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); + addWellContributionToMassBalanceEq(cq_s); + addWellControlEq(state, well_state, aliveWells); } @@ -790,13 +874,35 @@ namespace detail { template - void BlackoilModelBase::addWellEq(const SolutionState& state, - WellState& xw, - V& aliveWells) + void + BlackoilModelBase::addWellContributionToMassBalanceEq(const std::vector& cq_s) + { + // Add well contributions to mass balance equations + const int nc = Opm::AutoDiffGrid::numCells(grid_); + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + const int np = wells().number_of_phases; + const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); + for (int phase = 0; phase < np; ++phase) { + residual_.material_balance_eq[phase] -= superset(cq_s[phase], well_cells, nc); + } + } + + + + + + 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) { if( ! wellsActive() ) return ; - const int nc = Opm::AutoDiffGrid::numCells(grid_); const int np = wells().number_of_phases; const int nw = wells().number_of_wells; const int nperf = wells().well_connpos[nw]; @@ -806,17 +912,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); - } + const ADB& rv_perfcells = subset(state.rv, well_cells); + const ADB& rs_perfcells = subset(state.rs, well_cells); // Perforation pressure const ADB perfpressure = (wops_.w2p * state.bhp) + cdp; @@ -841,7 +940,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) { @@ -858,7 +956,6 @@ namespace detail { } // HANDLE FLOW OUT FROM WELLBORE - // Using total mobilities ADB total_mob = mob_perfcells[0]; for (int phase = 1; phase < np; ++phase) { @@ -883,7 +980,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()); @@ -913,16 +1009,10 @@ namespace detail { ADB cqt_is = cqt_i/volumeRatio; // connection phase volumerates at standard conditions - std::vector cq_s(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; } - // Add well contributions to mass balance equations - for (int phase = 0; phase < np; ++phase) { - residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc); - } - // WELL EQUATIONS ADB qs = state.qs; for (int phase = 0; phase < np; ++phase) { @@ -1048,6 +1138,7 @@ namespace detail { + template void BlackoilModelBase::updateWellControls(WellState& xw) const { @@ -1124,6 +1215,91 @@ 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()); + std::vector indices = variableWellStateIndices(); + SolutionState state0 = state; + asImpl().makeConstantState(state0); + + 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()); + } + + int it = 0; + bool converged; + do { + // bhp and Q for the wells + std::vector vars0; + vars0.reserve(2); + variableWellStateInitials(well_state, vars0); + std::vector vars = ADB::variables(vars0); + + SolutionState 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); + converged = getWellConvergence(it); + + if (converged) { + break; + } + + ++it; + 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()); + assert(dx.size() == (well_state.numWells() * (well_state.numPhases()+1))); + updateWellState(dx.array(), well_state); + updateWellControls(well_state); + + } while (it < 15); + + 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); + } + + } + + + + + template void BlackoilModelBase::addWellControlEq(const SolutionState& state, const WellState& xw, @@ -1285,10 +1461,10 @@ namespace detail { const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null; varstart += dxvar.size(); - const V dqs = subset(dx, Span(np*nw, 1, varstart)); - varstart += dqs.size(); - const V dbhp = subset(dx, Span(nw, 1, varstart)); - varstart += dbhp.size(); + // Extract well parts np phase rates + bhp + const V dwells = subset(dx, Span((np+1)*nw, 1, varstart)); + varstart += dwells.size(); + assert(varstart == dx.size()); // Pressure update. @@ -1476,8 +1652,38 @@ namespace detail { std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin()); } + + updateWellState(dwells,well_state); + + // Update phase conditions used for property calculations. + updatePhaseCondFromPrimalVariable(); + } + + + + + + template + void + BlackoilModelBase::updateWellState(const V& dwells, + WellState& well_state) + { + if( wellsActive() ) { + const int np = wells().number_of_phases; + const int nw = wellsActive() ? wells().number_of_wells : 0; + + // Extract parts of dwells corresponding to each part. + int varstart = 0; + const V dqs = subset(dwells, Span(np*nw, 1, varstart)); + varstart += dqs.size(); + const V dbhp = subset(dwells, Span(nw, 1, varstart)); + varstart += dbhp.size(); + assert(varstart == dwells.size()); + const double dpmaxrel = dpMaxRel(); + + // Qs update. // Since we need to update the wellrates, that are ordered by wells, // from dqs which are ordered by phase, the simplest is to compute @@ -1494,9 +1700,6 @@ namespace detail { const V bhp = bhp_old - dbhp_limited; std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); } - - // Update phase conditions used for property calculations. - updatePhaseCondFromPrimalVariable(); } @@ -1654,6 +1857,7 @@ namespace detail { + template std::vector BlackoilModelBase::computeResidualNorms() const @@ -1696,6 +1900,7 @@ namespace detail { + template double BlackoilModelBase::convergenceReduction(const Eigen::Array& B, @@ -1780,6 +1985,7 @@ namespace detail { + template bool BlackoilModelBase::getConvergence(const double dt, const int iteration) @@ -1832,7 +2038,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); } @@ -1854,7 +2060,7 @@ namespace detail { std::isnan(well_flux_residual[Gas]) || well_flux_residual[Gas] > maxResidualAllowed() || std::isnan(residualWell) || residualWell > maxResidualAllowed() ) { - OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or to large!"); + OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or too large!"); } if ( terminal_output_ ) @@ -1885,6 +2091,84 @@ namespace detail { + + 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 too 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; + } + + + + + template ADB BlackoilModelBase::fluidViscosity(const int phase, @@ -1990,6 +2274,9 @@ namespace detail { } + + + template V BlackoilModelBase::fluidRvSat(const V& p, @@ -2014,6 +2301,8 @@ namespace detail { + + template ADB BlackoilModelBase::poroMult(const ADB& p) const diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index 7a18c2fed..5c5530d78 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..15e7042db 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 ); diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index d5b6ea344..ecd4de5dd 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -153,7 +153,7 @@ namespace Opm auto solver = asImpl().createSolver(wells); // If sub stepping is enabled allow the solver to sub cycle - // in case the report steps are to large for the solver to converge + // in case the report steps are too large for the solver to converge // // \Note: The report steps are met in any case // \Note: The sub stepping will require a copy of the state variables