moving computeWellConnectionPressures to StandardWells

the results look okay, while the running for flow_solvent needs further
investigation even the results with flow_solvent actually look okay.

With two different version of
computePropertiesForWellConnectionPressures, flow_solvent produces the
same results. This is something needs further investigation.

The current implementation requires a copy of
computeWellConnectionPressure in StandardWells and StandardWellsSolvent.
That means probably we need to introduce the asImpl() for the Wells.

opm/autodiff/BlackoilModelBase.hpp

@@ -380,9 +380,6 @@ namespace Opm {
         computeAccum(const SolutionState& state,
                      const int            aix  );
 
-        void computeWellConnectionPressures(const SolutionState& state,
-                                            const WellState& xw);
-
         void
         assembleMassBalanceEq(const SolutionState& state);
 
@@ -410,6 +407,8 @@ namespace Opm {
                               const std::vector<ADB>& mob_perfcells,
                               const std::vector<ADB>& b_perfcells,
                               WellState& well_state);
+
+        void
         addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
                                            const SolutionState& state,
                                            const WellState& xw);

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -798,41 +798,6 @@ namespace detail {
 
 
 
-    template <class Grid, class WellModel, class Implementation>
-    void 
-    BlackoilModelBase<Grid, WellModel, Implementation>::
-    computeWellConnectionPressures(const SolutionState& state,
-                                   const WellState& xw)
-    {
-        if( ! localWellsActive() ) return ;
-
-        using namespace Opm::AutoDiffGrid;
-        // 1. Compute properties required by computeConnectionPressureDelta().
-        //    Note that some of the complexity of this part is due to the function
-        //    taking std::vector<double> arguments, and not Eigen objects.
-        std::vector<double> b_perf;
-        std::vector<double> rsmax_perf;
-        std::vector<double> rvmax_perf;
-        std::vector<double> surf_dens_perf;
-        asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
-
-        // Extract well connection depths.
-        const typename WellModel::Vector depth = cellCentroidsZToEigen(grid_);
-        const typename WellModel::Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells);
-        const int nperf = wells().well_connpos[wells().number_of_wells];
-        const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
-
-        // Gravity
-        const double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
-
-        asImpl().stdWells().computeWellConnectionDensitesPressures(xw, fluid_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, grav);
-
-    }
-
-
-
-
-
     template <class Grid, class WellModel, class Implementation>
     void
     BlackoilModelBase<Grid, WellModel, Implementation>::
@@ -842,6 +807,9 @@ namespace detail {
     {
         using namespace Opm::AutoDiffGrid;
 
+        const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
+        const V depth = cellCentroidsZToEigen(grid_);
+
         // If we have VFP tables, we need the well connection
         // pressures for the "simple" hydrostatic correction
         // between well depth and vfp table depth.
@@ -849,14 +817,15 @@ namespace detail {
             SolutionState state = asImpl().variableState(reservoir_state, well_state);
             SolutionState state0 = state;
             asImpl().makeConstantState(state0);
-            asImpl().computeWellConnectionPressures(state0, well_state);
+            // asImpl().computeWellConnectionPressures(state0, well_state);
+            // Extract well connection depths.
+            asImpl().stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity);
         }
 
         // Possibly switch well controls and updating well state to
         // get reasonable initial conditions for the wells
         // asImpl().updateWellControls(well_state);
         // asImpl().stdWells().updateWellControls(well_state);
-        const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
         asImpl().stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state);
 
         // Create the primary variables.
@@ -869,7 +838,8 @@ namespace detail {
             // Compute initial accumulation contributions
             // and well connection pressures.
             asImpl().computeAccum(state0, 0);
-            asImpl().computeWellConnectionPressures(state0, well_state);
+            // asImpl().computeWellConnectionPressures(state0, well_state);
+            asImpl().stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity);
         }
 
         // OPM_AD_DISKVAL(state.pressure);
@@ -1133,6 +1103,9 @@ namespace detail {
             }
         }
 
