From 62cefb3a3e91300b701c7c65300cf7c7d4db61f8 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Tue, 7 Oct 2014 10:00:38 +0200 Subject: [PATCH 1/2] print warning when off-diagonal element is found in Schur complement instead of OPM_THROW. --- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index f7d3b2281..0a8b2cd01 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -241,7 +241,8 @@ namespace Opm // std::cout << "++++++++++++++++++++++++++++++++++++++++++++\n" // << D // << "++++++++++++++++++++++++++++++++++++++++++++\n" << std::endl; - OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); + std::cerr << "WARNING (ignored): Cannot do Schur complement with respect to non-diagonal block." << std::endl; + //OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); } V diag = D.diagonal(); Eigen::DiagonalMatrix invD = (1.0 / diag).matrix().asDiagonal(); @@ -296,7 +297,8 @@ namespace Opm // Find inv(D). const M& D = equation.derivative()[n]; if (!isDiagonal(D)) { - OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); + std::cerr << "WARNING (ignored): Cannot do Schur complement with respect to non-diagonal block." << std::endl; + //OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); } V diag = D.diagonal(); Eigen::DiagonalMatrix invD = (1.0 / diag).matrix().asDiagonal(); From 9f58ad54764d402974e6e5aa7fbdea86b00fe45a Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 9 Oct 2014 14:06:02 +0200 Subject: [PATCH 2/2] bugfix, Mismatched free() / delete / delete [] in DuneMatrix due to use of new instead of the matrix internal allocators. This fix also avoid the copying of the BCRSMatrix by providing a contructor that creates the DuneMatrix for a given Eigen SparseMatrix. --- opm/autodiff/DuneMatrix.hpp | 35 +++++++++++++++++---- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 20 ++---------- 2 files changed, 31 insertions(+), 24 deletions(-) diff --git a/opm/autodiff/DuneMatrix.hpp b/opm/autodiff/DuneMatrix.hpp index 3bcfd6a79..5c572ab97 100644 --- a/opm/autodiff/DuneMatrix.hpp +++ b/opm/autodiff/DuneMatrix.hpp @@ -1,5 +1,6 @@ /* Copyright 2014 SINTEF ICT, Applied Mathematics. + Copyright 2014 IRIS AS This file is part of the Open Porous Media project (OPM). @@ -26,6 +27,9 @@ #include +#include +#include + #include #include @@ -43,9 +47,26 @@ namespace Opm { 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) + { + // create BCRSMatrix from given CSR storage + init( rows, cols, ia, ja, sa ); + } + + /// \brief create an ISTL BCRSMatrix from a Eigen::SparseMatrix + DuneMatrix( const Eigen::SparseMatrix& matrix ) + { + // Create ISTL matrix. + const int rows = matrix.rows(); + const int cols = matrix.cols(); + const int* ia = matrix.outerIndexPtr(); + const int* ja = matrix.innerIndexPtr(); + const double* sa = matrix.valuePtr(); + // create BCRSMatrix from Eigen matrix + init( rows, cols, ia, ja, sa ); + } + + protected: + void init(const int rows, const int cols, const int* ia, const int* ja, const double* sa) { typedef Dune::BCRSMatrix< Dune::FieldMatrix > Super; typedef Super::block_type block_type; @@ -61,12 +82,14 @@ namespace Opm this->overflowsize = -1.0; #endif - this->a = new block_type[this->nnz]; + // make sure to use the allocators of this matrix + // because the same allocators are used to deallocate the data + this->a = this->allocator_.allocate(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(this->a)); - this->j.reset(new Super::size_type[this->nnz]); + this->j.reset(this->sizeAllocator_.allocate(this->nnz)); std::copy(ja, ja +this-> nnz, this->j.get()); - this->r = new Super::row_type[rows]; + this->r = rowAllocator_.allocate(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]); } diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index b3e3dbd67..41a415e3e 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -168,10 +168,10 @@ namespace Opm SolutionVector dx(SolutionVector::Zero(b.size())); // Create ISTL matrix. - Mat istlA = makeIstlMatrix(A); + DuneMatrix istlA( A ); // Create ISTL matrix for elliptic part. - Mat istlAe = makeIstlMatrix(A.topLeftCorner(nc, nc)); + DuneMatrix istlAe( A.topLeftCorner(nc, nc) ); // Construct operator, scalar product and vectors needed. typedef Dune::MatrixAdapter Operator; @@ -457,22 +457,6 @@ namespace Opm b = L * total_residual.value().matrix(); } - - - - - Mat makeIstlMatrix(const Eigen::SparseMatrix& matrix) - { - // Create ISTL matrix. - const int size = matrix.rows(); - const int* ia = matrix.outerIndexPtr(); - const int* ja = matrix.innerIndexPtr(); - const double* sa = matrix.valuePtr(); - return Opm::DuneMatrix(size, size, ia, ja, sa); - } - - - } // anonymous namespace