diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index 56641126d..a41dea00a 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -261,7 +261,7 @@ assemblePressureEq(const int seg, bool gfrac) const { MultisegmentWellEquationAccess eqns(eqns1); - eqns.residual()[seg][SPres] = pressure_equation.value(); + eqns.residual()[seg][SPres] += pressure_equation.value(); eqns.D()[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); eqns.D()[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq); if (wfrac) { diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index bb4a9f5eb..d15d1126c 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -218,7 +218,7 @@ assembleDefaultPressureEq(const int seg, { assert(seg != 0); // not top segment const int seg_upwind = segments_.upwinding_segment(seg); - const bool reverseFlow = seg != seg_upwind; + const bool reverseFlow = seg != seg_upwind; // special treatment for reverse flow // for top segment, the well control equation will be used. EvalWell pressure_equation = primary_variables_.getSegmentPressure(seg); @@ -231,12 +231,24 @@ assembleDefaultPressureEq(const int seg, auto& ws = well_state.well(baseif_.indexOfWell()); auto& segments = ws.segments; - //pressure_equation -= hydro_pressure_drop; + + // We approximate the hydrostatic pressure loss over the segment by using the mean of the densities + // at the current segment node and its outlet node. Since density derivatives are organized + // differently than what is required for assemblePressureEq, this part needs to be assembled separately. + const int seg_outlet = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); + const auto hydro_pressure_drop_seg = segments_.getHalfHydroPressureLoss(seg, seg); + const auto hydro_pressure_drop_outlet = segments_.getHalfHydroPressureLoss(seg, seg_outlet); + MultisegmentWellAssemble(baseif_). + assembleHydroPressureLoss(seg, seg_outlet, hydro_pressure_drop_seg, hydro_pressure_drop_outlet, linSys_); + segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop_seg.value() + hydro_pressure_drop_outlet.value(); if (this->frictionalPressureLossConsidered()) { const auto friction_pressure_drop = segments_.getFrictionPressureLoss(seg, false); if (reverseFlow){ + // call function once again to obtain/assemble remaining derivatives extra_derivatives = -segments_.getFrictionPressureLoss(seg, true); + MultisegmentWellAssemble(baseif_). + assemblePressureEqExtraDerivatives(seg, seg_upwind, extra_derivatives, linSys_); } pressure_equation -= friction_pressure_drop; segments.pressure_drop_friction[seg] = friction_pressure_drop.value(); @@ -246,28 +258,13 @@ assembleDefaultPressureEq(const int seg, const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); const EvalWell outlet_pressure = primary_variables_.getSegmentPressure(outlet_segment_index); - //const int seg_upwind = segments_.upwinding_segment(seg); MultisegmentWellAssemble(baseif_). assemblePressureEq(seg, seg_upwind, outlet_segment_index, pressure_equation, outlet_pressure, linSys_); - if (reverseFlow){ - MultisegmentWellAssemble(baseif_). - assemblePressureEqExtraDerivatives(seg, seg_upwind, extra_derivatives, linSys_); - } - if (this->accelerationalPressureLossConsidered()) { handleAccelerationPressureLoss(seg, well_state); } - - // Do separately since different organization of derivatives - const int seg_outlet = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); - const auto hydro_pressure_drop1 = segments_.getHalfHydroPressureLoss(seg, seg); - const auto hydro_pressure_drop2 = segments_.getHalfHydroPressureLoss(seg, seg_outlet); - MultisegmentWellAssemble(baseif_). - assembleHydroPressureLoss(seg, seg_outlet, hydro_pressure_drop1, hydro_pressure_drop2, linSys_); - - segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop1.value() + hydro_pressure_drop2.value(); } template diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index ca068e898..248ec03f6 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -322,10 +322,9 @@ template typename MultisegmentWellSegments::EvalWell MultisegmentWellSegments:: getHalfHydroPressureLoss(const int seg, - const int seg_side) const + const int seg_density) const { - //const int seg_upwind = upwinding_segments_[seg]; - return 0.5*densities_[seg_side] * well_.gravity() * depth_diffs_[seg]; + return 0.5*densities_[seg_density] * well_.gravity() * depth_diffs_[seg]; } @@ -487,32 +486,37 @@ getSurfaceVolume(const EvalWell& temperature, template typename MultisegmentWellSegments::EvalWell MultisegmentWellSegments:: -getFrictionPressureLoss(const int seg, const bool return_upwind_derivatives) const +getFrictionPressureLoss(const int seg, const bool return_extra_derivatives) const { EvalWell mass_rate = mass_rates_[seg]; const int seg_upwind = upwinding_segments_[seg]; EvalWell density = densities_[seg_upwind]; EvalWell visc = viscosities_[seg_upwind]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. + // In the reverse flow case, we don't have enough slots for all derivatives, e.g., + // upwind pressure and flow. We amend this by a second function call optioin, where + // only these remaining derivatives are considered. + // For reference: the pressure equation assumes pressure/flow derivatives are given + // at segment node while fraction derivatives are given at upwind node. if (seg != seg_upwind) { - if (!return_upwind_derivatives){ - const int WQTotal = Indices::numEq + PrimaryVariables::WQTotal; - const int SPres = Indices::numEq + PrimaryVariables::SPres; + if (!return_extra_derivatives){ + constexpr int WQTotal = Indices::numEq + PrimaryVariables::WQTotal; + constexpr int SPres = Indices::numEq + PrimaryVariables::SPres; density.setDerivative(WQTotal, 0.0); density.setDerivative(SPres, 0.0); visc.setDerivative(WQTotal, 0.0); visc.setDerivative(SPres, 0.0); } else { - const int WFrac = Indices::numEq + PrimaryVariables::WFrac; - const int GFrac = Indices::numEq + PrimaryVariables::GFrac; - // if has ... - density.setDerivative(WFrac, 0.0); - density.setDerivative(GFrac, 0.0); - visc.setDerivative(WFrac, 0.0); - visc.setDerivative(GFrac, 0.0); + if (PrimaryVariables::has_water){ + constexpr int WFrac = Indices::numEq + PrimaryVariables::WFrac; + density.setDerivative(WFrac, 0.0); + visc.setDerivative(WFrac, 0.0); + } + if (PrimaryVariables::has_gas){ + constexpr int GFrac = Indices::numEq + PrimaryVariables::GFrac; + density.setDerivative(GFrac, 0.0); + visc.setDerivative(GFrac, 0.0); + } mass_rate.clearDerivatives(); } }