Merge remote-tracking branch 'upstream/master'

CMakeLists.txt

@@ -47,6 +47,7 @@ include (CMakeLists_files.cmake)
 
 macro (config_hook)
 	opm_need_version_of ("dune-common")
+	opm_need_version_of ("dune-istl")
 endmacro (config_hook)
 
 macro (prereqs_hook)

CMakeLists_files.cmake

@@ -96,6 +96,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/BlackoilPropsAdFromDeck.hpp
 	opm/autodiff/BlackoilPropsAdInterface.hpp
 	opm/autodiff/CPRPreconditioner.hpp
+	opm/autodiff/DuneMatrix.hpp
 	opm/autodiff/GeoProps.hpp
 	opm/autodiff/GridHelpers.hpp
 	opm/autodiff/ImpesTPFAAD.hpp

opm/autodiff/DuneMatrix.hpp

@@ -0,0 +1,78 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_DUNEMATRIX_HEADER_INCLUDED
+#define OPM_DUNEMATRIX_HEADER_INCLUDED
+
+#ifdef DUNE_BCRSMATRIX_HH
+#error This header must be included before any bcrsmatrix.hh is included (directly or indirectly)
+#endif
+
+#include <opm/core/utility/platform_dependent/disable_warnings.h>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/version.hh>
+
+// Include matrix header with hackery to make it possible to inherit.
+#define private protected
+#include <dune/istl/bcrsmatrix.hh>
+#undef private
+
+#include <opm/core/utility/platform_dependent/reenable_warnings.h>
+
+namespace Opm
+{
+
+    class DuneMatrix : public Dune::BCRSMatrix< Dune::FieldMatrix<double, 1, 1> >
+    {
+    public:
+        DuneMatrix(const int rows, const int cols, const int* ia, const int* ja, const double* sa)
+            // : build_mode(unknown), ready(built), n(rows), m(cols), nnz(ia[rows]),
+            //   allocationSize(nnz), r(0), a(0),
+            //   avg(0), overflowsize(-1.0)
+        {
+            typedef Dune::BCRSMatrix< Dune::FieldMatrix<double, 1, 1> > Super;
+            typedef Super::block_type block_type;
+            this->build_mode = Super::unknown;
+            this->ready = Super::built;
+            this->n = rows;
+            this->m = cols;
+            this->nnz = ia[rows];
+
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
+            this->allocationSize = this->nnz;
+            this->avg = 0;
+            this->overflowsize = -1.0;
+#endif
+
+            this->a = new block_type[this->nnz];
+            static_assert(sizeof(block_type) == sizeof(double), "This constructor requires a block type that is the same as a double.");
+            std::copy(sa, sa + this->nnz, reinterpret_cast<double*>(this->a));
+            this->j.reset(new Super::size_type[this->nnz]);
+            std::copy(ja, ja +this-> nnz, this->j.get());
+            this->r = new Super::row_type[rows];
+            for (int row = 0; row < rows; ++row) {
+                this->r[row].set(ia[row+1] - ia[row], this->a + ia[row], this->j.get() + ia[row]);
+            }
+        }
+    };
+
+} // namespace Opm
+
+#endif // OPM_DUNEMATRIX_HEADER_INCLUDED

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -580,6 +580,7 @@ namespace {
         const ADB&              rs    = state.rs;
         const ADB&              rv    = state.rv;
 
+        const std::vector<ADB> pressures = computePressures(state);
         const std::vector<PhasePresence> cond = phaseCondition();
 
         const ADB pv_mult = poroMult(press);
@@ -588,7 +589,7 @@ namespace {
         for (int phase = 0; phase < maxnp; ++phase) {
             if (active_[ phase ]) {
                 const int pos = pu.phase_pos[ phase ];
-                rq_[pos].b = fluidReciprocFVF(phase, press, rs, rv, cond, cells_);
+                rq_[pos].b = fluidReciprocFVF(phase, pressures[phase], rs, rv, cond, cells_);
                 rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
                 // DUMP(rq_[pos].b);
                 // DUMP(rq_[pos].accum[aix]);

opm/autodiff/NewtonIterationBlackoilCPR.cpp

@@ -19,6 +19,8 @@
 
 #include <config.h>
 
+#include <opm/autodiff/DuneMatrix.hpp>
+
 #include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
@@ -29,7 +31,7 @@
 #include <opm/core/utility/platform_dependent/disable_warnings.h>
 
 #include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
+// #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/operators.hh>
 #include <dune/istl/io.hh>
 #include <dune/istl/owneroverlapcopy.hh>
@@ -463,23 +465,10 @@ namespace Opm
         {
             // Create ISTL matrix.
             const int size = matrix.rows();
-            const int nonzeros = matrix.nonZeros();
             const int* ia = matrix.outerIndexPtr();
             const int* ja = matrix.innerIndexPtr();
             const double* sa = matrix.valuePtr();
-            Mat A(size, size, nonzeros, Mat::row_wise);
-            for (Mat::CreateIterator row = A.createbegin(); row != A.createend(); ++row) {
-                const int ri = row.index();
-                for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
-                    row.insert(ja[i]);
-                }
-            }
-            for (int ri = 0; ri < size; ++ri) {
-                for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
-                    A[ri][ja[i]] = sa[i];
-                }
-            }
-            return A;
+            return Opm::DuneMatrix(size, size, ia, ja, sa);
         }