From 96327164ccfc20db032358d172c58b805b437986 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Thu, 23 Aug 2012 08:57:48 +0200
Subject: [PATCH] Added computeSurfacevol() method. The method is called by the
 reordering transport solver after computing new saturations in order to
 update the surface volumes.

---
 .../TransportModelCompressibleTwophase.cpp    | 11 ++-----
 opm/core/utility/miscUtilitiesBlackoil.cpp    | 30 +++++++++++++++++++
 opm/core/utility/miscUtilitiesBlackoil.hpp    | 19 ++++++++++++
 3 files changed, 52 insertions(+), 8 deletions(-)

diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
index 98aeefdf5..b81a7e60a 100644
--- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
+++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
@@ -24,6 +24,7 @@
 #include <opm/core/transport/reorder/reordersequence.h>
 #include <opm/core/utility/RootFinders.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 #include <opm/core/pressure/tpfa/trans_tpfa.h>
 
 #include <fstream>
@@ -112,15 +113,9 @@ 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];
-            }
-        } 
+        computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]);
     }
 
     // Residual function r(s) for a single-cell implicit Euler transport
diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp
index af088e059..1c8bceee2 100644
--- a/opm/core/utility/miscUtilitiesBlackoil.cpp
+++ b/opm/core/utility/miscUtilitiesBlackoil.cpp
@@ -220,5 +220,35 @@ namespace Opm
         }
     }
 
+    /// Computes the surface volume densities from saturations by the formula
+    ///     z = A s
+    /// for a number of data points, where z is the surface volume density,
+    /// s is the saturation (both as column vectors) and A is the
+    /// phase-to-component relation matrix.
+    /// @param[in]  n            number of data points
+    /// @param[in]  np           number of phases, must be 2 or 3
+    /// @param[in]  A            array containing n square matrices of size num_phases^2,
+    ///                          in Fortran ordering, typically the output of a call
+    ///                          to the matrix() method of a BlackoilProperties* class.
+    /// @param[in]  saturation   concatenated saturation values (for all P phases)
+    /// @param[out] surfacevol   concatenated surface-volume values (for all P phases)
+    void computeSurfacevol(const int n,
+                           const int np,
+                           const double* A,
+                           const double* saturation,
+                           double* surfacevol)
+    {
+        // Note: since this is a simple matrix-vector product, it can
+        // be done by a BLAS call, but then we have to reorder the A
+        // matrix data.
+        std::fill(surfacevol, surfacevol + n*np, 0.0);
+        for (int i = 0; i < n; ++i) {
+            for (int row = 0; row < np; ++row) {
+                for (int col = 0; col < np; ++col) {
+                    surfacevol[i*np + row] += A[i*np*np + row*np + col] * saturation[i*np + col];
+                }
+            }
+        }
+    }
 
 } // namespace Opm
diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp
index 929d42ecd..4a3d6448a 100644
--- a/opm/core/utility/miscUtilitiesBlackoil.hpp
+++ b/opm/core/utility/miscUtilitiesBlackoil.hpp
@@ -112,6 +112,25 @@ namespace Opm
                                const std::vector<double>& s,
                                std::vector<double>& fractional_flows);
 
+
+    /// Computes the surface volume densities from saturations by the formula
+    ///     z = A s
+    /// for a number of data points, where z is the surface volume density,
+    /// s is the saturation (both as column vectors) and A is the
+    /// phase-to-component relation matrix.
+    /// @param[in]  n            number of data points
+    /// @param[in]  np           number of phases, must be 2 or 3
+    /// @param[in]  A            array containing n square matrices of size num_phases,
+    ///                          in Fortran ordering, typically the output of a call
+    ///                          to the matrix() method of a BlackoilProperties* class.
+    /// @param[in]  saturation   concatenated saturation values (for all P phases)
+    /// @param[out] surfacevol   concatenated surface-volume values (for all P phases)
+    void computeSurfacevol(const int n,
+                           const int np,
+                           const double* A,
+                           const double* saturation,
+                           double* surfacevol);
+
 } // namespace Opm
 
 #endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED