ECL peaceman well: rename "bottom depth" to "reference depth"

this is because the "bottom hole pressure" by default is the the
presssure at the center of the _topmost_ cell penetrated by the
well. (The terminology keeps getting more weird the longer you look at
it.)

this patch also changes the default way to calculate the reference
depth from the center of the deepest to the center of the highest cell
and makes the reference depth settable in the deck.

applications/ebos/eclpeacemanwell.hh

@@ -465,9 +465,8 @@ public:
      */
     void beginSpec()
     {
-        // this is going to be increased by any realistic grid. Shall we bet?
+        // this is going to be set to a real value by any realistic grid. Shall we bet?
-        bottomDepth_ = -1e100;
+        refDepth_ = 1e100;
-        bottomDofGlobalIdx_ = -1;
 
         // By default, take the bottom hole pressure as a given
         controlMode_ = ControlMode::BottomHolePressure;
@@ -610,10 +609,8 @@ public:
 
         // we assume that the z-coordinate represents depth (and not
         // height) here...
-        if (dofPos[2] > bottomDepth_) {
+        if (dofPos[2] < refDepth_)
-            bottomDofGlobalIdx_ = globalDofIdx;
+            refDepth_ = dofPos[2];
-            bottomDepth_ = dofPos[2];
-        }
     }
 
     /*!
@@ -624,7 +621,7 @@ public:
         const auto& comm = simulator_.gridView().comm();
 
         // determine the maximum depth of the well over all processes
-        bottomDepth_ = comm.max(bottomDepth_);
+        refDepth_ = comm.min(refDepth_);
 
         // the total volume of the well must also be summed over all processes
         wellTotalVolume_ = comm.sum(wellTotalVolume_);
@@ -690,10 +687,16 @@ public:
     { injectedPhaseIdx_ = injPhaseIdx; }
 
     /*!
-     * \brief The Z-coordinate of the well's deepest degree of freedom
+     * \brief Sets the reference depth for the bottom hole pressure [m]
      */
-    Scalar bottomDepth() const
+    void setReferenceDepth(Scalar value)
-    { return bottomDepth_; }
+    { refDepth_ = value; }
+
+    /*!
+     * \brief The reference depth for the bottom hole pressure [m]
+     */
+    Scalar referenceDepth() const
+    { return refDepth_; }
 
     /*!
      * \brief Set whether the well is open,closed or shut
@@ -756,7 +759,7 @@ public:
         // warning: this is a bit hacky...
         Scalar rho = 650; // kg/m^3
         Scalar g = 9.81; // m/s^2
-        return actualBottomHolePressure_ + rho*bottomDepth_*g;
+        return actualBottomHolePressure_ + rho*refDepth_*g;
     }
 
     /*!
@@ -865,7 +868,7 @@ public:
             // assume a density of 650 kg/m^3 for the bottom hole pressure
             // calculation
             Scalar rho = 650.0;
-            targetBottomHolePressure_ = thpLimit_ + rho*bottomDepth_;
+            targetBottomHolePressure_ = thpLimit_ + rho*refDepth_;
         }
         else if (controlMode_ == ControlMode::BottomHolePressure)
             targetBottomHolePressure_ = bhpLimit_;
@@ -1181,10 +1184,10 @@ protected:
             Valgrind::CheckDefined(rho);
             Valgrind::CheckDefined(lambda);
             Valgrind::CheckDefined(depth);
-            Valgrind::CheckDefined(bottomDepth_);
+            Valgrind::CheckDefined(refDepth_);
 
             // pressure in the borehole ("hole pressure") at the given location
-            Scalar ph = pbh + rho*g*(bottomDepth_ - depth);
+            Scalar ph = pbh + rho*g*(refDepth_ - depth);
 
             // volumetric flux of the phase from the well to the reservoir
             volRates[phaseIdx] = Twj*lambda*(ph - p);
@@ -1350,9 +1353,8 @@ protected:
     Scalar computeRateEquivalentBhp_() const
     {
         if (wellStatus() == Shut)
-            // there is no flow happening in the well, so the "BHP" is the pressure of
+            // there is no flow happening in the well, so we return 0...
-            // the well's lowest DOF!
+            return 0.0;
-            return dofVariables_.at(bottomDofGlobalIdx_).pressure[oilPhaseIdx];
 
         // initialize the bottom hole pressure which we would like to calculate
         Scalar bhp = actualBottomHolePressure_;
@@ -1543,12 +1545,9 @@ protected:
 
     int injectedPhaseIdx_;
 
-    // the depth of the deepest DOF. (actually, the center of this
+    // the reference depth for the bottom hole pressure. if not specified otherwise, this
-    // DOF, but the difference should be minimal.)
+    // is the position of the _highest_ DOF in the well.
-    Scalar bottomDepth_;
+    Scalar refDepth_;
-
-    // global index of the DOF at the bottom of the well
-    int bottomDofGlobalIdx_;
 };
 } // namespace Ewoms
 

applications/ebos/eclwellmanager.hh

@@ -119,7 +119,7 @@ public:
 
         // set those parameters of the wells which do not change the topology of the
         // linearized system of equations
-        updateWellCompletions_(wellCompMap);
+        updateWellParameters_(episodeIdx, wellCompMap);
 
         const std::vector<Opm::WellConstPtr>& deckWells = deckSchedule->getWells(episodeIdx);
         // set the injection data for the respective wells.
@@ -626,8 +626,21 @@ protected:
         }
     }
 
-    void updateWellCompletions_(const WellCompletionsMap& wellCompletions)
+    void updateWellParameters_(int reportStepIdx, const WellCompletionsMap& wellCompletions)
     {
+        auto eclStatePtr = simulator_.gridManager().eclState();
+        auto deckSchedule = eclStatePtr->getSchedule();
+        const std::vector<Opm::WellConstPtr>& deckWells = deckSchedule->getWells(reportStepIdx);
+
+        // set the reference depth for all wells
+        for (size_t deckWellIdx = 0; deckWellIdx < deckWells.size(); ++deckWellIdx) {
+            Opm::WellConstPtr deckWell = deckWells[deckWellIdx];
+            const std::string& wellName = deckWell->name();
+
+            if (!deckWell->getRefDepthDefaulted())
+                wells_[wellIndex(wellName)]->setReferenceDepth(deckWell->getRefDepth());
+        }
+
         // associate the well completions with grid cells and register them in the
         // Peaceman well object
         const GridView gridView = simulator_.gridManager().gridView();