From a816f4c06cb6f047057fa2dc047ae08ec7d9693b Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 17 Aug 2015 13:04:07 +0200 Subject: [PATCH 1/4] optimizing the formInterleavedSystem(). The current implementation avoids the using of formEllipticSystem() and vercatCollapseJacs(), which take a significant amount of computing time during the non-linear solutions. --- .../NewtonIterationBlackoilInterleaved.cpp | 25 +++++++++++++++---- .../NewtonIterationBlackoilInterleaved.hpp | 1 - 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 122c5dc64..2a5e52fab 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -107,12 +107,26 @@ namespace Opm // Form modified system. Eigen::SparseMatrix A; - V b; - formEllipticSystem(np, eqs, A, b); + // calculating the size for b + int size_b = 0; + for (int elem = 0; elem < np; ++elem) { + const int loc_size = eqs[elem].size(); + size_b += loc_size; + } + + V b(size_b); + + int pos = 0; + for (int elem = 0; elem < np; ++elem) { + const int loc_size = eqs[elem].size(); + b.segment(pos, loc_size) = eqs[elem].value(); + pos += loc_size; + } + assert(pos == size); // Create ISTL matrix with interleaved rows and columns (block structured). Mat istlA; - formInterleavedSystem(eqs, A, istlA); + formInterleavedSystem(eqs, istlA); // Solve reduced system. SolutionVector dx(SolutionVector::Zero(b.size())); @@ -185,7 +199,6 @@ namespace Opm void NewtonIterationBlackoilInterleaved::formInterleavedSystem(const std::vector& eqs, - const Eigen::SparseMatrix& A, Mat& istlA) const { const int np = eqs.size(); @@ -197,6 +210,7 @@ namespace Opm for (int phase = 0; phase < np; ++phase) { structure += eqs[phase].derivative()[0]; } + Eigen::SparseMatrix s = structure; // Create ISTL matrix with interleaved rows and columns (block structured). @@ -211,6 +225,7 @@ namespace Opm row.insert(ja[i]); } } + const int size = s.rows(); Span span[3] = { Span(size, 1, 0), Span(size, 1, size), @@ -221,7 +236,7 @@ namespace Opm MatrixBlockType block; for (int p1 = 0; p1 < np; ++p1) { for (int p2 = 0; p2 < np; ++p2) { - block[p1][p2] = A.coeff(span[p1][row], span[p2][col]); + block[p1][p2] = eqs[p1].derivative()[p2].coeff(row, col); } } istlA[row][col] = block; diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index 37b021983..45e7aa69d 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -156,7 +156,6 @@ namespace Opm } void formInterleavedSystem(const std::vector& eqs, - const Eigen::SparseMatrix& A, Mat& istlA) const; mutable int iterations_; From 7d7b05a12657aab8b48461cab4c5617f391330ca Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 17 Aug 2015 15:04:31 +0200 Subject: [PATCH 2/4] removing unused Span variable. --- opm/autodiff/NewtonIterationBlackoilInterleaved.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 2a5e52fab..edf990522 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -227,9 +227,6 @@ namespace Opm } const int size = s.rows(); - Span span[3] = { Span(size, 1, 0), - Span(size, 1, size), - Span(size, 1, 2*size) }; for (int row = 0; row < size; ++row) { for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) { const int col = ja[col_ix]; From ce4b25f0b8fb6f184e3039f81aea5fa367e9e16b Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 17 Aug 2015 15:12:13 +0200 Subject: [PATCH 3/4] removing the unused SparseMatrix A --- opm/autodiff/NewtonIterationBlackoilInterleaved.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index edf990522..c5bb58a5a 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -105,8 +105,6 @@ namespace Opm eqs[phase] = eqs[phase] * matbalscale[phase]; } - // Form modified system. - Eigen::SparseMatrix A; // calculating the size for b int size_b = 0; for (int elem = 0; elem < np; ++elem) { From 997616658b1dc8b593fd01aa803d0602ad041cd2 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 17 Aug 2015 15:33:33 +0200 Subject: [PATCH 4/4] adding 0th phase only once when computing sparse pattern --- opm/autodiff/NewtonIterationBlackoilInterleaved.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index c5bb58a5a..34de4d10f 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -205,7 +205,7 @@ namespace Opm // corresponding to the jacobians with respect to pressure. // Use addition to get to the union structure. Eigen::SparseMatrix structure = eqs[0].derivative()[0]; - for (int phase = 0; phase < np; ++phase) { + for (int phase = 1; phase < np; ++phase) { structure += eqs[phase].derivative()[0]; }