diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 03e750c84..8a7368977 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -90,25 +90,7 @@ enum WellVariablePositions { StandardWellsDense(const Wells* wells_arg, WellCollection* well_collection, const ModelParameters& param, - const bool terminal_output) - : wells_active_(wells_arg!=nullptr) - , wells_(wells_arg) - , well_collection_(well_collection) - , param_(param) - , terminal_output_(terminal_output) - , well_perforation_efficiency_factors_((wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0), 1.0) - , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) - , well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) - , wellVariables_( wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0) - , F0_(wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0 ) - { - if( wells_ ) - { - invDuneD_.setBuildMode( Mat::row_wise ); - duneC_.setBuildMode( Mat::row_wise ); - duneB_.setBuildMode( Mat::row_wise ); - } - } + const bool terminal_output); void init(const PhaseUsage phase_usage_arg, const std::vector& active_arg, @@ -116,1276 +98,141 @@ enum WellVariablePositions { const double gravity_arg, const std::vector& depth_arg, const std::vector& pv_arg, - const RateConverterType* rate_converter) - { - - if ( ! localWellsActive() ) { - return; - } - - phase_usage_ = phase_usage_arg; - active_ = active_arg; - vfp_properties_ = vfp_properties_arg; - gravity_ = gravity_arg; - cell_depths_ = extractPerfData(depth_arg); - pv_ = pv_arg; - rate_converter_ = rate_converter; - - calculateEfficiencyFactors(); - - // setup sparsity pattern for the matrices - //[A B^T [x = [ res - // C D] x_well] res_well] - - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - const int nc = numCells(); - -#ifndef NDEBUG - const auto pu = phase_usage_; - const int np = pu.num_phases; - - // assumes the gas fractions are stored after water fractions - // WellVariablePositions needs to be changed for 2p runs - assert (np == 3 || (np == 2 && !pu.phase_used[Gas]) ); -#endif - - // set invDuneD - invDuneD_.setSize( nw, nw, nw ); - - // set duneC - duneC_.setSize( nw, nc, nperf ); - - // set duneB - duneB_.setSize( nw, nc, nperf ); - - for(auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) { - // Add nonzeros for diagonal - row.insert(row.index()); - } - - for(auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) { - // Add nonzeros for diagonal - for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { - const int cell_idx = wells().well_cells[perf]; - row.insert(cell_idx); - } - } - - // make the B^T matrix - for(auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) { - for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { - const int cell_idx = wells().well_cells[perf]; - row.insert(cell_idx); - } - } - - resWell_.resize( nw ); - - // resize temporary class variables - Cx_.resize( duneC_.N() ); - invDrw_.resize( invDuneD_.N() ); - } + const RateConverterType* rate_converter); template SimulatorReport assemble(Simulator& ebosSimulator, const int iterationIdx, const double dt, - WellState& well_state) { - - SimulatorReport report; - if ( ! localWellsActive() ) { - return report; - } - - if (param_.compute_well_potentials_) { - computeWellPotentials(ebosSimulator, well_state); - } - - resetWellControlFromState(well_state); - updateWellControls(well_state); - // Set the primary variables for the wells - setWellVariables(well_state); - - if (iterationIdx == 0) { - computeWellConnectionPressures(ebosSimulator, well_state); - computeAccumWells(); - } - - if (param_.solve_welleq_initially_ && iterationIdx == 0) { - // solve the well equations as a pre-processing step - report = solveWellEq(ebosSimulator, dt, well_state); - } - assembleWellEq(ebosSimulator, dt, well_state, false); - - report.converged = true; - return report; - } + WellState& well_state); template void assembleWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state, - bool only_wells) { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - - // clear all entries - duneB_ = 0.0; - duneC_ = 0.0; - invDuneD_ = 0.0; - resWell_ = 0.0; - - auto& ebosJac = ebosSimulator.model().linearizer().matrix(); - auto& ebosResid = ebosSimulator.model().linearizer().residual(); - - const double volume = 0.002831684659200; // 0.1 cu ft; - for (int w = 0; w < nw; ++w) { - bool allow_cf = allow_cross_flow(w, ebosSimulator); - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - - const int cell_idx = wells().well_cells[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - std::vector cq_s(np,0.0); - const EvalWell bhp = getBhp(w); - computeWellFlux(w, wells().WI[perf], intQuants, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); - - for (int p1 = 0; p1 < np; ++p1) { - - // the cq_s entering mass balance equations need to consider the efficiency factors. - const EvalWell cq_s_effective = cq_s[p1] * well_perforation_efficiency_factors_[perf]; - - if (!only_wells) { - // subtract sum of phase fluxes in the reservoir equation. - // need to consider the efficiency factor - ebosResid[cell_idx][flowPhaseToEbosCompIdx(p1)] -= cq_s_effective.value(); - } - - // subtract sum of phase fluxes in the well equations. - resWell_[w][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].value(); - - // assemble the jacobians - for (int p2 = 0; p2 < np; ++p2) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s_effective.derivative(p2); - duneB_[w][cell_idx][flowToEbosPvIdx(p2)][flowPhaseToEbosCompIdx(p1)] -= cq_s_effective.derivative(p2+blocksize); // intput in transformed matrix - duneC_[w][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s_effective.derivative(p2); - } - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2+blocksize); - } - - // add trivial equation for 2p cases (Only support water + oil) - if (np == 2) { - assert(!active_[ Gas ]); - invDuneD_[w][w][flowPhaseToEbosCompIdx(Gas)][flowToEbosPvIdx(Gas)] = 1.0; - } - - // Store the perforation phase flux for later usage. - well_state.perfPhaseRates()[perf*np + p1] = cq_s[p1].value(); - } - // Store the perforation pressure for later usage. - well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf]; - } - - // add vol * dF/dt + Q to the well equations; - for (int p1 = 0; p1 < np; ++p1) { - EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; - resWell_loc += getQs(w, p1); - for (int p2 = 0; p2 < np; ++p2) { - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivative(p2+blocksize); - } - resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value(); - } - } - - // do the local inversion of D. - localInvert( invDuneD_ ); - - } + bool only_wells); template - bool allow_cross_flow(const int w, Simulator& ebosSimulator) const { + bool allow_cross_flow(const int w, Simulator& ebosSimulator) const; - if (wells().allow_cf[w]) { - return true; - } - // check for special case where all perforations have cross flow - // then the wells must allow for cross flow - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - const int cell_idx = wells().well_cells[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& fs = intQuants.fluidState(); - EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); - EvalWell bhp = getBhp(w); - // Pressure drawdown (also used to determine direction of flow) - EvalWell well_pressure = bhp + wellPerforationPressureDiffs()[perf]; - EvalWell drawdown = pressure - well_pressure; + void localInvert(Mat& istlA) const; - if ( drawdown.value() < 0 && wells().type[w] == INJECTOR) { - return false; - } - if ( drawdown.value() > 0 && wells().type[w] == PRODUCER) { - return false; - } - } - return true; - } - - void localInvert(Mat& istlA) const { - for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { - for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { - //std::cout << (*col) << std::endl; - (*col).invert(); - } - } - } - - void print(Mat& istlA) const { - for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { - for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { - std::cout << row.index() << " " << col.index() << "/n \n"<<(*col) << std::endl; - } - } - } + void print(Mat& istlA) const; // substract Binv(D)rw from r; - void apply( BVector& r) const { - if ( ! localWellsActive() ) { - return; - } - - assert( invDrw_.size() == invDuneD_.N() ); - - invDuneD_.mv(resWell_,invDrw_); - duneB_.mmtv(invDrw_, r); - } + void apply( BVector& r) const; // subtract B*inv(D)*C * x from A*x - void apply(const BVector& x, BVector& Ax) { - if ( ! localWellsActive() ) { - return; - } - assert( Cx_.size() == duneC_.N() ); - - BVector& invDCx = invDrw_; - assert( invDCx.size() == invDuneD_.N()); - - duneC_.mv(x, Cx_); - invDuneD_.mv(Cx_, invDCx); - duneB_.mmtv(invDCx,Ax); - } + void apply(const BVector& x, BVector& Ax); // apply well model with scaling of alpha - void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) - { - if ( ! localWellsActive() ) { - return; - } - - if( scaleAddRes_.size() != Ax.size() ) { - scaleAddRes_.resize( Ax.size() ); - } - - scaleAddRes_ = 0.0; - apply( x, scaleAddRes_ ); - Ax.axpy( alpha, scaleAddRes_ ); - } + void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax); // xw = inv(D)*(rw - C*x) - void recoverVariable(const BVector& x, BVector& xw) const { - if ( ! localWellsActive() ) { - return; - } - BVector resWell = resWell_; - duneC_.mmv(x, resWell); - invDuneD_.mv(resWell, xw); - } + void recoverVariable(const BVector& x, BVector& xw) const; + int flowPhaseToEbosCompIdx( const int phaseIdx ) const; - int flowPhaseToEbosCompIdx( const int phaseIdx ) const - { - const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; - return phaseToComp[ phaseIdx ]; - } - - int flowToEbosPvIdx( const int flowPv ) const - { - const int flowToEbos[ 3 ] = { - BlackoilIndices::pressureSwitchIdx, - BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx - }; - return flowToEbos[ flowPv ]; - } - - int ebosCompToFlowPhaseIdx( const int compIdx ) const - { - const int compToPhase[ 3 ] = { Oil, Water, Gas }; - return compToPhase[ compIdx ]; - } + int flowToEbosPvIdx( const int flowPv ) const; + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; + int ebosCompToFlowPhaseIdx( const int compIdx ) const; std::vector - extractPerfData(const std::vector& in) const { - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - std::vector out(nperf); - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - const int well_idx = wells().well_cells[perf]; - out[perf] = in[well_idx]; - } - } - return out; - } + extractPerfData(const std::vector& in) const; - int numPhases() const { return wells().number_of_phases; } + int numPhases() const; - int numCells() const { return pv_.size();} + int numCells() const; - template - void resetWellControlFromState(WellState xw) { - const int nw = wells_->number_of_wells; - for (int w = 0; w < nw; ++w) { - WellControls* wc = wells_->ctrls[w]; - well_controls_set_current( wc, xw.currentControls()[w]); - } - } + void resetWellControlFromState(WellState xw) const; - const Wells& wells() const - { - assert(wells_ != 0); - return *(wells_); - } + const Wells& wells() const; - const Wells* wellsPointer() const - { - return wells_; - } + const Wells* wellsPointer() const; /// return true if wells are available in the reservoir - bool wellsActive() const - { - return wells_active_; - } - void setWellsActive(const bool wells_active) - { - wells_active_ = wells_active; - } + bool wellsActive() const; + + void setWellsActive(const bool wells_active); + /// return true if wells are available on this process - bool localWellsActive() const - { - return wells_ ? (wells_->number_of_wells > 0 ) : false; - } + bool localWellsActive() 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; - } + int numWellVars() const; /// Density of each well perforation - std::vector wellPerforationDensities() const { - return well_perforation_densities_; - } + const std::vector& wellPerforationDensities() const; /// Diff to bhp for each well perforation. - std::vector wellPerforationPressureDiffs() const - { - return well_perforation_pressure_diffs_; - } - + const std::vector& wellPerforationPressureDiffs() const; typedef DenseAd::Evaluation Eval; - EvalWell extendEval(Eval in) const { - EvalWell out = 0.0; - out.setValue(in.value()); - for(int i = 0; i < blocksize;++i) { - out.setDerivative(i, in.derivative(flowToEbosPvIdx(i))); - } - return out; - } - 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].setValue(xw.wellSolutions()[w + nw* phaseIdx]); - wellVariables_[w + nw*phaseIdx].setDerivative(blocksize + phaseIdx, 1.0); - } - } - } + EvalWell extendEval(Eval in) const; - void print(EvalWell in) const { - std::cout << in.value() << std::endl; - for (int i = 0; i < in.size; ++i) { - std::cout << in.derivative(i) << std::endl; - } - } - - 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(); - } - } - } + void setWellVariables(const WellState& xw); + void print(EvalWell in) const; + void computeAccumWells(); template void - computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants, const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector& cq_s) const - { - const Opm::PhaseUsage& pu = phase_usage_; - const int np = wells().number_of_phases; - std::vector cmix_s(np,0.0); - for (int phase = 0; phase < np; ++phase) { - //int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - cmix_s[phase] = wellVolumeFraction(w,phase); - } - const auto& fs = intQuants.fluidState(); - 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; - EvalWell drawdown = pressure - well_pressure; - - // injection perforations - if ( drawdown.value() > 0 ) { - - //Do nothing if crossflow is not allowed - if (!allow_cf && wells().type[w] == INJECTOR) - return; - // 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 * (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[phase] = cq_ps[phase]; - } - - } else { - //Do nothing if crossflow is not allowed - if (!allow_cf && wells().type[w] == PRODUCER) - return; - - // 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 * (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.setValue(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_s) const; template SimulatorReport solveWellEq(Simulator& ebosSimulator, const double dt, - WellState& well_state) - { - const int nw = wells().number_of_wells; - WellState well_state0 = well_state; + WellState& well_state); - int it = 0; - bool converged; - do { - assembleWellEq(ebosSimulator, dt, well_state, true); - converged = getWellConvergence(ebosSimulator, it); - - // checking whether the group targets are converged - if (wellCollection()->groupControlActive()) { - converged = converged && wellCollection()->groupTargetConverged(well_state.wellRates()); - } - - if (converged) { - break; - } - - ++it; - if( localWellsActive() ) - { - BVector dx_well (nw); - invDuneD_.mv(resWell_, dx_well); - - updateWellState(dx_well, well_state); - updateWellControls(well_state); - setWellVariables(well_state); - } - } while (it < 15); - - if (!converged) { - well_state = well_state0; - } - - SimulatorReport report; - report.converged = converged; - report.total_well_iterations = it; - return report; - } - - void printIf(int c, double x, double y, double eps, std::string type) { - if (std::abs(x-y) > eps) { - std::cout << type << " " << c << ": "< residual() { - if( ! wellsActive() ) - { - return std::vector(); - } - - const int np = numPhases(); - const int nw = wells().number_of_wells; - std::vector res(np*nw); - for( int p=0; p residual() const; template bool getWellConvergence(Simulator& ebosSimulator, - const int iteration) - { - typedef std::vector< double > Vector; - const int np = numPhases(); - const int nc = numCells(); - const double tol_wells = param_.tolerance_wells_; - const double maxResidualAllowed = param_.max_residual_allowed_; - - Vector R_sum(np); - Vector B_avg(np); - Vector maxCoeff(np); - Vector maxNormWell(np); - - std::vector< Vector > B( np, Vector( nc ) ); - std::vector< Vector > R2( np, Vector( nc ) ); - std::vector< Vector > tempV( np, Vector( nc ) ); - - for ( int idx = 0; idx < np; ++idx ) - { - Vector& B_idx = B[ idx ]; - const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); - - 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(); - - B_idx [cell_idx] = 1 / fs.invB(ebosPhaseIdx).value(); - } - } - - detail::convergenceReduction(B, tempV, R2, - R_sum, maxCoeff, B_avg, maxNormWell, - nc, np, pv_, residual() ); - - - 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; - } - - + const int iteration) const; 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); - computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, cell_depths_, gravity_); + const WellState& xw); - } - - template + 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 PhaseUsage& pu = phase_usage_; - const int np = phase_usage_.num_phases; - b_perf.resize(nperf*np); - surf_dens_perf.resize(nperf*np); + std::vector& surf_dens_perf) const; - //rs and rv are only used if both oil and gas is present - if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) { - rsmax_perf.resize(nperf); - rvmax_perf.resize(nperf); - } - - // Compute the average pressure in each well block - 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 = wells().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] : xw.perfPress()[perf - 1]; - const double p_avg = (xw.perfPress()[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; - int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; - - if (pu.phase_used[BlackoilPhases::Liquid]) { - int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; - const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers - rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); - if (oilrate > 0) { - const double gasrate = std::abs(xw.wellRates()[gaspos_well]); - double rv = 0.0; - if (gasrate > 0) { - rv = oilrate / gasrate; - } - rv = std::min(rv, rvmax_perf[perf]); - - b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); - } - else { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - - } else { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } - - if (pu.phase_used[BlackoilPhases::Liquid]) { - int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; - int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; - if (pu.phase_used[BlackoilPhases::Vapour]) { - rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); - int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; - const double gasrate = std::abs(xw.wellRates()[gaspos_well]); - if (gasrate > 0) { - const double oilrate = std::abs(xw.wellRates()[oilpos_well]); - double rs = 0.0; - if (oilrate > 0) { - rs = gasrate / oilrate; - } - rs = std::min(rs, rsmax_perf[perf]); - b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); - } else { - b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } else { - b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } - - // Surface density. - for (int p = 0; p < pu.num_phases; ++p) { - surf_dens_perf[np*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); - } - } - } - } - - template void updateWellState(const BVector& dwells, - WellState& well_state) - { - if( localWellsActive() ) - { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - double dFLimit = dWellFractionMax(); - double dBHPLimit = dbhpMaxRel(); - std::vector xvar_well_old = well_state.wellSolutions(); - - for (int w = 0; w < nw; ++w) { - - // update the second and third well variable (The flux fractions) - std::vector F(np,0.0); - if (active_[ Water ]) { - const int sign2 = dwells[w][flowPhaseToEbosCompIdx(WFrac)] > 0 ? 1: -1; - const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(WFrac)]),dFLimit); - well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited; - } - - if (active_[ Gas ]) { - const int sign3 = dwells[w][flowPhaseToEbosCompIdx(GFrac)] > 0 ? 1: -1; - const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(GFrac)]),dFLimit); - well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited; - } - - assert(active_[ Oil ]); - F[Oil] = 1.0; - if (active_[ Water ]) { - F[Water] = well_state.wellSolutions()[WFrac*nw + w]; - F[Oil] -= F[Water]; - } - - if (active_[ Gas ]) { - F[Gas] = well_state.wellSolutions()[GFrac*nw + w]; - F[Oil] -= F[Gas]; - } - - if (active_[ Water ]) { - if (F[Water] < 0.0) { - if (active_[ Gas ]) { - F[Gas] /= (1.0 - F[Water]); - } - F[Oil] /= (1.0 - F[Water]); - F[Water] = 0.0; - } - } - if (active_[ Gas ]) { - if (F[Gas] < 0.0) { - if (active_[ Water ]) { - F[Water] /= (1.0 - F[Gas]); - } - F[Oil] /= (1.0 - F[Gas]); - F[Gas] = 0.0; - } - } - if (F[Oil] < 0.0) { - if (active_[ Water ]) { - F[Water] /= (1.0 - F[Oil]); - } - if (active_[ Gas ]) { - F[Gas] /= (1.0 - F[Oil]); - } - F[Oil] = 0.0; - } - - if (active_[ Water ]) { - well_state.wellSolutions()[WFrac*nw + w] = F[Water]; - } - if (active_[ Gas ]) { - well_state.wellSolutions()[GFrac*nw + w] = F[Gas]; - } - - // The interpretation of the first well variable depends on the well control - 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 = well_state.currentControls()[w]; - const double target_rate = well_controls_iget_target(wc, current); - - 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 p = 0; p < np; ++p) { - F[p] /= distr[p]; - } - } else { - for (int p = 0; p < np; ++p) { - F[p] /= g[p]; - } - } - - switch (well_controls_iget_type(wc, current)) { - case THP: // The BHP and THP both uses the total rate as first well variable. - case BHP: - { - well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][flowPhaseToEbosCompIdx(XvarWell)]; - - switch (wells().type[w]) { - case INJECTOR: - for (int p = 0; p < np; ++p) { - const double comp_frac = wells().comp_frac[np*w + p]; - well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[nw*XvarWell + w]; - } - break; - case PRODUCER: - for (int p = 0; p < np; ++p) { - well_state.wellRates()[w*np + p] = well_state.wellSolutions()[nw*XvarWell + w] * F[p]; - } - break; - } - - if (well_controls_iget_type(wc, current) == THP) { - - // Calculate bhp from thp control and well rates - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - const Opm::PhaseUsage& pu = phase_usage_; - - 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(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]; - // pick the density in the top layer - const int perf = wells().well_connpos[w]; - const double rho = well_perforation_densities_[perf]; - - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - rho, gravity_); - - well_state.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(), - rho, gravity_); - - well_state.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); - } - } - - } - break; - case SURFACE_RATE: // Both rate controls use bhp as first well variable - case RESERVOIR_RATE: - { - const int sign1 = dwells[w][flowPhaseToEbosCompIdx(XvarWell)] > 0 ? 1: -1; - const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(XvarWell)]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit); - well_state.wellSolutions()[nw*XvarWell + w] = std::max(xvar_well_old[nw*XvarWell + w] - dx1_limited,1e5); - well_state.bhp()[w] = well_state.wellSolutions()[nw*XvarWell + w]; - - if (well_controls_iget_type(wc, current) == SURFACE_RATE) { - if (wells().type[w]==PRODUCER) { - - const double* distr = well_controls_iget_distr(wc, current); - - double F_target = 0.0; - for (int p = 0; p < np; ++p) { - F_target += distr[p] * F[p]; - } - for (int p = 0; p < np; ++p) { - well_state.wellRates()[np*w + p] = F[p] * target_rate / F_target; - } - } else { - - for (int p = 0; p < np; ++p) { - well_state.wellRates()[w*np + p] = wells().comp_frac[np*w + p] * target_rate; - } - } - } else { // RESERVOIR_RATE - for (int p = 0; p < np; ++p) { - well_state.wellRates()[np*w + p] = F[p] * target_rate; - } - } - } - break; - } - } - } - } + WellState& well_state) const; - template - void updateWellControls(WellState& xw) - { - if( !localWellsActive() ) return ; - - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - // keeping a copy of the current controls, to see whether control changes later. - std::vector old_control_index(nw, 0); - for (int w = 0; w < nw; ++w) { - old_control_index[w] = xw.currentControls()[w]; - } - - // Find, for each well, if any constraints are broken. If so, - // switch control to first broken constraint. - #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. - xw.currentControls()[w] = ctrl_index; - current = xw.currentControls()[w]; - well_controls_set_current( wc, current); - } - - // update whether well is under group control - if (wellCollection()->groupControlActive()) { - // get well node in the well collection - WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w])); - - // update whehter the well is under group control or individual control - if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) { - // under group control - well_node.setIndividualControl(false); - } else { - // individual control - well_node.setIndividualControl(true); - } - } - } - - // upate the well targets following group controls - if (wellCollection()->groupControlActive()) { - applyVREPGroupControl(xw); - wellCollection()->updateWellTargets(xw.wellRates()); - } - - // the new well control indices after all the related updates, - std::vector updated_control_index(nw, 0); - for (int w = 0; w < nw; ++w) { - updated_control_index[w] = xw.currentControls()[w]; - } - - // checking whether control changed - wellhelpers::WellSwitchingLogger logger; - for (int w = 0; w < nw; ++w) { - if (updated_control_index[w] != old_control_index[w]) { - WellControls* wc = wells().ctrls[w]; - logger.wellSwitched(wells().name[w], - well_controls_iget_type(wc, old_control_index[w]), - well_controls_iget_type(wc, updated_control_index[w])); - updateWellStateWithTarget(wc, updated_control_index[w], w, xw); - } - } - } - - - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const - { - const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; - return flowToEbos[ phaseIdx ]; - } + void updateWellControls(WellState& xw) const; /// upate the dynamic lists related to economic limits - template void updateListEconLimited(const Schedule& schedule, const int current_step, const Wells* wells_struct, const WellState& well_state, - DynamicListEconLimited& list_econ_limited) const - { - // With no wells (on process) wells_struct is a null pointer - const int nw = (wells_struct)? wells_struct->number_of_wells : 0; + DynamicListEconLimited& list_econ_limited) 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; - } - - // 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, 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 - well_perforation_densities_ = - WellDensitySegmented::computeConnectionDensities( - wells(), xw, phase_usage_, - b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - // Compute pressure deltas - well_perforation_pressure_diffs_ = - WellDensitySegmented::computeConnectionPressureDelta( - wells(), depth_perf, well_perforation_densities_, grav); - - - - } - - - + const double grav); // TODO: Later we might want to change the function to only handle one well, @@ -1394,220 +241,20 @@ enum WellVariablePositions { template void computeWellPotentials(const Simulator& ebosSimulator, - WellState& well_state) const - { - - // number of wells and phases - const int nw = wells().number_of_wells; - const int np = wells().number_of_phases; - - for (int w = 0; w < nw; ++w) { - // bhp needs to be determined for the well potential calculation - double bhp = 0.; - - const WellControls* well_control = wells().ctrls[w]; - // The number of the well controls - const int nwc = well_controls_get_num(well_control); - - // Finding a BHP control or a THP control - // IF we find a THP control, we calculate the BHP value. - // TODO: there is option to ignore the THP limit when calculating well potentials, - // we are not handling it for the moment. - for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(well_control, ctrl_index) == BHP) { - // set bhp to the bhp value - bhp = well_controls_iget_target(well_control, ctrl_index); - } - - - if (well_controls_iget_type(well_control, ctrl_index) == THP) { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - const Opm::PhaseUsage& pu = phase_usage_; - - 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(well_control, ctrl_index); - const double& thp = well_controls_iget_target(well_control, ctrl_index); - const double& alq = well_controls_iget_alq(well_control, ctrl_index); - - // Calculating the BHP value based on THP - const WellType& well_type = wells().type[w]; - const int first_perf = wells().well_connpos[w]; //first perforation - - if (well_type == INJECTOR) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities()[first_perf], gravity_); - const double bhp_calculated = 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_calculated < bhp) { - bhp = bhp_calculated; - } - } - else if (well_type == PRODUCER) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities()[first_perf], gravity_); - const double bhp_calculated = 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_calculated > bhp) { - bhp = bhp_calculated; - } - } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - } - } - - assert(bhp != 0.0); - - // Should we consider crossflow when calculating well potentionals? - const bool allow_cf = allow_cross_flow(w, ebosSimulator); - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const int cell_index = wells().well_cells[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0)); - std::vector well_potentials(np, 0.0); - computeWellFlux(w, wells().WI[perf], intQuants, bhp, wellPerforationPressureDiffs()[perf], allow_cf, well_potentials); - for(int p = 0; p < np; ++p) { - well_state.wellPotentials()[perf * np + p] = well_potentials[p].value(); - } - } - } - } - - - - - - WellCollection* wellCollection() const - { - return well_collection_; - } - - - + WellState& well_state) const; + WellCollection* wellCollection() const; const std::vector& - wellPerfEfficiencyFactors() const - { - return well_perforation_efficiency_factors_; - } - - - - - - void calculateEfficiencyFactors() - { - if ( !localWellsActive() ) { - return; - } - - const int nw = wells().number_of_wells; - - for (int w = 0; w < nw; ++w) { - const std::string well_name = wells().name[w]; - const WellNode& well_node = wellCollection()->findWellNode(well_name); - - const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor(); - - // assign the efficiency factor to each perforation related. - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w + 1]; ++perf) { - well_perforation_efficiency_factors_[perf] = well_efficiency_factor; - } - } - } - - + wellPerfEfficiencyFactors() const; + void calculateEfficiencyFactors(); void computeWellVoidageRates(const WellState& well_state, std::vector& well_voidage_rates, - std::vector& voidage_conversion_coeffs) const - { - if ( !localWellsActive() ) { - return; - } - // TODO: for now, we store the voidage rates for all the production wells. - // For injection wells, the rates are stored as zero. - // We only store the conversion coefficients for all the injection wells. - // Later, more delicate model will be implemented here. - // And for the moment, group control can only work for serial running. - const int nw = well_state.numWells(); - const int np = well_state.numPhases(); - - // we calculate the voidage rate for each well, that means the sum of all the phases. - well_voidage_rates.resize(nw, 0); - // store the conversion coefficients, while only for the use of injection wells. - voidage_conversion_coeffs.resize(nw * np, 1.0); - - std::vector well_rates(np, 0.0); - std::vector convert_coeff(np, 1.0); - - for (int w = 0; w < nw; ++w) { - const bool is_producer = wells().type[w] == PRODUCER; - - // not sure necessary to change all the value to be positive - if (is_producer) { - std::transform(well_state.wellRates().begin() + np * w, - well_state.wellRates().begin() + np * (w + 1), - well_rates.begin(), std::negate()); - - // the average hydrocarbon conditions of the whole field will be used - const int fipreg = 0; // Not considering FIP for the moment. - - rate_converter_->calcCoeff(well_rates, fipreg, convert_coeff); - well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(), - convert_coeff.begin(), 0.0); - } else { - // TODO: Not sure whether will encounter situation with all zero rates - // and whether it will cause problem here. - std::copy(well_state.wellRates().begin() + np * w, - well_state.wellRates().begin() + np * (w + 1), - well_rates.begin()); - // the average hydrocarbon conditions of the whole field will be used - const int fipreg = 0; // Not considering FIP for the moment. - rate_converter_->calcCoeff(well_rates, fipreg, convert_coeff); - std::copy(convert_coeff.begin(), convert_coeff.end(), - voidage_conversion_coeffs.begin() + np * w); - } - } - } - - - - - - void applyVREPGroupControl(WellState& well_state) const - { - if ( wellCollection()->havingVREPGroups() ) { - std::vector well_voidage_rates; - std::vector voidage_conversion_coeffs; - computeWellVoidageRates(well_state, well_voidage_rates, voidage_conversion_coeffs); - wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs); - - // for the wells under group control, update the currentControls for the well_state - for (const WellNode* well_node : wellCollection()->getLeafNodes()) { - if (well_node->isInjector() && !well_node->individualControl()) { - const int well_index = well_node->selfIndex(); - well_state.currentControls()[well_index] = well_node->groupControlIndex(); - } - } - } - } + std::vector& voidage_conversion_coeffs) const; + void applyVREPGroupControl(WellState& well_state) const; protected: @@ -1655,210 +302,17 @@ enum WellVariablePositions { double dWellFractionMax() const {return param_.dwell_fraction_max_; } // protected methods - 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.setValue(target_rate); - return bhp; - } else if (well_controls_get_current_type(wc) == THP) { - const int control = well_controls_get_current(wc); - const double thp = well_controls_get_current_target(wc); - const double alq = well_controls_iget_alq(wc, control); - const int table_id = well_controls_iget_vfp(wc, control); - EvalWell aqua = 0.0; - EvalWell liquid = 0.0; - EvalWell vapour = 0.0; - EvalWell bhp = 0.0; - double vfp_ref_depth = 0.0; + EvalWell getBhp(const int wellIdx) const; - const Opm::PhaseUsage& pu = phase_usage_; + EvalWell getQs(const int wellIdx, const int phaseIdx) const; - if (active_[ Water ]) { - aqua = getQs(wellIdx, pu.phase_pos[ Water]); - } - if (active_[ Oil ]) { - liquid = getQs(wellIdx, pu.phase_pos[ Oil ]); - } - if (active_[ Gas ]) { - vapour = getQs(wellIdx, pu.phase_pos[ Gas ]); - } - if (wells().type[wellIdx] == INJECTOR) { - bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp); - vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); - } else { - bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq); - vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); - } + EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const; - // pick the density in the top layer - const int perf = wells().well_connpos[wellIdx]; - const double rho = well_perforation_densities_[perf]; - const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_); - bhp -= dp; - return bhp; + EvalWell wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const; - } - const int nw = wells().number_of_wells; - return wellVariables_[nw*XvarWell + wellIdx]; - } - - EvalWell getQs(const int wellIdx, const int phaseIdx) const { - EvalWell qs = 0.0; - const WellControls* wc = wells().ctrls[wellIdx]; - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - 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 || well_controls_get_current_type(wc) == THP) { - return wellVariables_[nw*XvarWell + wellIdx]; - } - qs.setValue(target_rate); - return qs; - } - - // Producers - if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { - return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx); - } - - if (well_controls_get_current_type(wc) == SURFACE_RATE) { - // checking how many phases are included in the rate control - // to decide wheter it is a single phase rate control or not - const double* distr = well_controls_get_current_distr(wc); - int num_phases_under_rate_control = 0; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - num_phases_under_rate_control += 1; - } - } - - // there should be at least one phase involved - assert(num_phases_under_rate_control > 0); - - // when it is a single phase rate limit - if (num_phases_under_rate_control == 1) { - if (distr[phaseIdx] == 1.0) { - qs.setValue(target_rate); - return qs; - } - - int currentControlIdx = 0; - for (int i = 0; i < np; ++i) { - currentControlIdx += wells().comp_frac[np*wellIdx + i] * i; - } - - const double eps = 1e-6; - if (wellVolumeFractionScaled(wellIdx,currentControlIdx) < eps) { - return qs; - } - return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / wellVolumeFractionScaled(wellIdx,currentControlIdx)); - } - - // when it is a combined two phase rate limit, such like LRAT - // we neec to calculate the rate for the certain phase - if (num_phases_under_rate_control == 2) { - EvalWell combined_volume_fraction = 0.; - for (int p = 0; p < np; ++p) { - if (distr[p] == 1.0) { - combined_volume_fraction += wellVolumeFractionScaled(wellIdx, p); - } - } - return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / combined_volume_fraction); - } - - // suppose three phase combined limit is the same with RESV - // not tested yet. - } - // ReservoirRate - return target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx); - } - - EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const { - const int nw = wells().number_of_wells; - if (phaseIdx == Water) { - return wellVariables_[WFrac * nw + wellIdx]; - } - - if (phaseIdx == Gas) { - return wellVariables_[GFrac * nw + wellIdx]; - } - - // Oil fraction - EvalWell well_fraction = 1.0; - if (active_[Water]) { - well_fraction -= wellVariables_[WFrac * nw + wellIdx]; - } - - if (active_[Gas]) { - well_fraction -= wellVariables_[GFrac * nw + wellIdx]; - } - return well_fraction; - } - - 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 Opm::PhaseUsage& pu = phase_usage_; - 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; - } - + const int well_number) const; using WellMapType = typename WellState::WellMapType; using WellMapEntryType = typename WellState::mapentry_t; @@ -1877,300 +331,23 @@ enum WellVariablePositions { }; - template RatioCheckTuple 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. + const WellMapEntryType& map_entry) const; - 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 - { - bool water_cut_limit_violated = false; - int worst_offending_connection = INVALIDCONNECTION; - bool last_connection = false; - double violation_extent = -1.0; + const WellMapEntryType& map_entry) const; - const int np = well_state.numPhases(); - const Opm::PhaseUsage& pu = phase_usage_; - 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); - } - - - - - - template void updateWellStateWithTarget(const WellControls* wc, const int current, const int well_index, - WellState& xw) const - { - // number of phases - const int np = wells().number_of_phases; - // Updating well state and primary variables. - // 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()[well_index] = target; - break; - - case THP: { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - const Opm::PhaseUsage& pu = phase_usage_; - - if (active_[ Water ]) { - aqua = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ]; - } - if (active_[ Oil ]) { - liquid = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ]; - } - if (active_[ Gas ]) { - vapour = xw.wellRates()[well_index*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[well_index]; - - // pick the density in the top layer - const int perf = wells().well_connpos[well_index]; - const double rho = well_perforation_densities_[perf]; - - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - rho, gravity_); - - xw.bhp()[well_index] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - } - else if (well_type == PRODUCER) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - rho, gravity_); - - xw.bhp()[well_index] = 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[well_index]; - if (well_type == INJECTOR) { - for (int phase = 0; phase < np; ++phase) { - const double& compi = wells().comp_frac[np * well_index + phase]; - // TODO: it was commented out from the master branch already. - //if (compi > 0.0) { - xw.wellRates()[np*well_index + 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*well_index + phase] = target * distr[phase]; - } - } - } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - - break; - } // end of switch - - - std::vector g = {1.0, 1.0, 0.01}; - if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { - for (int phase = 0; phase < np; ++phase) { - g[phase] = distr[phase]; - } - } - - // the number of wells - const int nw = wells().number_of_wells; - - switch (well_controls_iget_type(wc, current)) { - case THP: - case BHP: { - const WellType& well_type = wells().type[well_index]; - xw.wellSolutions()[nw*XvarWell + well_index] = 0.0; - if (well_type == INJECTOR) { - for (int p = 0; p < np; ++p) { - xw.wellSolutions()[nw*XvarWell + well_index] += xw.wellRates()[np*well_index + p] * wells().comp_frac[np*well_index + p]; - } - } else { - for (int p = 0; p < np; ++p) { - xw.wellSolutions()[nw*XvarWell + well_index] += g[p] * xw.wellRates()[np*well_index + p]; - } - } - break; - } - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - xw.wellSolutions()[nw*XvarWell + well_index] = xw.bhp()[well_index]; - break; - } // end of switch - - double tot_well_rate = 0.0; - for (int p = 0; p < np; ++p) { - tot_well_rate += g[p] * xw.wellRates()[np*well_index + p]; - } - if(std::abs(tot_well_rate) > 0) { - if (active_[ Water ]) { - xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate; - } - if (active_[ Gas ]) { - xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; - } - } else { - if (active_[ Water ]) { - xw.wellSolutions()[WFrac*nw + well_index] = wells().comp_frac[np*well_index + Water]; - } - - if (active_[ Gas ]) { - xw.wellSolutions()[GFrac*nw + well_index] = wells().comp_frac[np*well_index + Gas]; - } - } - - } + WellState& xw) const; }; } // namespace Opm + +#include "StandardWellsDense_impl.hpp" #endif diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp new file mode 100644 index 000000000..8c6aae79d --- /dev/null +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -0,0 +1,2370 @@ + + + +namespace Opm { + + template + StandardWellsDense:: + StandardWellsDense(const Wells* wells_arg, + WellCollection* well_collection, + const ModelParameters& param, + const bool terminal_output) + : wells_active_(wells_arg!=nullptr) + , wells_(wells_arg) + , well_collection_(well_collection) + , param_(param) + , terminal_output_(terminal_output) + , well_perforation_efficiency_factors_((wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0), 1.0) + , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) + , well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) + , wellVariables_( wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0) + , F0_(wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0 ) + { + if( wells_ ) + { + invDuneD_.setBuildMode( Mat::row_wise ); + duneC_.setBuildMode( Mat::row_wise ); + duneB_.setBuildMode( Mat::row_wise ); + } + } + + + + + + template + void + StandardWellsDense:: + init(const PhaseUsage phase_usage_arg, + const std::vector& active_arg, + const VFPProperties* vfp_properties_arg, + const double gravity_arg, + const std::vector& depth_arg, + const std::vector& pv_arg, + const RateConverterType* rate_converter) + { + if ( ! localWellsActive() ) { + return; + } + + phase_usage_ = phase_usage_arg; + active_ = active_arg; + vfp_properties_ = vfp_properties_arg; + gravity_ = gravity_arg; + cell_depths_ = extractPerfData(depth_arg); + pv_ = pv_arg; + rate_converter_ = rate_converter; + + calculateEfficiencyFactors(); + + // setup sparsity pattern for the matrices + //[A B^T [x = [ res + // C D] x_well] res_well] + + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + const int nc = numCells(); + +#ifndef NDEBUG + const auto pu = phase_usage_; + const int np = pu.num_phases; + + // assumes the gas fractions are stored after water fractions + // WellVariablePositions needs to be changed for 2p runs + assert (np == 3 || (np == 2 && !pu.phase_used[Gas]) ); +#endif + + // set invDuneD + invDuneD_.setSize( nw, nw, nw ); + + // set duneC + duneC_.setSize( nw, nc, nperf ); + + // set duneB + duneB_.setSize( nw, nc, nperf ); + + for (auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) { + // Add nonzeros for diagonal + row.insert(row.index()); + } + + for (auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) { + // Add nonzeros for diagonal + for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { + const int cell_idx = wells().well_cells[perf]; + row.insert(cell_idx); + } + } + + // make the B^T matrix + for (auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) { + for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { + const int cell_idx = wells().well_cells[perf]; + row.insert(cell_idx); + } + } + + resWell_.resize( nw ); + + // resize temporary class variables + Cx_.resize( duneC_.N() ); + invDrw_.resize( invDuneD_.N() ); + } + + + + + + + template + template + SimulatorReport + StandardWellsDense:: + assemble(Simulator& ebosSimulator, + const int iterationIdx, + const double dt, + WellState& well_state) + { + SimulatorReport report; + if ( ! localWellsActive() ) { + return report; + } + + if (param_.compute_well_potentials_) { + computeWellPotentials(ebosSimulator, well_state); + } + + resetWellControlFromState(well_state); + updateWellControls(well_state); + // Set the primary variables for the wells + setWellVariables(well_state); + + if (iterationIdx == 0) { + computeWellConnectionPressures(ebosSimulator, well_state); + computeAccumWells(); + } + + if (param_.solve_welleq_initially_ && iterationIdx == 0) { + // solve the well equations as a pre-processing step + report = solveWellEq(ebosSimulator, dt, well_state); + } + assembleWellEq(ebosSimulator, dt, well_state, false); + + report.converged = true; + return report; + } + + + + + + template + template + void + StandardWellsDense:: + assembleWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells) + { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + // clear all entries + duneB_ = 0.0; + duneC_ = 0.0; + invDuneD_ = 0.0; + resWell_ = 0.0; + + auto& ebosJac = ebosSimulator.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + const double volume = 0.002831684659200; // 0.1 cu ft; + for (int w = 0; w < nw; ++w) { + bool allow_cf = allow_cross_flow(w, ebosSimulator); + const EvalWell bhp = getBhp(w); + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + + const int cell_idx = wells().well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + std::vector cq_s(np,0.0); + computeWellFlux(w, wells().WI[perf], intQuants, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); + + for (int p1 = 0; p1 < np; ++p1) { + + // the cq_s entering mass balance equations need to consider the efficiency factors. + const EvalWell cq_s_effective = cq_s[p1] * well_perforation_efficiency_factors_[perf]; + + if (!only_wells) { + // subtract sum of phase fluxes in the reservoir equation. + // need to consider the efficiency factor + ebosResid[cell_idx][flowPhaseToEbosCompIdx(p1)] -= cq_s_effective.value(); + } + + // subtract sum of phase fluxes in the well equations. + resWell_[w][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].value(); + + // assemble the jacobians + for (int p2 = 0; p2 < np; ++p2) { + if (!only_wells) { + // also need to consider the efficiency factor when manipulating the jacobians. + ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s_effective.derivative(p2); + duneB_[w][cell_idx][flowToEbosPvIdx(p2)][flowPhaseToEbosCompIdx(p1)] -= cq_s_effective.derivative(p2+blocksize); // intput in transformed matrix + duneC_[w][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s_effective.derivative(p2); + } + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2+blocksize); + } + + // add trivial equation for 2p cases (Only support water + oil) + if (np == 2) { + assert(!active_[ Gas ]); + invDuneD_[w][w][flowPhaseToEbosCompIdx(Gas)][flowToEbosPvIdx(Gas)] = 1.0; + } + + // Store the perforation phase flux for later usage. + well_state.perfPhaseRates()[perf*np + p1] = cq_s[p1].value(); + } + + // Store the perforation pressure for later usage. + well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf]; + } + + // add vol * dF/dt + Q to the well equations; + for (int p1 = 0; p1 < np; ++p1) { + EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; + resWell_loc += getQs(w, p1); + for (int p2 = 0; p2 < np; ++p2) { + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivative(p2+blocksize); + } + resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value(); + } + } + + // do the local inversion of D. + localInvert( invDuneD_ ); + } + + + + + + template + template + bool + StandardWellsDense:: + allow_cross_flow(const int w, Simulator& ebosSimulator) const + { + if (wells().allow_cf[w]) { + return true; + } + + // check for special case where all perforations have cross flow + // then the wells must allow for cross flow + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + const int cell_idx = wells().well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + EvalWell bhp = getBhp(w); + + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + wellPerforationPressureDiffs()[perf]; + EvalWell drawdown = pressure - well_pressure; + + if (drawdown.value() < 0 && wells().type[w] == INJECTOR) { + return false; + } + + if (drawdown.value() > 0 && wells().type[w] == PRODUCER) { + return false; + } + } + return true; + } + + + + + + template + void + StandardWellsDense:: + localInvert(Mat& istlA) const + { + for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { + for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { + //std::cout << (*col) << std::endl; + (*col).invert(); + } + } + } + + + + + + template + void + StandardWellsDense:: + print(Mat& istlA) const + { + for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { + for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { + std::cout << row.index() << " " << col.index() << "/n \n"<<(*col) << std::endl; + } + } + } + + + + + + template + void + StandardWellsDense:: + apply( BVector& r) const + { + if ( ! localWellsActive() ) { + return; + } + + assert( invDrw_.size() == invDuneD_.N() ); + + invDuneD_.mv(resWell_,invDrw_); + duneB_.mmtv(invDrw_, r); + } + + + + + + template + void + StandardWellsDense:: + apply(const BVector& x, BVector& Ax) + { + if ( ! localWellsActive() ) { + return; + } + + assert( Cx_.size() == duneC_.N() ); + + BVector& invDCx = invDrw_; + assert( invDCx.size() == invDuneD_.N()); + + duneC_.mv(x, Cx_); + invDuneD_.mv(Cx_, invDCx); + duneB_.mmtv(invDCx,Ax); + } + + + + + + template + void + StandardWellsDense:: + applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) + { + if ( ! localWellsActive() ) { + return; + } + + if( scaleAddRes_.size() != Ax.size() ) { + scaleAddRes_.resize( Ax.size() ); + } + + scaleAddRes_ = 0.0; + apply( x, scaleAddRes_ ); + Ax.axpy( alpha, scaleAddRes_ ); + } + + + + + + template + void + StandardWellsDense:: + recoverVariable(const BVector& x, BVector& xw) const + { + if ( ! localWellsActive() ) { + return; + } + BVector resWell = resWell_; + duneC_.mmv(x, resWell); + invDuneD_.mv(resWell, xw); + } + + + + + + template + int + StandardWellsDense:: + flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; + return phaseToComp[ phaseIdx ]; + } + + + + + + template + int + StandardWellsDense:: + flowToEbosPvIdx( const int flowPv ) const + { + const int flowToEbos[ 3 ] = { + BlackoilIndices::pressureSwitchIdx, + BlackoilIndices::waterSaturationIdx, + BlackoilIndices::compositionSwitchIdx + }; + return flowToEbos[ flowPv ]; + } + + + + + + template + int + StandardWellsDense:: + flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; + return flowToEbos[ phaseIdx ]; + } + + + + + + template + int + StandardWellsDense:: + ebosCompToFlowPhaseIdx( const int compIdx ) const + { + const int compToPhase[ 3 ] = { Oil, Water, Gas }; + return compToPhase[ compIdx ]; + } + + + + + + template + std::vector + StandardWellsDense:: + extractPerfData(const std::vector& in) const + { + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + std::vector out(nperf); + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + const int well_idx = wells().well_cells[perf]; + out[perf] = in[well_idx]; + } + } + return out; + } + + + + + + template + int + StandardWellsDense:: + numPhases() const + { + return wells().number_of_phases; + } + + + + + + template + int + StandardWellsDense:: + numCells() const + { + return pv_.size(); + } + + + + + + template + void + StandardWellsDense:: + resetWellControlFromState(WellState xw) const + { + const int nw = wells_->number_of_wells; + for (int w = 0; w < nw; ++w) { + WellControls* wc = wells_->ctrls[w]; + well_controls_set_current( wc, xw.currentControls()[w]); + } + } + + + + + + template + const Wells& + StandardWellsDense:: + wells() const + { + assert(wells_ != 0); + return *(wells_); + } + + + + + + template + const Wells* + StandardWellsDense:: + wellsPointer() const + { + return wells_; + } + + + + + + template + bool + StandardWellsDense:: + wellsActive() const + { + return wells_active_; + } + + + + + + template + void + StandardWellsDense:: + setWellsActive(const bool wells_active) + { + wells_active_ = wells_active; + } + + + + + + template + bool + StandardWellsDense:: + localWellsActive() const + { + return wells_ ? (wells_->number_of_wells > 0 ) : false; + } + + + + + + template + 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; + } + + + + + + template + const std::vector& + StandardWellsDense:: + wellPerforationDensities() const + { + return well_perforation_densities_; + } + + + + + + template + const std::vector& + StandardWellsDense:: + wellPerforationPressureDiffs() const + { + return well_perforation_pressure_diffs_; + } + + + + + + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + extendEval(Eval in) const { + EvalWell out = 0.0; + out.setValue(in.value()); + for(int i = 0; i < blocksize;++i) { + out.setDerivative(i, in.derivative(flowToEbosPvIdx(i))); + } + return out; + } + + + + + + template + void + StandardWellsDense:: + 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].setValue(xw.wellSolutions()[w + nw* phaseIdx]); + wellVariables_[w + nw*phaseIdx].setDerivative(blocksize + phaseIdx, 1.0); + } + } + } + + + + + + template + void + StandardWellsDense:: + print(EvalWell in) const + { + std::cout << in.value() << std::endl; + for (int i = 0; i < in.size; ++i) { + std::cout << in.derivative(i) << std::endl; + } + } + + + + + + template + void + StandardWellsDense:: + 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 + template + void + StandardWellsDense:: + computeWellFlux(const int& w, const double& Tw, + const intensiveQuants& intQuants, + const EvalWell& bhp, const double& cdp, + const bool& allow_cf, std::vector& cq_s) const + { + const Opm::PhaseUsage& pu = phase_usage_; + const int np = wells().number_of_phases; + std::vector cmix_s(np,0.0); + for (int phase = 0; phase < np; ++phase) { + //int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + cmix_s[phase] = wellVolumeFraction(w,phase); + } + + const auto& fs = intQuants.fluidState(); + 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; + EvalWell drawdown = pressure - well_pressure; + + // injection perforations + if ( drawdown.value() > 0 ) { + //Do nothing if crossflow is not allowed + if (!allow_cf && wells().type[w] == INJECTOR) { + return; + } + + // 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 * (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[phase] = cq_ps[phase]; + } + } else { + //Do nothing if crossflow is not allowed + if (!allow_cf && wells().type[w] == PRODUCER) { + return; + } + + // 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 * (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.setValue(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 " < + template + SimulatorReport + StandardWellsDense:: + solveWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state) + { + const int nw = wells().number_of_wells; + WellState well_state0 = well_state; + + int it = 0; + bool converged; + do { + assembleWellEq(ebosSimulator, dt, well_state, true); + converged = getWellConvergence(ebosSimulator, it); + + // checking whether the group targets are converged + if (wellCollection()->groupControlActive()) { + converged = converged && wellCollection()->groupTargetConverged(well_state.wellRates()); + } + + if (converged) { + break; + } + + ++it; + if( localWellsActive() ) + { + BVector dx_well (nw); + invDuneD_.mv(resWell_, dx_well); + + updateWellState(dx_well, well_state); + updateWellControls(well_state); + setWellVariables(well_state); + } + } while (it < 15); + + if (!converged) { + well_state = well_state0; + } + + SimulatorReport report; + report.converged = converged; + report.total_well_iterations = it; + return report; + } + + + + + + template + void + StandardWellsDense:: + printIf(const int c, const double x, const double y, const double eps, const std::string type) const + { + if (std::abs(x-y) > eps) { + std::cout << type << " " << c << ": "< + std::vector + StandardWellsDense:: + residual() const + { + if( ! wellsActive() ) + { + return std::vector(); + } + + const int np = numPhases(); + const int nw = wells().number_of_wells; + std::vector res(np*nw); + for( int p=0; p + template + bool + StandardWellsDense:: + getWellConvergence(Simulator& ebosSimulator, + const int iteration) const + { + typedef std::vector< double > Vector; + + const int np = numPhases(); + const int nc = numCells(); + const double tol_wells = param_.tolerance_wells_; + const double maxResidualAllowed = param_.max_residual_allowed_; + + Vector R_sum(np); + Vector B_avg(np); + Vector maxCoeff(np); + Vector maxNormWell(np); + + std::vector< Vector > B( np, Vector( nc ) ); + std::vector< Vector > R2( np, Vector( nc ) ); + std::vector< Vector > tempV( np, Vector( nc ) ); + + for ( int idx = 0; idx < np; ++idx ) + { + Vector& B_idx = B[ idx ]; + const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + + 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(); + + B_idx [cell_idx] = 1 / fs.invB(ebosPhaseIdx).value(); + } + } + + detail::convergenceReduction(B, tempV, R2, + R_sum, maxCoeff, B_avg, maxNormWell, + nc, np, pv_, residual() ); + + 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 + template + void + StandardWellsDense:: + 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); + computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, cell_depths_, gravity_); + } + + + + + + template + template + void + StandardWellsDense:: + 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 + { + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + const PhaseUsage& pu = phase_usage_; + const int np = phase_usage_.num_phases; + b_perf.resize(nperf*np); + surf_dens_perf.resize(nperf*np); + + //rs and rv are only used if both oil and gas is present + if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) { + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + } + + // Compute the average pressure in each well block + 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 = wells().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] : xw.perfPress()[perf - 1]; + const double p_avg = (xw.perfPress()[perf] + p_above)/2; + const 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]) { + const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); + if (oilrate > 0) { + const double gasrate = std::abs(xw.wellRates()[gaspos_well]); + double rv = 0.0; + if (gasrate > 0) { + rv = oilrate / gasrate; + } + rv = std::min(rv, rvmax_perf[perf]); + + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); + } + else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + } else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + if (pu.phase_used[BlackoilPhases::Vapour]) { + rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + const double gasrate = std::abs(xw.wellRates()[gaspos_well]); + if (gasrate > 0) { + const double oilrate = std::abs(xw.wellRates()[oilpos_well]); + double rs = 0.0; + if (oilrate > 0) { + rs = gasrate / oilrate; + } + rs = std::min(rs, rsmax_perf[perf]); + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } + + // Surface density. + for (int p = 0; p < pu.num_phases; ++p) { + surf_dens_perf[np*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); + } + } + } + } + + + + + + template + void + StandardWellsDense:: + updateWellState(const BVector& dwells, + WellState& well_state) const + { + if( !localWellsActive() ) return; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + double dFLimit = dWellFractionMax(); + double dBHPLimit = dbhpMaxRel(); + std::vector xvar_well_old = well_state.wellSolutions(); + + for (int w = 0; w < nw; ++w) { + + // update the second and third well variable (The flux fractions) + std::vector F(np,0.0); + if (active_[ Water ]) { + const int sign2 = dwells[w][flowPhaseToEbosCompIdx(WFrac)] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(WFrac)]),dFLimit); + well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited; + } + + if (active_[ Gas ]) { + const int sign3 = dwells[w][flowPhaseToEbosCompIdx(GFrac)] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(GFrac)]),dFLimit); + well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited; + } + + assert(active_[ Oil ]); + F[Oil] = 1.0; + if (active_[ Water ]) { + F[Water] = well_state.wellSolutions()[WFrac*nw + w]; + F[Oil] -= F[Water]; + } + + if (active_[ Gas ]) { + F[Gas] = well_state.wellSolutions()[GFrac*nw + w]; + F[Oil] -= F[Gas]; + } + + if (active_[ Water ]) { + if (F[Water] < 0.0) { + if (active_[ Gas ]) { + F[Gas] /= (1.0 - F[Water]); + } + F[Oil] /= (1.0 - F[Water]); + F[Water] = 0.0; + } + } + if (active_[ Gas ]) { + if (F[Gas] < 0.0) { + if (active_[ Water ]) { + F[Water] /= (1.0 - F[Gas]); + } + F[Oil] /= (1.0 - F[Gas]); + F[Gas] = 0.0; + } + } + if (F[Oil] < 0.0) { + if (active_[ Water ]) { + F[Water] /= (1.0 - F[Oil]); + } + if (active_[ Gas ]) { + F[Gas] /= (1.0 - F[Oil]); + } + F[Oil] = 0.0; + } + + if (active_[ Water ]) { + well_state.wellSolutions()[WFrac*nw + w] = F[Water]; + } + if (active_[ Gas ]) { + well_state.wellSolutions()[GFrac*nw + w] = F[Gas]; + } + + // The interpretation of the first well variable depends on the well control + 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 = well_state.currentControls()[w]; + const double target_rate = well_controls_iget_target(wc, current); + + 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 p = 0; p < np; ++p) { + F[p] /= distr[p]; + } + } else { + for (int p = 0; p < np; ++p) { + F[p] /= g[p]; + } + } + + switch (well_controls_iget_type(wc, current)) { + case THP: // The BHP and THP both uses the total rate as first well variable. + case BHP: + { + well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][flowPhaseToEbosCompIdx(XvarWell)]; + + switch (wells().type[w]) { + case INJECTOR: + for (int p = 0; p < np; ++p) { + const double comp_frac = wells().comp_frac[np*w + p]; + well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[nw*XvarWell + w]; + } + break; + case PRODUCER: + for (int p = 0; p < np; ++p) { + well_state.wellRates()[w*np + p] = well_state.wellSolutions()[nw*XvarWell + w] * F[p]; + } + break; + } + + if (well_controls_iget_type(wc, current) == THP) { + + // Calculate bhp from thp control and well rates + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = phase_usage_; + + 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(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]; + // pick the density in the top layer + const int perf = wells().well_connpos[w]; + const double rho = well_perforation_densities_[perf]; + + if (well_type == INJECTOR) { + const double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + well_state.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + const double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + well_state.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + } + } + break; + case SURFACE_RATE: // Both rate controls use bhp as first well variable + case RESERVOIR_RATE: + { + const int sign1 = dwells[w][flowPhaseToEbosCompIdx(XvarWell)] > 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(XvarWell)]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit); + well_state.wellSolutions()[nw*XvarWell + w] = std::max(xvar_well_old[nw*XvarWell + w] - dx1_limited,1e5); + well_state.bhp()[w] = well_state.wellSolutions()[nw*XvarWell + w]; + + if (well_controls_iget_type(wc, current) == SURFACE_RATE) { + if (wells().type[w]==PRODUCER) { + + const double* distr = well_controls_iget_distr(wc, current); + + double F_target = 0.0; + for (int p = 0; p < np; ++p) { + F_target += distr[p] * F[p]; + } + for (int p = 0; p < np; ++p) { + well_state.wellRates()[np*w + p] = F[p] * target_rate / F_target; + } + } else { + + for (int p = 0; p < np; ++p) { + well_state.wellRates()[w*np + p] = wells().comp_frac[np*w + p] * target_rate; + } + } + } else { // RESERVOIR_RATE + for (int p = 0; p < np; ++p) { + well_state.wellRates()[np*w + p] = F[p] * target_rate; + } + } + } + break; + } // end of switch (well_controls_iget_type(wc, current)) + } // end of for (int w = 0; w < nw; ++w) + } + + + + + + template + void + StandardWellsDense:: + updateWellControls(WellState& xw) const + { + if( !localWellsActive() ) return ; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + // keeping a copy of the current controls, to see whether control changes later. + std::vector old_control_index(nw, 0); + for (int w = 0; w < nw; ++w) { + old_control_index[w] = xw.currentControls()[w]; + } + + // Find, for each well, if any constraints are broken. If so, + // switch control to first broken constraint. +#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. + xw.currentControls()[w] = ctrl_index; + current = xw.currentControls()[w]; + well_controls_set_current( wc, current); + } + + // update whether well is under group control + if (wellCollection()->groupControlActive()) { + // get well node in the well collection + WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w])); + + // update whehter the well is under group control or individual control + if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) { + // under group control + well_node.setIndividualControl(false); + } else { + // individual control + well_node.setIndividualControl(true); + } + } + } + + // upate the well targets following group controls + if (wellCollection()->groupControlActive()) { + applyVREPGroupControl(xw); + wellCollection()->updateWellTargets(xw.wellRates()); + } + + // the new well control indices after all the related updates, + std::vector updated_control_index(nw, 0); + for (int w = 0; w < nw; ++w) { + updated_control_index[w] = xw.currentControls()[w]; + } + + // checking whether control changed + wellhelpers::WellSwitchingLogger logger; + for (int w = 0; w < nw; ++w) { + if (updated_control_index[w] != old_control_index[w]) { + WellControls* wc = wells().ctrls[w]; + logger.wellSwitched(wells().name[w], + well_controls_iget_type(wc, old_control_index[w]), + well_controls_iget_type(wc, updated_control_index[w])); + updateWellStateWithTarget(wc, updated_control_index[w], w, xw); + } + } + } + + + + + + template + void + StandardWellsDense:: + updateListEconLimited(const Schedule& schedule, + const int current_step, + const Wells* wells_struct, + const WellState& well_state, + DynamicListEconLimited& list_econ_limited) const + { + // With no wells (on process) wells_struct is a null pointer + const int nw = (wells_struct)? wells_struct->number_of_wells : 0; + + 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); + } + } + + } // for (int w = 0; w < nw; ++w) + } + + + + + + 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 + well_perforation_densities_ = + WellDensitySegmented::computeConnectionDensities( + wells(), xw, phase_usage_, + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + // Compute pressure deltas + well_perforation_pressure_diffs_ = + WellDensitySegmented::computeConnectionPressureDelta( + wells(), depth_perf, well_perforation_densities_, grav); + } + + + + + + template + template + void + StandardWellsDense:: + computeWellPotentials(const Simulator& ebosSimulator, + WellState& well_state) const + { + + // number of wells and phases + const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; + + for (int w = 0; w < nw; ++w) { + // bhp needs to be determined for the well potential calculation + double bhp = 0.; + + const WellControls* well_control = wells().ctrls[w]; + // The number of the well controls + const int nwc = well_controls_get_num(well_control); + + // Finding a BHP control or a THP control + // IF we find a THP control, we calculate the BHP value. + // TODO: there is option to ignore the THP limit when calculating well potentials, + // we are not handling it for the moment. + for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(well_control, ctrl_index) == BHP) { + // set bhp to the bhp value + bhp = well_controls_iget_target(well_control, ctrl_index); + } + + if (well_controls_iget_type(well_control, ctrl_index) == THP) { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = phase_usage_; + + 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(well_control, ctrl_index); + const double& thp = well_controls_iget_target(well_control, ctrl_index); + const double& alq = well_controls_iget_alq(well_control, ctrl_index); + + // Calculating the BHP value based on THP + const WellType& well_type = wells().type[w]; + const int first_perf = wells().well_connpos[w]; //first perforation + + if (well_type == INJECTOR) { + const double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities()[first_perf], gravity_); + const double bhp_calculated = 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_calculated < bhp) { + bhp = bhp_calculated; + } + } + else if (well_type == PRODUCER) { + const double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities()[first_perf], gravity_); + const double bhp_calculated = 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_calculated > bhp) { + bhp = bhp_calculated; + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + } + } + + assert(bhp != 0.0); + + // Should we consider crossflow when calculating well potentionals? + const bool allow_cf = allow_cross_flow(w, ebosSimulator); + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const int cell_index = wells().well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0)); + std::vector well_potentials(np, 0.0); + computeWellFlux(w, wells().WI[perf], intQuants, bhp, wellPerforationPressureDiffs()[perf], allow_cf, well_potentials); + for(int p = 0; p < np; ++p) { + well_state.wellPotentials()[perf * np + p] = well_potentials[p].value(); + } + } + } // end of for (int w = 0; w < nw; ++w) + } + + + + + + template + WellCollection* + StandardWellsDense:: + wellCollection() const + { + return well_collection_; + } + + + + + template + const std::vector& + StandardWellsDense:: + wellPerfEfficiencyFactors() const + { + return well_perforation_efficiency_factors_; + } + + + + + + template + void + StandardWellsDense:: + calculateEfficiencyFactors() + { + if ( !localWellsActive() ) { + return; + } + + const int nw = wells().number_of_wells; + + for (int w = 0; w < nw; ++w) { + const std::string well_name = wells().name[w]; + const WellNode& well_node = wellCollection()->findWellNode(well_name); + + const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor(); + + // assign the efficiency factor to each perforation related. + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w + 1]; ++perf) { + well_perforation_efficiency_factors_[perf] = well_efficiency_factor; + } + } + } + + + + + + template + void + StandardWellsDense:: + computeWellVoidageRates(const WellState& well_state, + std::vector& well_voidage_rates, + std::vector& voidage_conversion_coeffs) const + { + if ( !localWellsActive() ) { + return; + } + // TODO: for now, we store the voidage rates for all the production wells. + // For injection wells, the rates are stored as zero. + // We only store the conversion coefficients for all the injection wells. + // Later, more delicate model will be implemented here. + // And for the moment, group control can only work for serial running. + const int nw = well_state.numWells(); + const int np = well_state.numPhases(); + + // we calculate the voidage rate for each well, that means the sum of all the phases. + well_voidage_rates.resize(nw, 0); + // store the conversion coefficients, while only for the use of injection wells. + voidage_conversion_coeffs.resize(nw * np, 1.0); + + std::vector well_rates(np, 0.0); + std::vector convert_coeff(np, 1.0); + + for (int w = 0; w < nw; ++w) { + const bool is_producer = wells().type[w] == PRODUCER; + + // not sure necessary to change all the value to be positive + if (is_producer) { + std::transform(well_state.wellRates().begin() + np * w, + well_state.wellRates().begin() + np * (w + 1), + well_rates.begin(), std::negate()); + + // the average hydrocarbon conditions of the whole field will be used + const int fipreg = 0; // Not considering FIP for the moment. + + rate_converter_->calcCoeff(well_rates, fipreg, convert_coeff); + well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(), + convert_coeff.begin(), 0.0); + } else { + // TODO: Not sure whether will encounter situation with all zero rates + // and whether it will cause problem here. + std::copy(well_state.wellRates().begin() + np * w, + well_state.wellRates().begin() + np * (w + 1), + well_rates.begin()); + // the average hydrocarbon conditions of the whole field will be used + const int fipreg = 0; // Not considering FIP for the moment. + rate_converter_->calcCoeff(well_rates, fipreg, convert_coeff); + std::copy(convert_coeff.begin(), convert_coeff.end(), + voidage_conversion_coeffs.begin() + np * w); + } + } + } + + + + + + template + void + StandardWellsDense:: + applyVREPGroupControl(WellState& well_state) const + { + if ( wellCollection()->havingVREPGroups() ) { + std::vector well_voidage_rates; + std::vector voidage_conversion_coeffs; + computeWellVoidageRates(well_state, well_voidage_rates, voidage_conversion_coeffs); + wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs); + + // for the wells under group control, update the currentControls for the well_state + for (const WellNode* well_node : wellCollection()->getLeafNodes()) { + if (well_node->isInjector() && !well_node->individualControl()) { + const int well_index = well_node->selfIndex(); + well_state.currentControls()[well_index] = well_node->groupControlIndex(); + } + } + } + } + + + + + + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + 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.setValue(target_rate); + return bhp; + } else if (well_controls_get_current_type(wc) == THP) { + const int control = well_controls_get_current(wc); + const double thp = well_controls_get_current_target(wc); + const double alq = well_controls_iget_alq(wc, control); + const int table_id = well_controls_iget_vfp(wc, control); + EvalWell aqua = 0.0; + EvalWell liquid = 0.0; + EvalWell vapour = 0.0; + EvalWell bhp = 0.0; + double vfp_ref_depth = 0.0; + + const Opm::PhaseUsage& pu = phase_usage_; + + if (active_[ Water ]) { + aqua = getQs(wellIdx, pu.phase_pos[ Water]); + } + if (active_[ Oil ]) { + liquid = getQs(wellIdx, pu.phase_pos[ Oil ]); + } + if (active_[ Gas ]) { + vapour = getQs(wellIdx, pu.phase_pos[ Gas ]); + } + if (wells().type[wellIdx] == INJECTOR) { + bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp); + vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); + } else { + bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq); + vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + } + + // pick the density in the top layer + const int perf = wells().well_connpos[wellIdx]; + const double rho = well_perforation_densities_[perf]; + const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_); + bhp -= dp; + return bhp; + } + + const int nw = wells().number_of_wells; + return wellVariables_[nw*XvarWell + wellIdx]; + } + + + + + + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + getQs(const int wellIdx, const int phaseIdx) const + { + EvalWell qs = 0.0; + const WellControls* wc = wells().ctrls[wellIdx]; + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + 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 || well_controls_get_current_type(wc) == THP) { + return wellVariables_[nw*XvarWell + wellIdx]; + } + qs.setValue(target_rate); + return qs; + } + + // Producers + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { + return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx); + } + + if (well_controls_get_current_type(wc) == SURFACE_RATE) { + // checking how many phases are included in the rate control + // to decide wheter it is a single phase rate control or not + const double* distr = well_controls_get_current_distr(wc); + int num_phases_under_rate_control = 0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + num_phases_under_rate_control += 1; + } + } + + // there should be at least one phase involved + assert(num_phases_under_rate_control > 0); + + // when it is a single phase rate limit + if (num_phases_under_rate_control == 1) { + if (distr[phaseIdx] == 1.0) { + qs.setValue(target_rate); + return qs; + } + + int currentControlIdx = 0; + for (int i = 0; i < np; ++i) { + currentControlIdx += wells().comp_frac[np*wellIdx + i] * i; + } + + const double eps = 1e-6; + if (wellVolumeFractionScaled(wellIdx,currentControlIdx) < eps) { + return qs; + } + return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / wellVolumeFractionScaled(wellIdx,currentControlIdx)); + } + + // when it is a combined two phase rate limit, such like LRAT + // we neec to calculate the rate for the certain phase + if (num_phases_under_rate_control == 2) { + EvalWell combined_volume_fraction = 0.; + for (int p = 0; p < np; ++p) { + if (distr[p] == 1.0) { + combined_volume_fraction += wellVolumeFractionScaled(wellIdx, p); + } + } + return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / combined_volume_fraction); + } + + // suppose three phase combined limit is the same with RESV + // not tested yet. + } + + // ReservoirRate + return target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx); + } + + + + + + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + wellVolumeFraction(const int wellIdx, const int phaseIdx) const + { + const int nw = wells().number_of_wells; + if (phaseIdx == Water) { + return wellVariables_[WFrac * nw + wellIdx]; + } + + if (phaseIdx == Gas) { + return wellVariables_[GFrac * nw + wellIdx]; + } + + // Oil fraction + EvalWell well_fraction = 1.0; + if (active_[Water]) { + well_fraction -= wellVariables_[WFrac * nw + wellIdx]; + } + + if (active_[Gas]) { + well_fraction -= wellVariables_[GFrac * nw + wellIdx]; + } + return well_fraction; + } + + + + + + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + 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 + StandardWellsDense:: + checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const int well_number) const + { + const Opm::PhaseUsage& pu = phase_usage_; + 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 + typename 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 + typename 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 = phase_usage_; + 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); + } + + + + + + template + void + StandardWellsDense:: + updateWellStateWithTarget(const WellControls* wc, + const int current, + const int well_index, + WellState& xw) const + { + // number of phases + const int np = wells().number_of_phases; + // Updating well state and primary variables. + // 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()[well_index] = target; + break; + + case THP: { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = phase_usage_; + + if (active_[ Water ]) { + aqua = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ]; + } + if (active_[ Oil ]) { + liquid = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ]; + } + if (active_[ Gas ]) { + vapour = xw.wellRates()[well_index*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[well_index]; + + // pick the density in the top layer + const int perf = wells().well_connpos[well_index]; + const double rho = well_perforation_densities_[perf]; + + if (well_type == INJECTOR) { + const double dp = wellhelpers::computeHydrostaticCorrection( + wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + xw.bhp()[well_index] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + const double dp = wellhelpers::computeHydrostaticCorrection( + wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + xw.bhp()[well_index] = 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[well_index]; + if (well_type == INJECTOR) { + for (int phase = 0; phase < np; ++phase) { + const double& compi = wells().comp_frac[np * well_index + phase]; + // TODO: it was commented out from the master branch already. + //if (compi > 0.0) { + xw.wellRates()[np*well_index + 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*well_index + phase] = target * distr[phase]; + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + + break; + } // end of switch + + + std::vector g = {1.0, 1.0, 0.01}; + if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { + for (int phase = 0; phase < np; ++phase) { + g[phase] = distr[phase]; + } + } + + // the number of wells + const int nw = wells().number_of_wells; + + switch (well_controls_iget_type(wc, current)) { + case THP: + case BHP: { + const WellType& well_type = wells().type[well_index]; + xw.wellSolutions()[nw*XvarWell + well_index] = 0.0; + if (well_type == INJECTOR) { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[nw*XvarWell + well_index] += xw.wellRates()[np*well_index + p] * wells().comp_frac[np*well_index + p]; + } + } else { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[nw*XvarWell + well_index] += g[p] * xw.wellRates()[np*well_index + p]; + } + } + break; + } + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + xw.wellSolutions()[nw*XvarWell + well_index] = xw.bhp()[well_index]; + break; + } // end of switch + + double tot_well_rate = 0.0; + for (int p = 0; p < np; ++p) { + tot_well_rate += g[p] * xw.wellRates()[np*well_index + p]; + } + if(std::abs(tot_well_rate) > 0) { + if (active_[ Water ]) { + xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate; + } + if (active_[ Gas ]) { + xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; + } + } else { + if (active_[ Water ]) { + xw.wellSolutions()[WFrac*nw + well_index] = wells().comp_frac[np*well_index + Water]; + } + + if (active_[ Gas ]) { + xw.wellSolutions()[GFrac*nw + well_index] = wells().comp_frac[np*well_index + Gas]; + } + } + + } + + +} // namespace Opm