Add missing reverse flow derivatives for valve/icd/aicd

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -283,19 +283,31 @@ assembleICDPressureEq(const int seg,
     // the pressure equation is something like
     // p_seg - deltaP - p_outlet = 0.
     // the major part is how to calculate the deltaP
+    const int seg_upwind = segments_.upwinding_segment(seg);
+    const bool reverseFlow = seg != seg_upwind; // special treatment for reverse flow
 
     EvalWell pressure_equation = primary_variables_.getSegmentPressure(seg);
 
     EvalWell icd_pressure_drop;
+    EvalWell extra_derivatives;
     switch(this->segmentSet()[seg].segmentType()) {
         case Segment::SegmentType::SICD :
-            icd_pressure_drop = segments_.pressureDropSpiralICD(seg);
+            icd_pressure_drop = segments_.pressureDropSpiralICD(seg, /*extra derivatives*/false);
+            if (reverseFlow){
+                extra_derivatives = -segments_.pressureDropSpiralICD(seg, /*extra derivatives*/true);
+            }
             break;
         case Segment::SegmentType::AICD :
-            icd_pressure_drop = segments_.pressureDropAutoICD(seg, unit_system);
+            icd_pressure_drop = segments_.pressureDropAutoICD(seg, unit_system, /*extra derivatives*/false);
+            if (reverseFlow){
+                extra_derivatives = -segments_.pressureDropAutoICD(seg, unit_system, /*extra derivatives*/true);
+            }
             break;
         case Segment::SegmentType::VALVE :
-            icd_pressure_drop = segments_.pressureDropValve(seg);
+            icd_pressure_drop = segments_.pressureDropValve(seg, /*extra derivatives*/false);
+            if (reverseFlow){
+                extra_derivatives = -segments_.pressureDropValve(seg, /*extra derivatives*/true);
+            }
             break;
         default: {
             OPM_DEFLOG_THROW(std::runtime_error,
@@ -305,6 +317,12 @@ assembleICDPressureEq(const int seg,
                              deferred_logger);
         }
     }
+    if (reverseFlow){
+        MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
+            assemblePressureEqExtraDerivatives(seg, seg_upwind, extra_derivatives, linSys_);
+
+    }
+
     pressure_equation = pressure_equation - icd_pressure_drop;
     auto& ws = well_state.well(baseif_.indexOfWell());
     ws.segments.pressure_drop_friction[seg] = icd_pressure_drop.value();
@@ -313,7 +331,6 @@ assembleICDPressureEq(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,

opm/simulators/wells/MultisegmentWellSegments.cpp

@@ -549,7 +549,8 @@ getFrictionPressureLoss(const int seg, const bool return_extra_derivatives) cons
 template<class FluidSystem, class Indices, class Scalar>
 typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
-pressureDropSpiralICD(const int seg) const
+pressureDropSpiralICD(const int seg,
+                      const bool return_extra_derivatives) const
 {
     const auto& segment_set = well_.wellEcl().getSegments();
     const SICD& sicd = segment_set[seg].spiralICD();
@@ -583,17 +584,39 @@ pressureDropSpiralICD(const int seg) const
     }
 
     EvalWell density = densities_[seg_upwind];
-    // WARNING
+    EvalWell mass_rate = mass_rates_[seg];
-    // We disregard the derivatives from the upwind density to make sure derivatives
+    // In the reverse flow case, we don't have enough slots for all derivatives, e.g.,
-    // wrt. to different segments dont get mixed.
+    // upwind pressure and flow. We amend this by a second function call option, 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) {
-        water_fraction.clearDerivatives();
+        constexpr int nvar = FluidSystem::numPhases + 1;
-        water_viscosity.clearDerivatives();
+        std::vector<bool> zero_mask(nvar, false);
-        oil_fraction.clearDerivatives();
+        if (!return_extra_derivatives){
-        oil_viscosity.clearDerivatives();
+            zero_mask[PrimaryVariables::WQTotal] = true;
-        gas_fraction.clearDerivatives();
+            zero_mask[PrimaryVariables::SPres] = true;
-        gas_viscosity.clearDerivatives();
+        } else {
-        density.clearDerivatives();
+            if (PrimaryVariables::has_water){
+                zero_mask[PrimaryVariables::WFrac] = true;
+            }
+            if (PrimaryVariables::has_gas){
+                zero_mask[PrimaryVariables::GFrac] = true;
+            }
+            // mass_rate has no extra derivatives (they are organized as in equations)
+            mass_rate.clearDerivatives();
+        }
+        for (int ii = 0; ii < nvar; ++ii) {
+            if (zero_mask[ii]) {
+                water_fraction.setDerivative(Indices::numEq + ii, 0.0);
+                water_viscosity.setDerivative(Indices::numEq + ii, 0.0);
+                oil_fraction.setDerivative(Indices::numEq + ii, 0.0);
+                oil_viscosity.setDerivative(Indices::numEq + ii, 0.0);
+                gas_fraction.setDerivative(Indices::numEq + ii, 0.0);
+                gas_viscosity.setDerivative(Indices::numEq + ii, 0.0);
+                density.setDerivative(Indices::numEq + ii, 0.0);
+            }
+        }
     }
 
     const EvalWell liquid_fraction = water_fraction + oil_fraction;
@@ -604,7 +627,7 @@ pressureDropSpiralICD(const int seg) const
 
     const EvalWell mixture_viscosity = liquid_viscosity_fraction + gas_fraction * gas_viscosity;
 
-    const EvalWell reservoir_rate = mass_rates_[seg] / density;
+    const EvalWell reservoir_rate = mass_rate / density;
 
     const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor();
 
@@ -630,7 +653,8 @@ template<class FluidSystem, class Indices, class Scalar>
 typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
 pressureDropAutoICD(const int seg,
-                    const UnitSystem& unit_system) const
+                    const UnitSystem& unit_system, 
+                    const bool return_extra_derivatives) const
 {
     const auto& segment_set = well_.wellEcl().getSegments();
     const AutoICD& aicd = segment_set[seg].autoICD();
@@ -671,20 +695,42 @@ pressureDropAutoICD(const int seg,
     }
 
     EvalWell density = densities_[seg_upwind];
-    // WARNING
+    EvalWell mass_rate = mass_rates_[seg];
-    // We disregard the derivatives from the upwind density to make sure derivatives
+    // In the reverse flow case, we don't have enough slots for all derivatives, e.g.,
-    // wrt. to different segments dont get mixed.
+    // upwind pressure and flow. We amend this by a second function call option, 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) {
-        water_fraction.clearDerivatives();
+        constexpr int nvar = FluidSystem::numPhases + 1;
-        water_viscosity.clearDerivatives();
+        std::vector<bool> zero_mask(nvar, false);
-        water_density.clearDerivatives();
+        if (!return_extra_derivatives){
-        oil_fraction.clearDerivatives();
+            zero_mask[PrimaryVariables::WQTotal] = true;
-        oil_viscosity.clearDerivatives();
+            zero_mask[PrimaryVariables::SPres] = true;
-        oil_density.clearDerivatives();
+        } else {
-        gas_fraction.clearDerivatives();
+            if (PrimaryVariables::has_water){
-        gas_viscosity.clearDerivatives();
+                zero_mask[PrimaryVariables::WFrac] = true;
-        gas_density.clearDerivatives();
+            }
-        density.clearDerivatives();
+            if (PrimaryVariables::has_gas){
+                zero_mask[PrimaryVariables::GFrac] = true;
+            }
+            // mass_rate has no extra derivatives (they are organized as in equations)
+            mass_rate.clearDerivatives();
+        }
+        for (int ii = 0; ii < nvar; ++ii) {
+            if (zero_mask[ii]) {
+                water_fraction.setDerivative(Indices::numEq + ii, 0.0);
+                water_viscosity.setDerivative(Indices::numEq + ii, 0.0);
+                water_density.setDerivative(Indices::numEq + ii, 0.0);
+                oil_fraction.setDerivative(Indices::numEq + ii, 0.0);
+                oil_viscosity.setDerivative(Indices::numEq + ii, 0.0);
+                oil_density.setDerivative(Indices::numEq + ii, 0.0);
+                gas_fraction.setDerivative(Indices::numEq + ii, 0.0);
+                gas_viscosity.setDerivative(Indices::numEq + ii, 0.0);
+                gas_density.setDerivative(Indices::numEq + ii, 0.0);
+                density.setDerivative(Indices::numEq + ii, 0.0);
+            }
+        }
     }
 
     using MathTool = MathToolbox<EvalWell>;
@@ -698,7 +744,7 @@ pressureDropAutoICD(const int seg,
 
     const double rho_reference = aicd.densityCalibration();
     const double visc_reference = aicd.viscosityCalibration();
-    const auto volume_rate_icd = mass_rates_[seg] * aicd.scalingFactor() / mixture_density;
+    const auto volume_rate_icd = mass_rate * aicd.scalingFactor() / mixture_density;
     const double sign = volume_rate_icd <= 0. ? 1.0 : -1.0;
     // convert 1 unit volume rate
     using M  = UnitSystem::measure;
@@ -715,21 +761,42 @@ pressureDropAutoICD(const int seg,
 template<class FluidSystem, class Indices, class Scalar>
 typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
-pressureDropValve(const int seg) const
+pressureDropValve(const int seg,
+                  const bool return_extra_derivatives) const
 {
     const auto& segment_set = well_.wellEcl().getSegments();
     const Valve& valve = segment_set[seg].valve();
 
-    const EvalWell& mass_rate = mass_rates_[seg];
+    EvalWell mass_rate = mass_rates_[seg];
     const int seg_upwind = upwinding_segments_[seg];
     EvalWell visc = viscosities_[seg_upwind];
     EvalWell density = densities_[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) {
-        visc.clearDerivatives();
+        if (!return_extra_derivatives){
-        density.clearDerivatives();
+            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 {
+            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();
+        }
     }
 
     const double additional_length = valve.pipeAdditionalLength();

opm/simulators/wells/MultisegmentWellSegments.hpp

@@ -74,15 +74,18 @@ public:
     EvalWell getFrictionPressureLoss(const int seg, const bool return_upwind_derivatives) const;
 
     // pressure drop for Spiral ICD segment (WSEGSICD)
-    EvalWell pressureDropSpiralICD(const int seg) const;
+    EvalWell pressureDropSpiralICD(const int seg,
+                                   const bool return_extra_derivatives) const;
 
     // pressure drop for Autonomous ICD segment (WSEGAICD)
     EvalWell pressureDropAutoICD(const int seg,
-                                 const UnitSystem& unit_system) const;
+                                 const UnitSystem& unit_system,
+                                 const bool return_extra_derivatives) const;
 
     // pressure drop for sub-critical valve (WSEGVALV)
-    EvalWell pressureDropValve(const int seg) const;
+    EvalWell pressureDropValve(const int seg, 
-
+                               const bool return_extra_derivatives) const;
+                               
     // pressure loss due to acceleration
     EvalWell accelerationPressureLoss(const int seg) const;