Revert to original segment density computations

(keep average version as option)

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -164,6 +164,10 @@ template<class TypeTag, class MyTypeTag>
 struct MaximumNumberOfWellSwitches {
     using type = UndefinedProperty;
 };
+template<class TypeTag, class MyTypeTag>
+struct UseAverageDensityMsWells {
+    using type = UndefinedProperty;
+};
 
 template<class TypeTag>
 struct DbhpMaxRel<TypeTag, TTag::FlowModelParameters> {
@@ -308,6 +312,14 @@ template<class TypeTag>
 struct MaximumNumberOfWellSwitches<TypeTag, TTag::FlowModelParameters> {
     static constexpr int value = 3;
 };
+template<class TypeTag>
+struct UseAverageDensityMsWells<TypeTag, TTag::FlowModelParameters> {
+    static constexpr bool value = false;
+};
+
+
+
+
 
 // if openMP is available, determine the number threads per process automatically.
 #if _OPENMP
@@ -421,6 +433,9 @@ namespace Opm
         /// Maximum number of times a well can switch to the same controt
         int max_number_of_well_switches_;
 
+        /// Whether to approximate segment densities by averaging over segment and its outlet 
+        bool use_average_density_ms_wells_;
+
 
 
         /// Construct from user parameters or defaults.
@@ -457,6 +472,7 @@ namespace Opm
             check_well_operability_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
             check_well_operability_iter_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter);
             max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches);
+            use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells);
             deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
         }
 
@@ -495,6 +511,7 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheck, "Enable the well operability checking");
             EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter, "Enable the well operability checking during iterations");
             EWOMS_REGISTER_PARAM(TypeTag, int, MaximumNumberOfWellSwitches, "Maximum number of times a well can switch to the same control");
+            EWOMS_REGISTER_PARAM(TypeTag, bool, UseAverageDensityMsWells, "Approximate segment densitities by averaging over segment and its outlet");
         }
     };
 } // namespace Opm

opm/simulators/flow/FlowMainEbos.hpp

@@ -220,6 +220,9 @@ void handleExtraConvergenceOutput(SimulatorReport& report,
             EWOMS_HIDE_PARAM(TypeTag, VtkWriteTortuosities);
             EWOMS_HIDE_PARAM(TypeTag, VtkWriteDiffusionCoefficients);
             EWOMS_HIDE_PARAM(TypeTag, VtkWriteEffectiveDiffusionCoefficients);
+            
+            // hide average density option 
+            EWOMS_HIDE_PARAM(TypeTag, UseAverageDensityMsWells);
 
             EWOMS_END_PARAM_REGISTRATION(TypeTag);
 

opm/simulators/wells/MultisegmentWellAssemble.cpp

@@ -219,19 +219,16 @@ assemblePressureLoss(const int seg,
 template<class FluidSystem, class Indices, class Scalar>
 void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
 assembleHydroPressureLoss(const int seg, 
-                          const int seg_outlet, 
+                          const int seg_density, 
                           const EvalWell& hydro_pressure_drop_seg, 
-                          const EvalWell& hydro_pressure_drop_outlet, 
                           Equations& eqns1) const
 {
     MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
-    eqns.residual()[seg][SPres] -= (hydro_pressure_drop_seg.value() + hydro_pressure_drop_outlet.value());
+    eqns.residual()[seg][SPres] -= hydro_pressure_drop_seg.value();
     for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
-        eqns.D()[seg][seg][SPres][pv_idx] -= hydro_pressure_drop_seg.derivative(pv_idx + Indices::numEq);
+        eqns.D()[seg][seg_density][SPres][pv_idx] -= hydro_pressure_drop_seg.derivative(pv_idx + Indices::numEq);
-    }
-    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
-        eqns.D()[seg][seg_outlet][SPres][pv_idx] -= hydro_pressure_drop_outlet.derivative(pv_idx + Indices::numEq);
     }
+
 }
 
 template<class FluidSystem, class Indices, class Scalar>

opm/simulators/wells/MultisegmentWellAssemble.hpp