+        // gravity
+        const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
+
         int it  = 0;
         bool converged;
         do {
@@ -1147,7 +1120,6 @@ namespace detail {
             asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
             asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
             asImpl().stdWells().addWellFluxEq(cq_s, wellSolutionState, residual_);
-            const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
             asImpl().stdWells().addWellControlEq(wellSolutionState, well_state, aliveWells,
                                                  active_, vfp_properties_, gravity, residual_);
             converged = getWellConvergence(it);
@@ -1201,7 +1173,9 @@ namespace detail {
                 std::vector<ADB::M> old_derivs = state.qs.derivative();
                 state.qs = ADB::function(std::move(new_qs), std::move(old_derivs));
             }
-            asImpl().computeWellConnectionPressures(state, well_state);
+            // asImpl().computeWellConnectionPressures(state, well_state);
+            const ADB::V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_);
+            asImpl().stdWells().computeWellConnectionPressures(state, well_state, fluid_, active_, phaseCondition(), depth, gravity);
         }
 
         if (!converged) {

opm/autodiff/BlackoilSolventModel.hpp

@@ -146,7 +146,7 @@ namespace Opm {
         using Base::drMaxRel;
         using Base::maxResidualAllowed;
         // using Base::updateWellControls;
-        using Base::computeWellConnectionPressures;
+        // using Base::computeWellConnectionPressures;
         // using Base::addWellControlEq;
         // using Base::computePropertiesForWellConnectionPressures;
 

opm/autodiff/StandardWells.hpp

@@ -157,6 +157,16 @@ namespace Opm {
                                   const double gravity,
                                   LinearisedBlackoilResidual& residual);
 
+            template <class SolutionState, class WellState>
+            void computeWellConnectionPressures(const SolutionState& state,
+                                                const WellState& xw,
+                                                const BlackoilPropsAdInterface& fluid,
+                                                const std::vector<bool>& active,
+                                                const std::vector<PhasePresence>& phaseCondition,
+                                                const Vector& depth,
+                                                const double gravity);
+
+
 
 
         protected:

opm/autodiff/StandardWellsSolvent.hpp

@@ -34,6 +34,7 @@ namespace Opm {
         public:
 
             using Base = StandardWells;
+            using Base::computeWellConnectionDensitesPressures;
 
             // ---------  Public methods  ---------
             explicit StandardWellsSolvent(const Wells* wells);
@@ -63,6 +64,16 @@ namespace Opm {
                                            const std::vector<bool>& active,
                                            std::vector<ADB>& mob_perfcells,
                                            std::vector<ADB>& b_perfcells) const;
+
+            template <class SolutionState, class WellState>
+            void computeWellConnectionPressures(const SolutionState& state,
+                                                const WellState& xw,
+                                                const BlackoilPropsAdInterface& fluid,
+                                                const std::vector<bool>& active,
+                                                const std::vector<PhasePresence>& phaseCondition,
+                                                const Vector& depth,
+                                                const double gravity);
+
         protected:
             const SolventPropsAdFromDeck* solvent_props_;
             int solvent_pos_;

opm/autodiff/StandardWellsSolvent_impl.hpp

@@ -197,6 +197,42 @@ namespace Opm
 
 
 
+    template <class SolutionState, class WellState>
+    void
+    StandardWellsSolvent::
+    computeWellConnectionPressures(const SolutionState& state,
+                                   const WellState& xw,
+                                   const BlackoilPropsAdInterface& fluid,
+                                   const std::vector<bool>& active,
+                                   const std::vector<PhasePresence>& phaseCondition,
+                                   const Vector& depth,
+                                   const double gravity)
+    {
+        if( ! localWellsActive() ) return ;
+        // 1. Compute properties required by computeConnectionPressureDelta().
+        //    Note that some of the complexity of this part is due to the function
+        //    taking std::vector<double> arguments, and not Eigen objects.
+        std::vector<double> b_perf;
+        std::vector<double> rsmax_perf;
+        std::vector<double> rvmax_perf;
+        std::vector<double> surf_dens_perf;
+        computePropertiesForWellConnectionPressures(state, xw, fluid, active, phaseCondition,
+                                                    b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
+
+        // Extract well connection depths
+        // TODO: depth_perf should be a member of the StandardWells class
+        const Vector pdepth = subset(depth, wellOps().well_cells);
+        const int nperf = wells().well_connpos[wells().number_of_wells];
+        const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
+
+        computeWellConnectionDensitesPressures(xw, fluid, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity);
+
+    }
+
+
+
+
+
     template <class ReservoirResidualQuant, class SolutionState>
     void
     StandardWellsSolvent::

opm/autodiff/StandardWells_impl.hpp

@@ -289,6 +289,42 @@ namespace Opm
 
 
 
+    template <class SolutionState, class WellState>
+    void
+    StandardWells::
+    computeWellConnectionPressures(const SolutionState& state,
+                                   const WellState& xw,
+                                   const BlackoilPropsAdInterface& fluid,
+                                   const std::vector<bool>& active,
+                                   const std::vector<PhasePresence>& phaseCondition,
+                                   const Vector& depth,
+                                   const double gravity)
+    {
+        if( ! localWellsActive() ) return ;
+        // 1. Compute properties required by computeConnectionPressureDelta().
+        //    Note that some of the complexity of this part is due to the function
+        //    taking std::vector<double> arguments, and not Eigen objects.
+        std::vector<double> b_perf;
+        std::vector<double> rsmax_perf;
+        std::vector<double> rvmax_perf;
+        std::vector<double> surf_dens_perf;
+        computePropertiesForWellConnectionPressures(state, xw, fluid, active, phaseCondition,
+                                                    b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
+
+        // Extract well connection depths
+        // TODO: depth_perf should be a member of the StandardWells class
+        const Vector pdepth = subset(depth, wellOps().well_cells);
+        const int nperf = wells().well_connpos[wells().number_of_wells];
+        const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
+
+        computeWellConnectionDensitesPressures(xw, fluid, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity);
+
+    }
+
+
+
+
+
     template <class ReservoirResidualQuant, class SolutionState>
     void
     StandardWells::

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -202,7 +202,7 @@ namespace Opm {
         using Base::maxResidualAllowed;
 
         // using Base::updateWellControls;
-        using Base::computeWellConnectionPressures;
+        // using Base::computeWellConnectionPressures;
         // using Base::addWellControlEq;
         using Base::computeRelPerm;
 

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -499,6 +499,7 @@ namespace Opm {
         // Possibly switch well controls and updating well state to
         // get reasonable initial conditions for the wells
         const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
+        const V depth = cellCentroidsZToEigen(grid_);
         // updateWellControls(well_state);
         stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state);
 
@@ -512,7 +513,8 @@ namespace Opm {
             // Compute initial accumulation contributions
             // and well connection pressures.
             computeAccum(state0, 0);
-            computeWellConnectionPressures(state0, well_state);
+            // computeWellConnectionPressures(state0, well_state);
+            stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity);
         }
 
         // OPM_AD_DISKVAL(state.pressure);