From 9ed0923279fbd33355c6ba95cee5edef148cfdb6 Mon Sep 17 00:00:00 2001
From: Xavier Raynaud <xavier.raynaud@sintef.no>
Date: Wed, 22 Aug 2012 10:11:32 +0200
Subject: [PATCH] Update surface volume in transport solver.

---
 .../TransportModelCompressibleTwophase.cpp        | 15 ++++++++++++---
 .../TransportModelCompressibleTwophase.hpp        |  4 ++--
 2 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
index 531aca4f5..06708da70 100644
--- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
+++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
@@ -75,15 +75,15 @@ namespace Opm
 
     void TransportModelCompressibleTwophase::solve(const double* darcyflux,
                                                    const double* pressure,
-                                                   const double* surfacevol0,
                                                    const double* porevolume0,
                                                    const double* porevolume,
                                                    const double* source,
                                                    const double dt,
-                                                   std::vector<double>& saturation)
+                                                   std::vector<double>& saturation,
+                                                   std::vector<double>& surfacevol)
     {
         darcyflux_ = darcyflux;
-        surfacevol0_ = surfacevol0;
+        surfacevol0_ = &surfacevol[0];
         porevolume0_ = porevolume0;
         porevolume_ = porevolume;
         source_ = source;
@@ -107,6 +107,15 @@ namespace Opm
                                &ia_downw_[0], &ja_downw_[0]);
         reorderAndTransport(grid_, darcyflux);
         toBothSat(saturation_, saturation);
+        
+        // Compute surface volume as a postprocessing step from saturation and A_
+        surfacevol = saturation;
+        const int np = props_.numPhases();
+        for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
+            for (int phase = 0; phase < np; ++phase) {
+                surfacevol[np*cell + phase] *= A_[np*np*cell + np*phase + phase];
+            }
+        } 
     }
 
     // Residual function r(s) for a single-cell implicit Euler transport
diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
index 19c44b22d..542363cd3 100644
--- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
+++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
@@ -54,12 +54,12 @@ namespace Opm
         /// \param[in, out] saturation   Phase saturations.
         void solve(const double* darcyflux,
                    const double* pressure,
-                   const double* surfacevol0,
                    const double* porevolume0,
                    const double* porevolume,
                    const double* source,
                    const double dt,
-                   std::vector<double>& saturation);
+                   std::vector<double>& saturation,
+                   std::vector<double>& surfacevol);
 
         /// Initialise quantities needed by gravity solver.
         void initGravity();