Inluding all derivatives for reversed flow

opm/simulators/wells/MultisegmentWellAssemble.cpp

@@ -216,6 +216,35 @@ assemblePressureLoss(const int seg,
     }
 }
 
+template<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assembleHydroPressureLoss(const int seg, 
+                          const int seg_upwind, 
+                          const EvalWell& hydro_pressure_drop, 
+                          Equations& eqns1) const
+{
+    MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
+    eqns.residual()[seg][SPres] -= hydro_pressure_drop.value();
+    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+        eqns.D()[seg][seg_upwind][SPres][pv_idx] -= hydro_pressure_drop.derivative(pv_idx + Indices::numEq);
+    }
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assemblePressureEqExtraDerivatives(const int seg, 
+                                   const int seg_upwind, 
+                                   const EvalWell& extra_derivatives, 
+                                   Equations& eqns1) const
+{
+    MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
+    // diregard residual
+    // Frac - derivatives are zero (they belong to upwind^2)
+    eqns.D()[seg][seg_upwind][SPres][SPres] += extra_derivatives.derivative(SPres + Indices::numEq);
+    eqns.D()[seg][seg_upwind][SPres][WQTotal] += extra_derivatives.derivative(WQTotal + Indices::numEq);
+}                                   
+
+
 template<class FluidSystem, class Indices, class Scalar>
 void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
 assemblePressureEq(const int seg,

opm/simulators/wells/MultisegmentWellAssemble.hpp

@@ -87,6 +87,17 @@ public:
                               const EvalWell& accelerationPressureLoss,
                               Equations& eqns) const;
 
+
+    void assembleHydroPressureLoss(const int seg, 
+                                   const int seg_upwind, 
+                                   const EvalWell& hydro_pressure_drop, 
+                                   Equations& eqs1) const;
+
+    void assemblePressureEqExtraDerivatives(const int seg, 
+                                            const int seg_upwind, 
+                                            const EvalWell& extra_derivatives, 
+                                            Equations& eqns1) const;
+
     //! \brief Assemble pressure terms.
     void assemblePressureEq(const int seg,
                             const int seg_upwind,

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -217,22 +217,27 @@ assembleDefaultPressureEq(const int seg,
                           WellState& well_state)
 {
     assert(seg != 0); // not top segment
+    const int seg_upwind = segments_.upwinding_segment(seg);
+    const bool reverseFlow = seg != seg_upwind;
 
     // for top segment, the well control equation will be used.
     EvalWell pressure_equation = primary_variables_.getSegmentPressure(seg);
+    EvalWell extra_derivatives;
 
     // we need to handle the pressure difference between the two segments
     // we only consider the hydrostatic pressure loss first
     // TODO: we might be able to add member variables to store these values, then we update well state
     // after converged
-    const auto hydro_pressure_drop = segments_.getHydroPressureLoss(seg);
+
     auto& ws = well_state.well(baseif_.indexOfWell());
     auto& segments = ws.segments;
-    segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop.value();
+    //pressure_equation -= hydro_pressure_drop;
-    pressure_equation -= hydro_pressure_drop;
 
     if (this->frictionalPressureLossConsidered()) {
-        const auto friction_pressure_drop = segments_.getFrictionPressureLoss(seg);
+        const auto friction_pressure_drop = segments_.getFrictionPressureLoss(seg, false);
+        if (reverseFlow){
+            extra_derivatives = -segments_.getFrictionPressureLoss(seg, true);
+        }
         pressure_equation -= friction_pressure_drop;
         segments.pressure_drop_friction[seg] = friction_pressure_drop.value();
     }
@@ -241,14 +246,25 @@ 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);
+    //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 auto hydro_pressure_drop = segments_.getHydroPressureLoss(seg);
+    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
+        assembleHydroPressureLoss(seg, seg_upwind, hydro_pressure_drop, linSys_);
+    segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop.value();
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>

opm/simulators/wells/MultisegmentWellSegments.cpp

@@ -322,10 +322,12 @@ template<class FluidSystem, class Indices, class Scalar>
 typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
 getHydroPressureLoss(const int seg) const
-{
+{   
-    return densities_[seg] * well_.gravity() * depth_diffs_[seg];
+    const int seg_upwind = upwinding_segments_[seg];
+    return densities_[seg_upwind] * well_.gravity() * depth_diffs_[seg];
 }
 
+
 template<class FluidSystem, class Indices, class Scalar>
 Scalar MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
 getPressureDiffSegPerf(const int seg,
@@ -484,19 +486,36 @@ 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
+getFrictionPressureLoss(const int seg, const bool return_upwind_derivatives) const
 {
-    const EvalWell mass_rate = mass_rates_[seg];
+    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.
+    
     if (seg != seg_upwind) {
-        density.clearDerivatives();
+        if (!return_upwind_derivatives){
-        visc.clearDerivatives();
+            const int WQTotal = Indices::numEq + PrimaryVariables::WQTotal;
+            const 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);
+            mass_rate.clearDerivatives();
+        }
     }
+    
     const auto& segment_set = well_.wellEcl().getSegments();
     const int outlet_segment_index = segment_set.segmentNumberToIndex(segment_set[seg].outletSegment());
     const double length = segment_set[seg].totalLength() - segment_set[outlet_segment_index].totalLength();

opm/simulators/wells/MultisegmentWellSegments.hpp

@@ -63,7 +63,7 @@ public:
                               const int pvt_region_index,
                               const int seg_idx) const;
 
-    EvalWell getFrictionPressureLoss(const int seg) const;
+    EvalWell getFrictionPressureLoss(const int seg, const bool return_upwind_derivatives) const;
 
     // pressure drop for Spiral ICD segment (WSEGSICD)
     EvalWell pressureDropSpiralICD(const int seg) const;