From ce74f956e7fa268ff957905c05b9e742a6c627ad Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 20 Apr 2012 14:00:38 +0200
Subject: [PATCH] Accounting for (constant) formation volume factor in
 incompressible fluids.

---
 opm/core/fluid/IncompPropertiesFromDeck.cpp    |  2 +-
 opm/core/fluid/PvtPropertiesIncompFromDeck.cpp | 18 +++++++++++++++---
 opm/core/fluid/PvtPropertiesIncompFromDeck.hpp | 18 ++++++++++++++----
 3 files changed, 30 insertions(+), 8 deletions(-)

diff --git a/opm/core/fluid/IncompPropertiesFromDeck.cpp b/opm/core/fluid/IncompPropertiesFromDeck.cpp
index 9febc0857..2789fe920 100644
--- a/opm/core/fluid/IncompPropertiesFromDeck.cpp
+++ b/opm/core/fluid/IncompPropertiesFromDeck.cpp
@@ -87,7 +87,7 @@ namespace Opm
     /// \return Array of P density values.
     const double* IncompPropertiesFromDeck::density() const
     {
-	return pvt_.surfaceDensities();
+	return pvt_.reservoirDensities();
     }
 
     /// \param[in]  n      Number of data points.
diff --git a/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp b/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
index dc5e80436..62cc1b231 100644
--- a/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
+++ b/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
@@ -51,18 +51,23 @@ namespace Opm
 	if (deck.hasField("DENSITY")) {
 	    const std::vector<double>& d = deck.getDENSITY().densities_[region_number];
 	    enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
-	    density_[phase_usage.phase_pos[PhaseUsage::Aqua]]   = d[ECL_water];
-	    density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
+	    surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]]   = d[ECL_water];
+	    surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
 	} else {
 	    THROW("Input is missing DENSITY\n");
 	}
 
+        // Make reservoir densities the same as surface densities initially.
+        // We will modify them with formation volume factors if found.
+        reservoir_density_ = surface_density_;
+
         // Water viscosity.
 	if (deck.hasField("PVTW")) {
 	    const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number];
 	    if (pvtw[2] != 0.0 || pvtw[4] != 0.0) {
 		MESSAGE("Compressibility effects in PVTW are ignored.");
 	    }
+            reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1];
 	    viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3];
 	} else {
 	    // Eclipse 100 default.
@@ -76,6 +81,7 @@ namespace Opm
 	    if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) {
 		MESSAGE("Compressibility effects in PVCDO are ignored.");
 	    }
+            reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1];
 	    viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3];
 	} else {
 	    THROW("Input is missing PVCDO\n");
@@ -84,7 +90,13 @@ namespace Opm
 
     const double* PvtPropertiesIncompFromDeck::surfaceDensities() const
     {
-        return density_.data();
+        return surface_density_.data();
+    }
+
+
+    const double* PvtPropertiesIncompFromDeck::reservoirDensities() const
+    {
+        return reservoir_density_.data();
     }
 
 
diff --git a/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp b/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp
index 84f143740..acbabaa49 100644
--- a/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp
+++ b/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp
@@ -31,9 +31,6 @@ namespace Opm
     /// eclipse input (keywords DENSITY, PVTW, PVCDO).
     ///
     /// All phases are incompressible and have constant viscosities.
-    /// For all the methods, the following apply: p and z are unused.
-    /// Output arrays shall be of size n*numPhases(), and must be valid
-    /// before calling the method.
     /// NOTE: This class is intentionally similar to BlackoilPvtProperties.
     class PvtPropertiesIncompFromDeck
     {
@@ -51,11 +48,24 @@ namespace Opm
         /// \return  Array of size numPhases().
 	const double* surfaceDensities() const;
 
+        /// Densities of stock components at reservoir conditions.
+        /// Note: a reasonable question to ask is why there can be
+        /// different densities at surface and reservoir conditions,
+        /// when the phases are assumed incompressible. The answer is
+        /// that even if we approximate the phases as being
+        /// incompressible during simulation, the density difference
+        /// between surface and reservoir may be larger. For accurate
+        /// reporting and using data given in terms of surface values,
+        /// we need to handle this difference.
+        /// \return  Array of size numPhases().
+	const double* reservoirDensities() const;
+
         /// Viscosities.
         const double* viscosity() const;
 
     private:
-	std::tr1::array<double, 2> density_;
+	std::tr1::array<double, 2> surface_density_;
+	std::tr1::array<double, 2> reservoir_density_;
 	std::tr1::array<double, 2> viscosity_;
     };