mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
third part in separating the StandardWellsDense.hpp implementations.
This commit is contained in:
Kai Bao
parent
8de7795629
commit
2a3fe58ac2
@ -179,705 +179,42 @@ enum WellVariablePositions {
|
||||
|
||||
template<typename intensiveQuants>
|
||||
void
|
||||
computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants, const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector<EvalWell>& cq_s) const
|
||||
{
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
const int np = wells().number_of_phases;
|
||||
std::vector<EvalWell> 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<EvalWell> b_perfcells_dense(np, 0.0);
|
||||
std::vector<EvalWell> 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<EvalWell> 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 " <<tmp_oil << std::endl;
|
||||
volumeRatio += tmp_oil / b_perfcells_dense[oilpos];
|
||||
|
||||
const EvalWell tmp_gas = (cmix_s[gaspos] - rsPerf * cmix_s[oilpos]) / d;
|
||||
//std::cout << "tmp_gas " <<tmp_gas << std::endl;
|
||||
volumeRatio += tmp_gas / b_perfcells_dense[gaspos];
|
||||
}
|
||||
else {
|
||||
if (active_[Oil]) {
|
||||
const int oilpos = pu.phase_pos[Oil];
|
||||
volumeRatio += cmix_s[oilpos] / b_perfcells_dense[oilpos];
|
||||
}
|
||||
if (active_[Gas]) {
|
||||
const int gaspos = pu.phase_pos[Gas];
|
||||
volumeRatio += cmix_s[gaspos] / b_perfcells_dense[gaspos];
|
||||
}
|
||||
}
|
||||
// injecting connections total volumerates at standard conditions
|
||||
EvalWell cqt_is = cqt_i/volumeRatio;
|
||||
//std::cout << "volrat " << volumeRatio << " " << volrat_perf_[perf] << std::endl;
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
cq_s[phase] = cmix_s[phase] * cqt_is; // * b_perfcells_dense[phase];
|
||||
}
|
||||
}
|
||||
}
|
||||
computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants,
|
||||
const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector<EvalWell>& cq_s) const;
|
||||
|
||||
template <typename Simulator>
|
||||
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 << ": "<<x << " " << y << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
std::vector<double> residual() {
|
||||
if( ! wellsActive() )
|
||||
{
|
||||
return std::vector<double>();
|
||||
}
|
||||
|
||||
const int np = numPhases();
|
||||
const int nw = wells().number_of_wells;
|
||||
std::vector<double> res(np*nw);
|
||||
for( int p=0; p<np; ++p) {
|
||||
const int ebosCompIdx = flowPhaseToEbosCompIdx(p);
|
||||
for (int i = 0; i < nw; ++i) {
|
||||
int idx = i + nw*p;
|
||||
res[idx] = resWell_[ i ][ ebosCompIdx ];
|
||||
}
|
||||
}
|
||||
return res;
|
||||
}
|
||||
void printIf(const int c, const double x, const double y, const double eps, const std::string type) const;
|
||||
|
||||
std::vector<double> residual() const;
|
||||
|
||||
template <typename Simulator>
|
||||
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<typename Simulator>
|
||||
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<double> arguments, and not Eigen objects.
|
||||
std::vector<double> b_perf;
|
||||
std::vector<double> rsmax_perf;
|
||||
std::vector<double> rvmax_perf;
|
||||
std::vector<double> 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<typename Simulator, class WellState>
|
||||
template<typename Simulator>
|
||||
void
|
||||
computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
|
||||
const WellState& xw,
|
||||
std::vector<double>& b_perf,
|
||||
std::vector<double>& rsmax_perf,
|
||||
std::vector<double>& rvmax_perf,
|
||||
std::vector<double>& 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<double>& 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 <class WellState>
|
||||
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<double> 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<double> 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<double> 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 <class WellState>
|
||||
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<int> 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<int> 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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void updateWellControls(WellState& xw) const;
|
||||
|
||||
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
|
||||
{
|
||||
|
||||
@ -182,12 +182,12 @@ namespace Opm {
|
||||
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<EvalWell> 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) {
|
||||
@ -659,4 +659,767 @@ namespace Opm {
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
template<typename intensiveQuants>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
computeWellFlux(const int& w, const double& Tw,
|
||||
const intensiveQuants& intQuants,
|
||||
const EvalWell& bhp, const double& cdp,
|
||||
const bool& allow_cf, std::vector<EvalWell>& cq_s) const
|
||||
{
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
const int np = wells().number_of_phases;
|
||||
std::vector<EvalWell> 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<EvalWell> b_perfcells_dense(np, 0.0);
|
||||
std::vector<EvalWell> 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<EvalWell> 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 " <<tmp_oil << std::endl;
|
||||
volumeRatio += tmp_oil / b_perfcells_dense[oilpos];
|
||||
|
||||
const EvalWell tmp_gas = (cmix_s[gaspos] - rsPerf * cmix_s[oilpos]) / d;
|
||||
//std::cout << "tmp_gas " <<tmp_gas << std::endl;
|
||||
volumeRatio += tmp_gas / b_perfcells_dense[gaspos];
|
||||
}
|
||||
else {
|
||||
if (active_[Oil]) {
|
||||
const int oilpos = pu.phase_pos[Oil];
|
||||
volumeRatio += cmix_s[oilpos] / b_perfcells_dense[oilpos];
|
||||
}
|
||||
if (active_[Gas]) {
|
||||
const int gaspos = pu.phase_pos[Gas];
|
||||
volumeRatio += cmix_s[gaspos] / b_perfcells_dense[gaspos];
|
||||
}
|
||||
}
|
||||
|
||||
// injecting connections total volumerates at standard conditions
|
||||
EvalWell cqt_is = cqt_i/volumeRatio;
|
||||
//std::cout << "volrat " << volumeRatio << " " << volrat_perf_[perf] << std::endl;
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
cq_s[phase] = cmix_s[phase] * cqt_is; // * b_perfcells_dense[phase];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
template <typename Simulator>
|
||||
SimulatorReport
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
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<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
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 << ": "<<x << " " << y << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
std::vector<double>
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
residual() const
|
||||
{
|
||||
if( ! wellsActive() )
|
||||
{
|
||||
return std::vector<double>();
|
||||
}
|
||||
|
||||
const int np = numPhases();
|
||||
const int nw = wells().number_of_wells;
|
||||
std::vector<double> res(np*nw);
|
||||
for( int p=0; p<np; ++p) {
|
||||
const int ebosCompIdx = flowPhaseToEbosCompIdx(p);
|
||||
for (int i = 0; i < nw; ++i) {
|
||||
int idx = i + nw*p;
|
||||
res[idx] = resWell_[ i ][ ebosCompIdx ];
|
||||
}
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
template <typename Simulator>
|
||||
bool
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
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<typename FluidSystem, typename BlackoilIndices>
|
||||
template <typename Simulator>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
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<double> arguments, and not Eigen objects.
|
||||
std::vector<double> b_perf;
|
||||
std::vector<double> rsmax_perf;
|
||||
std::vector<double> rvmax_perf;
|
||||
std::vector<double> 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<typename FluidSystem, typename BlackoilIndices>
|
||||
template <typename Simulator>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
|
||||
const WellState& xw,
|
||||
std::vector<double>& b_perf,
|
||||
std::vector<double>& rsmax_perf,
|
||||
std::vector<double>& rvmax_perf,
|
||||
std::vector<double>& 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<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
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<double> 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<double> 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<double> 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<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
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<int> 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<int> 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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
Loading…
Reference in New Issue
Block a user