From 4bdf74eb11ec23afc007a8cab4206ab39013a104 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 23 Aug 2016 09:45:37 +0200 Subject: [PATCH] Reorganize the wellModel - All well-related stuff is moved to StandardWellsDense.hpp - Some unused functions are removed from BlackoilModelEbos.hpp --- opm/autodiff/BlackoilDetails.hpp | 102 + opm/autodiff/BlackoilModelEbos.hpp | 1039 +------- .../SimulatorFullyImplicitBlackoilEbos.hpp | 10 +- opm/autodiff/StandardWellsDense.hpp | 1854 +++++++++++++- opm/autodiff/StandardWellsDense_impl.hpp | 2215 ----------------- 5 files changed, 1888 insertions(+), 3332 deletions(-) diff --git a/opm/autodiff/BlackoilDetails.hpp b/opm/autodiff/BlackoilDetails.hpp index c1c36d26e..7f72af998 100644 --- a/opm/autodiff/BlackoilDetails.hpp +++ b/opm/autodiff/BlackoilDetails.hpp @@ -188,6 +188,108 @@ namespace detail { return product; } } + /// \brief Compute the reduction within the convergence check. + /// \param[in] B A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. B.col(i) contains the values + /// for phase i. + /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. tempV.col(i) contains the + /// values + /// for phase i. + /// \param[in] R A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. B.col(i) contains the values + /// for phase i. + /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum + /// of R for the phase i. + /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the + /// maximum of tempV for the phase i. + /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average + /// of B for the phase i. + /// \param[out] maxNormWell The maximum of the well flux equations for each phase. + /// \param[in] nc The number of cells of the local grid. + /// \return The total pore volume over all cells. + inline + double + convergenceReduction(const Eigen::Array& B, + const Eigen::Array& tempV, + const Eigen::Array& R, + std::vector& R_sum, + std::vector& maxCoeff, + std::vector& B_avg, + std::vector& maxNormWell, + int nc, + int np, + const std::vector pv, + LinearisedBlackoilResidual residual) + { + const int nw = residual.well_flux_eq.size() / np; + assert(nw * np == int(residual.well_flux_eq.size())); + + // Do the global reductions +#if HAVE_MPI + if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast(linsolver_.parallelInformation()); + + // Compute the global number of cells and porevolume + std::vector v(nc, 1); + auto nc_and_pv = std::tuple(0, 0.0); + auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), + Opm::Reduction::makeGlobalSumFunctor()); + auto nc_and_pv_containers = std::make_tuple(v, pv); + info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv); + + for ( int idx = 0; idx < np; ++idx ) + { + auto values = std::tuple(0.0 ,0.0 ,0.0); + auto containers = std::make_tuple(B.col(idx), + tempV.col(idx), + R.col(idx)); + auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), + Opm::Reduction::makeGlobalMaxFunctor(), + Opm::Reduction::makeGlobalSumFunctor()); + info.computeReduction(containers, operators, values); + B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv); + maxCoeff[idx] = std::get<1>(values); + R_sum[idx] = std::get<2>(values); + assert(np >= np); + if (idx < np) { + maxNormWell[idx] = 0.0; + for ( int w = 0; w < nw; ++w ) { + maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w])); + } + } + } + info.communicator().max(maxNormWell.data(), np); + // Compute pore volume + return std::get<1>(nc_and_pv); + } + else +#endif + { + B_avg.resize(np); + maxCoeff.resize(np); + R_sum.resize(np); + maxNormWell.resize(np); + for ( int idx = 0; idx < np; ++idx ) + { + B_avg[idx] = B.col(idx).sum()/nc; + maxCoeff[idx] = tempV.col(idx).maxCoeff(); + R_sum[idx] = R.col(idx).sum(); + + assert(np >= np); + if (idx < np) { + maxNormWell[idx] = 0.0; + for ( int w = 0; w < nw; ++w ) { + maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual.well_flux_eq.value()[nw*idx + w])); + } + } + } + // Compute total pore volume + return std::accumulate(pv.begin(), pv.end(), 0.0); + } + } /// \brief Compute the L-infinity norm of a vector representing a well equation. /// \param a The container to compute the infinity norm on. diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 1145a2320..c4a4b2ee5 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -138,7 +138,7 @@ namespace Opm { const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , const RockCompressibility* rock_comp_props, - const StandardWellsDense& well_model, + const StandardWellsDense& well_model, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output) : ebosSimulator_(ebosSimulator) @@ -172,8 +172,6 @@ namespace Opm { , terminal_output_ (terminal_output) , current_relaxation_(1.0) , isBeginReportStep_(false) - , wellVariables_(wells().number_of_wells * wells().number_of_phases) - , F0_(wells().number_of_wells * wells().number_of_phases) { const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); @@ -283,574 +281,16 @@ namespace Opm { { using namespace Opm::AutoDiffGrid; - //wellModel().assembleAndUpdate(well_state); - // Possibly switch well controls and updating well state to - // get reasonable initial conditions for the wells - wellModel().resetWellControlFromState(well_state); - wellModel().updateWellControls(well_state); - - // Set the primary variables for the wells - setWellVariables(well_state); - - //SolutionState state(/*numPhases=*/3); - //setupLegacyState(state, reservoir_state, well_state); - // -------- Mass balance equations -------- assembleMassBalanceEq(timer, iterationIdx, reservoir_state); // -------- Well equations ---------- - - if (iterationIdx == 0) { - // Create the (constant, derivativeless) initial state. - // Compute initial accumulation contributions - // and well connection pressures. - computeWellConnectionPressures(well_state); - - computeAccumWells(); - // Create the (constant, derivativeless) initial state. - //SolutionState state0 = state; - //makeConstantState(state0); - // Compute initial accumulation contributions - // and well connection pressures. - //computeWellConnectionPressures(state0, well_state); - //wellModel().computeAccumWells(state0); - } - - IterationReport iter_report = {false, false, 0, 0}; - if ( ! wellsActive() ) { - return iter_report; - } - - double dt = timer.currentStepLength(); - std::vector mob_perfcells; - std::vector b_perfcells; - //wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); - if (param_.solve_welleq_initially_ && iterationIdx == 0) { - // solve the well equations as a pre-processing step - iter_report = solveWellEq(mob_perfcells, b_perfcells, dt, well_state); - } - std::vector cq_s; - computeWellFluxDense(cq_s, 4); - updatePerfPhaseRatesAndPressures(cq_s, well_state); + IterationReport iter_report = wellModel().assemble(ebosSimulator_, iterationIdx, dt, well_state, residual_); - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - for (int p = 0; p < np; ++p) { - for (int w = 0; w < nw; ++w) { - //std::cout << p << " " << w << " " << getQs(w,p).value< Eval; - EvalWell extendEval(Eval in) const { - EvalWell out = 0.0; - out.value = in.value; - for(int i = 0;i<3;++i) { - out.derivatives[i] = in.derivatives[flowToEbosPvIdx(i)]; - } - return out; - } - - void print(EvalWell in) const { - std::cout << in.value << std::endl; - for (int i = 0; i < in.derivatives.size(); ++i) { - std::cout << in.derivatives[i] << std::endl; - } - } - - void - setWellVariables(const WellState& xw) { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - for (int w = 0; w < nw; ++w) { - wellVariables_[w + nw*phaseIdx] = 0.0; - wellVariables_[w + nw*phaseIdx].value = xw.wellSolutions()[w + nw* phaseIdx]; - wellVariables_[w + nw*phaseIdx].derivatives[np + phaseIdx] = 1.0; - } - } - } - - void - computeAccumWells() { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - for (int w = 0; w < nw; ++w) { - F0_[w + nw * phaseIdx] = wellVolumeFraction(w,phaseIdx).value; - } - } - } - - - template - void - updatePerfPhaseRatesAndPressures(const std::vector& cq_s, - WellState& xw) const - { - if ( !localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } - - // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - - V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); - for (int phase = 1; phase < np; ++phase) { - cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); - } - xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); - - // Update the perforation pressures. - const V& cdp = wellModel().wellPerforationPressureDiffs(); - for (int w = 0; w < nw; ++w ) { - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - xw.perfPress()[perf] = cdp[perf] + xw.bhp()[w]; - } - } - - } - - void - addWellFluxEq(std::vector cq_s, - const double dt, - const int numBlocks) - { - if( !localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - double volume = 0.002831684659200; // 0.1 cu ft; - //std::vector F = wellVolumeFractions(state); - //std::cout << F0_[0] << std::endl; - //std::cout << F[0] << std::endl; - //std::cout << "før Ebos" < res_vec(nw); - for (int w = 0; w < nw; ++w) { - - EvalWell res = (wellVolumeFraction(w, p) - F0_[w + nw*p]) * volume / dt; - res += getQs(w, p); - //for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - // res -= cq_s[perf*np + p]; - //} - res_vec[w] = res; - - } - ADB tmp = convertToADBWell(res_vec, numBlocks); - qs += superset(tmp,Span(nw,1,p*nw), nw*np); - } - //std::cout << residual_.well_flux_eq << std::endl; - - - //wellModel().convertToADB(res_vec, well_cells, nc, well_id, nw, numBlocks); - //ADB qs = state.qs; - for (int phase = 0; phase < np; ++phase) { - qs -= superset(wellModel().wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); - //qs += superset((F[phase]-F0_[phase]) * vol_dt, Span(nw,1,phase*nw), nw*np); - } - - residual_.well_flux_eq = qs; - - //std::cout << "etter Ebos" << residual_.well_flux_eq << std::endl; - - } - const AutoDiffBlock convertToADBWell(const std::vector& local, const int numVars) const - { - typedef typename ADB::M M; - const int nLocal = local.size(); - typename ADB::V value( nLocal ); - //const int numVars = 5; - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - Eigen::SparseMatrix matFlux(nLocal,np*nw); - - for( int i=0; i jacs( numVars ); - if (numVars == 4) { - for( int d=0; d(deri[d]); - //jacs[ d ] = M(mat[d]); - } - - jacs[3] = M(matFlux); - //jacs[4] = M(matBHP); - } - else if (numVars == 1) { - jacs[0] = M(matFlux); - //jacs[1] = M(matBHP); - } - //std::cout << numVars << std::endl; - - return ADB::function( std::move( value ), std::move( jacs )); - } - - - void - computeWellFluxDense(std::vector& cq_s, const int numBlocks) const - { - if( ! localWellsActive() ) return ; - 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(wellModel().wells().WI, nperf); - const std::vector& well_cells = wellModel().wellOps().well_cells; - std::vector well_id(nperf); - std::vector> cq_s_dense(np, std::vector(nperf,0.0)); - - - // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = wellModel().wellPerforationPressureDiffs(); - - //std::vector> cq_s_dense(np, std::vector(nperf,0.0)); - //std::vector cmix_s_ADB = wellModel().wellVolumeFractions(state); - - for (int w = 0; w < nw; ++w) { - - EvalWell bhp = getBhp(w); //wellModel().extractDenseADWell(state.bhp,w); -// std::cout << "well " << w << std::endl; -// std::cout << "bhpF " << std::endl; -// print(bhp); -// std::cout << "state.bhp " << std::endl; -// print(wellModel().extractDenseADWell(state.bhp,w)); - - - // TODO: fix for 2-phase case - std::vector cmix_s(np,0.0); - for (int phase = 0; phase < np; ++phase) { - //int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - cmix_s[phase] = wellVolumeFraction(w,phase); - } - - //std::cout <<"cmix gas "<< w<< " "< b_perfcells_dense(np, 0.0); - std::vector mob_perfcells_dense(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); - mob_perfcells_dense[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); - } - - // Pressure drawdown (also used to determine direction of flow) - EvalWell well_pressure = bhp + cdp[perf]; - EvalWell drawdown = pressure - well_pressure; - - // injection perforations - if ( drawdown.value > 0 ) { - - //Do nothing if crossflow is not allowed - if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) - continue; - // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - const EvalWell cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); - cq_ps[phase] = b_perfcells_dense[phase] * cq_p; - } - - - if ((active_)[Oil] && (active_)[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const EvalWell cq_psOil = cq_ps[oilpos]; - const EvalWell cq_psGas = cq_ps[gaspos]; - cq_ps[gaspos] += rs * cq_psOil; - cq_ps[oilpos] += rv * cq_psGas; - } - - - // map to ADB - for (int phase = 0; phase < np; ++phase) { - cq_s_dense[phase][perf] = cq_ps[phase]; - - } - - } else { - //Do nothing if crossflow is not allowed - if (!wells().allow_cf[w] && wells().type[w] == PRODUCER) - continue; - - // Using total mobilities - EvalWell total_mob_dense = mob_perfcells_dense[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob_dense += mob_perfcells_dense[phase]; - } - // injection perforations total volume rates - const EvalWell cqt_i = - Tw[perf] * (total_mob_dense * drawdown); - - // compute volume ratio between connection at standard conditions - EvalWell volumeRatio = 0.0; - if ((active_)[Water]) { - const int watpos = pu.phase_pos[Water]; - volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; - } - - if ((active_)[Oil] && (active_)[Gas]) { - EvalWell well_temperature = extendEval(fs.temperature(FluidSystem::oilPhaseIdx)); - EvalWell rsSatEval = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); - EvalWell rvSatEval = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); - - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - EvalWell rvPerf = 0.0; - if (cmix_s[gaspos] > 0) - rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; - - if (rvPerf.value > rvSatEval.value) { - rvPerf = rvSatEval; - //rvPerf.value = rvSatEval.value; - } - - EvalWell rsPerf = 0.0; - if (cmix_s[oilpos] > 0) - rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; - - if (rsPerf.value > rsSatEval.value) { - //rsPerf = 0.0; - rsPerf= rsSatEval; - } - - // Incorporate RS/RV factors if both oil and gas active - const EvalWell d = 1.0 - rvPerf * rsPerf; - - const EvalWell tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; - //std::cout << "tmp_oil " < cq_s2; - //Vector aliveWells; - //computeWellFlux(state,mob_perfcells,b_perfcells, aliveWells,cq_s2); - - //for (int phase = 0; phase < np; ++phase) { - //if( !(((cq_s[phase].value() - cq_s2[phase].value()).abs()<1e-10).all()) ) { - //std::cout << "phase " << phase << std::endl; - //std::cout << cq_s2[phase].value() << std::endl; - //std::cout << cq_s[phase].value() << std::endl; - //} - //} - } - void - computeWellConnectionPressures(const WellState& xw) - { - if( ! localWellsActive() ) return ; - // 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. - std::vector b_perf; - std::vector rsmax_perf; - std::vector rvmax_perf; - std::vector surf_dens_perf; - computePropertiesForWellConnectionPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); - const V& pdepth = subset(depth, wellModel().wellOps().well_cells); - const int nperf = wells().well_connpos[wells().number_of_wells]; - const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); - - wellModel().computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity); - - } - - template - void - computePropertiesForWellConnectionPressures(const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf) - { - const int nperf = wells().well_connpos[wells().number_of_wells]; - const int nw = wells().number_of_wells; - const std::vector& well_cells = wellModel().wellOps().well_cells; - const PhaseUsage& pu = fluid_.phaseUsage(); - const int np = fluid_.numPhases(); - b_perf.resize(nperf*np); - rsmax_perf.resize(nperf); - rvmax_perf.resize(nperf); - - std::vector perf_cond(nperf); - for (int perf = 0; perf < nperf; ++perf) { - perf_cond[perf] = phaseCondition_[well_cells[perf]]; - } - - // Compute the average pressure in each well block - const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); - //V avg_press = perf_press*0; - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const int cell_idx = well_cells[perf]; - const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& fs = intQuants.fluidState(); - - const double p_above = perf == wells().well_connpos[w] ? xw.bhp()[w] : perf_press[perf - 1]; - const double p_avg = (perf_press[perf] + p_above)/2; - double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value; - - if (pu.phase_used[BlackoilPhases::Aqua]) { - b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * pu.num_phases] = - FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - - - if (pu.phase_used[BlackoilPhases::Vapour]) { - int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; - if (perf_cond[perf].hasFreeOil()) { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - else { - double rv = fs.Rv().value; - b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); - } - } - - if (pu.phase_used[BlackoilPhases::Liquid]) { - int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; - if (perf_cond[perf].hasFreeGas()) { - b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - else { - double rs = fs.Rs().value; - b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); - } - } - - if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { - rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); - rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - - } - } - - -// // Use cell values for the temperature as the wells don't knows its temperature yet. -// const ADB perf_temp = subset(state.temperature, well_cells); - -// // Compute b, rsmax, rvmax values for perforations. -// // Evaluate the properties using average well block pressures -// // and cell values for rs, rv, phase condition and temperature. -// const ADB avg_press_ad = ADB::constant(avg_press); -// // const std::vector& pc = phaseCondition(); - -// DataBlock b(nperf, pu.num_phases); -// if (pu.phase_used[BlackoilPhases::Aqua]) { -// const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); -// b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; -// } -// assert((*active_)[Oil]); -// const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); -// if (pu.phase_used[BlackoilPhases::Liquid]) { -// const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null(); -// const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); -// b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; -// } -// if (pu.phase_used[BlackoilPhases::Vapour]) { -// const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null(); -// const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); -// b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; -// } -// if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { -// const V rssat = fluid_.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); -// //rsmax_perf.assign(rssat.data(), rssat.data() + nperf); - -// const V rvsat = fluid_.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); -// //rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); -// } - -// // b is row major, so can just copy data. -// //b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); - - // Surface density. - // The compute density segment wants the surface densities as - // an np * number of wells cells array - V rho = superset(fluid_.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); - for (int phase = 1; phase < pu.num_phases; ++phase) { - rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); - } - surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); - - } - - /// \brief Compute the residual norms of the mass balance for each phase, @@ -1188,7 +628,7 @@ namespace Opm { } - wellModel().updateWellState(dwells, dpMaxRel(), well_state); + wellModel().updateWellState(dwells, well_state); // Update phase conditions used for property calculations. updatePhaseCondFromPrimalVariable(reservoir_state); @@ -1250,9 +690,10 @@ namespace Opm { tempV.col(idx) = R.col(idx).abs()/pv; } - const double pvSum = convergenceReduction(B, tempV, R, + std::vector pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); + const double pvSum = detail::convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, - nc); + nc, np, pv_vector, residual_); std::vector CNV(np); std::vector mass_balance_residual(np); @@ -1430,7 +871,7 @@ namespace Opm { std::vector phaseCondition_; // Well Model - StandardWellsDense well_model_; + StandardWellsDense well_model_; V isRs_; V isRv_; @@ -1447,113 +888,16 @@ namespace Opm { double current_relaxation_; V dx_old_; - std::vector wellVariables_; - std::vector F0_; + // --------- Protected methods --------- public: - EvalWell getBhp(const int wellIdx) const { - const WellControls* wc = wells().ctrls[wellIdx]; - if (well_controls_get_current_type(wc) == BHP) { - EvalWell bhp = 0.0; - const double target_rate = well_controls_get_current_target(wc); - bhp.value = target_rate; - return bhp; - } - return wellVariables_[wellIdx]; - } - - EvalWell getQs(const int wellIdx, const int phaseIdx) const { - EvalWell qs = 0.0; - const WellControls* wc = wells().ctrls[wellIdx]; - const int np = fluid_.numPhases(); - const double target_rate = well_controls_get_current_target(wc); -// std::cout << "well info " << std::endl; -// std::cout << wellIdx << " " << wells().type[wellIdx] << " " < g = {1,1,0.01}; - return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]); - } - - - /// return the StandardWells object - StandardWellsDense& wellModel() { return well_model_; } - const StandardWellsDense& wellModel() const { return well_model_; } + StandardWellsDense& wellModel() { return well_model_; } + const StandardWellsDense& wellModel() const { return well_model_; } /// return the Well struct in the StandardWells const Wells& wells() const { return well_model_.wells(); } @@ -1564,87 +908,6 @@ namespace Opm { /// return true if wells are available on this process bool localWellsActive() const { return well_model_.localWellsActive(); } - void - makeConstantState(SolutionState& state) const - { - // HACK: throw away the derivatives. this may not be the most - // performant way to do things, but it will make the state - // automatically consistent with variableState() (and doing - // things automatically is all the rage in this module ;) - state.pressure = ADB::constant(state.pressure.value()); - state.temperature = ADB::constant(state.temperature.value()); - state.rs = ADB::constant(state.rs.value()); - state.rv = ADB::constant(state.rv.value()); - const int num_phases = state.saturation.size(); - for (int phaseIdx = 0; phaseIdx < num_phases; ++ phaseIdx) { - state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value()); - } - state.qs = ADB::constant(state.qs.value()); - state.bhp = ADB::constant(state.bhp.value()); - assert(state.canonical_phase_pressures.size() == static_cast(Opm::BlackoilPhases::MaxNumPhases)); - for (int canphase = 0; canphase < Opm::BlackoilPhases::MaxNumPhases; ++canphase) { - ADB& pp = state.canonical_phase_pressures[canphase]; - pp = ADB::constant(pp.value()); - } - state.wellVariables = ADB::constant(state.wellVariables.value()); - } - - void setupLegacyState(SolutionState& state, - const ReservoirState& x, - const WellState& xw) const - { - const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int np = x.numPhases(); - - std::vector vars0(np, V::Ones(nc, 1)); - wellModel().variableWellStateInitials(xw, vars0); - std::vector vars = ADB::variables(vars0); - - std::vector indices = {{Pressure, Sw, Xvar}}; - int foo = indices.size(); - indices.resize(4); - wellModel().variableStateWellIndices(indices, foo); - - const ADB& ones = ADB::constant(V::Ones(nc, 1)); - - - // temperature cannot be a variable at this time (only constant). - state.temperature = ones; - - // saturations - state.saturation[Water] = std::move(vars[indices[Sw]]); - - const ADB& sw = state.saturation[Water]; - const ADB& xvar = vars[indices[Xvar]]; - state.saturation[Gas] = xvar; - state.saturation[Oil] = sw + xvar; - - // pressures - state.pressure = std::move(vars[indices[Pressure]]); - const ADB& po = state.pressure; - - const ADB& tmp = po + sw + xvar; - state.canonical_phase_pressures[Gas] = tmp; - state.canonical_phase_pressures[Water] = tmp; - state.canonical_phase_pressures[Oil] = po; - - if (has_disgas_) { - state.rs = po + xvar; - } else { - state.rs = po; - } - if (has_vapoil_) { - state.rv = po + xvar; - } else { - state.rv = po; - } - - // Note that so is never a primary variable. - //state.saturation[Oil] = std::move(so); - - // wells - wellModel().variableStateExtractWellsVars(indices, vars, state, xw); - } V fluidRsSat(const V& p, @@ -1767,11 +1030,7 @@ namespace Opm { return phaseToComp[ phaseIdx ]; } - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const - { - const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; - return flowToEbos[ phaseIdx ]; - } + private: void convertResults(const Simulator& simulator) @@ -1878,6 +1137,12 @@ namespace Opm { } } + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; + return flowToEbos[ phaseIdx ]; + } + void updateLegacyState(const Simulator& simulator, SolutionState& legacyState) { const int numPhases = 3; @@ -2157,7 +1422,6 @@ namespace Opm { prevEpisodeIdx = ebosSimulator_.episodeIndex(); convertResults(ebosSimulator_); - //updateLegacyState(ebosSimulator_, state); if (param_.update_equations_scaling_) { updateEquationsScaling(); @@ -2165,173 +1429,13 @@ namespace Opm { } - IterationReport solveWellEq(const std::vector& mob_perfcells, - const std::vector& b_perfcells, - const double dt, - WellState& well_state) - { - 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]; - std::vector cq_s(np, ADB::null()); - WellState well_state0 = well_state; - - int it = 0; - bool converged; - do { - // bhp and Q for the wells - computeWellFluxDense(cq_s, 1); - updatePerfPhaseRatesAndPressures(cq_s, well_state); - addWellFluxEq(cq_s, dt, 1); - converged = getWellConvergence(it); - - if (converged) { - break; - } - - ++it; - if( localWellsActive() ) - { - std::vector eqs; - eqs.reserve(1); - 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(); - typedef Eigen::SparseMatrix Sp; - Sp Jn0; - Jn[0].toSparse(Jn0); - const Eigen::SparseLU< Sp > solver(Jn0); - ADB::V total_residual_v = total_residual.value(); - const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); - assert(dx.size() == total_residual_v.size()); - wellModel().updateWellState(dx.array(), dpMaxRel(), well_state); - wellModel().updateWellControls(well_state); - setWellVariables(well_state); - } - } while (it < 15); - - if (!converged) { - well_state = well_state0; - } - - const bool failed = false; // Not needed in this method. - const int linear_iters = 0; // Not needed in this method - return IterationReport{failed, converged, linear_iters, it}; - } - void - addWellContributionToMassBalanceEq(const std::vector& cq_s) - { - if ( !localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } - - // Add well contributions to mass balance equations - const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int np = numPhases(); - for (int phase = 0; phase < np; ++phase) { - residual_.material_balance_eq[phase] -= superset(cq_s[phase], wellModel().wellOps().well_cells, nc); - } - } - bool getWellConvergence(const int iteration) - { - const double tol_wells = param_.tolerance_wells_; - const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int np = numPhases(); - const V& pv = geo_.poreVolume(); - std::vector R_sum(np); - std::vector B_avg(np); - std::vector maxCoeff(np); - std::vector maxNormWell(np); - Eigen::Array B(nc, np); - Eigen::Array R(nc, np); - Eigen::Array tempV(nc, np); - for ( int idx = 0; idx < np; ++idx ) - { - V b(nc); - for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { - const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& fs = intQuants.fluidState(); - - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); - - b[cell_idx] = 1 / fs.invB(ebosPhaseIdx).value; - } - B.col(idx) = b; - } - - for ( int idx = 0; idx < np; ++idx ) - { - //const ADB& tempB = rq_[idx].b; - //B.col(idx) = 1./tempB.value(); - R.col(idx) = residual_.material_balance_eq[idx].value(); - tempV.col(idx) = R.col(idx).abs()/pv; - } - - convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc); - - std::vector well_flux_residual(np); - bool converged_Well = true; - // Finish computation - for ( int idx = 0; idx < np; ++idx ) - { - well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx]; - converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); - } - - const double residualWell = detail::infinityNormWell(residual_.well_eq, - linsolver_.parallelInformation()); - converged_Well = converged_Well && (residualWell < Opm::unit::barsa); - const bool converged = converged_Well; - - // if one of the residuals is NaN, throw exception, so that the solver can be restarted - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); - - if (std::isnan(well_flux_residual[phaseIdx])) { - OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName); - } - if (well_flux_residual[phaseIdx] > maxResidualAllowed()) { - OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName); - } - } - - if ( terminal_output_ ) - { - // Only rank 0 does print to std::cout - if (iteration == 0) { - std::string msg; - msg = "Iter"; - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); - msg += " W-FLUX(" + phaseName + ")"; - } - OpmLog::note(msg); - } - std::ostringstream ss; - const std::streamsize oprec = ss.precision(3); - const std::ios::fmtflags oflags = ss.setf(std::ios::scientific); - ss << std::setw(4) << iteration; - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - ss << std::setw(11) << well_flux_residual[phaseIdx]; - } - ss.precision(oprec); - ss.flags(oflags); - OpmLog::note(ss.str()); - } - return converged; - } std::vector computePressures(const ADB& po, @@ -2432,115 +1536,6 @@ namespace Opm { } - // TODO: added since the interfaces of the function are different - // TODO: for StandardWells and MultisegmentWells -// void -// computeWellConnectionPressures(const SolutionState& state, -// const WellState& well_state) -// { -// wellModel().computeWellConnectionPressures(state, well_state); -// } - - - /// \brief Compute the reduction within the convergence check. - /// \param[in] B A matrix with MaxNumPhases columns and the same number rows - /// as the number of cells of the grid. B.col(i) contains the values - /// for phase i. - /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows - /// as the number of cells of the grid. tempV.col(i) contains the - /// values - /// for phase i. - /// \param[in] R A matrix with MaxNumPhases columns and the same number rows - /// as the number of cells of the grid. B.col(i) contains the values - /// for phase i. - /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum - /// of R for the phase i. - /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the - /// maximum of tempV for the phase i. - /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average - /// of B for the phase i. - /// \param[out] maxNormWell The maximum of the well flux equations for each phase. - /// \param[in] nc The number of cells of the local grid. - /// \return The total pore volume over all cells. - double - convergenceReduction(const Eigen::Array& B, - const Eigen::Array& tempV, - const Eigen::Array& R, - std::vector& R_sum, - std::vector& maxCoeff, - std::vector& B_avg, - std::vector& maxNormWell, - int nc) const - { - const int np = numPhases(); - const int nw = residual_.well_flux_eq.size() / np; - assert(nw * np == int(residual_.well_flux_eq.size())); - - // Do the global reductions -#if HAVE_MPI - if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& info = - boost::any_cast(linsolver_.parallelInformation()); - - // Compute the global number of cells and porevolume - std::vector v(nc, 1); - auto nc_and_pv = std::tuple(0, 0.0); - auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), - Opm::Reduction::makeGlobalSumFunctor()); - auto nc_and_pv_containers = std::make_tuple(v, geo_.poreVolume()); - info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv); - - for ( int idx = 0; idx < np; ++idx ) - { - auto values = std::tuple(0.0 ,0.0 ,0.0); - auto containers = std::make_tuple(B.col(idx), - tempV.col(idx), - R.col(idx)); - auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), - Opm::Reduction::makeGlobalMaxFunctor(), - Opm::Reduction::makeGlobalSumFunctor()); - info.computeReduction(containers, operators, values); - B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv); - maxCoeff[idx] = std::get<1>(values); - R_sum[idx] = std::get<2>(values); - assert(np >= np); - if (idx < np) { - maxNormWell[idx] = 0.0; - for ( int w = 0; w < nw; ++w ) { - maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w])); - } - } - } - info.communicator().max(maxNormWell.data(), np); - // Compute pore volume - return std::get<1>(nc_and_pv); - } - else -#endif - { - B_avg.resize(np); - maxCoeff.resize(np); - R_sum.resize(np); - maxNormWell.resize(np); - for ( int idx = 0; idx < np; ++idx ) - { - B_avg[idx] = B.col(idx).sum()/nc; - maxCoeff[idx] = tempV.col(idx).maxCoeff(); - R_sum[idx] = R.col(idx).sum(); - - assert(np >= np); - if (idx < np) { - maxNormWell[idx] = 0.0; - for ( int w = 0; w < nw; ++w ) { - maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w])); - } - } - } - // Compute total pore volume - return geo_.poreVolume().sum(); - } - } double dpMaxRel() const { return param_.dp_max_rel_; } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index e52b573e3..98cfc7fb4 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -33,7 +33,7 @@ namespace Opm { class SimulatorFullyImplicitBlackoilEbos; -class StandardWellsDense; +//class StandardWellsDense; /// a simulator for the blackoil model class SimulatorFullyImplicitBlackoilEbos @@ -42,13 +42,16 @@ public: typedef typename TTAG(EclFlowProblem) TypeTag; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + typedef WellStateFullyImplicitBlackoil WellState; typedef BlackoilState ReservoirState; typedef BlackoilOutputWriter OutputWriter; typedef BlackoilModelEbos Model; typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; - typedef StandardWellsDense WellModel; + typedef StandardWellsDense WellModel; /// Initialise from parameters and objects to observe. @@ -229,7 +232,8 @@ public: // Run a multiple steps of the solver depending on the time step control. solver_timer.start(); - const WellModel well_model(wells); + const std::vector pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); + const WellModel well_model(wells, model_param_, terminal_output_, pv); auto solver = createSolver(well_model); diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 71c9052a7..d1b2a1f41 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -40,6 +40,16 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include + + + #include #include @@ -48,10 +58,42 @@ namespace Opm { /// Class for handling the standard well model. + template class StandardWellsDense { public: struct WellOps { - explicit WellOps(const Wells* wells); + WellOps(const Wells* wells) + : w2p(), + p2w(), + well_cells() + { + if( wells ) + { + w2p = Eigen::SparseMatrix(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); + p2w = Eigen::SparseMatrix(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); + + const int nw = wells->number_of_wells; + const int* const wpos = wells->well_connpos; + + typedef Eigen::Triplet Tri; + + 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()); + + well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]); + } + } Eigen::SparseMatrix w2p; // well -> perf (scatter) Eigen::SparseMatrix p2w; // perf -> well (gather) std::vector well_cells; // the set of perforated cells @@ -61,9 +103,13 @@ namespace Opm { using ADB = AutoDiffBlock; typedef DenseAd::Evaluation EvalWell; + typedef WellStateFullyImplicitBlackoil WellState; + typedef BlackoilModelParameters ModelParameters; + //typedef AutoDiffBlock ADB; using Vector = ADB::V; using V = ADB::V; + typedef ADB::M M; // copied from BlackoilModelBase @@ -73,18 +119,92 @@ namespace Opm { Eigen::Dynamic, Eigen::RowMajor>; // --------- Public methods --------- - explicit StandardWellsDense(const Wells* wells_arg); + StandardWellsDense(const Wells* wells_arg, + const ModelParameters& param, + const bool terminal_output, + const std::vector& pv) + : wells_active_(wells_arg!=nullptr) + , wells_(wells_arg) + , wops_(wells_arg) + , fluid_(nullptr) + , active_(nullptr) + , phase_condition_(nullptr) + , vfp_properties_(nullptr) + , well_perforation_densities_(Vector()) + , well_perforation_pressure_diffs_(Vector()) + , store_well_perforation_fluxes_(false) + , wellVariables_(wells_arg->number_of_wells * wells_arg->number_of_phases) + , F0_(wells_arg->number_of_wells * wells_arg->number_of_phases) + , param_(param) + , terminal_output_(terminal_output) + , pv_(pv) + { + } + + + template + IterationReport assemble(const Simulator& ebosSimulator, + const int iterationIdx, + const double dt, + WellState& well_state, + LinearisedBlackoilResidual& residual) { + resetWellControlFromState(well_state); + updateWellControls(well_state); + // Set the primary variables for the wells + setWellVariables(well_state); + + if (iterationIdx == 0) { + computeWellConnectionPressures(ebosSimulator, well_state); + computeAccumWells(); + } + + IterationReport iter_report = {false, false, 0, 0}; + if ( ! wellsActive() ) { + return iter_report; + } + //wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); + if (param_.solve_welleq_initially_ && iterationIdx == 0) { + // solve the well equations as a pre-processing step + iter_report = solveWellEq(ebosSimulator, dt, well_state); + } + std::vector cq_s; + computeWellFluxDense(ebosSimulator, cq_s, 4); + updatePerfPhaseRatesAndPressures(cq_s, well_state); + addWellFluxEq(cq_s, dt, 4, residual); + addWellContributionToMassBalanceEq(cq_s, residual); + + if (param_.compute_well_potentials_) { + //wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state); + } + return iter_report; + } + + void init(const BlackoilPropsAdInterface* fluid_arg, const std::vector* active_arg, const std::vector* pc_arg, const VFPProperties* vfp_properties_arg, const double gravity_arg, - const Vector& depth_arg); + const Vector& depth_arg) + { + fluid_ = fluid_arg; + active_ = active_arg; + phase_condition_ = pc_arg; + vfp_properties_ = vfp_properties_arg; + gravity_ = gravity_arg; + perf_cell_depth_ = subset(depth_arg, wellOps().well_cells); + } - const WellOps& wellOps() const; + const WellOps& wellOps() const + { + return wops_; + } - int numPhases() const { return wells().number_of_phases; }; + + int numPhases() const { return wells().number_of_phases; } + + int numCells() const { return pv_.size();} template void resetWellControlFromState(WellState xw) { @@ -95,82 +215,1206 @@ namespace Opm { } } - const Wells& wells() const; + const Wells& wells() const + { + assert(wells_ != 0); + return *(wells_); + } - const Wells* wellsPointer() const; + const Wells* wellsPointer() const + { + return wells_; + } /// return true if wells are available in the reservoir - bool wellsActive() const; - void setWellsActive(const bool wells_active); + bool wellsActive() const + { + return wells_active_; + } + void setWellsActive(const bool wells_active) + { + wells_active_ = wells_active; + } /// return true if wells are available on this process - bool localWellsActive() const; + bool localWellsActive() const + { + return wells_ ? (wells_->number_of_wells > 0 ) : false; + } - int numWellVars() const; + int numWellVars() const + { + if ( !localWellsActive() ) + { + return 0; + } + + // For each well, we have a bhp variable, and one flux per phase. + const int nw = wells().number_of_wells; + return (numPhases() + 1) * nw; + } /// Density of each well perforation - Vector& wellPerforationDensities(); // mutable version kept for BlackoilMultisegmentModel - const Vector& wellPerforationDensities() const; + Vector& wellPerforationDensities() // mutable version kept for BlackoilMultisegmentModel + { + return well_perforation_densities_; + } + + const Vector& wellPerforationDensities() const { + return well_perforation_densities_; + } /// Diff to bhp for each well perforation. - Vector& wellPerforationPressureDiffs(); // mutable version kept for BlackoilMultisegmentModel - const Vector& wellPerforationPressureDiffs() const; + Vector& wellPerforationPressureDiffs() { // mutable version kept for BlackoilMultisegmentModel + return well_perforation_pressure_diffs_; + } + const Vector& wellPerforationPressureDiffs() const + { + return well_perforation_pressure_diffs_; + } - template - void extractWellPerfProperties(const SolutionState& state, - const std::vector& rq, - std::vector& mob_perfcells, - std::vector& b_perfcells) const; + typedef DenseAd::Evaluation Eval; + EvalWell extendEval(Eval in) const { + EvalWell out = 0.0; + out.value = in.value; + for(int i = 0;i<3;++i) { + out.derivatives[i] = in.derivatives[flowToEbosPvIdx(i)]; + } + return out; + } - template - void computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - Vector& aliveWells, - std::vector& cq_s) const; - - template void - computeWellFluxDense(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - std::vector& cq_s) const; + setWellVariables(const WellState& xw) { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + for (int w = 0; w < nw; ++w) { + wellVariables_[w + nw*phaseIdx] = 0.0; + wellVariables_[w + nw*phaseIdx].value = xw.wellSolutions()[w + nw* phaseIdx]; + wellVariables_[w + nw*phaseIdx].derivatives[np + phaseIdx] = 1.0; + } + } + } - EvalWell extractDenseAD(const ADB& data, int i, int j) const; - EvalWell extractDenseADWell(const ADB& data, int i) const; + void print(EvalWell in) const { + std::cout << in.value << std::endl; + for (int i = 0; i < in.derivatives.size(); ++i) { + std::cout << in.derivatives[i] << std::endl; + } + } - const ADB convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const; + void + computeAccumWells() { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + for (int w = 0; w < nw; ++w) { + F0_[w + nw * phaseIdx] = wellVolumeFraction(w,phaseIdx).value; + } + } + } + + EvalWell extractDenseAD(const ADB& data, int i, int j) const + { + EvalWell output = 0.0; + output.value = data.value()[i]; + const int np = wells().number_of_phases; + const std::vector& jac = data.derivative(); + //std::cout << jac.size() << std::endl; + int numblocs = jac.size(); + for (int b = 0; b < numblocs; ++b) { + if (b < np) { // don't copy well blocks) + //std::cout << jac[b].coeff(i,j) << std::endl; + output.derivatives[b] = jac[b].coeff(i,j); + } + } + return output; + } + + EvalWell extractDenseADWell(const ADB& data, int i) const + { + EvalWell output = 0.0; + output.value = data.value()[i]; + const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; + const std::vector& jac = data.derivative(); + //std::cout << jac.size() << std::endl; + int numblocs = jac.size(); + for (int b = 0; b < np; ++b) { + output.derivatives[b+np] = jac[numblocs-1].coeff(i, b*nw + i); + } + return output; + } + + const ADB convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const + { + typedef typename ADB::M M; + const int nLocal = local.size(); + typename ADB::V value( nLocal ); + //const int numVars = 5; + const int np = wells().number_of_phases; + + std::vector> mat(np, Eigen::SparseMatrix(nLocal,nc)); + Eigen::SparseMatrix matFlux(nLocal,np*nw); + Eigen::SparseMatrix matBHP(nLocal,nw); + + for( int i=0; i jacs( numVars ); + if (numVars == 4) { + for( int d=0; d(deri[d]); + jacs[ d ] = M(mat[d]); + } + + jacs[3] = M(matFlux); + //jacs[4] = M(matBHP); + } + else if (numVars == 1) { + jacs[0] = M(matFlux); + //jacs[1] = M(matBHP); + } + //std::cout << numVars << std::endl; + + return ADB::function( std::move( value ), std::move( jacs )); + } - template + template void updatePerfPhaseRatesAndPressures(const std::vector& cq_s, - const SolutionState& state, - WellState& xw) const; + WellState& xw) const + { + if ( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + + V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); + } + xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); + + // Update the perforation pressures. + const V& cdp = wellPerforationPressureDiffs(); + for (int w = 0; w < nw; ++w ) { + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + xw.perfPress()[perf] = cdp[perf] + xw.bhp()[w]; + } + } + + } + + void + addWellFluxEq(std::vector cq_s, + const double dt, + const int numBlocks, + LinearisedBlackoilResidual& residual) + { + if( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + double volume = 0.002831684659200; // 0.1 cu ft; + //std::vector F = wellVolumeFractions(state); + //std::cout << F0_[0] << std::endl; + //std::cout << F[0] << std::endl; + //std::cout << "før Ebos" < res_vec(nw); + for (int w = 0; w < nw; ++w) { + + EvalWell res = (wellVolumeFraction(w, p) - F0_[w + nw*p]) * volume / dt; + res += getQs(w, p); + //for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + // res -= cq_s[perf*np + p]; + //} + res_vec[w] = res; + + } + ADB tmp = convertToADBWell(res_vec, numBlocks); + qs += superset(tmp,Span(nw,1,p*nw), nw*np); + } + //std::cout << residual_.well_flux_eq << std::endl; + + + //wellModel().convertToADB(res_vec, well_cells, nc, well_id, nw, numBlocks); + //ADB qs = state.qs; + for (int phase = 0; phase < np; ++phase) { + qs -= superset(wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); + //qs += superset((F[phase]-F0_[phase]) * vol_dt, Span(nw,1,phase*nw), nw*np); + } + + residual.well_flux_eq = qs; + + //std::cout << "etter Ebos" << residual_.well_flux_eq << std::endl; + + } + const AutoDiffBlock convertToADBWell(const std::vector& local, const int numVars) const + { + typedef typename ADB::M M; + const int nLocal = local.size(); + typename ADB::V value( nLocal ); + //const int numVars = 5; + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + Eigen::SparseMatrix matFlux(nLocal,np*nw); + + for( int i=0; i jacs( numVars ); + if (numVars == 4) { + for( int d=0; d(deri[d]); + //jacs[ d ] = M(mat[d]); + } + + jacs[3] = M(matFlux); + //jacs[4] = M(matBHP); + } + else if (numVars == 1) { + jacs[0] = M(matFlux); + //jacs[1] = M(matBHP); + } + //std::cout << numVars << std::endl; + + return ADB::function( std::move( value ), std::move( jacs )); + } + + + template + void + computeWellFluxDense(const Simulator& ebosSimulator, + std::vector& cq_s, + const int numBlocks) const + { + if( ! localWellsActive() ) return ; + 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 = wellOps().well_cells; + std::vector well_id(nperf); + std::vector> cq_s_dense(np, std::vector(nperf,0.0)); + + + // pressure diffs computed already (once per step, not changing per iteration) + const V& cdp = wellPerforationPressureDiffs(); + + for (int w = 0; w < nw; ++w) { + + EvalWell bhp = getBhp(w); + + // TODO: fix for 2-phase case + std::vector cmix_s(np,0.0); + for (int phase = 0; phase < np; ++phase) { + //int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + cmix_s[phase] = wellVolumeFraction(w,phase); + } + + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + const int cell_idx = well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + well_id[perf] = w; + EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + EvalWell rs = extendEval(fs.Rs()); + EvalWell rv = extendEval(fs.Rv()); + std::vector b_perfcells_dense(np, 0.0); + std::vector mob_perfcells_dense(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); + mob_perfcells_dense[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); + } + + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + cdp[perf]; + EvalWell drawdown = pressure - well_pressure; + + // injection perforations + if ( drawdown.value > 0 ) { + + //Do nothing if crossflow is not allowed + if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) + continue; + // compute phase volumetric rates at standard conditions + std::vector cq_ps(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + const EvalWell cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); + cq_ps[phase] = b_perfcells_dense[phase] * cq_p; + } + + + if ((*active_)[Oil] && (*active_)[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const EvalWell cq_psOil = cq_ps[oilpos]; + const EvalWell cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs * cq_psOil; + cq_ps[oilpos] += rv * cq_psGas; + } + + + // map to ADB + for (int phase = 0; phase < np; ++phase) { + cq_s_dense[phase][perf] = cq_ps[phase]; + + } + + } else { + //Do nothing if crossflow is not allowed + if (!wells().allow_cf[w] && wells().type[w] == PRODUCER) + continue; + + // Using total mobilities + EvalWell total_mob_dense = mob_perfcells_dense[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob_dense += mob_perfcells_dense[phase]; + } + // injection perforations total volume rates + const EvalWell cqt_i = - Tw[perf] * (total_mob_dense * drawdown); + + // compute volume ratio between connection at standard conditions + EvalWell volumeRatio = 0.0; + if ((*active_)[Water]) { + const int watpos = pu.phase_pos[Water]; + volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; + } + + if ((*active_)[Oil] && (*active_)[Gas]) { + EvalWell well_temperature = extendEval(fs.temperature(FluidSystem::oilPhaseIdx)); + EvalWell rsSatEval = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); + EvalWell rvSatEval = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); + + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + EvalWell rvPerf = 0.0; + if (cmix_s[gaspos] > 0) + rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; + + if (rvPerf.value > rvSatEval.value) { + rvPerf = rvSatEval; + //rvPerf.value = rvSatEval.value; + } + + EvalWell rsPerf = 0.0; + if (cmix_s[oilpos] > 0) + rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; + + if (rsPerf.value > rsSatEval.value) { + //rsPerf = 0.0; + rsPerf= rsSatEval; + } + + // Incorporate RS/RV factors if both oil and gas active + const EvalWell d = 1.0 - rvPerf * rsPerf; + + const EvalWell tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; + //std::cout << "tmp_oil " < + IterationReport solveWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state) + { + const int np = wells().number_of_phases; + std::vector cq_s(np, ADB::null()); + WellState well_state0 = well_state; + + LinearisedBlackoilResidual residual( { std::vector(3, ADB::null()), + ADB::null(), + ADB::null(), + { 1.1169, 1.0031, 0.0031 }, // the default magic numbers + false } ); + + int it = 0; + bool converged; + do { + // bhp and Q for the wells + computeWellFluxDense(ebosSimulator, cq_s, 1); + updatePerfPhaseRatesAndPressures(cq_s, well_state); + addWellFluxEq(cq_s, dt, 1, residual); + converged = getWellConvergence(ebosSimulator, it, residual); + + if (converged) { + break; + } + + ++it; + if( localWellsActive() ) + { + std::vector eqs; + eqs.reserve(1); + 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(); + typedef Eigen::SparseMatrix Sp; + Sp Jn0; + Jn[0].toSparse(Jn0); + const Eigen::SparseLU< Sp > solver(Jn0); + ADB::V total_residual_v = total_residual.value(); + const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); + assert(dx.size() == total_residual_v.size()); + updateWellState(dx.array(), well_state); + updateWellControls(well_state); + setWellVariables(well_state); + } + } while (it < 15); + + if (!converged) { + well_state = well_state0; + } + + const bool failed = false; // Not needed in this method. + const int linear_iters = 0; // Not needed in this method + return IterationReport{failed, converged, linear_iters, it}; + } + + + template + bool getWellConvergence(Simulator& ebosSimulator, + const int iteration, + const LinearisedBlackoilResidual& residual) + { + const int np = numPhases(); + const int nc = numCells(); + const double tol_wells = param_.tolerance_wells_; + const double maxResidualAllowed = param_.max_residual_allowed_; + + std::vector R_sum(np); + std::vector B_avg(np); + std::vector maxCoeff(np); + std::vector maxNormWell(np); + Eigen::Array B(nc, np); + Eigen::Array R(nc, np); + Eigen::Array tempV(nc, np); + + for ( int idx = 0; idx < np; ++idx ) + { + V b(nc); + for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + + b[cell_idx] = 1 / fs.invB(ebosPhaseIdx).value; + } + B.col(idx) = b; + } + + detail::convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc, np, pv_, residual); + + std::vector well_flux_residual(np); + bool converged_Well = true; + // Finish computation + for ( int idx = 0; idx < np; ++idx ) + { + well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx]; + converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); + } + + // if one of the residuals is NaN, throw exception, so that the solver can be restarted + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + + if (std::isnan(well_flux_residual[phaseIdx])) { + OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName); + } + if (well_flux_residual[phaseIdx] > maxResidualAllowed) { + OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName); + } + } + + if ( terminal_output_ ) + { + // Only rank 0 does print to std::cout + if (iteration == 0) { + std::string msg; + msg = "Iter"; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + msg += " W-FLUX(" + phaseName + ")"; + } + OpmLog::note(msg); + } + std::ostringstream ss; + const std::streamsize oprec = ss.precision(3); + const std::ios::fmtflags oflags = ss.setf(std::ios::scientific); + ss << std::setw(4) << iteration; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + ss << std::setw(11) << well_flux_residual[phaseIdx]; + } + ss.precision(oprec); + ss.flags(oflags); + OpmLog::note(ss.str()); + } + return converged_Well; + } + + + + template + void + computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw) + { + if( ! localWellsActive() ) return ; + // 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. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(ebosSimulator, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const V& pdepth = perf_cell_depth_; + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_); + + } + + template + void + computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + const std::vector& well_cells = wellOps().well_cells; + const PhaseUsage& pu = fluid_->phaseUsage(); + const int np = fluid_->numPhases(); + b_perf.resize(nperf*np); + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + + std::vector perf_cond(nperf); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = (*phase_condition_)[well_cells[perf]]; + } + + // Compute the average pressure in each well block + const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); + //V avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const int cell_idx = well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + const double p_above = perf == wells().well_connpos[w] ? xw.bhp()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * pu.num_phases] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; + if (perf_cond[perf].hasFreeOil()) { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + else { + double rv = fs.Rv().value; + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; + if (perf_cond[perf].hasFreeGas()) { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + else { + double rs = fs.Rs().value; + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { + rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + } + } + + // Surface density. + // The compute density segment wants the surface densities as + // an np * number of wells cells array + V rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); + + } + + + void + addWellContributionToMassBalanceEq(const std::vector& cq_s, + LinearisedBlackoilResidual& residual) + { + if ( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + // Add well contributions to mass balance equations + const int np = numPhases(); + for (int phase = 0; phase < np; ++phase) { + residual.material_balance_eq[phase] -= superset(cq_s[phase], wellOps().well_cells, numCells()); + } + } template void updateWellState(const Vector& dwells, - const double dpmaxrel, - WellState& well_state); + WellState& well_state) + { + if( localWellsActive() ) + { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + // Extract parts of dwells corresponding to each part. + int varstart = 0; + const Vector dxvar_well = subset(dwells, Span(np*nw, 1, varstart)); + //const Vector dqs = subset(dwells, Span(np*nw, 1, varstart)); + varstart += dxvar_well.size(); + //const Vector dbhp = subset(dwells, Span(nw, 1, varstart)); + //varstart += dbhp.size(); + assert(varstart == dwells.size()); + + // 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 Vector xvar_well_old = Eigen::Map(&well_state.wellSolutions()[0], nw*np); + + double dFLimit = 0.2; + double dBHPLimit = 2; + double dTotalRateLimit = 0.5; + //std::cout << "dxvar_well "< F(np,0.0); + const int sign2 = dxvar_well[nw + w] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dxvar_well[nw + w]),dFLimit); + well_state.wellSolutions()[nw + w] = xvar_well_old[nw + w] - dx2_limited; + const int sign3 = dxvar_well[2*nw + w] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dxvar_well[2*nw + w]),dFLimit); + well_state.wellSolutions()[2*nw + w] = xvar_well_old[2*nw + w] - dx3_limited; + F[Water] = well_state.wellSolutions()[nw + w]; + F[Gas] = well_state.wellSolutions()[2*nw + w]; + F[Oil] = 1.0 - F[Water] - F[Gas]; + + // const double dFw = dxvar_well[nw + w]; + // const double dFg = dxvar_well[nw*2 + w]; + // const double dFo = - dFw - dFg; + // //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; + // double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit; + // step = std::min(step, 1.0); + // //std::cout << step << std::endl; + // F[Water] = xvar_well_old[nw + w] - step*dFw; + // F[Gas] = xvar_well_old[2*nw + w] - step*dFg; + // F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo; + + if (F[Water] < 0.0) { + F[Gas] /= (1.0 - F[Water]); + F[Oil] /= (1.0 - F[Water]); + F[Water] = 0.0; + } + if (F[Gas] < 0.0) { + F[Water] /= (1.0 - F[Gas]); + F[Oil] /= (1.0 - F[Gas]); + F[Gas] = 0.0; + } + if (F[Oil] < 0.0) { + F[Water] /= (1.0 - F[Oil]); + F[Gas] /= (1.0 - F[Oil]); + F[Oil] = 0.0; + } + well_state.wellSolutions()[nw + w] = F[Water]; + well_state.wellSolutions()[2*nw + w] = F[Gas]; + + //std::cout << wells().name[w] << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; + + + std::vector g = {1,1,0.01}; + + if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { + for (int p = 0; p < np; ++p) { + F[p] /= distr[p]; + } + } else { + for (int p = 0; p < np; ++p) { + F[p] /= g[p]; + } + } + //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; + + + + + + + // const double dFw = dxvar_well[nw + w]; + // const double dFg = dxvar_well[nw*2 + w]; + // const double dFo = - dFw - dFg; + //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; + // double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit; + // step = std::min(step, 1.0); + // std::cout << step << std::endl; + // F[Water] = xvar_well_old[nw + w] - step*dFw; + // F[Gas] = xvar_well_old[2*nw + w] - step*dFg; + // F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo; + // double sumF = F[Water]+F[Gas]+F[Oil]; + // F[Water] /= sumF; + // F[Gas] /= sumF; + // F[Oil] /= sumF; + // well_state.wellSolutions()[nw + w] = F[Water]; + // well_state.wellSolutions()[2 * nw + w] = F[Gas]; + + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + //const int sign1 = dxvar_well[w] > 0 ? 1: -1; + //const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dTotalRateLimit); + well_state.wellSolutions()[w] = xvar_well_old[w] - dxvar_well[w]; + + switch (wells().type[w]) { + case INJECTOR: + for (int p = 0; p < np; ++p) { + const double comp_frac = wells().comp_frac[np*w + p]; + //if (comp_frac > 0) { + well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[w]; + //} + + } + break; + case PRODUCER: + for (int p = 0; p < np; ++p) { + well_state.wellRates()[w*np + p] = well_state.wellSolutions()[w] * F[p]; + } + break; + } + } + break; + case SURFACE_RATE: + { + const int sign1 = dxvar_well[w] > 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dBHPLimit); + well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited; + //const int sign = (dxvar_well1[w] < 0) ? -1 : 1; + //well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ; + well_state.bhp()[w] = well_state.wellSolutions()[w]; + + if (wells().type[w]==PRODUCER) { + + double F_target = 0.0; + for (int p = 0; p < np; ++p) { + F_target += wells().comp_frac[np*w + p] * F[p]; + } + for (int p = 0; p < np; ++p) { + //std::cout << F[p] << std::endl; + well_state.wellRates()[np*w + p] = F[p] * target_rate /F_target; + } + } else { + + for (int p = 0; p < np; ++p) { + //std::cout << wells().comp_frac[np*w + p] << " " < 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dBHPLimit); + well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited; + //const int sign = (dxvar_well1[w] < 0) ? -1 : 1; + //well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ; + well_state.bhp()[w] = well_state.wellSolutions()[w]; + for (int p = 0; p < np; ++p) { + well_state.wellRates()[np*w + p] = F[p] * target_rate; + } + } + break; + + } + } + + + + + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + //Loop over all wells + #pragma omp parallel for schedule(static) + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + const int nwc = well_controls_get_num(wc); + //Loop over all controls until we find a THP control + //that specifies what we need... + //Will only update THP for wells with THP control + for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(wc, ctrl_index) == THP) { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + if ((*active_)[ Water ]) { + aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + double alq = well_controls_iget_alq(wc, ctrl_index); + int table_id = well_controls_iget_vfp(wc, ctrl_index); + + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp); + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp, alq); + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + + //Assume only one THP control specified for each well + break; + } + } + } + } + } template - void updateWellControls(WellState& xw); + void updateWellControls(WellState& xw) + { + if( !localWellsActive() ) return ; - // TODO: should LinearisedBlackoilResidual also be a template class? - template - void addWellFluxEq(const std::vector& cq_s, - const SolutionState& state, - const double dt, - LinearisedBlackoilResidual& residual); + std::string modestring[4] = { "BHP", "THP", "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; + #pragma omp parallel for schedule(dynamic) + for (int w = 0; w < nw; ++w) { + 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. + int current = xw.currentControls()[w]; + // Loop over all controls except the current one, and also + // skip any RESERVOIR_RATE controls, since we cannot + // handle those. + const int nwc = well_controls_get_num(wc); + int ctrl_index = 0; + for (; ctrl_index < nwc; ++ctrl_index) { + if (ctrl_index == current) { + // This is the currently used control, so it is + // used as an equation. So this is not used as an + // inequality constraint, and therefore skipped. + continue; + } + if (wellhelpers::constraintBroken( + xw.bhp(), xw.thp(), xw.wellRates(), + 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. + // We disregard terminal_ouput here as with it only messages + // for wells on one process will be printed. + std::ostringstream ss; + ss << "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; + OpmLog::info(ss.str()); + xw.currentControls()[w] = ctrl_index; + current = xw.currentControls()[w]; + well_controls_set_current( wc, current); - // TODO: some parameters, like gravity, maybe it is better to put in the member list - template - void addWellControlEq(const SolutionState& state, - const WellState& xw, - const Vector& aliveWells, - LinearisedBlackoilResidual& residual); - template - void computeWellConnectionPressures(const SolutionState& state, - const WellState& xw); + + // Updating well state and primary variables if constraint is broken + + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(wc, current); + const double* distr = well_controls_iget_distr(wc, current); + switch (well_controls_iget_type(wc, current)) { + case BHP: + xw.bhp()[w] = target; + break; + + case THP: { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + if ((*active_)[ Water ]) { + aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + const int vfp = well_controls_iget_vfp(wc, current); + const double& thp = well_controls_iget_target(wc, current); + const double& alq = well_controls_iget_alq(wc, current); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + break; + } + + case RESERVOIR_RATE: + // No direct change to any observable quantity at + // surface condition. In this case, use existing + // flow rates as initial conditions as reservoir + // rate acts only in aggregate. + break; + + case SURFACE_RATE: + // assign target value as initial guess for injectors and + // single phase producers (orat, grat, wrat) + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + for (int phase = 0; phase < np; ++phase) { + const double& compi = wells().comp_frac[np * w + phase]; + //if (compi > 0.0) { + xw.wellRates()[np*w + phase] = target * compi; + //} + } + } else if (well_type == PRODUCER) { + + // only set target as initial rates for single phase + // producers. (orat, grat and wrat, and not lrat) + // lrat will result in numPhasesWithTargetsUnderThisControl == 2 + int numPhasesWithTargetsUnderThisControl = 0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + numPhasesWithTargetsUnderThisControl += 1; + } + } + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { + xw.wellRates()[np*w + phase] = target * distr[phase]; + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + + + break; + } + std::vector g = {1,1,0.01}; + if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { + const double* distr = well_controls_iget_distr(wc, current); + for (int phase = 0; phase < np; ++phase) { + g[phase] = distr[phase]; + } + } + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + const WellType& well_type = wells().type[w]; + xw.wellSolutions()[w] = 0.0; + if (well_type == INJECTOR) { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p]; + } + } else { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[w] += g[p] * xw.wellRates()[np*w + p]; + } + } + } + break; + + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + xw.wellSolutions()[w] = xw.bhp()[w]; + } + break; + } + + double tot_well_rate = 0.0; + for (int p = 0; p < np; ++p) { + tot_well_rate += g[p] * xw.wellRates()[np*w + p]; + } + if(std::abs(tot_well_rate) > 0) { + xw.wellSolutions()[nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; //wells->comp_frac[np*w + Water]; // Water; + xw.wellSolutions()[2*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; //wells->comp_frac[np*w + Gas]; //Gas + } else { + //xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; + //xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; + } + } + } + } + + + + + + + template + void computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw) + { + if( ! localWellsActive() ) return ; + // 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. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(ebosSimulator, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const Vector& pdepth = perf_cell_depth_; + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_); + + } // state0 is non-constant, while it will not be used outside of the function template @@ -178,52 +1422,229 @@ namespace Opm { computeWellPotentials(const std::vector& mob_perfcells, const std::vector& b_perfcells, SolutionState& state0, - WellState& well_state); + WellState& well_state) + { + const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - template - void - variableStateExtractWellsVars(const std::vector& indices, - std::vector& vars, - SolutionState& state, - WellState& well_state) const; + Vector bhps = Vector::Zero(nw); + for (int w = 0; w < nw; ++w) { + const WellControls* ctrl = wells().ctrls[w]; + const int nwc = well_controls_get_num(ctrl); + //Loop over all controls until we find a BHP control + //or a THP control that specifies what we need. + //Pick the value that gives the most restrictive flow + for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(ctrl, ctrl_index) == BHP) { + bhps[w] = well_controls_iget_target(ctrl, ctrl_index); + } - void - variableStateWellIndices(std::vector& indices, - int& next) const; + if (well_controls_iget_type(ctrl, ctrl_index) == THP) { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; - std::vector - variableWellStateIndices() const; + if ((*active_)[ Water ]) { + aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + const int vfp = well_controls_iget_vfp(ctrl, ctrl_index); + const double& thp = well_controls_iget_target(ctrl, ctrl_index); + const double& alq = well_controls_iget_alq(ctrl, ctrl_index); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + // apply the strictest of the bhp controlls i.e. smallest bhp for injectors + if ( bhp < bhps[w]) { + bhps[w] = bhp; + } + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + // apply the strictest of the bhp controlls i.e. largest bhp for producers + if ( bhp > bhps[w]) { + bhps[w] = bhp; + } + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + } + } + + } + + // use bhp limit from control + state0.bhp = ADB::constant(bhps); + + // compute well potentials + Vector aliveWells; + std::vector well_potentials; + computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, well_potentials); + + // store well potentials in the well state + // transform to a single vector instead of separate vectors pr phase + const int nperf = wells().well_connpos[nw]; + Vector cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np); + } + well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np); + } - template - void - variableWellStateInitials(const WellState& xw, - std::vector& vars0) const; /// If set, computeWellFlux() will additionally store the /// total reservoir volume perforation fluxes. - void setStoreWellPerforationFluxesFlag(const bool store_fluxes); + void setStoreWellPerforationFluxesFlag(const bool store_fluxes) + { + store_well_perforation_fluxes_ = store_fluxes; + } /// Retrieves the stored fluxes. It is an error to call this /// unless setStoreWellPerforationFluxesFlag(true) has been /// called. - const Vector& getStoredWellPerforationFluxes() const; + const Vector& getStoredWellPerforationFluxes() const + { + assert(store_well_perforation_fluxes_); + return well_perforation_fluxes_; + } + + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; + return flowToEbos[ phaseIdx ]; + } + + int flowToEbosPvIdx( const int flowPv ) const + { + const int flowToEbos[ 3 ] = { + BlackoilIndices::pressureSwitchIdx, + BlackoilIndices::waterSaturationIdx, + BlackoilIndices::compositionSwitchIdx + }; + return flowToEbos[ flowPv ]; + } /// upate the dynamic lists related to economic limits template void updateListEconLimited(ScheduleConstPtr schedule, const int current_step, - const Wells* wells, + const Wells* wells_struct, const WellState& well_state, - DynamicListEconLimited& list_econ_limited) const; + DynamicListEconLimited& list_econ_limited) const + { + const int nw = wells_struct->number_of_wells; - template - std::vector wellVolumeFractions(const SolutionState& state) const; + for (int w = 0; w < nw; ++w) { + // flag to check if the mim oil/gas rate limit is violated + bool rate_limit_violated = false; + const std::string& well_name = wells_struct->name[w]; + const Well* well_ecl = schedule->getWell(well_name); + const WellEconProductionLimits& econ_production_limits = well_ecl->getEconProductionLimits(current_step); + // economic limits only apply for production wells. + if (wells_struct->type[w] != PRODUCER) { + continue; + } - template - void computeAccumWells(const SolutionState& state); + // if no limit is effective here, then continue to the next well + if ( !econ_production_limits.onAnyEffectiveLimit() ) { + continue; + } + // for the moment, we only handle rate limits, not handling potential limits + // the potential limits should not be difficult to add + const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit(); + if (quantity_limit == WellEcon::POTN) { + const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n") + + std::string("All the limits will be evaluated based on RATE. "); + OpmLog::warning("NOT_SUPPORTING_POTN", msg); + } + + const WellMapType& well_map = well_state.wellMap(); + const typename WellMapType::const_iterator i_well = well_map.find(well_name); + assert(i_well != well_map.end()); // should always be found? + const WellMapEntryType& map_entry = i_well->second; + const int well_number = map_entry[0]; + + if (econ_production_limits.onAnyRateLimit()) { + rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state, well_number); + } + + if (rate_limit_violated) { + if (econ_production_limits.endRun()) { + const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n") + + std::string("the program will keep running after ") + well_name + std::string(" is closed"); + OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message); + } + + if (econ_production_limits.validFollowonWell()) { + OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet"); + } + + if (well_ecl->getAutomaticShutIn()) { + list_econ_limited.addShutWell(well_name); + const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit"); + OpmLog::info(msg); + } else { + list_econ_limited.addStoppedWell(well_name); + const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit"); + OpmLog::info(msg); + } + // the well is closed, not need to check other limits + continue; + } + + // checking for ratio related limits, mostly all kinds of ratio. + bool ratio_limits_violated = false; + RatioCheckTuple ratio_check_return; + + if (econ_production_limits.onAnyRatioLimit()) { + ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state, map_entry); + ratio_limits_violated = std::get<0>(ratio_check_return); + } + + if (ratio_limits_violated) { + const bool last_connection = std::get<1>(ratio_check_return); + const int worst_offending_connection = std::get<2>(ratio_check_return); + + const int perf_start = map_entry[1]; + + assert((worst_offending_connection >= 0) && (worst_offending_connection < map_entry[2])); + + const int cell_worst_offending_connection = wells_struct->well_cells[perf_start + worst_offending_connection]; + list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection); + const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ") + + well_name + std::string(" will be closed due to economic limit"); + OpmLog::info(msg); + + if (last_connection) { + list_econ_limited.addShutWell(well_name); + const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed"); + OpmLog::info(msg2); + } + } + + } + } template void computeWellConnectionDensitesPressures(const WellState& xw, @@ -232,12 +1653,32 @@ namespace Opm { const std::vector& rvmax_perf, const std::vector& surf_dens_perf, const std::vector& depth_perf, - const double grav); + const double grav) { + // Compute densities + std::vector cd = + WellDensitySegmented::computeConnectionDensities( + wells(), xw, fluid_->phaseUsage(), + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const int nperf = wells().well_connpos[wells().number_of_wells]; + + // Compute pressure deltas + std::vector cdp = + WellDensitySegmented::computeConnectionPressureDelta( + wells(), depth_perf, cd, grav); + + // Store the results + well_perforation_densities_ = Eigen::Map(cd.data(), nperf); + well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); + + } protected: bool wells_active_; const Wells* wells_; const WellOps wops_; + ModelParameters param_; + bool terminal_output_; const BlackoilPropsAdInterface* fluid_; const std::vector* active_; @@ -254,24 +1695,137 @@ namespace Opm { bool store_well_perforation_fluxes_; Vector well_perforation_fluxes_; - std::vector F0_; + std::vector wellVariables_; + std::vector F0_; + const std::vector& pv_; // protected methods - template - void computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf); + EvalWell getBhp(const int wellIdx) const { + const WellControls* wc = wells().ctrls[wellIdx]; + if (well_controls_get_current_type(wc) == BHP) { + EvalWell bhp = 0.0; + const double target_rate = well_controls_get_current_target(wc); + bhp.value = target_rate; + return bhp; + } + return wellVariables_[wellIdx]; + } + EvalWell getQs(const int wellIdx, const int phaseIdx) const { + EvalWell qs = 0.0; + const WellControls* wc = wells().ctrls[wellIdx]; + const int np = fluid_->numPhases(); + const double target_rate = well_controls_get_current_target(wc); + + if (wells().type[wellIdx] == INJECTOR) { + const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx]; + if (comp_frac == 0.0) + return qs; + + if (well_controls_get_current_type(wc) == BHP) { + return wellVariables_[wellIdx]; + } + qs.value = target_rate; + return qs; + } + + // Producers + if (well_controls_get_current_type(wc) == BHP) { + return wellVariables_[wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx); + } + if (well_controls_get_current_type(wc) == SURFACE_RATE) { + const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx]; + + if (comp_frac == 1.0) { + qs.value = target_rate; + return qs; + } + int currentControlIdx = 0; + for (int i = 0; i < np; ++i) { + currentControlIdx += wells().comp_frac[np*wellIdx + i] * i; + } + + if (wellVolumeFractionScaled(wellIdx,currentControlIdx) == 0) { + return qs; + } + return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / wellVolumeFractionScaled(wellIdx,currentControlIdx)); + } + // ReservoirRate + return target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx); + } + + EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const { + assert(fluid_.numPhases() == 3); + const int nw = wells().number_of_wells; + if (phaseIdx == Water) { + return wellVariables_[nw + wellIdx]; + } + + if (phaseIdx == Gas) { + return wellVariables_[2*nw + wellIdx]; + } + + // Oil + return 1.0 - wellVariables_[nw + wellIdx] - wellVariables_[2 * nw + wellIdx]; + } + + EvalWell wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const { + const WellControls* wc = wells().ctrls[wellIdx]; + if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + const double* distr = well_controls_get_current_distr(wc); + return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx]; + } + std::vector g = {1,1,0.01}; + return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]); + } template bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, - const int well_number) const; + const int well_number) const + { + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + const int np = well_state.numPhases(); + + if (econ_production_limits.onMinOilRate()) { + assert((*active_)[Oil]); + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double min_oil_rate = econ_production_limits.minOilRate(); + if (std::abs(oil_rate) < min_oil_rate) { + return true; + } + } + + if (econ_production_limits.onMinGasRate() ) { + assert((*active_)[Gas]); + const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ]; + const double min_gas_rate = econ_production_limits.minGasRate(); + if (std::abs(gas_rate) < min_gas_rate) { + return true; + } + } + + if (econ_production_limits.onMinLiquidRate() ) { + assert((*active_)[Oil]); + assert((*active_)[Water]); + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; + const double liquid_rate = oil_rate + water_rate; + const double min_liquid_rate = econ_production_limits.minLiquidRate(); + if (std::abs(liquid_rate) < min_liquid_rate) { + return true; + } + } + + if (econ_production_limits.onMinReservoirFluidRate()) { + OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet"); + } + + return false; + } + using WellMapType = typename WellState::WellMapType; using WellMapEntryType = typename WellState::mapentry_t; @@ -293,12 +1847,128 @@ namespace Opm { template RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, - const WellMapEntryType& map_entry) const; + const WellMapEntryType& map_entry) const + { + // TODO: not sure how to define the worst-offending connection when more than one + // ratio related limit is violated. + // The defintion used here is that we define the violation extent based on the + // ratio between the value and the corresponding limit. + // For each violated limit, we decide the worst-offending connection separately. + // Among the worst-offending connections, we use the one has the biggest violation + // extent. + + bool any_limit_violated = false; + bool last_connection = false; + int worst_offending_connection = INVALIDCONNECTION; + double violation_extent = -1.0; + + if (econ_production_limits.onMaxWaterCut()) { + const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry); + bool water_cut_violated = std::get<0>(water_cut_return); + if (water_cut_violated) { + any_limit_violated = true; + const double violation_extent_water_cut = std::get<3>(water_cut_return); + if (violation_extent_water_cut > violation_extent) { + violation_extent = violation_extent_water_cut; + worst_offending_connection = std::get<2>(water_cut_return); + last_connection = std::get<1>(water_cut_return); + } + } + } + + if (econ_production_limits.onMaxGasOilRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!"); + } + + if (econ_production_limits.onMaxWaterGasRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!"); + } + + if (econ_production_limits.onMaxGasLiquidRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!"); + } + + if (any_limit_violated) { + assert(worst_offending_connection >=0); + assert(violation_extent > 1.); + } + + return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent); + } template RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, - const WellMapEntryType& map_entry) const; + const WellMapEntryType& map_entry) const + { + bool water_cut_limit_violated = false; + int worst_offending_connection = INVALIDCONNECTION; + bool last_connection = false; + double violation_extent = -1.0; + + const int np = well_state.numPhases(); + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + const int well_number = map_entry[0]; + + assert((*active_)[Oil]); + assert((*active_)[Water]); + + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; + const double liquid_rate = oil_rate + water_rate; + double water_cut; + if (std::abs(liquid_rate) != 0.) { + water_cut = water_rate / liquid_rate; + } else { + water_cut = 0.0; + } + + const double max_water_cut_limit = econ_production_limits.maxWaterCut(); + if (water_cut > max_water_cut_limit) { + water_cut_limit_violated = true; + } + + if (water_cut_limit_violated) { + // need to handle the worst_offending_connection + const int perf_start = map_entry[1]; + const int perf_number = map_entry[2]; + + std::vector water_cut_perf(perf_number); + for (int perf = 0; perf < perf_number; ++perf) { + const int i_perf = perf_start + perf; + const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ]; + const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ]; + const double liquid_perf_rate = oil_perf_rate + water_perf_rate; + if (std::abs(liquid_perf_rate) != 0.) { + water_cut_perf[perf] = water_perf_rate / liquid_perf_rate; + } else { + water_cut_perf[perf] = 0.; + } + } + + last_connection = (perf_number == 1); + if (last_connection) { + worst_offending_connection = 0; + violation_extent = water_cut_perf[0] / max_water_cut_limit; + return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + double max_water_cut_perf = 0.; + for (int perf = 0; perf < perf_number; ++perf) { + if (water_cut_perf[perf] > max_water_cut_perf) { + worst_offending_connection = perf; + max_water_cut_perf = water_cut_perf[perf]; + } + } + + assert(max_water_cut_perf != 0.); + assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number)); + + violation_extent = max_water_cut_perf / max_water_cut_limit; + } + + return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); + } }; diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index cedcf233b..2d104747f 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -30,2218 +30,3 @@ #include - - - -namespace Opm -{ - - - StandardWellsDense:: - WellOps::WellOps(const Wells* wells) - : w2p(), - p2w(), - well_cells() - { - if( wells ) - { - w2p = Eigen::SparseMatrix(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); - p2w = Eigen::SparseMatrix(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); - - const int nw = wells->number_of_wells; - const int* const wpos = wells->well_connpos; - - typedef Eigen::Triplet Tri; - - 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()); - - well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]); - } - } - - - - - StandardWellsDense::StandardWellsDense(const Wells* wells_arg) - : wells_active_(wells_arg!=nullptr) - , wells_(wells_arg) - , wops_(wells_arg) - , fluid_(nullptr) - , active_(nullptr) - , phase_condition_(nullptr) - , vfp_properties_(nullptr) - , well_perforation_densities_(Vector()) - , well_perforation_pressure_diffs_(Vector()) - , store_well_perforation_fluxes_(false) - { - } - - - - - - void - StandardWellsDense::init(const BlackoilPropsAdInterface* fluid_arg, - const std::vector* active_arg, - const std::vector* pc_arg, - const VFPProperties* vfp_properties_arg, - const double gravity_arg, - const Vector& depth_arg) - { - fluid_ = fluid_arg; - active_ = active_arg; - phase_condition_ = pc_arg; - vfp_properties_ = vfp_properties_arg; - gravity_ = gravity_arg; - perf_cell_depth_ = subset(depth_arg, wellOps().well_cells);; - } - - - - - - const Wells& StandardWellsDense::wells() const - { - assert(wells_ != 0); - return *(wells_); - } - - - const Wells* StandardWellsDense::wellsPointer() const - { - return wells_; - } - - - - bool StandardWellsDense::wellsActive() const - { - return wells_active_; - } - - - - - - void StandardWellsDense::setWellsActive(const bool wells_active) - { - wells_active_ = wells_active; - } - - - - - - bool StandardWellsDense::localWellsActive() const - { - return wells_ ? (wells_->number_of_wells > 0 ) : false; - } - - - - - - int - StandardWellsDense::numWellVars() const - { - if ( !localWellsActive() ) - { - return 0; - } - - // For each well, we have a bhp variable, and one flux per phase. - const int nw = wells().number_of_wells; - return (numPhases() + 1) * nw; - } - - - - - - const StandardWellsDense::WellOps& - StandardWellsDense::wellOps() const - { - return wops_; - } - - - - - - StandardWellsDense::Vector& StandardWellsDense::wellPerforationDensities() - { - return well_perforation_densities_; - } - - - - - - const StandardWellsDense::Vector& - StandardWellsDense::wellPerforationDensities() const - { - return well_perforation_densities_; - } - - - - - - StandardWellsDense::Vector& - StandardWellsDense::wellPerforationPressureDiffs() - { - return well_perforation_pressure_diffs_; - } - - - - - - const StandardWellsDense::Vector& - StandardWellsDense::wellPerforationPressureDiffs() const - { - return well_perforation_pressure_diffs_; - } - - - - - template - void - StandardWellsDense:: - computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf) - { - const int nperf = wells().well_connpos[wells().number_of_wells]; - const int nw = wells().number_of_wells; - - // Compute the average pressure in each well block - const Vector perf_press = Eigen::Map(xw.perfPress().data(), nperf); - Vector avg_press = perf_press*0; - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; - const double p_avg = (perf_press[perf] + p_above)/2; - avg_press[perf] = p_avg; - } - } - - const std::vector& well_cells = wellOps().well_cells; - - // Use cell values for the temperature as the wells don't knows its temperature yet. - const ADB perf_temp = subset(state.temperature, well_cells); - - // Compute b, rsmax, rvmax values for perforations. - // Evaluate the properties using average well block pressures - // and cell values for rs, rv, phase condition and temperature. - const ADB avg_press_ad = ADB::constant(avg_press); - std::vector perf_cond(nperf); - // const std::vector& pc = phaseCondition(); - for (int perf = 0; perf < nperf; ++perf) { - perf_cond[perf] = (*phase_condition_)[well_cells[perf]]; - } - const PhaseUsage& pu = fluid_->phaseUsage(); - DataBlock b(nperf, pu.num_phases); - if (pu.phase_used[BlackoilPhases::Aqua]) { - const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; - } - assert((*active_)[Oil]); - const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); - if (pu.phase_used[BlackoilPhases::Liquid]) { - const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null(); - const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; - } - if (pu.phase_used[BlackoilPhases::Vapour]) { - const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null(); - const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; - } - if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { - const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); - rsmax_perf.assign(rssat.data(), rssat.data() + nperf); - - const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); - rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); - } - - // b is row major, so can just copy data. - b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); - - // Surface density. - // The compute density segment wants the surface densities as - // an np * number of wells cells array - Vector rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); - for (int phase = 1; phase < pu.num_phases; ++phase) { - rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); - } - surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); - - } - - - - - - template - void - StandardWellsDense:: - computeWellConnectionDensitesPressures(const WellState& xw, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& surf_dens_perf, - const std::vector& depth_perf, - const double grav) - { - // Compute densities - std::vector cd = - WellDensitySegmented::computeConnectionDensities( - wells(), xw, fluid_->phaseUsage(), - b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - const int nperf = wells().well_connpos[wells().number_of_wells]; - - // Compute pressure deltas - std::vector cdp = - WellDensitySegmented::computeConnectionPressureDelta( - wells(), depth_perf, cd, grav); - - // Store the results - well_perforation_densities_ = Eigen::Map(cd.data(), nperf); - well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); - } - - - - - - template - void - StandardWellsDense:: - computeWellConnectionPressures(const SolutionState& state, - const WellState& xw) - { - if( ! localWellsActive() ) return ; - // 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. - std::vector b_perf; - std::vector rsmax_perf; - std::vector rvmax_perf; - std::vector surf_dens_perf; - computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - const Vector& pdepth = perf_cell_depth_; - const int nperf = wells().well_connpos[wells().number_of_wells]; - const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); - - computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_); - - } - - - - - - template - void - StandardWellsDense:: - extractWellPerfProperties(const SolutionState& /* state */, - const std::vector& rq, - std::vector& mob_perfcells, - std::vector& b_perfcells) const - { - // If we have wells, extract the mobilities and b-factors for - // the well-perforated cells. - if ( !localWellsActive() ) { - mob_perfcells.clear(); - b_perfcells.clear(); - return; - } else { - const std::vector& well_cells = wellOps().well_cells; - const int num_phases = wells().number_of_phases; - mob_perfcells.resize(num_phases, ADB::null()); - b_perfcells.resize(num_phases, ADB::null()); - for (int phase = 0; phase < num_phases; ++phase) { - mob_perfcells[phase] = subset(rq[phase].mob, well_cells); - b_perfcells[phase] = subset(rq[phase].b, well_cells); - } - } - } - - template - void - StandardWellsDense:: - computeWellFluxDense(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - std::vector& cq_s) const - { - if( ! localWellsActive() ) return ; - - 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(); - Vector Tw = Eigen::Map(wells().WI, nperf); - const std::vector& well_cells = wellOps().well_cells; - std::vector well_id(nperf); - - // pressure diffs computed already (once per step, not changing per iteration) - const Vector& cdp = wellPerforationPressureDiffs(); - - std::vector> cq_s_dense(np, std::vector(nperf,0.0)); - std::vector cmix_s_ADB = wellVolumeFractions(state); - - const int oilpos = pu.phase_pos[Oil]; - ADB perfpressure = (wellOps().w2p * state.bhp) + cdp; - const ADB rsSat = fluid_->rsSat(perfpressure, cmix_s_ADB[oilpos], well_cells); - const ADB rvSat = fluid_->rvSat(perfpressure, cmix_s_ADB[oilpos], well_cells); - - for (int w = 0; w < nw; ++w) { - - EvalWell bhp = extractDenseADWell(state.bhp,w); - - - // TODO: fix for 2-phase case - std::vector cmix_s(np,0.0); - for (int phase = 0; phase < np; ++phase) { - cmix_s[phase] = extractDenseADWell(cmix_s_ADB[phase],w); - } - - //std::cout <<"cmix gas "<< w<< " "< b_perfcells_dense(np, 0.0); - std::vector mob_perfcells_dense(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - b_perfcells_dense[phase] = extractDenseAD(b_perfcells[phase], perf, cell_idx); - mob_perfcells_dense[phase] = extractDenseAD(mob_perfcells[phase], perf, cell_idx); - - } - - // Pressure drawdown (also used to determine direction of flow) - EvalWell drawdown = pressure - bhp - cdp[perf]; - - // injection perforations - if ( drawdown.value > 0 ) { - - //Do nothing if crossflow is not allowed - if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) - continue; - // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - const EvalWell cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); - cq_ps[phase] = b_perfcells_dense[phase] * cq_p; - } - - - if ((*active_)[Oil] && (*active_)[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const EvalWell cq_psOil = cq_ps[oilpos]; - const EvalWell cq_psGas = cq_ps[gaspos]; - cq_ps[gaspos] += rs * cq_psOil; - cq_ps[oilpos] += rv * cq_psGas; - } - - - // map to ADB - for (int phase = 0; phase < np; ++phase) { - cq_s_dense[phase][perf] = cq_ps[phase]; - } - - } else { - //Do nothing if crossflow is not allowed - if (!wells().allow_cf[w] && wells().type[w] == PRODUCER) - continue; - - // Using total mobilities - EvalWell total_mob_dense = mob_perfcells_dense[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob_dense += mob_perfcells_dense[phase]; - } - // injection perforations total volume rates - const EvalWell cqt_i = - Tw[perf] * (total_mob_dense * drawdown); - - // compute volume ratio between connection at standard conditions - EvalWell volumeRatio = 0.0; - if ((*active_)[Water]) { - const int watpos = pu.phase_pos[Water]; - volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; - } - - if ((*active_)[Oil] && (*active_)[Gas]) { - - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - EvalWell rvPerf = 0.0; - if (cmix_s[gaspos] > 0) - rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; - - if (rvPerf.value > rvSat.value()[w]) { - rvPerf = 0.0; - rvPerf.value = rvSat.value()[w]; - } - - EvalWell rsPerf = 0.0; - if (cmix_s[oilpos] > 0) - rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; - - if (rsPerf.value > rsSat.value()[w]) { - rsPerf = 0.0; - rsPerf.value = rsSat.value()[w]; - } - - // Incorporate RS/RV factors if both oil and gas active - const EvalWell d = 1.0 - rvPerf * rsPerf; - - const EvalWell tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; - //std::cout << "tmp_oil " < cq_s2; - //Vector aliveWells; - //computeWellFlux(state,mob_perfcells,b_perfcells, aliveWells,cq_s2); - - //for (int phase = 0; phase < np; ++phase) { - //if( !(((cq_s[phase].value() - cq_s2[phase].value()).abs()<1e-10).all()) ) { - //std::cout << "phase " << phase << std::endl; - //std::cout << cq_s2[phase].value() << std::endl; - //std::cout << cq_s[phase].value() << std::endl; - //} - //} - } - - - - - - - template - void - StandardWellsDense:: - computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - Vector& aliveWells, - std::vector& cq_s) const - { - if( ! localWellsActive() ) return ; - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - Vector Tw = Eigen::Map(wells().WI, nperf); - const std::vector& well_cells = wellOps().well_cells; - - // pressure diffs computed already (once per step, not changing per iteration) - const Vector& cdp = wellPerforationPressureDiffs(); - // Extract needed quantities for the perforation cells - const ADB& p_perfcells = subset(state.pressure, well_cells); - - // Perforation pressure - const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp; - - // Pressure drawdown (also used to determine direction of flow) - const ADB drawdown = p_perfcells - perfpressure; - - // Compute vectors with zero and ones that - // selects the wanted quantities. - - // selects injection perforations - Vector selectInjectingPerforations = Vector::Zero(nperf); - // selects producing perforations - Vector selectProducingPerforations = Vector::Zero(nperf); - for (int c = 0; c < nperf; ++c){ - if (drawdown.value()[c] < 0) - selectInjectingPerforations[c] = 1; - else - selectProducingPerforations[c] = 1; - } - std::vector distri(np); - for (int p = 0; p < np; ++p) { - distri[p] = V::Zero(nw); - } - - - Vector isResv = Vector::Zero(nw); - for (int w = 0; w < nw; ++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. - const double* distr = well_controls_get_current_distr(wc); - if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { - isResv[w] = 1; - } - for (int p = 0; p < np; ++p) { - distri[p][w] = distr[p]; - //std::cout << "distr " << distr[p] << " comp_frac " << comp_frac << std::endl; - } - - } - - - // Handle cross flow - const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); - const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); - for (int w = 0; w < nw; ++w) { - if (!wells().allow_cf[w]) { - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - // Crossflow is not allowed; reverse flow is prevented. - // At least one of the perforation must be open in order to have a meeningful - // equation to solve. For the special case where all perforations have reverse flow, - // and the target rate is non-zero all of the perforations are keept open. - if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { - selectProducingPerforations[perf] = 0.0; - } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ - selectInjectingPerforations[perf] = 0.0; - } - } - } - } - - // HANDLE FLOW INTO WELLBORE - // compute phase volumetric rates at standard conditions - std::vector cq_p(np, ADB::null()); - std::vector cq_ps(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); - cq_ps[phase] = b_perfcells[phase] * cq_p[phase]; - } - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - if ((*active_)[Oil] && (*active_)[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const ADB cq_psOil = cq_ps[oilpos]; - const ADB cq_psGas = cq_ps[gaspos]; - const ADB& rv_perfcells = subset(state.rv, well_cells); - const ADB& rs_perfcells = subset(state.rs, well_cells); - cq_ps[gaspos] += rs_perfcells * cq_psOil; - cq_ps[oilpos] += rv_perfcells * cq_psGas; - } - - // HANDLE FLOW OUT FROM WELLBORE - // Using total mobilities - ADB total_mob = mob_perfcells[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob += mob_perfcells[phase]; - } - // injection perforations total volume rates - const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); - - // Store well perforation total fluxes (reservor volumes) if requested. - if (store_well_perforation_fluxes_) { - // Ugly const-cast, but unappealing alternatives. - Vector& wf = const_cast(well_perforation_fluxes_); - wf = cqt_i.value(); - for (int phase = 0; phase < np; ++phase) { - wf += cq_p[phase].value(); - } - } - - // compute wellbore mixture for injecting perforations - // The wellbore mixture depends on the inflow from the reservoar - // and the well injection rates. - - // compute avg. and total wellbore phase volumetric rates at standard conds -// const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); -// std::vector wbq(np, ADB::null()); -// ADB wbqt = ADB::constant(Vector::Zero(nw)); -// for (int phase = 0; phase < np; ++phase) { -// const ADB& q_ps = wellOps().p2w * cq_ps[phase]; -// const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); -// Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); -// const int pos = pu.phase_pos[phase]; -// wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps; -// wbqt += wbq[phase]; -// } - // compute wellbore mixture at standard conditions. - - - std::vector cmix_s(np, ADB::null()); - Vector ones = Vector::Constant(nperf,1.0); - cmix_s[Water] = wellOps().w2p * subset(state.wellVariables, Span(nw, 1, 1*nw)); - cmix_s[Gas] = wellOps().w2p * subset(state.wellVariables, Span(nw, 1, 2*nw)); - cmix_s[Oil] = (ones - cmix_s[Water] - cmix_s[Gas]); -// std::cout << cmix_s[Gas].value() << std::endl; - -// std::vector g(np, Vector::Constant(nw,1.0)); -// g[Gas] = Vector::Constant(nw,0.01); -// const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); -// std::vector wbq(np, ADB::null()); -// ADB wbqt = ADB::constant(Vector::Zero(nw)); -// for (int phase = 0; phase < np; ++phase) { -// const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); -// wbq[phase] = q_s * ( (Vector::Constant(nw,1.0) - isResv) * g[phase] + isResv * distri[phase]); -// wbqt += wbq[phase]; -// } -// Selector notDeadWells_selector(wbqt.value(), Selector::Zero); -// std::vector cmix_s(np, ADB::null()); -// for (int phase = 0; phase < np; ++phase) { -// const int pos = pu.phase_pos[phase]; -// cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); -// } - - - // compute volume ratio between connection at standard conditions - ADB volumeRatio = ADB::constant(Vector::Zero(nperf)); - - if ((*active_)[Water]) { - const int watpos = pu.phase_pos[Water]; - volumeRatio += cmix_s[watpos] / b_perfcells[watpos]; - } - - - if ((*active_)[Oil] && (*active_)[Gas]) { - // Incorporate RS/RV factors if both oil and gas active -// const ADB& rv = subset(state.rv, well_cells); -// const ADB& rs = subset(state.rs, well_cells); -// const ADB d = Vector::Constant(nperf,1.0) - rv* rs; - - - - - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - - Selector noGas_selector(cmix_s[gaspos].value(), Selector::GreaterZero); - Selector noOil_selector(cmix_s[oilpos].value(), Selector::GreaterZero); - ADB rv = noGas_selector.select(cmix_s[oilpos]/cmix_s[gaspos], ADB::constant(Vector::Constant(nperf,0.0))); - ADB rs = noOil_selector.select(cmix_s[gaspos]/cmix_s[oilpos], ADB::constant(Vector::Constant(nperf,0.0))); - - const ADB rsSat = fluid_->rsSat(perfpressure, cmix_s[oilpos], well_cells); - const ADB rvSat = fluid_->rvSat(perfpressure, cmix_s[oilpos], well_cells); - - Selector maxRs_selector(rs.value() - rsSat.value(), Selector::GreaterZero); - Selector maxRv_selector(rv.value() - rvSat.value(), Selector::GreaterZero); - rv = maxRv_selector.select(rvSat, rv); - rs = maxRs_selector.select(rsSat, rs); - - const ADB d = Vector::Constant(nperf,1.0) - rv * rs; - - - const ADB tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d; - volumeRatio += tmp_oil / b_perfcells[oilpos]; - - const ADB tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d; - volumeRatio += tmp_gas / b_perfcells[gaspos]; - } - else { - if ((*active_)[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos]; - } - if ((*active_)[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos]; - } - } - - - // injecting connections total volumerates at standard conditions - //std::cout << volumeRatio.value() << std::endl; - ADB cqt_is = cqt_i/volumeRatio; - - // connection phase volumerates at standard conditions - cq_s.resize(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; //*b_perfcells[phase]; - } - - // check for dead wells (used in the well controll equations) - aliveWells = Vector::Constant(nw, 1.0); -// for (int w = 0; w < nw; ++w) { -// if (wbqt.value()[w] == 0) { -// aliveWells[w] = 0.0; -// } -// } - } - - typedef DenseAd::Evaluation EvalWell; - EvalWell - StandardWellsDense:: - extractDenseAD(const ADB& data, int i, int j) const - { - EvalWell output = 0.0; - output.value = data.value()[i]; - const int np = wells().number_of_phases; - const std::vector& jac = data.derivative(); - //std::cout << jac.size() << std::endl; - int numblocs = jac.size(); - for (int b = 0; b < numblocs; ++b) { - if (b < np) { // don't copy well blocks) - //std::cout << jac[b].coeff(i,j) << std::endl; - output.derivatives[b] = jac[b].coeff(i,j); - } - } - return output; - } - - typedef DenseAd::Evaluation EvalWell; - EvalWell - StandardWellsDense:: - extractDenseADWell(const ADB& data, int i) const - { - EvalWell output = 0.0; - output.value = data.value()[i]; - const int nw = wells().number_of_wells; - const int np = wells().number_of_phases; - const std::vector& jac = data.derivative(); - //std::cout << jac.size() << std::endl; - int numblocs = jac.size(); - for (int b = 0; b < np; ++b) { - output.derivatives[b+np] = jac[numblocs-1].coeff(i, b*nw + i); - } - return output; - } - - const AutoDiffBlock StandardWellsDense::convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const - { - typedef typename ADB::M M; - const int nLocal = local.size(); - typename ADB::V value( nLocal ); - //const int numVars = 5; - const int np = wells().number_of_phases; - - std::vector> mat(np, Eigen::SparseMatrix(nLocal,nc)); - Eigen::SparseMatrix matFlux(nLocal,np*nw); - Eigen::SparseMatrix matBHP(nLocal,nw); - - for( int i=0; i jacs( numVars ); - if (numVars == 4) { - for( int d=0; d(deri[d]); - jacs[ d ] = M(mat[d]); - } - - jacs[3] = M(matFlux); - //jacs[4] = M(matBHP); - } - else if (numVars == 1) { - jacs[0] = M(matFlux); - //jacs[1] = M(matBHP); - } - //std::cout << numVars << std::endl; - - return ADB::function( std::move( value ), std::move( jacs )); - } - - - - - template - void - StandardWellsDense:: - updatePerfPhaseRatesAndPressures(const std::vector& cq_s, - const SolutionState& state, - WellState& xw) const - { - if ( !localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } - - // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - - Vector cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); - for (int phase = 1; phase < np; ++phase) { - cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); - } - xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); - - // Update the perforation pressures. - const Vector& cdp = wellPerforationPressureDiffs(); - const Vector perfpressure = (wellOps().w2p * state.bhp.value().matrix()).array() + cdp; - xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); - } - - - - - - - - template - void - StandardWellsDense:: - updateWellState(const Vector& dwells, - const double dpmaxrel, - WellState& well_state) - { - if( localWellsActive() ) - { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - // Extract parts of dwells corresponding to each part. - int varstart = 0; - const Vector dxvar_well = subset(dwells, Span(np*nw, 1, varstart)); - //const Vector dqs = subset(dwells, Span(np*nw, 1, varstart)); - varstart += dxvar_well.size(); - //const Vector dbhp = subset(dwells, Span(nw, 1, varstart)); - //varstart += dbhp.size(); - assert(varstart == dwells.size()); - - // 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 Vector xvar_well_old = Eigen::Map(&well_state.wellSolutions()[0], nw*np); - - double dFLimit = 0.2; - double dBHPLimit = 2; - double dTotalRateLimit = 0.5; - //std::cout << "dxvar_well "< F(np,0.0); - const int sign2 = dxvar_well[nw + w] > 0 ? 1: -1; - const double dx2_limited = sign2 * std::min(std::abs(dxvar_well[nw + w]),dFLimit); - well_state.wellSolutions()[nw + w] = xvar_well_old[nw + w] - dx2_limited; - const int sign3 = dxvar_well[2*nw + w] > 0 ? 1: -1; - const double dx3_limited = sign3 * std::min(std::abs(dxvar_well[2*nw + w]),dFLimit); - well_state.wellSolutions()[2*nw + w] = xvar_well_old[2*nw + w] - dx3_limited; - F[Water] = well_state.wellSolutions()[nw + w]; - F[Gas] = well_state.wellSolutions()[2*nw + w]; - F[Oil] = 1.0 - F[Water] - F[Gas]; - -// const double dFw = dxvar_well[nw + w]; -// const double dFg = dxvar_well[nw*2 + w]; -// const double dFo = - dFw - dFg; -// //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; -// double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit; -// step = std::min(step, 1.0); -// //std::cout << step << std::endl; -// F[Water] = xvar_well_old[nw + w] - step*dFw; -// F[Gas] = xvar_well_old[2*nw + w] - step*dFg; -// F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo; - - if (F[Water] < 0.0) { - F[Gas] /= (1.0 - F[Water]); - F[Oil] /= (1.0 - F[Water]); - F[Water] = 0.0; - } - if (F[Gas] < 0.0) { - F[Water] /= (1.0 - F[Gas]); - F[Oil] /= (1.0 - F[Gas]); - F[Gas] = 0.0; - } - if (F[Oil] < 0.0) { - F[Water] /= (1.0 - F[Oil]); - F[Gas] /= (1.0 - F[Oil]); - F[Oil] = 0.0; - } - well_state.wellSolutions()[nw + w] = F[Water]; - well_state.wellSolutions()[2*nw + w] = F[Gas]; - - //std::cout << wells().name[w] << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; - - - std::vector g = {1,1,0.01}; - - if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { - for (int p = 0; p < np; ++p) { - F[p] /= distr[p]; - } - } else { - for (int p = 0; p < np; ++p) { - F[p] /= g[p]; - } - } - //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; - - - - - - -// const double dFw = dxvar_well[nw + w]; -// const double dFg = dxvar_well[nw*2 + w]; -// const double dFo = - dFw - dFg; - //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; -// double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit; -// step = std::min(step, 1.0); -// std::cout << step << std::endl; -// F[Water] = xvar_well_old[nw + w] - step*dFw; -// F[Gas] = xvar_well_old[2*nw + w] - step*dFg; -// F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo; -// double sumF = F[Water]+F[Gas]+F[Oil]; -// F[Water] /= sumF; -// F[Gas] /= sumF; -// F[Oil] /= sumF; -// well_state.wellSolutions()[nw + w] = F[Water]; -// well_state.wellSolutions()[2 * nw + w] = F[Gas]; - - switch (well_controls_iget_type(wc, current)) { - case BHP: - { - //const int sign1 = dxvar_well[w] > 0 ? 1: -1; - //const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dTotalRateLimit); - well_state.wellSolutions()[w] = xvar_well_old[w] - dxvar_well[w]; - - switch (wells().type[w]) { - case INJECTOR: - for (int p = 0; p < np; ++p) { - const double comp_frac = wells().comp_frac[np*w + p]; - //if (comp_frac > 0) { - well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[w]; - //} - - } - break; - case PRODUCER: - for (int p = 0; p < np; ++p) { - well_state.wellRates()[w*np + p] = well_state.wellSolutions()[w] * F[p]; - } - break; - } - } - break; - case SURFACE_RATE: - { - const int sign1 = dxvar_well[w] > 0 ? 1: -1; - const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dBHPLimit); - well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited; - //const int sign = (dxvar_well1[w] < 0) ? -1 : 1; - //well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ; - well_state.bhp()[w] = well_state.wellSolutions()[w]; - - if (wells().type[w]==PRODUCER) { - - double F_target = 0.0; - for (int p = 0; p < np; ++p) { - F_target += wells().comp_frac[np*w + p] * F[p]; - } - for (int p = 0; p < np; ++p) { - //std::cout << F[p] << std::endl; - well_state.wellRates()[np*w + p] = F[p] * target_rate /F_target; - } - } else { - - for (int p = 0; p < np; ++p) { - //std::cout << wells().comp_frac[np*w + p] << " " < 0 ? 1: -1; - const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dBHPLimit); - well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited; - //const int sign = (dxvar_well1[w] < 0) ? -1 : 1; - //well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ; - well_state.bhp()[w] = well_state.wellSolutions()[w]; - for (int p = 0; p < np; ++p) { - well_state.wellRates()[np*w + p] = F[p] * target_rate; - } - } - break; - - } - } - - - - - - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - //Loop over all wells -#pragma omp parallel for schedule(static) - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; - const int nwc = well_controls_get_num(wc); - //Loop over all controls until we find a THP control - //that specifies what we need... - //Will only update THP for wells with THP control - for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(wc, ctrl_index) == THP) { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - if ((*active_)[ Water ]) { - aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if ((*active_)[ Oil ]) { - liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if ((*active_)[ Gas ]) { - vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; - } - - double alq = well_controls_iget_alq(wc, ctrl_index); - int table_id = well_controls_iget_vfp(wc, ctrl_index); - - const WellType& well_type = wells().type[w]; - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(), - wellPerforationDensities(), gravity_); - - well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp); - } - else if (well_type == PRODUCER) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(), - wellPerforationDensities(), gravity_); - - well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp, alq); - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); - } - - //Assume only one THP control specified for each well - break; - } - } - } - } - } - - - - - - template - void - StandardWellsDense:: - updateWellControls(WellState& xw) - { - if( !localWellsActive() ) return ; - - std::string modestring[4] = { "BHP", "THP", "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; -#pragma omp parallel for schedule(dynamic) - for (int w = 0; w < nw; ++w) { - 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. - int current = xw.currentControls()[w]; - // Loop over all controls except the current one, and also - // skip any RESERVOIR_RATE controls, since we cannot - // handle those. - const int nwc = well_controls_get_num(wc); - int ctrl_index = 0; - for (; ctrl_index < nwc; ++ctrl_index) { - if (ctrl_index == current) { - // This is the currently used control, so it is - // used as an equation. So this is not used as an - // inequality constraint, and therefore skipped. - continue; - } - if (wellhelpers::constraintBroken( - xw.bhp(), xw.thp(), xw.wellRates(), - 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. - // We disregard terminal_ouput here as with it only messages - // for wells on one process will be printed. - std::ostringstream ss; - ss << "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; - OpmLog::info(ss.str()); - xw.currentControls()[w] = ctrl_index; - current = xw.currentControls()[w]; - well_controls_set_current( wc, current); - - - - // Updating well state and primary variables if constraint is broken - - // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE - const double target = well_controls_iget_target(wc, current); - const double* distr = well_controls_iget_distr(wc, current); - switch (well_controls_iget_type(wc, current)) { - case BHP: - xw.bhp()[w] = target; - break; - - case THP: { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - - if ((*active_)[ Water ]) { - aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if ((*active_)[ Oil ]) { - liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if ((*active_)[ Gas ]) { - vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; - } - - const int vfp = well_controls_iget_vfp(wc, current); - const double& thp = well_controls_iget_target(wc, current); - const double& alq = well_controls_iget_alq(wc, current); - - //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[w]; - - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); - - xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - } - else if (well_type == PRODUCER) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); - - xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - } - else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - break; - } - - case RESERVOIR_RATE: - // No direct change to any observable quantity at - // surface condition. In this case, use existing - // flow rates as initial conditions as reservoir - // rate acts only in aggregate. - break; - - case SURFACE_RATE: - // assign target value as initial guess for injectors and - // single phase producers (orat, grat, wrat) - const WellType& well_type = wells().type[w]; - if (well_type == INJECTOR) { - for (int phase = 0; phase < np; ++phase) { - const double& compi = wells().comp_frac[np * w + phase]; - //if (compi > 0.0) { - xw.wellRates()[np*w + phase] = target * compi; - //} - } - } else if (well_type == PRODUCER) { - - // only set target as initial rates for single phase - // producers. (orat, grat and wrat, and not lrat) - // lrat will result in numPhasesWithTargetsUnderThisControl == 2 - int numPhasesWithTargetsUnderThisControl = 0; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - numPhasesWithTargetsUnderThisControl += 1; - } - } - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { - xw.wellRates()[np*w + phase] = target * distr[phase]; - } - } - } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - - - break; - } - std::vector g = {1,1,0.01}; - if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { - const double* distr = well_controls_iget_distr(wc, current); - for (int phase = 0; phase < np; ++phase) { - g[phase] = distr[phase]; - } - } - switch (well_controls_iget_type(wc, current)) { - case BHP: - { - const WellType& well_type = wells().type[w]; - xw.wellSolutions()[w] = 0.0; - if (well_type == INJECTOR) { - for (int p = 0; p < np; ++p) { - xw.wellSolutions()[w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p]; - } - } else { - for (int p = 0; p < np; ++p) { - xw.wellSolutions()[w] += g[p] * xw.wellRates()[np*w + p]; - } - } - } - break; - - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - { - xw.wellSolutions()[w] = xw.bhp()[w]; - } - break; - } - - double tot_well_rate = 0.0; - for (int p = 0; p < np; ++p) { - tot_well_rate += g[p] * xw.wellRates()[np*w + p]; - } - if(std::abs(tot_well_rate) > 0) { - xw.wellSolutions()[nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; //wells->comp_frac[np*w + Water]; // Water; - xw.wellSolutions()[2*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; //wells->comp_frac[np*w + Gas]; //Gas - } else { - //xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; - //xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; - } - } - } - } - - - - - - template - void - StandardWellsDense:: - addWellFluxEq(const std::vector& cq_s, - const SolutionState& state, - const double dt, - LinearisedBlackoilResidual& residual) - { - if( !localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - - double volume = 0.002831684659200; // 0.1 cu ft; - const Vector vol_dt = Vector::Constant(nw,volume/dt); - std::vector F = wellVolumeFractions(state); - //std::cout << F0_[0] << std::endl; - //std::cout << F[0] << std::endl; - - ADB qs = state.qs; - for (int phase = 0; phase < np; ++phase) { - qs -= superset(wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); - qs += superset((F[phase]-F0_[phase]) * vol_dt, Span(nw,1,phase*nw), nw*np); - } - - - residual.well_flux_eq = qs; - //std::cout << "etter dense " < - void - StandardWellsDense::addWellControlEq(const SolutionState& state, - const WellState& xw, - const Vector& aliveWells, - LinearisedBlackoilResidual& residual) - { - if( ! localWellsActive() ) return; - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - ADB aqua = ADB::constant(Vector::Zero(nw)); - ADB liquid = ADB::constant(Vector::Zero(nw)); - ADB vapour = ADB::constant(Vector::Zero(nw)); - - if ((*active_)[Water]) { - aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw)); - } - if ((*active_)[Oil]) { - liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); - } - if ((*active_)[Gas]) { - vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); - } - - //THP calculation variables - std::vector inj_table_id(nw, -1); - std::vector prod_table_id(nw, -1); - Vector thp_inj_target_v = Vector::Zero(nw); - Vector thp_prod_target_v = Vector::Zero(nw); - Vector alq_v = Vector::Zero(nw); - - //Hydrostatic correction variables - Vector rho_v = Vector::Zero(nw); - Vector vfp_ref_depth_v = Vector::Zero(nw); - - //Target vars - Vector bhp_targets = Vector::Zero(nw); - Vector rate_targets = Vector::Zero(nw); - Eigen::SparseMatrix rate_distr(nw, np*nw); - - //Selection variables - std::vector bhp_elems; - std::vector thp_inj_elems; - std::vector thp_prod_elems; - std::vector rate_elems; - - //Run through all wells to calculate BHP/RATE targets - //and gather info about current control - for (int w = 0; w < nw; ++w) { - auto 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. - const int current = xw.currentControls()[w]; - - switch (well_controls_iget_type(wc, current)) { - case BHP: - { - bhp_elems.push_back(w); - bhp_targets(w) = well_controls_iget_target(wc, current); - rate_targets(w) = -1e100; - } - break; - - case THP: - { - const int perf = wells().well_connpos[w]; - rho_v[w] = wellPerforationDensities()[perf]; - - const int table_id = well_controls_iget_vfp(wc, current); - const double target = well_controls_iget_target(wc, current); - - const WellType& well_type = wells().type[w]; - if (well_type == INJECTOR) { - inj_table_id[w] = table_id; - thp_inj_target_v[w] = target; - alq_v[w] = -1e100; - - vfp_ref_depth_v[w] = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); - - thp_inj_elems.push_back(w); - } - else if (well_type == PRODUCER) { - prod_table_id[w] = table_id; - thp_prod_target_v[w] = target; - alq_v[w] = well_controls_iget_alq(wc, current); - - vfp_ref_depth_v[w] = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); - - thp_prod_elems.push_back(w); - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER type well"); - } - bhp_targets(w) = -1e100; - rate_targets(w) = -1e100; - } - break; - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - { - rate_elems.push_back(w); - // RESERVOIR and SURFACE rates look the same, from a - // high-level point of view, in the system of - // simultaneous linear equations. - - const double* const distr = - well_controls_iget_distr(wc, current); - - for (int p = 0; p < np; ++p) { - rate_distr.insert(w, p*nw + w) = distr[p]; - } - - bhp_targets(w) = -1.0e100; - rate_targets(w) = well_controls_iget_target(wc, current); - } - break; - } - } - - //Calculate BHP target from THP - const ADB thp_inj_target = ADB::constant(thp_inj_target_v); - const ADB thp_prod_target = ADB::constant(thp_prod_target_v); - const ADB alq = ADB::constant(alq_v); - const ADB bhp_from_thp_inj = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target); - const ADB bhp_from_thp_prod = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq); - - //Perform hydrostatic correction to computed targets - const Vector dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, wellPerforationDensities(), gravity_); - const ADB dp = ADB::constant(dp_v); - const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw); - const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw); - - //Calculate residuals - const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj; - const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod; - const ADB bhp_residual = state.bhp - bhp_targets; - const ADB rate_residual = rate_distr * state.qs - rate_targets; - - //Select the right residual for each well - residual.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) + - superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) + - superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) + - superset(subset(rate_residual, rate_elems), rate_elems, nw); - - // For wells that are dead (not flowing), and therefore not communicating - // with the reservoir, we set the equation to be equal to the well's total - // flow. This will be a solution only if the target rate is also zero. - Eigen::SparseMatrix rate_summer(nw, np*nw); - for (int w = 0; w < nw; ++w) { - for (int phase = 0; phase < np; ++phase) { - rate_summer.insert(w, phase*nw + w) = 1.0; - } - } - Selector alive_selector(aliveWells, Selector::NotEqualZero); - residual.well_eq = alive_selector.select(residual.well_eq, rate_summer * state.qs); - // OPM_AD_DUMP(residual_.well_eq); - } - - - - - - template - void - StandardWellsDense::computeWellPotentials(const std::vector& mob_perfcells, - const std::vector& b_perfcells, - SolutionState& state0, - WellState& well_state) - { - const int nw = wells().number_of_wells; - const int np = wells().number_of_phases; - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - - Vector bhps = Vector::Zero(nw); - for (int w = 0; w < nw; ++w) { - const WellControls* ctrl = wells().ctrls[w]; - const int nwc = well_controls_get_num(ctrl); - //Loop over all controls until we find a BHP control - //or a THP control that specifies what we need. - //Pick the value that gives the most restrictive flow - for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { - - if (well_controls_iget_type(ctrl, ctrl_index) == BHP) { - bhps[w] = well_controls_iget_target(ctrl, ctrl_index); - } - - if (well_controls_iget_type(ctrl, ctrl_index) == THP) { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - if ((*active_)[ Water ]) { - aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if ((*active_)[ Oil ]) { - liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if ((*active_)[ Gas ]) { - vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; - } - - const int vfp = well_controls_iget_vfp(ctrl, ctrl_index); - const double& thp = well_controls_iget_target(ctrl, ctrl_index); - const double& alq = well_controls_iget_alq(ctrl, ctrl_index); - - //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[w]; - - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); - const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - // apply the strictest of the bhp controlls i.e. smallest bhp for injectors - if ( bhp < bhps[w]) { - bhps[w] = bhp; - } - } - else if (well_type == PRODUCER) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); - - const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - // apply the strictest of the bhp controlls i.e. largest bhp for producers - if ( bhp > bhps[w]) { - bhps[w] = bhp; - } - } - else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - } - } - - } - - // use bhp limit from control - state0.bhp = ADB::constant(bhps); - - // compute well potentials - Vector aliveWells; - std::vector well_potentials; - computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, well_potentials); - - // store well potentials in the well state - // transform to a single vector instead of separate vectors pr phase - const int nperf = wells().well_connpos[nw]; - Vector cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np); - for (int phase = 1; phase < np; ++phase) { - cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np); - } - well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np); - } - - - - - - void - StandardWellsDense::variableStateWellIndices(std::vector& indices, - int& next) const - { - indices[Qs] = next++; - //indices[Bhp] = next++; - } - - - - - template - void - StandardWellsDense:: - variableStateExtractWellsVars(const std::vector& indices, - std::vector& vars, - SolutionState& state, - WellState &xw) const - { - - state.wellVariables = std::move(vars[indices[Qs]]); - //std::cout << "state.wellVariables " << state.wellVariables.value() << std::endl; - - const int nw = wells().number_of_wells; - const int np = wells().number_of_phases; - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - - Vector isBHPControlled = Vector::Zero(nw); - Vector isReservoirRateControlled = Vector::Zero(nw); - Vector isSurfaceRateControlled = Vector::Zero(nw); - - Vector isInjector = Vector::Zero(nw); - for (int w = 0; w < nw; ++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. - const int current = xw.currentControls()[w]; - - switch (well_controls_iget_type(wc, current)) { - - case BHP: - isBHPControlled[w] = 1; - break; - - case RESERVOIR_RATE: - isReservoirRateControlled[w] = 1; - - break; - case SURFACE_RATE: - isSurfaceRateControlled[w] = 1; - break; - } - - if (wells().type[w] == INJECTOR) { - isInjector[w] = 1; - } - } - - Vector ones = Vector::Constant(nw,1); - - - const ADB& xvar_well1 = subset(state.wellVariables,Span(nw,1,0)); - //std::cout << xvar_well1.value() << std::endl; - const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); - - state.bhp = isBHPControlled * bhp + (1-isBHPControlled) * xvar_well1; - - const ADB Fw = subset(state.wellVariables,Span(nw,1,nw)); - const ADB Fg = subset(state.wellVariables,Span(nw,1,nw*2)); - const ADB Fo = ones - Fw - Fg; - - const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); - - V target_rates = V::Zero(nw); - V isNoFlow = V::Zero(nw); - std::vector distri(np); - for (int p = 0; p < np; ++p) { - distri[p] = V::Zero(nw); - } - - - for (int w = 0; w < nw; ++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. - const int current = xw.currentControls()[w]; - const double* distr = well_controls_iget_distr(wc, current); - - - target_rates[w] = well_controls_iget_target(wc, current); - if (target_rates[w]==0) - isNoFlow[w] = 1; - for (int p = 0; p < np; ++p) { - const double comp_frac = wells().comp_frac[np*w + p]; - distri[p][w] = distr[p]; - //std::cout << "distr " << distr[p] << " comp_frac " << comp_frac << std::endl; - } - - } - - std::vector F(np, ADB::null()); - F[Water] = Fw; - F[Oil] = Fo; - F[Gas] = Fg; // / Vector::Constant(nw,0.01); - - std::vector g(np,Vector::Constant(nw,1.0)); - g[Gas] = Vector::Constant(nw,0.01); -// g[Water] = Vector::Constant(nw,1.0); -// g[Oil] = Vector::Constant(nw,1.0); - for (int p = 0; p < np; ++p) { - const V tmp = (ones - isReservoirRateControlled) * g[p] + isReservoirRateControlled * distri[p]; - F[p] = F[p] / tmp; - } - - ADB targetF = compi.col(pu.phase_pos[0])*F[0]; - for (int p = 1; p < np; ++p) { - targetF += compi.col(pu.phase_pos[p])*F[p]; - } - - - ADB sumF = distri[0] * F[0]; - for (int p = 1; p < np; ++p) { - sumF += distri[p] * (F[p]); - } - Selector noTargetFSelector(targetF.value(), Selector::Zero); - Selector noSumFSelector(sumF.value(), Selector::Zero); - - - //std::cout << target_rates << std::endl; - ADB qs = ADB::constant(V::Zero(nw*np)); - for (int p = 0; p < np; ++p) { - - const int pos = pu.phase_pos[p]; - ADB Q = isSurfaceRateControlled * ADB::constant(compi.col(pos) * target_rates); - //DUMPVAL(Q); - Q += isBHPControlled * isInjector * xvar_well1 * ADB::constant(compi.col(pos)); - //DUMPVAL(Q); - Q += isBHPControlled * (ones - isInjector) * xvar_well1 * F[p]; - //DUMPVAL(Q); - Q += (ones - isInjector) * isSurfaceRateControlled * noTargetFSelector.select( ADB::constant(V::Zero(nw)), (ones - compi.col(pos)) * target_rates * F[p] / targetF); - - //Q += isNoFlow * -1e-24 * F[p]; - - //Q += isSurfaceRateControlled * isInjector * -1e-12 * F[p] * (ones - compi.col(pos)); - - Q += isReservoirRateControlled * F[p] * target_rates; - //DUMPVAL(Q); - qs += superset(Q, Span(nw, 1, p*nw), nw*np); - } - state.qs = qs; - //if (isNoFlow.sum() > 0) - //std::cout << state.qs << std::endl; - - // Bhp. - //state.bhp = - //state.qs = - } - - - - - - std::vector - StandardWellsDense::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; - - variableStateWellIndices(indices, next); - - assert(next == 1); - return indices; - } - - - - - - template - void - StandardWellsDense::variableWellStateInitials(const WellState& xw, - std::vector& vars0) const - { - // Initial well rates. - if ( localWellsActive() ) - { - // 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.wellSolutions()[0], nw, np).transpose(); - //const Vector qs = Eigen::Map(wrates.data(), nw*np); - const Vector qs = Eigen::Map(& xw.wellSolutions()[0], xw.wellSolutions().size()); - vars0.push_back(qs); - - // Initial well bottom-hole pressure. -// assert (not xw.bhp().empty()); -// const Vector 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()); - } - } - - template - std::vector> StandardWellsDense::wellVolumeFractions(const SolutionState& state) const - { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - std::vector F(np,ADB::null()); - F[Water] = subset(state.wellVariables,Span(nw,1,nw)); - F[Gas] = subset(state.wellVariables,Span(nw,1,2 * nw)); - F[Oil] = Vector::Constant(nw,1.0) - F[Water] - F[Gas]; - return F; - } - - - template - void StandardWellsDense::computeAccumWells(const SolutionState& state) { - F0_ = wellVolumeFractions(state); - } - - - - void - StandardWellsDense::setStoreWellPerforationFluxesFlag(const bool store_fluxes) - { - store_well_perforation_fluxes_ = store_fluxes; - } - - - - - - const StandardWellsDense::Vector& - StandardWellsDense::getStoredWellPerforationFluxes() const - { - assert(store_well_perforation_fluxes_); - return well_perforation_fluxes_; - } - - - - - - - template - void - StandardWellsDense:: - updateListEconLimited(ScheduleConstPtr schedule, - const int current_step, - const Wells* wells_struct, - const WellState& well_state, - DynamicListEconLimited& list_econ_limited) const - { - const int nw = wells_struct->number_of_wells; - - for (int w = 0; w < nw; ++w) { - // flag to check if the mim oil/gas rate limit is violated - bool rate_limit_violated = false; - const std::string& well_name = wells_struct->name[w]; - const Well* well_ecl = schedule->getWell(well_name); - const WellEconProductionLimits& econ_production_limits = well_ecl->getEconProductionLimits(current_step); - - // economic limits only apply for production wells. - if (wells_struct->type[w] != PRODUCER) { - continue; - } - - // if no limit is effective here, then continue to the next well - if ( !econ_production_limits.onAnyEffectiveLimit() ) { - continue; - } - // for the moment, we only handle rate limits, not handling potential limits - // the potential limits should not be difficult to add - const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit(); - if (quantity_limit == WellEcon::POTN) { - const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n") - + std::string("All the limits will be evaluated based on RATE. "); - OpmLog::warning("NOT_SUPPORTING_POTN", msg); - } - - const WellMapType& well_map = well_state.wellMap(); - const typename WellMapType::const_iterator i_well = well_map.find(well_name); - assert(i_well != well_map.end()); // should always be found? - const WellMapEntryType& map_entry = i_well->second; - const int well_number = map_entry[0]; - - if (econ_production_limits.onAnyRateLimit()) { - rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state, well_number); - } - - if (rate_limit_violated) { - if (econ_production_limits.endRun()) { - const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n") - + std::string("the program will keep running after ") + well_name + std::string(" is closed"); - OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message); - } - - if (econ_production_limits.validFollowonWell()) { - OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet"); - } - - if (well_ecl->getAutomaticShutIn()) { - list_econ_limited.addShutWell(well_name); - const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit"); - OpmLog::info(msg); - } else { - list_econ_limited.addStoppedWell(well_name); - const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit"); - OpmLog::info(msg); - } - // the well is closed, not need to check other limits - continue; - } - - // checking for ratio related limits, mostly all kinds of ratio. - bool ratio_limits_violated = false; - RatioCheckTuple ratio_check_return; - - if (econ_production_limits.onAnyRatioLimit()) { - ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state, map_entry); - ratio_limits_violated = std::get<0>(ratio_check_return); - } - - if (ratio_limits_violated) { - const bool last_connection = std::get<1>(ratio_check_return); - const int worst_offending_connection = std::get<2>(ratio_check_return); - - const int perf_start = map_entry[1]; - - assert((worst_offending_connection >= 0) && (worst_offending_connection < map_entry[2])); - - const int cell_worst_offending_connection = wells_struct->well_cells[perf_start + worst_offending_connection]; - list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection); - const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ") - + well_name + std::string(" will be closed due to economic limit"); - OpmLog::info(msg); - - if (last_connection) { - list_econ_limited.addShutWell(well_name); - const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed"); - OpmLog::info(msg2); - } - } - - } - } - - - - - - template - bool - StandardWellsDense:: - checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, - const WellState& well_state, - const int well_number) const - { - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - const int np = well_state.numPhases(); - - if (econ_production_limits.onMinOilRate()) { - assert((*active_)[Oil]); - const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; - const double min_oil_rate = econ_production_limits.minOilRate(); - if (std::abs(oil_rate) < min_oil_rate) { - return true; - } - } - - if (econ_production_limits.onMinGasRate() ) { - assert((*active_)[Gas]); - const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ]; - const double min_gas_rate = econ_production_limits.minGasRate(); - if (std::abs(gas_rate) < min_gas_rate) { - return true; - } - } - - if (econ_production_limits.onMinLiquidRate() ) { - assert((*active_)[Oil]); - assert((*active_)[Water]); - const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; - const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; - const double liquid_rate = oil_rate + water_rate; - const double min_liquid_rate = econ_production_limits.minLiquidRate(); - if (std::abs(liquid_rate) < min_liquid_rate) { - return true; - } - } - - if (econ_production_limits.onMinReservoirFluidRate()) { - OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet"); - } - - return false; - } - - - - - - template - StandardWellsDense::RatioCheckTuple - StandardWellsDense:: - checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, - const WellState& well_state, - const WellMapEntryType& map_entry) const - { - // TODO: not sure how to define the worst-offending connection when more than one - // ratio related limit is violated. - // The defintion used here is that we define the violation extent based on the - // ratio between the value and the corresponding limit. - // For each violated limit, we decide the worst-offending connection separately. - // Among the worst-offending connections, we use the one has the biggest violation - // extent. - - bool any_limit_violated = false; - bool last_connection = false; - int worst_offending_connection = INVALIDCONNECTION; - double violation_extent = -1.0; - - if (econ_production_limits.onMaxWaterCut()) { - const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry); - bool water_cut_violated = std::get<0>(water_cut_return); - if (water_cut_violated) { - any_limit_violated = true; - const double violation_extent_water_cut = std::get<3>(water_cut_return); - if (violation_extent_water_cut > violation_extent) { - violation_extent = violation_extent_water_cut; - worst_offending_connection = std::get<2>(water_cut_return); - last_connection = std::get<1>(water_cut_return); - } - } - } - - if (econ_production_limits.onMaxGasOilRatio()) { - OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!"); - } - - if (econ_production_limits.onMaxWaterGasRatio()) { - OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!"); - } - - if (econ_production_limits.onMaxGasLiquidRatio()) { - OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!"); - } - - if (any_limit_violated) { - assert(worst_offending_connection >=0); - assert(violation_extent > 1.); - } - - return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent); - } - - - - - - template - StandardWellsDense::RatioCheckTuple - StandardWellsDense:: - checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, - const WellState& well_state, - const WellMapEntryType& map_entry) const - { - bool water_cut_limit_violated = false; - int worst_offending_connection = INVALIDCONNECTION; - bool last_connection = false; - double violation_extent = -1.0; - - const int np = well_state.numPhases(); - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - const int well_number = map_entry[0]; - - assert((*active_)[Oil]); - assert((*active_)[Water]); - - const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; - const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; - const double liquid_rate = oil_rate + water_rate; - double water_cut; - if (std::abs(liquid_rate) != 0.) { - water_cut = water_rate / liquid_rate; - } else { - water_cut = 0.0; - } - - const double max_water_cut_limit = econ_production_limits.maxWaterCut(); - if (water_cut > max_water_cut_limit) { - water_cut_limit_violated = true; - } - - if (water_cut_limit_violated) { - // need to handle the worst_offending_connection - const int perf_start = map_entry[1]; - const int perf_number = map_entry[2]; - - std::vector water_cut_perf(perf_number); - for (int perf = 0; perf < perf_number; ++perf) { - const int i_perf = perf_start + perf; - const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ]; - const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ]; - const double liquid_perf_rate = oil_perf_rate + water_perf_rate; - if (std::abs(liquid_perf_rate) != 0.) { - water_cut_perf[perf] = water_perf_rate / liquid_perf_rate; - } else { - water_cut_perf[perf] = 0.; - } - } - - last_connection = (perf_number == 1); - if (last_connection) { - worst_offending_connection = 0; - violation_extent = water_cut_perf[0] / max_water_cut_limit; - return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); - } - - double max_water_cut_perf = 0.; - for (int perf = 0; perf < perf_number; ++perf) { - if (water_cut_perf[perf] > max_water_cut_perf) { - worst_offending_connection = perf; - max_water_cut_perf = water_cut_perf[perf]; - } - } - - assert(max_water_cut_perf != 0.); - assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number)); - - violation_extent = max_water_cut_perf / max_water_cut_limit; - } - - return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); - } - - -} // namespace Opm