From 17400cebafa579cc2335c68b2c9f8e43bccf6f95 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Mon, 28 Apr 2014 13:19:18 +0200
Subject: [PATCH] Fix bug in use of capillary pressure.

Pcow must be subtracted, not added to oil pressure to obtain water pressure.
---
 .../FullyImplicitBlackoilSolver_impl.hpp       | 18 ++++++++++++++++--
 1 file changed, 16 insertions(+), 2 deletions(-)

diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index 4f43a8f5a..ded3d6672 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -56,6 +56,12 @@
                   << foo.value() << std::endl;                          \
     } while (0)
 
+#define DISKVAL(foo)                                                    \
+    do {                                                                \
+        std::ofstream os(#foo);                                         \
+        os.precision(16);                                               \
+        os << foo.value() << std::endl;                                 \
+    } while (0)
 
 
 namespace Opm {
@@ -1573,9 +1579,17 @@ namespace {
             pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
         }
 
-        // add the total pressure to the capillary pressures
+        // Since pcow = po - pw, but pcog = pg - po,
+        // we have
+        //   pw = po - pcow
+        //   pg = po + pcgo
+        // This is an unfortunate inconsistency, but a convention we must handle.
         for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
-            pressure[phaseIdx] += state.pressure;
+            if (phaseIdx == BlackoilPhases::Aqua) {
+                pressure[phaseIdx] = state.pressure - pressure[phaseIdx];
+            } else {
+                pressure[phaseIdx] += state.pressure;
+            }
         }
 
         return pressure;