From 463e4bc5e332b39fdce9b7e955004988bd2a973f Mon Sep 17 00:00:00 2001 From: Robert K Date: Mon, 19 Jan 2015 14:14:18 +0100 Subject: [PATCH] BlockOilSimulator: allow to run without wells (mainly for testing and debugging). --- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 14 +- .../FullyImplicitBlackoilSolver_impl.hpp | 190 +++++++++++------- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 35 ++-- .../SimulatorFullyImplicitBlackoil_impl.hpp | 23 ++- .../WellStateFullyImplicitBlackoil.hpp | 4 +- 5 files changed, 159 insertions(+), 107 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 8c21f77de..868e3b17e 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -90,7 +90,7 @@ namespace Opm { const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , const RockCompressibility* rock_comp_props, - const Wells& wells, + const Wells* wells, const NewtonIterationBlackoilInterface& linsolver, const bool has_disgas, const bool has_vapoil ); @@ -147,11 +147,11 @@ namespace Opm { ADB rs; ADB rv; ADB qs; - ADB bhp; + ADB bhp; }; struct WellOps { - WellOps(const Wells& wells); + WellOps(const Wells* wells); M w2p; // well -> perf (scatter) M p2w; // perf -> well (gather) }; @@ -169,7 +169,7 @@ namespace Opm { const BlackoilPropsAdInterface& fluid_; const DerivedGeology& geo_; const RockCompressibility* rock_comp_props_; - const Wells& wells_; + const Wells* wells_; const NewtonIterationBlackoilInterface& linsolver_; // For each canonical phase -> true if active const std::vector active_; @@ -194,6 +194,12 @@ namespace Opm { std::vector primalVariable_; // Private methods. + + // return true if wells are available + bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; } + // return wells object + const Wells& wells () const { assert( wells_ ); return *wells_; } + SolutionState constantState(const BlackoilState& x, const WellStateFullyImplicitBlackoil& xw); diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 271c5a305..768ab4e02 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -211,7 +211,7 @@ namespace { const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , const RockCompressibility* rock_comp_props, - const Wells& wells, + const Wells* wells, const NewtonIterationBlackoilInterface& linsolver, const bool has_disgas, const bool has_vapoil) @@ -225,7 +225,7 @@ namespace { , canph_ (active2Canonical(fluid.phaseUsage())) , cells_ (buildAllCells(Opm::AutoDiffGrid::numCells(grid))) , ops_ (grid) - , wops_ (wells) + , wops_ (wells_) , has_disgas_(has_disgas) , has_vapoil_(has_vapoil) , param_( param ) @@ -372,30 +372,34 @@ namespace { template FullyImplicitBlackoilSolver:: - WellOps::WellOps(const Wells& wells) - : w2p(wells.well_connpos[ wells.number_of_wells ], - wells.number_of_wells) - , p2w(wells.number_of_wells, - wells.well_connpos[ wells.number_of_wells ]) + WellOps::WellOps(const Wells* wells) + : w2p(), + p2w() { - const int nw = wells.number_of_wells; - const int* const wpos = wells.well_connpos; + if( wells ) + { + w2p = M(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); + p2w = M(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); - typedef Eigen::Triplet Tri; + const int nw = wells->number_of_wells; + const int* const wpos = wells->well_connpos; - std::vector scatter, gather; - scatter.reserve(wpos[nw]); - gather .reserve(wpos[nw]); + typedef Eigen::Triplet Tri; - for (int w = 0, i = 0; w < nw; ++w) { - for (; i < wpos[ w + 1 ]; ++i) { - scatter.push_back(Tri(i, w, 1.0)); - gather .push_back(Tri(w, i, 1.0)); + std::vector scatter, gather; + scatter.reserve(wpos[nw]); + gather .reserve(wpos[nw]); + + for (int w = 0, i = 0; w < nw; ++w) { + for (; i < wpos[ w + 1 ]; ++i) { + scatter.push_back(Tri(i, w, 1.0)); + gather .push_back(Tri(w, i, 1.0)); + } } - } - w2p.setFromTriplets(scatter.begin(), scatter.end()); - p2w.setFromTriplets(gather .begin(), gather .end()); + w2p.setFromTriplets(scatter.begin(), scatter.end()); + p2w.setFromTriplets(gather .begin(), gather .end()); + } } @@ -490,21 +494,29 @@ namespace { } - // Initial well rates. - assert (not xw.wellRates().empty()); - // Need to reshuffle well rates, from phase running fastest - // to wells running fastest. - const int nw = wells_.number_of_wells; - // The transpose() below switches the ordering. - const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); - const V qs = Eigen::Map(wrates.data(), nw*np); - vars0.push_back(qs); + if( wellsActive() ) + { + // Need to reshuffle well rates, from phase running fastest + // to wells running fastest. + const int nw = wells().number_of_wells; - // Initial well bottom-hole pressure. - assert (not xw.bhp().empty()); - const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); - vars0.push_back(bhp); + // The transpose() below switches the ordering. + const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); + const V qs = Eigen::Map(wrates.data(), nw*np); + vars0.push_back(qs); + + // Initial well bottom-hole pressure. + assert (not xw.bhp().empty()); + const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); + vars0.push_back(bhp); + } + else + { + // push null states for qs and bhp + vars0.push_back(V()); + vars0.push_back(V()); + } std::vector vars = ADB::variables(vars0); @@ -564,6 +576,7 @@ namespace { state.saturation[pu.phase_pos[ Oil ]] = so; } } + // Qs. state.qs = vars[ nextvar++ ]; @@ -631,12 +644,14 @@ namespace { void FullyImplicitBlackoilSolver::computeWellConnectionPressures(const SolutionState& state, const WellStateFullyImplicitBlackoil& xw) { + if( ! wellsActive() ) return ; + using namespace Opm::AutoDiffGrid; // 1. Compute properties required by computeConnectionPressureDelta(). // Note that some of the complexity of this part is due to the function // taking std::vector arguments, and not Eigen objects. - const int nperf = wells_.well_connpos[wells_.number_of_wells]; - const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); // Compute b, rsmax, rvmax values for perforations. const std::vector pressures = computePressures(state); const ADB perf_temp = subset(state.temperature, well_cells); @@ -694,7 +709,7 @@ namespace { // 2. Compute pressure deltas, and store the results. std::vector cdp = WellDensitySegmented - ::computeConnectionPressureDelta(wells_, xw, fluid_.phaseUsage(), + ::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(), b_perf, rssat_perf, rvsat_perf, perf_depth, surf_dens, grav); well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); @@ -796,13 +811,15 @@ namespace { WellStateFullyImplicitBlackoil& xw, V& aliveWells) { + 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]; + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - V Tw = Eigen::Map(wells_.WI, nperf); - const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); + V Tw = Eigen::Map(wells().WI, nperf); + const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); // pressure diffs computed already (once per step, not changing per iteration) const V& cdp = well_perforation_pressure_diffs_; @@ -824,7 +841,7 @@ namespace { // injector == 1, producer == 0 V isInj = V::Zero(nw); for (int w = 0; w < nw; ++w) { - if (wells_.type[w] == INJECTOR) { + if (wells().type[w] == INJECTOR) { isInj[w] = 1; } } @@ -889,7 +906,7 @@ namespace { } // compute avg. and total wellbore phase volumetric rates at std. conds - const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np); + const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); std::vector wbq(np, ADB::null()); ADB wbqt = ADB::constant(V::Zero(nw), state.pressure.blockPattern()); for (int phase = 0; phase < np; ++phase) { @@ -1083,15 +1100,17 @@ namespace { ADB& well_phase_flow_rate, WellStateFullyImplicitBlackoil& xw) const { + if( ! wellsActive() ) return ; + std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" }; // Find, for each well, if any constraints are broken. If so, // switch control to first broken constraint. - const int np = wells_.number_of_phases; - const int nw = wells_.number_of_wells; + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; bool bhp_changed = false; bool rates_changed = false; for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells_.ctrls[w]; + const WellControls* wc = wells().ctrls[w]; // The current control in the well state overrides // the current control set in the Wells struct, which // is instead treated as a default. @@ -1108,14 +1127,14 @@ namespace { // inequality constraint, and therefore skipped. continue; } - if (constraintBroken(bhp, well_phase_flow_rate, w, np, wells_.type[w], wc, ctrl_index)) { + if (constraintBroken(bhp, well_phase_flow_rate, w, np, wells().type[w], wc, ctrl_index)) { // ctrl_index will be the index of the broken constraint after the loop. break; } } if (ctrl_index != nwc) { // Constraint number ctrl_index was broken, switch to it. - std::cout << "Switching control mode for well " << wells_.name[w] + std::cout << "Switching control mode for well " << wells().name[w] << " from " << modestring[well_controls_iget_type(wc, current)] << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; xw.currentControls()[w] = ctrl_index; @@ -1176,14 +1195,16 @@ namespace { const WellStateFullyImplicitBlackoil& xw, const V& aliveWells) { - const int np = wells_.number_of_phases; - const int nw = wells_.number_of_wells; + if( ! wellsActive() ) return; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; V bhp_targets = V::Zero(nw); V rate_targets = V::Zero(nw); M rate_distr(nw, np*nw); for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells_.ctrls[w]; + const WellControls* wc = wells().ctrls[w]; // The current control in the well state overrides // the current control set in the Wells struct, which // is instead treated as a default. @@ -1254,6 +1275,17 @@ namespace { struct Chop01 { double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); } }; + + double infinityNorm( const ADB& a ) + { + if( a.value().size() > 0 ) { + return a.value().matrix().template lpNorm (); + } + else { // this situation can occur when no wells are present + return 0.0; + } + } + } @@ -1268,7 +1300,7 @@ namespace { using namespace Opm::AutoDiffGrid; const int np = fluid_.numPhases(); const int nc = numCells(grid_); - const int nw = wells_.number_of_wells; + const int nw = wellsActive() ? wells().number_of_wells : 0; const V null; assert(null.size() == 0); const V zero = V::Zero(nc); @@ -1497,21 +1529,24 @@ namespace { std::copy(&rv[0], &rv[0] + nc, state.rv().begin()); } - // 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 - // dwr, which is the data from dqs but ordered by wells. - const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); - const V dwr = Eigen::Map(wwr.data(), nw*np); - const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); - const V wr = wr_old - dwr; - std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + if( wellsActive() ) + { + // 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 + // dwr, which is the data from dqs but ordered by wells. + const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); + const V dwr = Eigen::Map(wwr.data(), nw*np); + const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); + const V wr = wr_old - dwr; + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); - // Bhp update. - const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); - const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); - const V bhp = bhp_old - dbhp_limited; - std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + // Bhp update. + const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); + const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); + 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(); @@ -1615,8 +1650,8 @@ namespace { const DataBlock& well_s, const std::vector& well_cells) const { - const int nw = wells_.number_of_wells; - const int nperf = wells_.well_connpos[nw]; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; const std::vector& bpat = state.pressure.blockPattern(); const ADB null = ADB::constant(V::Zero(nperf), bpat); @@ -1753,7 +1788,7 @@ namespace { const std::vector::const_iterator endMassBalanceIt = residual_.material_balance_eq.end(); for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) { - const double massBalanceResid = (*massBalanceIt).value().matrix().template lpNorm(); + const double massBalanceResid = infinityNorm( (*massBalanceIt) ); if (!std::isfinite(massBalanceResid)) { OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual"); @@ -1762,14 +1797,14 @@ namespace { } // the following residuals are not used in the oscillation detection now - const double wellFluxResid = residual_.well_flux_eq.value().matrix().template lpNorm(); + const double wellFluxResid = infinityNorm( residual_.well_flux_eq ); if (!std::isfinite(wellFluxResid)) { - OPM_THROW(Opm::NumericalProblem, + OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual"); } residual.push_back(wellFluxResid); - const double wellResid = residual_.well_eq.value().matrix().template lpNorm(); + const double wellResid = infinityNorm( residual_.well_eq ); if (!std::isfinite(wellResid)) { OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual"); @@ -1917,16 +1952,15 @@ namespace { const double mass_balance_residual_gas = fabs(BG_avg*RG_sum) * dt / pvSum; bool converged_MB = (mass_balance_residual_water < tol_mb) - && (mass_balance_residual_oil< tol_mb) - && (mass_balance_residual_gas < tol_mb); + && (mass_balance_residual_oil < tol_mb) + && (mass_balance_residual_gas < tol_mb); bool converged_CNV = (CNVW < tol_cnv) && (CNVO < tol_cnv) && (CNVG < tol_cnv); - double residualWellFlux = residual_.well_flux_eq.value().matrix().template lpNorm(); - double residualWell = residual_.well_eq.value().matrix().template lpNorm(); + const double residualWellFlux = infinityNorm( residual_.well_flux_eq ); + const double residualWell = infinityNorm( residual_.well_eq ); bool converged_Well = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa); - bool converged = converged_MB && converged_CNV && converged_Well; if (iteration == 0) { @@ -2241,7 +2275,7 @@ namespace { // For oil only cells Rs is used as primal variable. For cells almost full of water // the default primal variable (Sg) is used. - if (has_disgas_) { + if (has_disgas_) { for (V::Index c = 0, e = sg.size(); c != e; ++c) { if ( !watOnly[c] && hasOil[c] && !hasGas[c] ) {primalVariable_[c] = PrimalVariables::RS; } } diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index 842c2650b..0eb295e8b 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -142,17 +142,23 @@ namespace Opm for (int phase = 0; phase < np; ++phase) { eqs.push_back(residual.material_balance_eq[phase]); } - eqs.push_back(residual.well_flux_eq); - eqs.push_back(residual.well_eq); - // Eliminate the well-related unknowns, and corresponding equations. + // check if wells are present + const bool hasWells = residual.well_flux_eq.size() > 0 ; std::vector elim_eqs; - elim_eqs.reserve(2); - elim_eqs.push_back(eqs[np]); - eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns. - elim_eqs.push_back(eqs[np]); - eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns. - assert(int(eqs.size()) == np); + if( hasWells ) + { + eqs.push_back(residual.well_flux_eq); + eqs.push_back(residual.well_eq); + + // Eliminate the well-related unknowns, and corresponding equations. + elim_eqs.reserve(2); + elim_eqs.push_back(eqs[np]); + eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns. + elim_eqs.push_back(eqs[np]); + eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns. + assert(int(eqs.size()) == np); + } // Scale material balance equations. const double matbalscale[3] = { 1.1169, 1.0031, 0.0031 }; // HACK hardcoded instead of computed. @@ -219,10 +225,13 @@ namespace Opm // Copy solver output to dx. std::copy(x.begin(), x.end(), dx.data()); - // Compute full solution using the eliminated equations. - // Recovery in inverse order of elimination. - dx = recoverVariable(elim_eqs[1], dx, np); - dx = recoverVariable(elim_eqs[0], dx, np); + if( hasWells ) + { + // Compute full solution using the eliminated equations. + // Recovery in inverse order of elimination. + dx = recoverVariable(elim_eqs[1], dx, np); + dx = recoverVariable(elim_eqs[0], dx, np); + } return dx; } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index dc5c995d0..0b1301edd 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -356,7 +356,7 @@ namespace Opm // Run a multiple steps of the solver depending on the time step control. solver_timer.start(); - FullyImplicitBlackoilSolver solver(solverParam, grid_, props_, geo_, rock_comp_props_, *wells, solver_, has_disgas_, has_vapoil_); + FullyImplicitBlackoilSolver solver(solverParam, grid_, props_, geo_, rock_comp_props_, wells, solver_, has_disgas_, has_vapoil_); if (!threshold_pressures_by_face_.empty()) { solver.setThresholdPressures(threshold_pressures_by_face_); } @@ -478,18 +478,20 @@ namespace Opm } inline std::vector - resvProducers(const Wells& wells, + resvProducers(const Wells* wells, const std::size_t step, const WellMap& wmap) { std::vector resv_prod; - - for (int w = 0, nw = wells.number_of_wells; w < nw; ++w) { - if (is_resv_prod(wells, w) || - ((wells.name[w] != 0) && - is_resv_prod(wmap, wells.name[w], step))) - { - resv_prod.push_back(w); + if( wells ) + { + for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) { + if (is_resv_prod(*wells, w) || + ((wells->name[w] != 0) && + is_resv_prod(wmap, wells->name[w], step))) + { + resv_prod.push_back(w); + } } } @@ -541,8 +543,7 @@ namespace Opm const std::vector& w_ecl = eclipse_state_->getSchedule()->getWells(step); const WellMap& wmap = SimFIBODetails::mapWells(w_ecl); - const std::vector& resv_prod = - SimFIBODetails::resvProducers(*wells, step, wmap); + const std::vector& resv_prod = SimFIBODetails::resvProducers(wells, step, wmap); if (! resv_prod.empty()) { const PhaseUsage& pu = props_.phaseUsage(); diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index 677a4aece..18e370ff2 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -62,12 +62,14 @@ namespace Opm return; } + const int nw = wells->number_of_wells; + if( nw == 0 ) return ; + // We use the WellState::init() function to do bhp and well rates init. // The alternative would be to copy that function wholesale. BaseType :: init(wells, state); // Initialize perfphaserates_, which must be done here. - const int nw = wells->number_of_wells; const int np = wells->number_of_phases; const int nperf = wells->well_connpos[nw]; perfphaserates_.resize(nperf * np, 0.0);