@@ -89,9 +89,8 @@ public:
 
 
     void assembleHydroPressureLoss(const int seg, 
-                                   const int seg_outlet, 
+                                   const int seg_density, 
                                    const EvalWell& hydro_pressure_drop_seg, 
-                                   const EvalWell& hydro_pressure_drop_outlet, 
                                    Equations& eqns1) const;                                   
 
     void assemblePressureEqExtraDerivatives(const int seg, 

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -214,7 +214,8 @@ template<typename FluidSystem, typename Indices, typename Scalar>
 void
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 assembleDefaultPressureEq(const int seg,
-                          WellState& well_state)
+                          WellState& well_state, 
+                          const bool use_average_density)
 {
     assert(seg != 0); // not top segment
     const int seg_upwind = segments_.upwinding_segment(seg);
@@ -256,15 +257,22 @@ assembleDefaultPressureEq(const int seg,
         handleAccelerationPressureLoss(seg, well_state);
     }
 
-    // We approximate the hydrostatic pressure loss over the segment by using the mean of the densities
+    // Since density derivatives are organized differently than what is required for assemblePressureEq, 
-    // at the current segment node and its outlet node. Since density derivatives are organized 
+    // this part needs to be assembled separately. Optionally use average density variant.  
-    // differently than what is required for assemblePressureEq, this part needs to be assembled separately.  
+    const auto hydro_pressure_drop_seg = segments_.getHydroPressureLoss(seg, seg);
-    const int seg_outlet = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
+    if (!use_average_density){
-    const auto hydro_pressure_drop_seg = segments_.getHalfHydroPressureLoss(seg, seg);
+        MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
-    const auto hydro_pressure_drop_outlet = segments_.getHalfHydroPressureLoss(seg, seg_outlet);
+            assembleHydroPressureLoss(seg, seg, hydro_pressure_drop_seg, linSys_);
-    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
+            segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop_seg.value();
-        assembleHydroPressureLoss(seg, seg_outlet, hydro_pressure_drop_seg, hydro_pressure_drop_outlet, linSys_);
+    } else {
-    segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop_seg.value() + hydro_pressure_drop_outlet.value();
+        const int seg_outlet = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
+        const auto hydro_pressure_drop_outlet = segments_.getHydroPressureLoss(seg, seg_outlet);
+        MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
+            assembleHydroPressureLoss(seg, seg, 0.5*hydro_pressure_drop_seg, linSys_);
+        MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
+            assembleHydroPressureLoss(seg, seg_outlet, 0.5*hydro_pressure_drop_outlet, linSys_);
+        segments.pressure_drop_hydrostatic[seg] = 0.5*hydro_pressure_drop_seg.value() + 0.5*hydro_pressure_drop_outlet.value();
+    }
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -338,6 +346,7 @@ MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 assemblePressureEq(const int seg,
                    const UnitSystem& unit_system,
                    WellState& well_state,
+                   const bool use_average_density,
                    DeferredLogger& deferred_logger)
 {
     switch(this->segmentSet()[seg].segmentType()) {
@@ -348,7 +357,7 @@ assemblePressureEq(const int seg,
             break;
         }
         default :
-            assembleDefaultPressureEq(seg, well_state);
+            assembleDefaultPressureEq(seg, well_state, use_average_density);
     }
 }
 

opm/simulators/wells/MultisegmentWellEval.hpp

@@ -76,7 +76,8 @@ protected:
     void initMatrixAndVectors(const int num_cells);
 
     void assembleDefaultPressureEq(const int seg,
-                                   WellState& well_state);
+                                   WellState& well_state,
+                                   const bool use_average_density);
 
     // assemble pressure equation for ICD segments
     void assembleICDPressureEq(const int seg,
@@ -88,6 +89,7 @@ protected:
     void assemblePressureEq(const int seg,
                             const UnitSystem& unit_system,
                             WellState& well_state,
+                            const bool use_average_density,
                             DeferredLogger& deferred_logger);
 
     /// check whether the well equations get converged for this well

opm/simulators/wells/MultisegmentWellSegments.cpp

@@ -321,10 +321,10 @@ updateUpwindingSegments(const PrimaryVariables& primary_variables)
 template<class FluidSystem, class Indices, class Scalar>
 typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
-getHalfHydroPressureLoss(const int seg,
+getHydroPressureLoss(const int seg,
                          const int seg_density) const
 {   
-    return 0.5*densities_[seg_density] * well_.gravity() * depth_diffs_[seg];
+    return densities_[seg_density] * well_.gravity() * depth_diffs_[seg];
 }
 
 

opm/simulators/wells/MultisegmentWellSegments.hpp

@@ -51,7 +51,7 @@ public:
     //! \brief Update upwinding segments.
     void updateUpwindingSegments(const PrimaryVariables& primary_variables);
 
-    EvalWell getHalfHydroPressureLoss(const int seg,
+    EvalWell getHydroPressureLoss(const int seg,
                                       const int seg_side) const;
 
     //! Pressure difference between segment and perforation.

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1589,7 +1589,7 @@ namespace Opm
                                         deferred_logger);
             } else {
                 const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
-                this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
+                this->assemblePressureEq(seg, unit_system, well_state, this->param_.use_average_density_ms_wells_, deferred_logger);
             }
         }