From acd58f5272270bbbb6ade15220228efac9f5b86c Mon Sep 17 00:00:00 2001 From: babrodtk Date: Thu, 27 Aug 2015 14:53:55 +0200 Subject: [PATCH] Fixed commented out functions and some warnings --- opm/autodiff/AutoDiffMatrix.hpp | 12 ++++---- .../NewtonIterationBlackoilSimple.cpp | 5 ++-- opm/autodiff/NewtonIterationUtilities.cpp | 30 +++++++++++-------- opm/autodiff/fastSparseProduct.hpp | 2 ++ 4 files changed, 28 insertions(+), 21 deletions(-) diff --git a/opm/autodiff/AutoDiffMatrix.hpp b/opm/autodiff/AutoDiffMatrix.hpp index 4f0cbb7a9..2c5c8c5fc 100644 --- a/opm/autodiff/AutoDiffMatrix.hpp +++ b/opm/autodiff/AutoDiffMatrix.hpp @@ -47,10 +47,10 @@ namespace Opm - AutoDiffMatrix(const int rows, const int cols) + AutoDiffMatrix(const int num_rows, const int num_cols) : type_(Z), - rows_(rows), - cols_(cols) + rows_(num_rows), + cols_(num_cols) { } @@ -59,10 +59,10 @@ namespace Opm enum CreationType { ZeroMatrix, IdentityMatrix }; - AutoDiffMatrix(const CreationType t, const int rows) + AutoDiffMatrix(const CreationType t, const int num_rows) : type_(t == ZeroMatrix ? Z : I), - rows_(rows), - cols_(rows) + rows_(num_rows), + cols_(num_rows) { } diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.cpp b/opm/autodiff/NewtonIterationBlackoilSimple.cpp index abdf31d9e..dffe989a3 100644 --- a/opm/autodiff/NewtonIterationBlackoilSimple.cpp +++ b/opm/autodiff/NewtonIterationBlackoilSimple.cpp @@ -47,7 +47,6 @@ namespace Opm NewtonIterationBlackoilSimple::SolutionVector NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const { - /* typedef LinearisedBlackoilResidual::ADB ADB; const int np = residual.material_balance_eq.size(); ADB mass_res = residual.material_balance_eq[0]; @@ -57,7 +56,8 @@ namespace Opm const ADB well_res = vertcat(residual.well_flux_eq, residual.well_eq); const ADB total_residual = collapseJacs(vertcat(mass_res, well_res)); - const Eigen::SparseMatrix matr = total_residual.derivative()[0]; + Eigen::SparseMatrix matr; + total_residual.derivative()[0].toSparse(matr); SolutionVector dx(SolutionVector::Zero(total_residual.size())); Opm::LinearSolverInterface::LinearSolverReport rep @@ -74,7 +74,6 @@ namespace Opm "Linear solver convergence failure."); } return dx; - */ } const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const diff --git a/opm/autodiff/NewtonIterationUtilities.cpp b/opm/autodiff/NewtonIterationUtilities.cpp index 8919a31af..62f88404a 100644 --- a/opm/autodiff/NewtonIterationUtilities.cpp +++ b/opm/autodiff/NewtonIterationUtilities.cpp @@ -39,6 +39,7 @@ namespace Opm typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; + typedef Eigen::SparseMatrix S; std::vector eliminateVariable(const std::vector& eqs, const int n) @@ -75,7 +76,8 @@ namespace Opm const Sp Di = solver.solve(id); // compute inv(D)*bn for the update of the right hand side - const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix()); + ADB::V eqs_n_v = eqs[n].value(); + const Eigen::VectorXd& Dibn = solver.solve(eqs_n_v.matrix()); std::vector vals(num_eq); // Number n will remain empty. std::vector> jacs(num_eq); // Number n will remain empty. @@ -198,7 +200,6 @@ namespace Opm Eigen::SparseMatrix& A, V& b) { - /* if (num_phases != 3) { OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases."); } @@ -215,14 +216,17 @@ namespace Opm // The l1 block indicates if the equation for a given cell and phase is // sufficiently strong on the diagonal. Block l1 = Block::Zero(n, num_phases); - for (int phase = 0; phase < num_phases; ++phase) { - const M& J = eqs[phase].derivative()[0]; - V dj = J.diagonal().cwiseAbs(); - V sod = V::Zero(n); - for (int elem = 0; elem < n; ++elem) { - sod(elem) = J.col(elem).cwiseAbs().sum() - dj(elem); + { + S J; + for (int phase = 0; phase < num_phases; ++phase) { + eqs[phase].derivative()[0].toSparse(J); + V dj = J.diagonal().cwiseAbs(); + V sod = V::Zero(n); + for (int elem = 0; elem < n; ++elem) { + sod(elem) = J.col(elem).cwiseAbs().sum() - dj(elem); + } + l1.col(phase) = (dj/sod > ratio_limit).cast(); } - l1.col(phase) = (dj/sod > ratio_limit).cast(); } // By default, replace first equation with sum of all phase equations. @@ -269,16 +273,18 @@ namespace Opm t.emplace_back(i3[ii], i1[ii], l31(ii)); t.emplace_back(i3[ii], i3[ii], l33(ii)); } - M L(3*n, 3*n); + S L(3*n, 3*n); L.setFromTriplets(t.begin(), t.end()); // Combine in single block. ADB total_residual = vertcatCollapseJacs(eqs); + S derivative; + total_residual.derivative()[0].toSparse(derivative); + // Create output as product of L with equations. - A = L * total_residual.derivative()[0]; + A = L * derivative; b = L * total_residual.value().matrix(); - */ } diff --git a/opm/autodiff/fastSparseProduct.hpp b/opm/autodiff/fastSparseProduct.hpp index c5ba0defe..4fa2f7527 100644 --- a/opm/autodiff/fastSparseProduct.hpp +++ b/opm/autodiff/fastSparseProduct.hpp @@ -151,6 +151,7 @@ inline void fastDiagSparseProduct(// const Eigen::DiagonalMatrix::InnerIterator It; for (It it(res, col); it; ++it) { @@ -171,6 +172,7 @@ inline void fastSparseDiagProduct(const Eigen::SparseMatrix& lhs, // Multiply columns by diagonal rhs. int n = res.cols(); +#pragma omp parallel for for (int col = 0; col < n; ++col) { typedef Eigen::SparseMatrix::InnerIterator It; for (It it(res, col); it; ++it) {