From 4b3a6f37d578a28197037a904c04d9c5b8c1aab1 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Mon, 19 Dec 2022 15:12:43 +0100
Subject: [PATCH] move calculation part of handleAccelerationPressureLoss to
 MultisegmentwellSegments

---
 opm/simulators/wells/MultisegmentWellEval.cpp | 36 ++--------------
 .../wells/MultisegmentWellSegments.cpp        | 41 +++++++++++++++++++
 .../wells/MultisegmentWellSegments.hpp        |  3 ++
 3 files changed, 48 insertions(+), 32 deletions(-)

diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp
index 26e7cfe67..7a9316951 100644
--- a/opm/simulators/wells/MultisegmentWellEval.cpp
+++ b/opm/simulators/wells/MultisegmentWellEval.cpp
@@ -185,43 +185,15 @@ MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 handleAccelerationPressureLoss(const int seg,
                                WellState& well_state)
 {
-    const double area = this->segmentSet()[seg].crossArea();
-    const EvalWell mass_rate = segments_.mass_rates_[seg];
-    const int seg_upwind = segments_.upwinding_segments_[seg];
-    EvalWell density = segments_.densities_[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();
-    }
-
-    EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density);
-    // handling the velocity head of intlet segments
-    for (const int inlet : this->segments_.inlets_[seg]) {
-        const int seg_upwind_inlet = segments_.upwinding_segments_[inlet];
-        const double inlet_area = this->segmentSet()[inlet].crossArea();
-        EvalWell inlet_density = this->segments_.densities_[seg_upwind_inlet];
-        // WARNING
-        // We disregard the derivatives from the upwind density to make sure derivatives
-        // wrt. to different segments dont get mixed.
-        if (inlet != seg_upwind_inlet) {
-            inlet_density.clearDerivatives();
-        }
-        const EvalWell inlet_mass_rate = segments_.mass_rates_[inlet];
-        accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density);
-    }
-
-    // We change the sign of the accelerationPressureLoss for injectors.
-    // Is this correct? Testing indicates that this is what the reference simulator does
-    const double sign = mass_rate < 0. ? 1.0 : - 1.0;
-    accelerationPressureLoss *= sign;
+    const EvalWell accelerationPressureLoss = segments_.accelerationPressureLoss(seg);
 
     auto& segments = well_state.well(baseif_.indexOfWell()).segments;
     segments.pressure_drop_accel[seg] = accelerationPressureLoss.value();
 
     MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
-        assemblePressureLoss(seg, seg_upwind, accelerationPressureLoss, linSys_);
+        assemblePressureLoss(seg,
+                             segments_.upwinding_segments_[seg],
+                             accelerationPressureLoss, linSys_);
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp
index 5ef5d9ea8..651b4f7cb 100644
--- a/opm/simulators/wells/MultisegmentWellSegments.cpp
+++ b/opm/simulators/wells/MultisegmentWellSegments.cpp
@@ -673,6 +673,47 @@ pressureDropValve(const int seg) const
     return sign * (friction_pressure_loss + constriction_pressure_loss);
 }
 
+template<class FluidSystem, class Indices, class Scalar>
+typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
+MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
+accelerationPressureLoss(const int seg) const
+{
+    const auto& segment_set = well_.wellEcl().getSegments();
+    const double area = segment_set[seg].crossArea();
+    const EvalWell mass_rate = mass_rates_[seg];
+    const int seg_upwind = upwinding_segments_[seg];
+    EvalWell density = densities_[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();
+    }
+
+    EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density);
+    // handling the velocity head of intlet segments
+    for (const int inlet : inlets_[seg]) {
+        const int seg_upwind_inlet = upwinding_segments_[inlet];
+        const double inlet_area = segment_set[inlet].crossArea();
+        EvalWell inlet_density = densities_[seg_upwind_inlet];
+        // WARNING
+        // We disregard the derivatives from the upwind density to make sure derivatives
+        // wrt. to different segments dont get mixed.
+        if (inlet != seg_upwind_inlet) {
+            inlet_density.clearDerivatives();
+        }
+        const EvalWell inlet_mass_rate = mass_rates_[inlet];
+        accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density);
+    }
+
+    // We change the sign of the accelerationPressureLoss for injectors.
+    // Is this correct? Testing indicates that this is what the reference simulator does
+    const double sign = mass_rate < 0. ? 1.0 : - 1.0;
+    accelerationPressureLoss *= sign;
+
+    return accelerationPressureLoss;
+}
+
 #define INSTANCE(...) \
 template class MultisegmentWellSegments<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;
 
diff --git a/opm/simulators/wells/MultisegmentWellSegments.hpp b/opm/simulators/wells/MultisegmentWellSegments.hpp
index a86bec795..0f728e69c 100644
--- a/opm/simulators/wells/MultisegmentWellSegments.hpp
+++ b/opm/simulators/wells/MultisegmentWellSegments.hpp
@@ -68,6 +68,9 @@ public:
     // pressure drop for sub-critical valve (WSEGVALV)
     EvalWell pressureDropValve(const int seg) const;
 
+    // pressure loss due to acceleration
+    EvalWell accelerationPressureLoss(const int seg) const;
+
     // TODO: trying to use the information from the Well opm-parser as much
     // as possible, it will possibly be re-implemented later for efficiency reason.