25 template<
typename TypeTag>
26 StandardWell<TypeTag>::
27 StandardWell(
const Well* well,
const int time_step,
const Wells* wells,
const ModelParameters& param)
28 : Base(well, time_step, wells, param)
29 , perf_densities_(number_of_perforations_)
30 , perf_pressure_diffs_(number_of_perforations_)
31 , primary_variables_(numWellEq, 0.0)
32 , primary_variables_evaluation_(numWellEq)
35 duneB_.setBuildMode( OffDiagMatWell::row_wise );
36 duneC_.setBuildMode( OffDiagMatWell::row_wise );
37 invDuneD_.setBuildMode( DiagMatWell::row_wise );
44 template<
typename TypeTag>
46 StandardWell<TypeTag>::
47 init(
const PhaseUsage* phase_usage_arg,
48 const std::vector<bool>* active_arg,
49 const std::vector<double>& depth_arg,
50 const double gravity_arg,
53 Base::init(phase_usage_arg, active_arg,
54 depth_arg, gravity_arg, num_cells);
56 perf_depth_.resize(number_of_perforations_, 0.);
57 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
58 const int cell_idx = well_cells_[perf];
59 perf_depth_[perf] = depth_arg[cell_idx];
66 invDuneD_.setSize(1, 1, 1);
67 duneB_.setSize(1, num_cells, number_of_perforations_);
68 duneC_.setSize(1, num_cells, number_of_perforations_);
70 for (
auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) {
72 row.insert(row.index());
75 for (
auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) {
76 for (
int perf = 0 ; perf < number_of_perforations_; ++perf) {
77 const int cell_idx = well_cells_[perf];
83 for (
auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) {
84 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
85 const int cell_idx = well_cells_[perf];
93 Bx_.resize( duneB_.N() );
94 invDrw_.resize( invDuneD_.N() );
101 template<
typename TypeTag>
102 void StandardWell<TypeTag>::
103 initPrimaryVariablesEvaluation()
const 108 for (
int eqIdx = 0; eqIdx < numComponents(); ++eqIdx) {
109 assert( (
size_t)eqIdx < primary_variables_.size() );
111 primary_variables_evaluation_[eqIdx] = 0.0;
112 primary_variables_evaluation_[eqIdx].setValue(primary_variables_[eqIdx]);
113 primary_variables_evaluation_[eqIdx].setDerivative(numEq + eqIdx, 1.0);
121 template<
typename TypeTag>
122 typename StandardWell<TypeTag>::EvalWell
123 StandardWell<TypeTag>::
126 const WellControls* wc = well_controls_;
127 if (well_controls_get_current_type(wc) == BHP) {
129 const double target_rate = well_controls_get_current_target(wc);
130 bhp.setValue(target_rate);
132 }
else if (well_controls_get_current_type(wc) == THP) {
133 const int control = well_controls_get_current(wc);
135 const Opm::PhaseUsage& pu = phaseUsage();
136 std::vector<EvalWell> rates(3, 0.0);
137 if (active()[ Water ]) {
138 rates[ Water ]= getQs(pu.phase_pos[ Water]);
140 if (active()[ Oil ]) {
141 rates[ Oil ] = getQs(pu.phase_pos[ Oil ]);
143 if (active()[ Gas ]) {
144 rates[ Gas ] = getQs(pu.phase_pos[ Gas ]);
146 return calculateBhpFromThp(rates, control);
149 return primary_variables_evaluation_[XvarWell];
156 template<
typename TypeTag>
157 typename StandardWell<TypeTag>::EvalWell
158 StandardWell<TypeTag>::
159 getQs(
const int comp_idx)
const 163 const WellControls* wc = well_controls_;
164 const int np = number_of_phases_;
165 const double target_rate = well_controls_get_current_target(wc);
167 assert(comp_idx < numComponents());
168 const auto pu = phaseUsage();
172 if (well_type_ == INJECTOR) {
176 double comp_frac = 0.0;
177 if (has_solvent && comp_idx == contiSolventEqIdx) {
178 comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent();
179 }
else if (comp_idx == pu.phase_pos[ Gas ]) {
180 comp_frac = comp_frac_[comp_idx] * (1.0 - wsolvent());
182 comp_frac = comp_frac_[comp_idx];
184 if (comp_frac == 0.0) {
188 if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
189 return comp_frac * primary_variables_evaluation_[XvarWell];
192 qs.setValue(comp_frac * target_rate);
196 const double comp_frac = comp_frac_[comp_idx];
197 if (comp_frac == 0.0) {
201 if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
202 return primary_variables_evaluation_[XvarWell];
204 qs.setValue(target_rate);
209 if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) {
210 return primary_variables_evaluation_[XvarWell] * wellVolumeFractionScaled(comp_idx);
213 if (well_controls_get_current_type(wc) == SURFACE_RATE) {
216 const double* distr = well_controls_get_current_distr(wc);
217 int num_phases_under_rate_control = 0;
218 for (
int phase = 0; phase < np; ++phase) {
219 if (distr[phase] > 0.0) {
220 num_phases_under_rate_control += 1;
225 assert(num_phases_under_rate_control > 0);
228 if (num_phases_under_rate_control == 1) {
231 int phase_under_control = -1;
232 for (
int phase = 0; phase < np; ++phase) {
233 if (distr[phase] > 0.0) {
234 phase_under_control = phase;
239 assert(phase_under_control >= 0);
241 EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(phase_under_control);
242 if (has_solvent && phase_under_control == Gas) {
244 wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(contiSolventEqIdx);
247 if (comp_idx == phase_under_control) {
248 if (has_solvent && phase_under_control == Gas) {
249 qs.setValue(target_rate * wellVolumeFractionScaled(Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() );
252 qs.setValue(target_rate);
257 const double eps = 1e-6;
258 if (wellVolumeFractionScaledPhaseUnderControl < eps) {
261 return (target_rate * wellVolumeFractionScaled(comp_idx) / wellVolumeFractionScaledPhaseUnderControl);
266 if (num_phases_under_rate_control == 2) {
267 EvalWell combined_volume_fraction = 0.;
268 for (
int p = 0; p < np; ++p) {
269 if (distr[p] == 1.0) {
270 combined_volume_fraction += wellVolumeFractionScaled(p);
273 return (target_rate * wellVolumeFractionScaled(comp_idx) / combined_volume_fraction);
277 if (num_phases_under_rate_control == 3) {
278 return target_rate * wellSurfaceVolumeFraction(comp_idx);
280 }
else if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
282 return target_rate * wellVolumeFractionScaled(comp_idx);
284 OPM_THROW(std::logic_error,
"Unknown control type for well " << name());
296 template<
typename TypeTag>
297 typename StandardWell<TypeTag>::EvalWell
298 StandardWell<TypeTag>::
299 wellVolumeFractionScaled(
const int compIdx)
const 302 const double scal = scalingFactor(compIdx);
304 return wellVolumeFraction(compIdx) / scal;
307 return wellVolumeFraction(compIdx);
314 template<
typename TypeTag>
315 typename StandardWell<TypeTag>::EvalWell
316 StandardWell<TypeTag>::
317 wellVolumeFraction(
const int compIdx)
const 319 const auto& pu = phaseUsage();
320 if (active()[Water] && compIdx == pu.phase_pos[Water]) {
321 return primary_variables_evaluation_[WFrac];
324 if (active()[Gas] && compIdx == pu.phase_pos[Gas]) {
325 return primary_variables_evaluation_[GFrac];
328 if (has_solvent && compIdx == contiSolventEqIdx) {
329 return primary_variables_evaluation_[SFrac];
333 EvalWell well_fraction = 1.0;
334 if (active()[Water]) {
335 well_fraction -= primary_variables_evaluation_[WFrac];
339 well_fraction -= primary_variables_evaluation_[GFrac];
342 well_fraction -= primary_variables_evaluation_[SFrac];
344 return well_fraction;
351 template<
typename TypeTag>
352 typename StandardWell<TypeTag>::EvalWell
353 StandardWell<TypeTag>::
354 wellSurfaceVolumeFraction(
const int compIdx)
const 356 EvalWell sum_volume_fraction_scaled = 0.;
357 const int numComp = numComponents();
358 for (
int idx = 0; idx < numComp; ++idx) {
359 sum_volume_fraction_scaled += wellVolumeFractionScaled(idx);
362 assert(sum_volume_fraction_scaled.value() != 0.);
364 return wellVolumeFractionScaled(compIdx) / sum_volume_fraction_scaled;
371 template<
typename TypeTag>
372 typename StandardWell<TypeTag>::EvalWell
373 StandardWell<TypeTag>::
374 extendEval(
const Eval& in)
const 377 out.setValue(in.value());
378 for(
int eqIdx = 0; eqIdx < numEq;++eqIdx) {
379 out.setDerivative(eqIdx, in.derivative(eqIdx));
388 template<
typename TypeTag>
390 StandardWell<TypeTag>::
391 computePerfRate(
const IntensiveQuantities& intQuants,
392 const std::vector<EvalWell>& mob_perfcells_dense,
393 const double Tw,
const EvalWell& bhp,
const double& cdp,
394 const bool& allow_cf, std::vector<EvalWell>& cq_s)
const 396 const Opm::PhaseUsage& pu = phaseUsage();
397 const int np = number_of_phases_;
398 const int numComp = numComponents();
399 std::vector<EvalWell> cmix_s(numComp,0.0);
400 for (
int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
401 cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
403 const auto& fs = intQuants.fluidState();
405 const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
406 const EvalWell rs = extendEval(fs.Rs());
407 const EvalWell rv = extendEval(fs.Rv());
408 std::vector<EvalWell> b_perfcells_dense(numComp, 0.0);
409 for (
int phase = 0; phase < np; ++phase) {
410 const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
411 b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx));
414 b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
418 const EvalWell well_pressure = bhp + cdp;
419 const EvalWell drawdown = pressure - well_pressure;
422 if ( drawdown.value() > 0 ) {
424 if (!allow_cf && well_type_ == INJECTOR) {
429 for (
int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
430 const EvalWell cq_p = - Tw * (mob_perfcells_dense[componentIdx] * drawdown);
431 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
434 if (active()[Oil] && active()[Gas]) {
435 const int oilpos = pu.phase_pos[Oil];
436 const int gaspos = pu.phase_pos[Gas];
437 const EvalWell cq_sOil = cq_s[oilpos];
438 const EvalWell cq_sGas = cq_s[gaspos];
439 cq_s[gaspos] += rs * cq_sOil;
440 cq_s[oilpos] += rv * cq_sGas;
445 if (!allow_cf && well_type_ == PRODUCER) {
450 EvalWell total_mob_dense = mob_perfcells_dense[0];
451 for (
int componentIdx = 1; componentIdx < numComp; ++componentIdx) {
452 total_mob_dense += mob_perfcells_dense[componentIdx];
456 const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown);
459 EvalWell volumeRatio = 0.0;
460 if (active()[Water]) {
461 const int watpos = pu.phase_pos[Water];
462 volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos];
466 volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx];
469 if (active()[Oil] && active()[Gas]) {
470 const int oilpos = pu.phase_pos[Oil];
471 const int gaspos = pu.phase_pos[Gas];
474 const EvalWell d = 1.0 - rv * rs;
476 if (d.value() == 0.0) {
477 OPM_THROW(Opm::NumericalProblem,
"Zero d value obtained for well " << name() <<
" during flux calcuation" 478 <<
" with rs " << rs <<
" and rv " << rv);
481 const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d;
483 volumeRatio += tmp_oil / b_perfcells_dense[oilpos];
485 const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d;
487 volumeRatio += tmp_gas / b_perfcells_dense[gaspos];
491 const int oilpos = pu.phase_pos[Oil];
492 volumeRatio += cmix_s[oilpos] / b_perfcells_dense[oilpos];
495 const int gaspos = pu.phase_pos[Gas];
496 volumeRatio += cmix_s[gaspos] / b_perfcells_dense[gaspos];
501 EvalWell cqt_is = cqt_i/volumeRatio;
503 for (
int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
504 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
513 template<
typename TypeTag>
515 StandardWell<TypeTag>::
516 assembleWellEq(Simulator& ebosSimulator,
518 WellState& well_state,
521 const int numComp = numComponents();
522 const int np = number_of_phases_;
532 auto& ebosJac = ebosSimulator.model().linearizer().matrix();
533 auto& ebosResid = ebosSimulator.model().linearizer().residual();
536 const double volume = 0.002831684659200;
538 const bool allow_cf = crossFlowAllowed(ebosSimulator);
540 const EvalWell& bhp = getBhp();
542 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
544 const int cell_idx = well_cells_[perf];
545 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
546 std::vector<EvalWell> cq_s(numComp,0.0);
547 std::vector<EvalWell> mob(numComp, 0.0);
548 getMobility(ebosSimulator, perf, mob);
549 computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
551 for (
int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
553 const EvalWell cq_s_effective = cq_s[componentIdx] * well_efficiency_factor_;
558 ebosResid[cell_idx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.value();
562 resWell_[0][componentIdx] -= cq_s_effective.value();
565 for (
int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
568 duneC_[0][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq);
570 invDuneD_[0][0][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx+numEq);
573 for (
int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
576 ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s_effective.derivative(pvIdx);
577 duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
582 if (has_solvent && componentIdx == contiSolventEqIdx) {
583 well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value();
585 well_state.perfPhaseRates()[(first_perf_ + perf) * np + componentIdx] = cq_s[componentIdx].value();
590 EvalWell cq_s_poly = cq_s[Water];
591 if (well_type_ == INJECTOR) {
592 cq_s_poly *= wpolymer();
594 cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
598 for (
int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
599 ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx);
601 ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value();
606 well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf];
610 for (
int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
611 EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;
612 resWell_loc += getQs(componentIdx) * well_efficiency_factor_;
613 for (
int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
614 invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq);
616 resWell_[0][componentIdx] += resWell_loc.value();
620 invDuneD_[0][0].invert();
627 template<
typename TypeTag>
629 StandardWell<TypeTag>::
630 crossFlowAllowed(
const Simulator& ebosSimulator)
const 632 if (getAllowCrossFlow()) {
640 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
641 const int cell_idx = well_cells_[perf];
642 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
643 const auto& fs = intQuants.fluidState();
644 EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
645 EvalWell bhp = getBhp();
648 EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
649 EvalWell drawdown = pressure - well_pressure;
651 if (drawdown.value() < 0 && well_type_ == INJECTOR) {
655 if (drawdown.value() > 0 && well_type_ == PRODUCER) {
666 template<
typename TypeTag>
668 StandardWell<TypeTag>::
669 getMobility(
const Simulator& ebosSimulator,
671 std::vector<EvalWell>& mob)
const 673 const int np = number_of_phases_;
674 const int cell_idx = well_cells_[perf];
675 assert (
int(mob.size()) == numComponents());
676 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
677 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
681 const int satid = saturation_table_number_[perf] - 1;
682 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
683 if( satid == satid_elem ) {
685 for (
int phase = 0; phase < np; ++phase) {
686 int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
687 mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx));
690 mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
694 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
695 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
696 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
699 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
702 for (
int phase = 0; phase < np; ++phase) {
703 int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
704 mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx));
709 OPM_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent");
716 EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
718 if (well_type_ == INJECTOR) {
719 const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex());
720 mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration,
true) );
723 if (PolymerModule::hasPlyshlog()) {
725 const int numComp = numComponents();
726 const bool allow_cf = crossFlowAllowed(ebosSimulator);
727 const EvalWell& bhp = getBhp();
728 std::vector<EvalWell> cq_s(numComp,0.0);
729 computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
731 double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf];
732 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
733 const auto& scaledDrainageInfo =
734 materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx);
735 const Scalar& Swcr = scaledDrainageInfo.Swcr;
736 const EvalWell poro = extendEval(intQuants.porosity());
737 const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water)));
739 const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12);
740 EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water)));
742 if (PolymerModule::hasShrate()) {
746 waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / bore_diameters_[perf];
748 EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
749 EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration,
750 intQuants.pvtRegionIndex(),
754 mob[ Water ] /= shearFactor;
763 template<
typename TypeTag>
765 StandardWell<TypeTag>::
766 updateWellState(
const BVectorWell& dwells,
767 WellState& well_state)
const 769 const int np = number_of_phases_;
770 const double dBHPLimit = param_.dbhp_max_rel_;
771 const double dFLimit = param_.dwell_fraction_max_;
772 const auto pu = phaseUsage();
774 const std::vector<double> xvar_well_old = primary_variables_;
777 std::vector<double> F(np,0.0);
778 if (active()[ Water ]) {
779 const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
780 const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac]),dFLimit);
781 primary_variables_[WFrac] = xvar_well_old[WFrac] - dx2_limited;
784 if (active()[ Gas ]) {
785 const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
786 const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac]),dFLimit);
787 primary_variables_[GFrac] = xvar_well_old[GFrac] - dx3_limited;
791 const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
792 const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]),dFLimit);
793 primary_variables_[SFrac] = xvar_well_old[SFrac] - dx4_limited;
796 assert(active()[ Oil ]);
797 F[pu.phase_pos[Oil]] = 1.0;
799 if (active()[ Water ]) {
800 F[pu.phase_pos[Water]] = primary_variables_[WFrac];
801 F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
804 if (active()[ Gas ]) {
805 F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
806 F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
809 double F_solvent = 0.0;
811 F_solvent = primary_variables_[SFrac];
812 F[pu.phase_pos[Oil]] -= F_solvent;
815 if (active()[ Water ]) {
816 if (F[Water] < 0.0) {
817 if (active()[ Gas ]) {
818 F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]);
821 F_solvent /= (1.0 - F[pu.phase_pos[Water]]);
823 F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]);
824 F[pu.phase_pos[Water]] = 0.0;
828 if (active()[ Gas ]) {
829 if (F[pu.phase_pos[Gas]] < 0.0) {
830 if (active()[ Water ]) {
831 F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]);
834 F_solvent /= (1.0 - F[pu.phase_pos[Gas]]);
836 F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]);
837 F[pu.phase_pos[Gas]] = 0.0;
841 if (F[pu.phase_pos[Oil]] < 0.0) {
842 if (active()[ Water ]) {
843 F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]);
845 if (active()[ Gas ]) {
846 F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]);
849 F_solvent /= (1.0 - F[pu.phase_pos[Oil]]);
851 F[pu.phase_pos[Oil]] = 0.0;
854 if (active()[ Water ]) {
855 primary_variables_[WFrac] = F[pu.phase_pos[Water]];
857 if (active()[ Gas ]) {
858 primary_variables_[GFrac] = F[pu.phase_pos[Gas]];
861 primary_variables_[SFrac] = F_solvent;
867 F[pu.phase_pos[Gas]] += F_solvent;
871 const WellControls* wc = well_controls_;
875 const int current = well_state.currentControls()[index_of_well_];
876 const double target_rate = well_controls_iget_target(wc, current);
878 for (
int p = 0; p < np; ++p) {
879 const double scal = scalingFactor(p);
887 switch (well_controls_iget_type(wc, current)) {
891 primary_variables_[XvarWell] = xvar_well_old[XvarWell] - dwells[0][XvarWell];
893 switch (well_type_) {
895 for (
int p = 0; p < np; ++p) {
896 const double comp_frac = comp_frac_[p];
897 well_state.wellRates()[index_of_well_ * np + p] = comp_frac * primary_variables_[XvarWell];
901 for (
int p = 0; p < np; ++p) {
902 well_state.wellRates()[index_of_well_ * np + p] = primary_variables_[XvarWell] * F[p];
907 if (well_controls_iget_type(wc, current) == THP) {
910 std::vector<double> rates(3, 0.0);
912 const Opm::PhaseUsage& pu = phaseUsage();
913 if (active()[ Water ]) {
914 rates[ Water ] = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ];
916 if (active()[ Oil ]) {
917 rates[ Oil ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
919 if (active()[ Gas ]) {
920 rates[ Gas ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ];
923 well_state.bhp()[index_of_well_] = calculateBhpFromThp(rates, current);
930 const int sign1 = dwells[0][XvarWell] > 0 ? 1: -1;
931 const double dx1_limited = sign1 * std::min(std::abs(dwells[0][XvarWell]),std::abs(xvar_well_old[XvarWell])*dBHPLimit);
932 primary_variables_[XvarWell] = std::max(xvar_well_old[XvarWell] - dx1_limited,1e5);
933 well_state.bhp()[index_of_well_] = primary_variables_[XvarWell];
935 if (well_controls_iget_type(wc, current) == SURFACE_RATE) {
936 if (well_type_ == PRODUCER) {
938 const double* distr = well_controls_iget_distr(wc, current);
940 double F_target = 0.0;
941 for (
int p = 0; p < np; ++p) {
942 F_target += distr[p] * F[p];
944 for (
int p = 0; p < np; ++p) {
945 well_state.wellRates()[np * index_of_well_ + p] = F[p] * target_rate / F_target;
949 for (
int p = 0; p < np; ++p) {
950 well_state.wellRates()[index_of_well_ * np + p] = comp_frac_[p] * target_rate;
954 for (
int p = 0; p < np; ++p) {
955 well_state.wellRates()[np * index_of_well_ + p] = F[p] * target_rate;
965 const int nwc = well_controls_get_num(wc);
968 for ( ; ctrl_index < nwc; ++ctrl_index) {
969 if (well_controls_iget_type(wc, ctrl_index) == THP) {
971 const int current = well_state.currentControls()[index_of_well_];
973 if (current == ctrl_index) {
974 const double thp_target = well_controls_iget_target(wc, current);
975 well_state.thp()[index_of_well_] = thp_target;
977 const Opm::PhaseUsage& pu = phaseUsage();
978 std::vector<double> rates(3, 0.0);
980 if (active()[ Water ]) {
981 rates[ Water ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Water ] ];
983 if (active()[ Oil ]) {
984 rates[ Oil ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Oil ] ];
986 if (active()[ Gas ]) {
987 rates[ Gas ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Gas ] ];
990 const double bhp = well_state.bhp()[index_of_well_];
992 well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, ctrl_index, bhp);
1001 if (ctrl_index == nwc) {
1002 well_state.thp()[index_of_well_] = 0.0;
1010 template<
typename TypeTag>
1017 const int np = number_of_phases_;
1018 const int well_index = index_of_well_;
1019 const WellControls* wc = well_controls_;
1022 const double target = well_controls_iget_target(wc, current);
1023 const double* distr = well_controls_iget_distr(wc, current);
1024 switch (well_controls_iget_type(wc, current)) {
1026 xw.bhp()[well_index] = target;
1032 xw.thp()[well_index] = target;
1034 const Opm::PhaseUsage& pu = phaseUsage();
1035 std::vector<double> rates(3, 0.0);
1036 if (active()[ Water ]) {
1037 rates[ Water ] = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
1039 if (active()[ Oil ]) {
1040 rates[ Oil ] = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
1042 if (active()[ Gas ]) {
1043 rates[ Gas ] = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
1046 xw.bhp()[well_index] = calculateBhpFromThp(rates, current);
1050 case RESERVOIR_RATE:
1053 int numPhasesWithTargetsUnderThisControl = 0;
1054 for (
int phase = 0; phase < np; ++phase) {
1055 if (distr[phase] > 0.0) {
1056 numPhasesWithTargetsUnderThisControl += 1;
1060 assert(numPhasesWithTargetsUnderThisControl > 0);
1062 if (well_type_ == INJECTOR) {
1065 assert(numPhasesWithTargetsUnderThisControl == 1);
1067 for (
int phase = 0; phase < np; ++phase) {
1068 if (distr[phase] > 0.) {
1069 xw.wellRates()[np*well_index + phase] = target / distr[phase];
1071 xw.wellRates()[np * well_index + phase] = 0.;
1074 }
else if (well_type_ == PRODUCER) {
1078 double original_rates_under_phase_control = 0.0;
1079 for (
int phase = 0; phase < np; ++phase) {
1080 if (distr[phase] > 0.0) {
1081 original_rates_under_phase_control += xw.wellRates()[np * well_index + phase] * distr[phase];
1085 if (original_rates_under_phase_control != 0.0 ) {
1086 double scaling_factor = target / original_rates_under_phase_control;
1088 for (
int phase = 0; phase < np; ++phase) {
1089 xw.wellRates()[np * well_index + phase] *= scaling_factor;
1093 const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl;
1094 for (
int phase = 0; phase < np; ++phase) {
1095 if (distr[phase] > 0.0) {
1096 xw.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase];
1099 xw.wellRates()[np * well_index + phase] = target_rate_divided;
1104 OPM_THROW(std::logic_error,
"Expected PRODUCER or INJECTOR type of well");
1110 updatePrimaryVariables(xw);
1117 template<
typename TypeTag>
1121 const WellState& xw,
1122 std::vector<double>& b_perf,
1123 std::vector<double>& rsmax_perf,
1124 std::vector<double>& rvmax_perf,
1125 std::vector<double>& surf_dens_perf)
const 1127 const int nperf = number_of_perforations_;
1129 const int numComp = numComponents();
1130 const PhaseUsage& pu = phaseUsage();
1131 b_perf.resize(nperf*numComp);
1132 surf_dens_perf.resize(nperf*numComp);
1133 const int w = index_of_well_;
1136 if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_used[BlackoilPhases::Liquid]) {
1137 rsmax_perf.resize(nperf);
1138 rvmax_perf.resize(nperf);
1142 for (
int perf = 0; perf < nperf; ++perf) {
1143 const int cell_idx = well_cells_[perf];
1144 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1145 const auto& fs = intQuants.fluidState();
1149 const double p_above = perf == 0 ? xw.bhp()[w] : xw.perfPress()[first_perf_ + perf - 1];
1150 const double p_avg = (xw.perfPress()[first_perf_ + perf] + p_above)/2;
1151 const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1153 if (pu.phase_used[BlackoilPhases::Aqua]) {
1154 b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * numComp] =
1155 FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1158 if (pu.phase_used[BlackoilPhases::Vapour]) {
1159 const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * numComp;
1160 const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases;
1162 if (pu.phase_used[BlackoilPhases::Liquid]) {
1163 const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases;
1164 const double oilrate = std::abs(xw.wellRates()[oilpos_well]);
1165 rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1167 const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
1170 rv = oilrate / gasrate;
1172 rv = std::min(rv, rvmax_perf[perf]);
1174 b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv);
1177 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1181 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1185 if (pu.phase_used[BlackoilPhases::Liquid]) {
1186 const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp;
1187 const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases;
1188 if (pu.phase_used[BlackoilPhases::Vapour]) {
1189 rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
1190 const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases;
1191 const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
1193 const double oilrate = std::abs(xw.wellRates()[oilpos_well]);
1196 rs = gasrate / oilrate;
1198 rs = std::min(rs, rsmax_perf[perf]);
1199 b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs);
1201 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1204 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1209 for (
int p = 0; p < pu.num_phases; ++p) {
1210 surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex());
1215 b_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
1216 surf_dens_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventRefDensity();
1225 template<
typename TypeTag>
1227 StandardWell<TypeTag>::
1228 computeConnectionDensities(
const std::vector<double>& perfComponentRates,
1229 const std::vector<double>& b_perf,
1230 const std::vector<double>& rsmax_perf,
1231 const std::vector<double>& rvmax_perf,
1232 const std::vector<double>& surf_dens_perf)
1235 const int np = number_of_phases_;
1236 const int nperf = number_of_perforations_;
1237 const int num_comp = numComponents();
1238 const PhaseUsage& phase_usage = phaseUsage();
1244 std::vector<double> q_out_perf(nperf*num_comp);
1248 for (
int perf = nperf - 1; perf >= 0; --perf) {
1249 for (
int component = 0; component < num_comp; ++component) {
1250 if (perf == nperf - 1) {
1252 q_out_perf[perf*num_comp+ component] = 0.0;
1255 q_out_perf[perf*num_comp + component] = q_out_perf[(perf+1)*num_comp + component];
1258 q_out_perf[perf*num_comp + component] -= perfComponentRates[perf*num_comp + component];
1266 const int gaspos = phase_usage.phase_pos[BlackoilPhases::Vapour];
1267 const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid];
1268 std::vector<double> mix(num_comp,0.0);
1269 std::vector<double> x(num_comp);
1270 std::vector<double> surf_dens(num_comp);
1271 std::vector<double> dens(nperf);
1273 for (
int perf = 0; perf < nperf; ++perf) {
1275 const double tot_surf_rate = std::accumulate(q_out_perf.begin() + num_comp*perf,
1276 q_out_perf.begin() + num_comp*(perf+1), 0.0);
1277 if (tot_surf_rate != 0.0) {
1278 for (
int component = 0; component < num_comp; ++component) {
1279 mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate);
1283 for (
int phase = 0; phase < np; ++phase) {
1284 mix[phase] = comp_frac_[phase];
1292 if (!rsmax_perf.empty() && mix[oilpos] > 0.0) {
1293 rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]);
1295 if (!rvmax_perf.empty() && mix[gaspos] > 0.0) {
1296 rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]);
1300 x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv);
1304 x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);;
1306 double volrat = 0.0;
1307 for (
int component = 0; component < num_comp; ++component) {
1308 volrat += x[component] / b_perf[perf*num_comp+ component];
1310 for (
int component = 0; component < num_comp; ++component) {
1311 surf_dens[component] = surf_dens_perf[perf*num_comp+ component];
1315 perf_densities_[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat;
1322 template<
typename TypeTag>
1324 StandardWell<TypeTag>::
1325 computeConnectionPressureDelta()
1341 const int nperf = number_of_perforations_;
1342 perf_pressure_diffs_.resize(nperf, 0.0);
1344 for (
int perf = 0; perf < nperf; ++perf) {
1345 const double z_above = perf == 0 ? ref_depth_ : perf_depth_[perf - 1];
1346 const double dz = perf_depth_[perf] - z_above;
1347 perf_pressure_diffs_[perf] = dz * perf_densities_[perf] * gravity_;
1354 const auto beg = perf_pressure_diffs_.begin();
1355 const auto end = perf_pressure_diffs_.end();
1356 std::partial_sum(beg, end, beg);
1363 template<
typename TypeTag>
1364 typename StandardWell<TypeTag>::ConvergenceReport
1368 typedef double Scalar;
1369 typedef std::vector< Scalar > Vector;
1371 const int np = number_of_phases_;
1372 const int numComp = numComponents();
1376 assert((
int(B_avg.size()) == numComp) || has_polymer);
1378 const double tol_wells = param_.tolerance_wells_;
1379 const double maxResidualAllowed = param_.max_residual_allowed_;
1383 std::vector<Scalar> res(numComp);
1384 for (
int comp = 0; comp < numComp; ++comp) {
1386 res[comp] = std::abs(resWell_[0][comp]);
1389 Vector well_flux_residual(numComp);
1392 for (
int compIdx = 0; compIdx < numComp; ++compIdx )
1394 well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx];
1400 for (
int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
1401 const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
1403 if (std::isnan(well_flux_residual[phaseIdx])) {
1404 report.nan_residual_found =
true;
1405 const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName};
1406 report.nan_residual_wells.push_back(problem_well);
1408 if (well_flux_residual[phaseIdx] > maxResidualAllowed) {
1409 report.too_large_residual_found =
true;
1410 const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName};
1411 report.too_large_residual_wells.push_back(problem_well);
1416 if ( !(report.nan_residual_found || report.too_large_residual_found) ) {
1418 for (
int compIdx = 0; compIdx < numComp; ++compIdx )
1420 report.converged = report.converged && (well_flux_residual[compIdx] < tol_wells);
1423 report.converged =
false;
1433 template<
typename TypeTag>
1437 const std::vector<double>& b_perf,
1438 const std::vector<double>& rsmax_perf,
1439 const std::vector<double>& rvmax_perf,
1440 const std::vector<double>& surf_dens_perf)
1443 const int nperf = number_of_perforations_;
1444 const int numComponent = numComponents();
1445 const int np = number_of_phases_;
1446 std::vector<double> perfRates(b_perf.size(),0.0);
1448 for (
int perf = 0; perf < nperf; ++perf) {
1449 for (
int phase = 0; phase < np; ++phase) {
1450 perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[(first_perf_ + perf) * np + phase];
1453 perfRates[perf*numComponent + contiSolventEqIdx] = xw.perfRateSolvent()[first_perf_ + perf];
1457 computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1459 computeConnectionPressureDelta();
1467 template<
typename TypeTag>
1469 StandardWell<TypeTag>::
1470 computeWellConnectionPressures(
const Simulator& ebosSimulator,
1471 const WellState& well_state)
1476 std::vector<double> b_perf;
1477 std::vector<double> rsmax_perf;
1478 std::vector<double> rvmax_perf;
1479 std::vector<double> surf_dens_perf;
1480 computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1481 computeWellConnectionDensitesPressures(well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1488 template<
typename TypeTag>
1490 StandardWell<TypeTag>::
1491 solveEqAndUpdateWellState(WellState& well_state)
1495 BVectorWell dx_well(1);
1496 invDuneD_.mv(resWell_, dx_well);
1498 updateWellState(dx_well, well_state);
1505 template<
typename TypeTag>
1507 StandardWell<TypeTag>::
1508 calculateExplicitQuantities(
const Simulator& ebosSimulator,
1509 const WellState& well_state)
1511 computeWellConnectionPressures(ebosSimulator, well_state);
1519 template<
typename TypeTag>
1521 StandardWell<TypeTag>::
1526 for (
int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
1527 F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value();
1535 template<
typename TypeTag>
1538 apply(
const BVector& x, BVector& Ax)
const 1540 assert( Bx_.size() == duneB_.N() );
1541 assert( invDrw_.size() == invDuneD_.N() );
1548 BVectorWell& invDBx = invDrw_;
1549 invDuneD_.mv(Bx_, invDBx);
1552 duneC_.mmtv(invDBx,Ax);
1558 template<
typename TypeTag>
1563 assert( invDrw_.size() == invDuneD_.N() );
1566 invDuneD_.mv(resWell_, invDrw_);
1568 duneC_.mmtv(invDrw_, r);
1575 template<
typename TypeTag>
1580 BVectorWell resWell = resWell_;
1582 duneB_.mmv(x, resWell);
1584 invDuneD_.mv(resWell, xw);
1591 template<
typename TypeTag>
1598 recoverSolutionWell(x, xw);
1599 updateWellState(xw, well_state);
1606 template<
typename TypeTag>
1610 const EvalWell& bhp,
1611 std::vector<double>& well_flux)
const 1613 const int np = number_of_phases_;
1614 const int numComp = numComponents();
1615 well_flux.resize(np, 0.0);
1617 const bool allow_cf = crossFlowAllowed(ebosSimulator);
1619 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
1620 const int cell_idx = well_cells_[perf];
1621 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1623 std::vector<EvalWell> cq_s(numComp, 0.0);
1624 std::vector<EvalWell> mob(numComp, 0.0);
1625 getMobility(ebosSimulator, perf, mob);
1626 computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
1628 for(
int p = 0; p < np; ++p) {
1629 well_flux[p] += cq_s[p].value();
1638 template<
typename TypeTag>
1640 StandardWell<TypeTag>::
1641 computeWellPotentialWithTHP(
const Simulator& ebosSimulator,
1642 const double initial_bhp,
1643 const std::vector<double>& initial_potential)
const 1647 const int np = number_of_phases_;
1649 assert( np ==
int(initial_potential.size()) );
1651 std::vector<double> potentials = initial_potential;
1652 std::vector<double> old_potentials = potentials;
1654 double bhp = initial_bhp;
1655 double old_bhp = bhp;
1657 bool converged =
false;
1658 const int max_iteration = 1000;
1659 const double bhp_tolerance = 1000.;
1663 while ( !converged && iteration < max_iteration ) {
1671 const int nwc = well_controls_get_num(well_controls_);
1673 for (
int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
1674 if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
1675 const Opm::PhaseUsage& pu = phaseUsage();
1677 std::vector<double> rates(3, 0.0);
1678 if (active()[ Water ]) {
1679 rates[ Water ] = potentials[pu.phase_pos[ Water ] ];
1681 if (active()[ Oil ]) {
1682 rates[ Oil ] = potentials[pu.phase_pos[ Oil ] ];
1684 if (active()[ Gas ]) {
1685 rates[ Gas ] = potentials[pu.phase_pos[ Gas ] ];
1688 const double bhp_calculated = calculateBhpFromThp(rates, ctrl_index);
1690 if (well_type_ == INJECTOR && bhp_calculated < bhp ) {
1691 bhp = bhp_calculated;
1694 if (well_type_ == PRODUCER && bhp_calculated > bhp) {
1695 bhp = bhp_calculated;
1701 if (std::isinf(bhp) || std::isnan(bhp)) {
1702 OPM_THROW(std::runtime_error,
"Unvalid bhp value obtained during the potential calculation for well " << name());
1705 converged = std::abs(old_bhp - bhp) < bhp_tolerance;
1707 computeWellRatesWithBhp(ebosSimulator, bhp, potentials);
1710 for (
const double value : potentials) {
1711 if (std::isinf(value) || std::isnan(value)) {
1712 OPM_THROW(std::runtime_error,
"Unvalid potential value obtained during the potential calculation for well " << name());
1718 for (
int p = 0; p < np; ++p) {
1721 const double potential_update_damping_factor = 0.001;
1722 potentials[p] = potential_update_damping_factor * potentials[p] + (1.0 - potential_update_damping_factor) * old_potentials[p];
1723 old_potentials[p] = potentials[p];
1731 OPM_THROW(std::runtime_error,
"Failed in getting converged for the potential calculation for well " << name());
1741 template<
typename TypeTag>
1746 std::vector<double>& well_potentials)
1748 updatePrimaryVariables(well_state);
1749 computeWellConnectionPressures(ebosSimulator, well_state);
1753 initPrimaryVariablesEvaluation();
1755 const int np = number_of_phases_;
1756 well_potentials.resize(np, 0.0);
1759 const double bhp = mostStrictBhpFromBhpLimits();
1762 if ( !wellHasTHPConstraints() ) {
1763 assert(std::abs(bhp) != std::numeric_limits<double>::max());
1765 computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials);
1769 if ( !well_state.isNewWell(index_of_well_) ) {
1770 for (
int p = 0; p < np; ++p) {
1773 well_potentials[p] = well_state.wellRates()[index_of_well_ * np + p];
1777 computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials);
1778 for (
double& value : well_potentials) {
1781 const double rate_safety_scaling_factor = 0.00001;
1782 value *= rate_safety_scaling_factor;
1786 well_potentials = computeWellPotentialWithTHP(ebosSimulator, bhp, well_potentials);
1794 template<
typename TypeTag>
1799 const int np = number_of_phases_;
1800 const int well_index = index_of_well_;
1801 const WellControls* wc = well_controls_;
1802 const double* distr = well_controls_get_current_distr(wc);
1803 const auto pu = phaseUsage();
1805 switch (well_controls_get_current_type(wc)) {
1808 primary_variables_[XvarWell] = 0.0;
1809 if (well_type_ == INJECTOR) {
1810 for (
int p = 0; p < np; ++p) {
1811 primary_variables_[XvarWell] += well_state.wellRates()[np*well_index + p] * comp_frac_[p];
1814 for (
int p = 0; p < np; ++p) {
1815 primary_variables_[XvarWell] += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
1820 case RESERVOIR_RATE:
1822 primary_variables_[XvarWell] = well_state.bhp()[well_index];
1826 double tot_well_rate = 0.0;
1827 for (
int p = 0; p < np; ++p) {
1828 tot_well_rate += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
1830 if(std::abs(tot_well_rate) > 0) {
1831 if (active()[ Water ]) {
1832 primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate;
1834 if (active()[ Gas ]) {
1835 primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ;
1838 primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / tot_well_rate ;
1841 if (well_type_ == INJECTOR) {
1843 if (active()[Water]) {
1844 if (distr[Water] > 0.0) {
1845 primary_variables_[WFrac] = 1.0;
1847 primary_variables_[WFrac] = 0.0;
1851 if (active()[Gas]) {
1852 if (distr[pu.phase_pos[Gas]] > 0.0) {
1853 primary_variables_[GFrac] = 1.0 - wsolvent();
1855 primary_variables_[SFrac] = wsolvent();
1858 primary_variables_[GFrac] = 0.0;
1865 }
else if (well_type_ == PRODUCER) {
1867 if (active()[Water]) {
1868 primary_variables_[WFrac] = 1.0 / np;
1870 if (active()[Gas]) {
1871 primary_variables_[GFrac] = 1.0 / np;
1874 OPM_THROW(std::logic_error,
"Expected PRODUCER or INJECTOR type of well");
1883 template<
typename TypeTag>
1884 template<
class ValueType>
1886 StandardWell<TypeTag>::
1887 calculateBhpFromThp(
const std::vector<ValueType>& rates,
1888 const int control_index)
const 1896 assert(
int(rates.size()) == 3);
1898 const ValueType aqua = rates[Water];
1899 const ValueType liquid = rates[Oil];
1900 const ValueType vapour = rates[Gas];
1902 const int vfp = well_controls_iget_vfp(well_controls_, control_index);
1903 const double& thp = well_controls_iget_target(well_controls_, control_index);
1904 const double& alq = well_controls_iget_alq(well_controls_, control_index);
1907 const double rho = perf_densities_[0];
1910 if (well_type_ == INJECTOR) {
1911 const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
1913 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
1915 bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
1917 else if (well_type_ == PRODUCER) {
1918 const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
1920 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
1922 bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
1925 OPM_THROW(std::logic_error,
"Expected INJECTOR or PRODUCER well");
1935 template<
typename TypeTag>
1937 StandardWell<TypeTag>::
1938 calculateThpFromBhp(
const std::vector<double>& rates,
1939 const int control_index,
1940 const double bhp)
const 1942 assert(
int(rates.size()) == 3);
1944 const double aqua = rates[Water];
1945 const double liquid = rates[Oil];
1946 const double vapour = rates[Gas];
1948 const int vfp = well_controls_iget_vfp(well_controls_, control_index);
1949 const double& alq = well_controls_iget_alq(well_controls_, control_index);
1952 const double rho = perf_densities_[0];
1955 if (well_type_ == INJECTOR) {
1956 const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
1958 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
1960 thp = vfp_properties_->getInj()->thp(vfp, aqua, liquid, vapour, bhp + dp);
1962 else if (well_type_ == PRODUCER) {
1963 const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
1965 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
1967 thp = vfp_properties_->getProd()->thp(vfp, aqua, liquid, vapour, bhp + dp, alq);
1970 OPM_THROW(std::logic_error,
"Expected INJECTOR or PRODUCER well");
1976 template<
typename TypeTag>
1978 StandardWell<TypeTag>::scalingFactor(
const int phaseIdx)
const 1980 const WellControls* wc = well_controls_;
1981 const double* distr = well_controls_get_current_distr(wc);
1983 if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
1984 if (has_solvent && phaseIdx == contiSolventEqIdx )
1985 OPM_THROW(std::runtime_error,
"RESERVOIR_RATE control in combination with solvent is not implemented");
1987 return distr[phaseIdx];
1989 const auto& pu = phaseUsage();
1990 if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
1992 if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
1994 if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
1996 if (has_solvent && phaseIdx == contiSolventEqIdx )
virtual void updateWellStateWithTarget(const int current, WellState &xw) const
updating the well state based the control mode specified with current
Definition: StandardWell_impl.hpp:1013
virtual void recoverWellSolutionAndUpdateWellState(const BVector &x, WellState &well_state) const
using the solution x to recover the solution xw for wells and applying xw to update Well State ...
Definition: StandardWell_impl.hpp:1594
Definition: StandardWell.hpp:33
a struct to collect information about the convergence checking
Definition: WellInterface.hpp:117
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
virtual ConvergenceReport getWellConvergence(const std::vector< double > &B_avg) const
check whether the well equations get converged for this well
Definition: StandardWell_impl.hpp:1366
virtual void computeWellPotentials(const Simulator &ebosSimulator, const WellState &well_state, std::vector< double > &well_potentials)
computing the well potentials for group control
Definition: StandardWell_impl.hpp:1744
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
virtual void apply(const BVector &x, BVector &Ax) const
Ax = Ax - C D^-1 B x.
Definition: StandardWell_impl.hpp:1538