From cb897b07d08ab806eccc7fa67f0232627d64e570 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 20 Oct 2016 20:10:07 +0200 Subject: [PATCH] Adding VREP injection support. As part of it, adding a function to calculate reservoir voidage rate. --- opm/autodiff/BlackoilModelBase.hpp | 9 ++ opm/autodiff/BlackoilModelBase_impl.hpp | 116 +++++++++++++++++++++++- 2 files changed, 124 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 24ba46ded..36b10d9d2 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -312,6 +312,14 @@ namespace Opm { computeFluidInPlace(const ReservoirState& x, const std::vector& fipnum); + /// Function to compute the resevoir voidage for the production wells. + /// TODO: it is just prototyping, and not sure where is the best place to + /// put this function yet. + void computeWellVoidageRates(const ReservoirState& reservoir_state, + const WellState& well_state, + std::vector& well_voidage_rates, + std::vector& voidage_conversion_coeffs); + protected: // --------- Types and enums --------- @@ -427,6 +435,7 @@ namespace Opm { IterationReport solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, + const ReservoirState& reservoir_state, SolutionState& state, WellState& well_state); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index f950dd36e..00d20453c 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -822,7 +822,7 @@ namespace detail { asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && initial_assembly) { // solve the well equations as a pre-processing step - iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); + iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state); } V aliveWells; std::vector cq_s; @@ -837,6 +837,7 @@ namespace detail { asImpl().makeConstantState(state0); asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state); } + return iter_report; } @@ -1020,6 +1021,7 @@ namespace detail { BlackoilModelBase:: solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, + const ReservoirState& reservoir_state, SolutionState& state, WellState& well_state) { @@ -1060,6 +1062,14 @@ namespace detail { asImpl().wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_); converged = getWellConvergence(it); + // Enforce the VREP control + if (it == 0 && asImpl().wellModel().wellCollection()->havingVREPGroups()) { + std::vector well_voidage_rates; + std::vector voidage_conversion_coeffs; + computeWellVoidageRates(reservoir_state, well_state, well_voidage_rates, voidage_conversion_coeffs); + asImpl().wellModel().wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs); + } + // When the well targets are just updated or need to be updated, we need at least one more iteration. if (asImpl().wellModel().wellCollection()->justUpdateWellTargets()) { converged = false; @@ -1091,6 +1101,14 @@ namespace detail { const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); assert(dx.size() == total_residual_v.size()); asImpl().wellModel().updateWellState(dx.array(), dpMaxRel(), well_state); + + // Enforce the VREP control + if (asImpl().wellModel().wellCollection()->havingVREPGroups()) { + std::vector well_voidage_rates; + std::vector voidage_conversion_coeffs; + computeWellVoidageRates(reservoir_state, well_state, well_voidage_rates, voidage_conversion_coeffs); + asImpl().wellModel().wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs); + } } // We have to update the well controls regardless whether there are local // wells active or not as parallel logging will take place that needs to @@ -1511,6 +1529,13 @@ namespace detail { asImpl().wellModel().updateWellState(dwells, dpMaxRel(), well_state); + if (asImpl().wellModel().wellCollection()->havingVREPGroups()) { + std::vector well_voidage_rates; + std::vector voidage_conversion_coeffs; + computeWellVoidageRates(reservoir_state, well_state, well_voidage_rates, voidage_conversion_coeffs); + asImpl().wellModel().wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs); + } + // Update phase conditions used for property calculations. updatePhaseCondFromPrimalVariable(reservoir_state); } @@ -2569,6 +2594,95 @@ namespace detail { return values; } + + + + + template + void + BlackoilModelBase:: + computeWellVoidageRates(const ReservoirState& reservoir_state, + const WellState& well_state, + std::vector& well_voidage_rates, + std::vector& voidage_conversion_coeffs) + { + // TODO: for now, we store the voidage rates for all the production wells. + // For injection wells, the rates are stored as zero. + // Later, more delicate model will be implemented here. + // And for the moment, group control can only work for serial running. + const int nw = well_state.numWells(); + const int np = numPhases(); + + const Wells* wells = asImpl().wellModel().wellsPointer(); + + // we calculate the voidage rate for each well, that means the sum of all the phases. + well_voidage_rates.resize(nw, 0); + // store the conversion coefficients, while only for the use of injection wells. + voidage_conversion_coeffs.resize(nw * np, 1.0); + + int global_number_wells = nw; + +#if HAVE_MPI + if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const auto& info = + boost::any_cast(linsolver_.parallelInformation()); + global_number_wells = info.communicator().sum(global_number_wells); + if ( global_number_wells ) + { + // At least one process has resv wells. Therefore rate converter needs + // to calculate averages over regions that might cross process + // borders. This needs to be done by all processes and therefore + // outside of the next if statement. + rate_converter_.defineState(reservoir_state, boost::any_cast(linsolver_.parallelInformation())); + } + } + else +#endif + { + if ( global_number_wells ) + { + rate_converter_.defineState(reservoir_state); + } + } + + std::vector well_rates(np); + std::vector convert_coeff(np); + + + if ( !well_voidage_rates.empty() ) { + for (int w = 0; w < nw; ++w) { + const bool is_producer = wells->type[w] == PRODUCER; + + // not sure necessary to change all the value to be positive + if (is_producer) { + std::transform(well_state.wellRates().begin() + np * w, + well_state.wellRates().begin() + np * (w + 1), + well_rates.begin(), std::negate()); + + const int fipreg = 0; // Not considering FIP for the moment. + // We will need convert_coeff later actually. + // They should all be the same, right? + rate_converter_.calcCoeff(well_rates, fipreg, convert_coeff); + well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(), + convert_coeff.begin(), 0.0); + } else { + // TODO: it is possible we should use the distribution coeffs from + // the well controls. + // It will be problem if the rates are all zero here. + // It also raises the question where we should call this function. + std::copy(well_state.wellRates().begin() + np * w, + well_state.wellRates().begin() + np * (w + 1), + well_rates.begin()); + const int fipreg = 0; // Not considering FIP for the moment. + rate_converter_.calcCoeff(well_rates, fipreg, convert_coeff); + std::copy(convert_coeff.begin(), convert_coeff.end(), + voidage_conversion_coeffs.begin() + np * w); + } + } + } + } + } // namespace Opm #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED