From 1a0e068f44f4e17bee50e26d85fe4ac685ae4125 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Thu, 5 Jan 2012 21:39:33 +0100
Subject: [PATCH] Implemented matrix() method.

---
 opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 57 +++++++++++++++++--
 opm/core/fluid/BlackoilPropertiesFromDeck.hpp |  4 ++
 2 files changed, 57 insertions(+), 4 deletions(-)

diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp
index 84c005839..a66c5f57a 100644
--- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp
+++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp
@@ -104,13 +104,62 @@ namespace Opm
     void BlackoilPropertiesFromDeck::matrix(const int n,
                                             const double* p,
                                             const double* z,
-                                            const int* cells,
+                                            const int* /*cells*/,
                                             double* A,
                                             double* dAdp) const
     {
-        THROW("BlackoilPropertiesFromDeck::matrix() not yet implemented.");
-    }
+        const int np = numPhases();
+        B_.resize(n*np);
+        R_.resize(n*np);
+        if (dAdp) {
+            dB_.resize(n*np);
+            dR_.resize(n*np);
+            pvt_.dBdp(n, p, z, &B_[0], &dB_[0]);
+            pvt_.dRdp(n, p, z, &R_[0], &dR_[0]);
+        } else {
+            pvt_.B(n, p, z, &B_[0]);
+            pvt_.R(n, p, z, &R_[0]);
+        }
+        const int* phase_pos = pvt_.phasePosition();
+        bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
+            pvt_.phaseUsed()[BlackoilPhases::Vapour];
+        const int o = phase_pos[BlackoilPhases::Liquid];
+        const int g = phase_pos[BlackoilPhases::Vapour];
 
+        // Compute A matrix
+#pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            double* m = A + i*np*np;
+            std::fill(m, m + np*np, 0.0);
+            // Diagonal entries.
+            for (int phase = 0; phase < np; ++phase) {
+                m[phase + phase*np] = 1.0/B_[i*np + phase];
+            }
+            // Off-diagonal entries.
+            if (oil_and_gas) {
+                m[o + g*np] = R_[i*np + g]/B_[i*np + g];
+                m[g + o*np] = R_[i*np + o]/B_[i*np + o];
+            }
+        }
+
+        // Derivative of A matrix.
+        if (dAdp) {
+#pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                double* m = dAdp + i*np*np;
+                std::fill(m, m + np*np, 0.0);
+                // Diagonal entries.
+                for (int phase = 0; phase < np; ++phase) {
+                    m[phase + phase*np] = -dB_[i*np + phase]/B_[i*np + phase]*B_[i*np + phase];
+                }
+                // Off-diagonal entries.
+                if (oil_and_gas) {
+                    m[o + g*np] = m[g + g*np]*R_[i*np + g] + dR_[i*np + g]/B_[i*np + g];
+                    m[g + o*np] = m[o + o*np]*R_[i*np + o] + dR_[i*np + o]/B_[i*np + o];
+                }
+            }
+        }
+    }
 
     /// \param[in]  n      Number of data points.
     /// \param[in]  A      Array of nP^2 values, where the P^2 values for a cell give the
@@ -122,7 +171,7 @@ namespace Opm
                                              const double* A,
                                              double* rho) const
     {
-        int np = numPhases();
+        const int np = numPhases();
         const double* sdens = pvt_.surfaceDensities();
 #pragma omp parallel for
         for (int i = 0; i < n; ++i) {
diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp
index 12913c6d6..731cbd286 100644
--- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp
+++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp
@@ -145,6 +145,10 @@ namespace Opm
         RockFromDeck rock_;
         BlackoilPvtProperties pvt_;
         SaturationPropsFromDeck satprops_;
+        mutable std::vector<double> B_;
+        mutable std::vector<double> dB_;
+        mutable std::vector<double> R_;
+        mutable std::vector<double> dR_;
     };