From 876c8864c7132e458548bd9340e7f0a9d3d2d774 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tor.harald.sandve@iris.no>
Date: Tue, 15 Sep 2015 10:44:07 +0200
Subject: [PATCH] Use cell thickness from EclipseGrid in well transmissibility
 calculations

The dz calculated in WellDetails::getCubeDim is not correct in cases
where the face centroid of the horizontal faces is located above or
below the face centroid or the vertical faces. The cell thickness in
EclipseGrid, calculated using the Z-coordinates, is therefore used
instead.
---
 opm/core/wells/WellsManager.hpp      |  1 +
 opm/core/wells/WellsManager_impl.hpp | 18 +++++++++++++++++-
 2 files changed, 18 insertions(+), 1 deletion(-)

diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp
index 58baaf7e..e29a1016 100644
--- a/opm/core/wells/WellsManager.hpp
+++ b/opm/core/wells/WellsManager.hpp
@@ -164,6 +164,7 @@ namespace Opm
                                    const int* cart_dims,
                                    FC begin_face_centroids, 
                                    int dimensions,
+                                   std::vector<double>& dz,
                                    std::vector<std::string>& well_names,
                                    std::vector<WellData>& well_data,
                                    std::map<std::string, int> & well_names_to_index,
diff --git a/opm/core/wells/WellsManager_impl.hpp b/opm/core/wells/WellsManager_impl.hpp
index b038bc0d..b48cfa93 100644
--- a/opm/core/wells/WellsManager_impl.hpp
+++ b/opm/core/wells/WellsManager_impl.hpp
@@ -107,6 +107,7 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
                                         const int* cart_dims,
                                         FC begin_face_centroids,
                                         int dimensions,
+                                        std::vector<double>& dz,
                                         std::vector<std::string>& well_names,
                                         std::vector<WellData>& well_data,
                                         std::map<std::string, int>& well_names_to_index,
@@ -178,9 +179,14 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
                                     OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
                                 }
 
-                                const std::array<double, 3> cubical =
+                                std::array<double, 3> cubical =
                                     WellsManagerDetail::getCubeDim<3>(c2f, begin_face_centroids, cell);
 
+                                // overwrite dz values calculated in getCubeDim.
+                                if (dz.size() > 0) {
+                                    cubical[2] = dz[cell];
+                                }
+
                                 const double* cell_perm = &permeability[dimensions*dimensions*cell];
                                 pd.well_index =
                                     WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm,
@@ -369,10 +375,20 @@ WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
     DoubleArray ntg_glob(eclipseState, "NTG", 1.0);
     NTGArray    ntg(ntg_glob, global_cell);
 
+    EclipseGridConstPtr eclGrid = eclipseState->getEclipseGrid();
+
+    // use cell thickness (dz) from eclGrid
+    // dz overwrites values calculated by WellDetails::getCubeDim
+    std::vector<double> dz(number_of_cells);
+    for (int cell = 0; cell < number_of_cells; ++cell) {
+        dz[cell] = eclGrid->getCellThicknes(global_cell[cell]);
+    }
+
     createWellsFromSpecs(wells, timeStep, cell_to_faces,
                          cart_dims,
                          begin_face_centroids,
                          dimensions,
+                         dz,
                          well_names, well_data, well_names_to_index,
                          pu, cartesian_to_compressed, permeability, ntg,
                          wells_on_proc);