Updates as requested PR review

opm/simulators/wells/MultisegmentWellAssemble.cpp

@@ -261,7 +261,7 @@ assemblePressureEq(const int seg,
                    bool gfrac) const
 {
     MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> 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) {

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<FluidSystem,Indices,Scalar>(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<FluidSystem,Indices,Scalar>(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<FluidSystem,Indices,Scalar>(baseif_).
         assemblePressureEq(seg, seg_upwind, outlet_segment_index,
                            pressure_equation, outlet_pressure, linSys_);
 
-    if (reverseFlow){
-        MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(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<FluidSystem,Indices,Scalar>(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<typename FluidSystem, typename Indices, typename Scalar>

opm/simulators/wells/MultisegmentWellSegments.cpp

@@ -322,10 +322,9 @@ template<class FluidSystem, class Indices, class Scalar>
 typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
 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_density] * well_.gravity() * depth_diffs_[seg];
-    return 0.5*densities_[seg_side] * well_.gravity() * depth_diffs_[seg];
 }
 
 
@@ -487,32 +486,37 @@ getSurfaceVolume(const EvalWell& temperature,
 template<class FluidSystem, class Indices, class Scalar>
 typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
-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
+    // In the reverse flow case, we don't have enough slots for all derivatives, e.g.,
-    // We disregard the derivatives from the upwind density to make sure derivatives
+    // upwind pressure and flow. We amend this by a second function call optioin, where 
-    // wrt. to different segments dont get mixed.
+    // 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){
+        if (!return_extra_derivatives){
-            const int WQTotal = Indices::numEq + PrimaryVariables::WQTotal;
+            constexpr int WQTotal = Indices::numEq + PrimaryVariables::WQTotal;
-            const int SPres = Indices::numEq + PrimaryVariables::SPres;
+            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;
+            if (PrimaryVariables::has_water){
-            const int GFrac = Indices::numEq + PrimaryVariables::GFrac;
+                constexpr int WFrac = Indices::numEq + PrimaryVariables::WFrac;
-            // if has ...
+                density.setDerivative(WFrac, 0.0);
-            density.setDerivative(WFrac, 0.0);
+                visc.setDerivative(WFrac, 0.0);
-            density.setDerivative(GFrac, 0.0);
+            }
-            visc.setDerivative(WFrac, 0.0);
+            if (PrimaryVariables::has_gas){
-            visc.setDerivative(GFrac, 0.0);
+                constexpr int GFrac = Indices::numEq + PrimaryVariables::GFrac;
+                density.setDerivative(GFrac, 0.0);
+                visc.setDerivative(GFrac, 0.0);
+            }
             mass_rate.clearDerivatives();
         }
     }