From 2f0d1a4c94cb7fbaa63c10b4c284a33f23a34ff3 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sat, 2 Sep 2017 10:48:13 +0200 Subject: [PATCH 01/10] [matrixconverter] Add first version of matrix converter * Convertes multitype block matrix to plain BCRS matrix * there is still a bug when numEq != 1 --- dumux/implicit/staggered/matrixconverter.hh | 188 +++++++++++++++++++ dumux/implicit/staggered/newtoncontroller.hh | 147 ++++++++------- 2 files changed, 265 insertions(+), 70 deletions(-) create mode 100644 dumux/implicit/staggered/matrixconverter.hh diff --git a/dumux/implicit/staggered/matrixconverter.hh b/dumux/implicit/staggered/matrixconverter.hh new file mode 100644 index 0000000000..c967cf8702 --- /dev/null +++ b/dumux/implicit/staggered/matrixconverter.hh @@ -0,0 +1,188 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program 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 2 of the License, or * + * (at your option) any later version. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief A 2p1cni specific controller for the newton solver. + * + * This controller 'knows' what a 'physically meaningful' solution is + * which allows the newton method to abort quicker if the solution is + * way out of bounds. + */ +#ifndef DUMUX_MATRIX_CONVERTER +#define DUMUX_MATRIX_CONVERTER + +#include "properties.hh" + +#include +#include +#include "newtonconvergencewriter.hh" +#include + +namespace Dumux { + +/*! + * \ingroup PNMModel + * \brief A PNM specific controller for the newton solver. + * + * This controller 'knows' what a 'physically meaningful' solution is + * which allows the newton method to abort quicker if the solution is + * way out of bounds. + */ + +template +class MatrixConverter +{ + using MatrixBlock = typename Dune::FieldMatrix; + using BCRSMatrix = typename Dune::BCRSMatrix; + +public: + + // copy the matrix and the vector to types the IterativeSolverBackend can handle + static auto multiTypeToBCRSMatrix(const MultiTypeBlockMatrix &A) + { + // get the number of rows for the converted matrix + const auto numRows = getNumRows_(A); + + auto M = BCRSMatrix(numRows, numRows, BCRSMatrix::random); + + setOccupationPattern_(M, A); + copyValues_(M, A); + + return M; + } + +private: + + /*! + * \brief Sets the occupation pattern (i.e. indices) for the converted matrix + * + * \param M The converted matrix + * \param A The original subMatrix of the multitype blockmatrix + */ + static void setOccupationPattern_(BCRSMatrix& M, const MultiTypeBlockMatrix& A) + { + // set the rowsizes + const auto numRows = M.N(); + Dune::MatrixIndexSet occupationPattern; + occupationPattern.resize(numRows, numRows); + + // lambda function to fill the occupation pattern + auto addIndices = [&occupationPattern](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol) + { + using BlockType = typename std::decay_t::block_type; + const auto blockSizeI = BlockType::rows; + const auto blockSizeJ = BlockType::cols; + for(auto row = subMatrix.begin(); row != subMatrix.end(); ++row) + for(auto col = row->begin(); col != row->end(); ++col) + for(std::size_t i = 0; i < blockSizeI; ++i) + for(std::size_t j = 0; j < blockSizeJ; ++j) + occupationPattern.add(startRow + row.index()*blockSizeI + i, startCol + col.index()*blockSizeJ + j); + }; + + int rowIndex = 0; + Dune::Hybrid::forEach(A, [&addIndices, &rowIndex, &occupationPattern, numRows](const auto& rowOfMultiTypeMatrix) + { + int colIndex = 0; + Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [&addIndices, &occupationPattern, &colIndex, &rowIndex, numRows](const auto& subMatrix) + { + addIndices(subMatrix, rowIndex, colIndex); + + const auto numEq = std::decay_t::block_type::cols; + colIndex += numEq*subMatrix.M(); + + // if we arrived at the right side of the matrix, increase the row index + if(colIndex == numRows) + rowIndex += numEq*subMatrix.N(); + }); + }); + + occupationPattern.exportIdx(M); + } + + /*! + * \brief Sets the occupation pattern (i.e. indices) for the converted matrix + * + * \param M The converted matrix + * \param A The original subMatrix of the multitype blockmatrix + */ + static void copyValues_(BCRSMatrix& M, const MultiTypeBlockMatrix& A) + { + // set the rowsizes + const auto numRows = M.N(); + + // lambda function to copy the values + auto copyValues = [&M](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol) + { + using BlockType = typename std::decay_t::block_type; + const auto blockSizeI = BlockType::rows; + const auto blockSizeJ = BlockType::cols; + for (auto row = subMatrix.begin(); row != subMatrix.end(); ++row) + for (auto col = row->begin(); col != row->end(); ++col) + for (std::size_t i = 0; i < blockSizeI; ++i) + for (std::size_t j = 0; j < blockSizeJ; ++j) + M[startRow + row.index()*blockSizeI + i][startCol + col.index()*blockSizeJ + j] = subMatrix[row.index()][col.index()][i][j]; + + }; + + int rowIndex = 0; + Dune::Hybrid::forEach(A, [©Values, &rowIndex, numRows](const auto& rowOfMultiTypeMatrix) + { + int colIndex = 0; + Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [©Values, &colIndex, &rowIndex, numRows](const auto& subMatrix) + { + copyValues(subMatrix, rowIndex, colIndex); + + const auto numEq = std::decay_t::block_type::cols; + colIndex += numEq*subMatrix.M(); + + // if we arrived at the right side of the matrix, increase the row index + if(colIndex == numRows) + rowIndex += numEq*subMatrix.N(); + }); + }); + } + + /*! + * \brief Calculates the total number of rows (== number of cols) for the converted matrix, assuming a block size of one + * + * \param A The original multitype blockmatrix + */ + static std::size_t getNumRows_(const MultiTypeBlockMatrix& A) + { + // problem-agnostic version + std::size_t numRows = 0; + Dune::Hybrid::forEach(A[Dune::Indices::_0], [&numRows](const auto& subMatrixInFirstRow) + { + // std::cout << "in matrix of first row:\n"; + // std::cout << "numRows: " << subMatrixInFirstRow.N() << ", numCols: " << subMatrixInFirstRow.M() << std::endl; + + // The number of rows of the individual submatrice's block equals the respective number of equations. + const auto numEq = std::decay_t::block_type::cols; // TODO: is this allways correct? + numRows += numEq * subMatrixInFirstRow.M(); + // std::cout << "numEq: " << numEq << std::endl; + // std::cout << "totalNumRows: " << numRows << std::endl; + }); + + return numRows; + } + +}; +} + +#endif diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh index edb72ac64e..6e964fc768 100644 --- a/dumux/implicit/staggered/newtoncontroller.hh +++ b/dumux/implicit/staggered/newtoncontroller.hh @@ -33,6 +33,10 @@ #include #include "newtonconvergencewriter.hh" +#include "matrixconverter.hh" + +#include + namespace Dumux { namespace Properties @@ -53,6 +57,7 @@ namespace Properties * which allows the newton method to abort quicker if the solution is * way out of bounds. */ + template class StaggeredNewtonController : public NewtonController { @@ -77,6 +82,75 @@ public: : ParentType(problem) {} + + auto createMatrix(JacobianMatrix &A) + { + // copy the matrix and the vector to types the IterativeSolverBackend can handle + + using MatrixBlock = typename Dune::FieldMatrix; + using SparseMatrix = typename Dune::BCRSMatrix; + + // problem-agnostic version + std::size_t numRows = 0; + Dune::Hybrid::forEach(A[Dune::Indices::_0], [&numRows](const auto& subMatrixInFirstRow) + { + std::cout << "in matrix of first row:\n"; + std::cout << "numRows: " << subMatrixInFirstRow.N() << ", numCols: " << subMatrixInFirstRow.M() << std::endl; + + // The number of rows of the individual submatrice's block equals the respective number of equations. + const auto numEq = std::decay_t::block_type::rows; // TODO: is this allways correct? + numRows += numEq * subMatrixInFirstRow.M(); + }); + + // /*const std::size_t*/ numRows = this->model_().numCellCenterDofs()*numEqCellCenter + this->model_().numFaceDofs()*numEqFace; + + + // set the rowsizes + Dune::MatrixIndexSet occupationPattern; + occupationPattern.resize(numRows, numRows); + + // lambda function to fill the occupation pattern + auto addIndices = [&occupationPattern](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol) + { + using BlockType = typename std::decay_t::block_type; + const auto blockSizeI = BlockType::rows; + const auto blockSizeJ = BlockType::cols; + for(auto row = subMatrix.begin(); row != subMatrix.end(); ++row) + for(auto col = row->begin(); col != row->end(); ++col) + for(std::size_t i = 0; i < blockSizeI; ++i) + for(std::size_t j = 0; j < blockSizeJ; ++j) + occupationPattern.add(startRow + row.index()*blockSizeI + i, startCol + col.index()*blockSizeJ + j); + + }; + + + int rowIndex = 0; + + Dune::Hybrid::forEach(A, [&rowIndex, &addIndices, numRows](const auto& rowOfMultiTypeMatrix) + { + int colIndex = 0; + + std::cout << "processing one row of multitypeblockmatrix " << std::endl; + Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [&colIndex, &addIndices, &rowIndex, numRows](const auto& subMatrix) + { + std::cout << "numRows: " << subMatrix.N() << ", numCols: " << subMatrix.M() << std::endl; + addIndices(subMatrix, rowIndex, colIndex); + + const auto numEq = std::decay_t::block_type::rows; + colIndex += numEq*subMatrix.M(); + + // if we arrived at the right side of the matrix, increase the row index + if(colIndex == numRows) + rowIndex += numEq*subMatrix.N(); + }); + }); + + auto M = SparseMatrix(numRows, numRows, SparseMatrix::random); + occupationPattern.exportIdx(M); + + return M; + } + /*! * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$. * @@ -96,6 +170,7 @@ public: SolutionVector &x, SolutionVector &b) { + // createMatrix(A); try { if (this->numSteps_ == 0) @@ -115,75 +190,7 @@ public: assert(numRows == numCols); // create the bcrs matrix the IterativeSolver backend can handle - auto M = SparseMatrix(numRows, numCols, SparseMatrix::random); - - // set the rowsizes - // A11 and A12 - for (auto row = A[cellCenterIdx][cellCenterIdx].begin(); row != A[cellCenterIdx][cellCenterIdx].end(); ++row) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - M.setrowsize(numEqCellCenter*row.index() + i, row->size()*numEqCellCenter); - for (auto row = A[cellCenterIdx][faceIdx].begin(); row != A[cellCenterIdx][faceIdx].end(); ++row) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - M.setrowsize(numEqCellCenter*row.index() + i, M.getrowsize(numEqCellCenter*row.index() + i) + row->size()*numEqFace); - // A21 and A22 - for (auto row = A[faceIdx][cellCenterIdx].begin(); row != A[faceIdx][cellCenterIdx].end(); ++row) - for (std::size_t i = 0; i < numEqFace; ++i) - M.setrowsize(numEqFace*row.index() + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter, row->size()*numEqCellCenter); - for (auto row = A[faceIdx][faceIdx].begin(); row != A[faceIdx][faceIdx].end(); ++row) - for (std::size_t i = 0; i < numEqFace; ++i) - M.setrowsize(numEqFace*row.index() + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter, M.getrowsize(numEqFace*row.index() + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter) + row->size()*numEqFace); - M.endrowsizes(); - - // set the indices - for (auto row = A[cellCenterIdx][cellCenterIdx].begin(); row != A[cellCenterIdx][cellCenterIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - M.addindex(row.index()*numEqCellCenter + i, col.index()*numEqCellCenter + j); - - for (auto row = A[cellCenterIdx][faceIdx].begin(); row != A[cellCenterIdx][faceIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - M.addindex(row.index()*numEqCellCenter + i, col.index()*numEqFace + j + A[cellCenterIdx][cellCenterIdx].M()*numEqCellCenter); - - for (auto row = A[faceIdx][cellCenterIdx].begin(); row != A[faceIdx][cellCenterIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqFace; ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - M.addindex(row.index()*numEqFace + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter, col.index()*numEqCellCenter + j); - - for (auto row = A[faceIdx][faceIdx].begin(); row != A[faceIdx][faceIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqFace; ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - M.addindex(row.index()*numEqFace + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter, col.index()*numEqFace + j + A[cellCenterIdx][cellCenterIdx].M()*numEqCellCenter); - M.endindices(); - - // copy values - for (auto row = A[cellCenterIdx][cellCenterIdx].begin(); row != A[cellCenterIdx][cellCenterIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - M[row.index()*numEqCellCenter + i][col.index()*numEqCellCenter + j] = A[cellCenterIdx][cellCenterIdx][row.index()][col.index()][i][j]; - - for (auto row = A[cellCenterIdx][faceIdx].begin(); row != A[cellCenterIdx][faceIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - M[row.index()*numEqCellCenter + i][col.index()*numEqFace + j + A[cellCenterIdx][cellCenterIdx].M()*numEqCellCenter] = A[cellCenterIdx][faceIdx][row.index()][col.index()][i][j]; - - for (auto row = A[faceIdx][cellCenterIdx].begin(); row != A[faceIdx][cellCenterIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqFace; ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - M[row.index()*numEqFace + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter][col.index()*numEqCellCenter + j] = A[faceIdx][cellCenterIdx][row.index()][col.index()][i][j]; - - for (auto row = A[faceIdx][faceIdx].begin(); row != A[faceIdx][faceIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqFace; ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - M[row.index()*numEqFace + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter][col.index()*numEqFace + j + A[cellCenterIdx][cellCenterIdx].M()*numEqCellCenter] = A[faceIdx][faceIdx][row.index()][col.index()][i][j]; + const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); // create the vector the IterativeSolver backend can handle using VectorBlock = typename Dune::FieldVector; @@ -202,7 +209,7 @@ public: // printmatrix(std::cout, M, "", ""); // solve - bool converged = this->linearSolver_.solve(M, y, bTmp); + const bool converged = this->linearSolver_.solve(M, y, bTmp); // copy back the result y into x for (std::size_t i = 0; i < x[cellCenterIdx].N(); ++i) -- GitLab From b37d3ac0fa6db76f02654bda8aa58841cd5c63cf Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sat, 2 Sep 2017 13:42:13 +0200 Subject: [PATCH 02/10] [matrixconverter] Fix bug and clean-up * use correct values of block size * remove unused code * clean-up docu --- dumux/implicit/staggered/matrixconverter.hh | 46 +++++++-------------- 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/dumux/implicit/staggered/matrixconverter.hh b/dumux/implicit/staggered/matrixconverter.hh index c967cf8702..d831874e40 100644 --- a/dumux/implicit/staggered/matrixconverter.hh +++ b/dumux/implicit/staggered/matrixconverter.hh @@ -18,31 +18,21 @@ *****************************************************************************/ /*! * \file - * \brief A 2p1cni specific controller for the newton solver. + * \brief A helper classe that converts a Dune::MultiTypeBlockMatrix into a plain Dune::BCRSMatrix * - * This controller 'knows' what a 'physically meaningful' solution is - * which allows the newton method to abort quicker if the solution is - * way out of bounds. */ #ifndef DUMUX_MATRIX_CONVERTER #define DUMUX_MATRIX_CONVERTER -#include "properties.hh" -#include -#include -#include "newtonconvergencewriter.hh" -#include +#include namespace Dumux { /*! - * \ingroup PNMModel - * \brief A PNM specific controller for the newton solver. + *A helper classe that converts a Dune::MultiTypeBlockMatrix into a plain Dune::BCRSMatrix + * TODO: allow block sizes for BCRSMatrix other than 1x1 ? * - * This controller 'knows' what a 'physically meaningful' solution is - * which allows the newton method to abort quicker if the solution is - * way out of bounds. */ template @@ -103,12 +93,13 @@ private: { addIndices(subMatrix, rowIndex, colIndex); - const auto numEq = std::decay_t::block_type::cols; - colIndex += numEq*subMatrix.M(); + using SubBlockType = typename std::decay_t::block_type; - // if we arrived at the right side of the matrix, increase the row index + colIndex += SubBlockType::cols * subMatrix.M(); + + // if we have arrived at the right side of the matrix, increase the row index if(colIndex == numRows) - rowIndex += numEq*subMatrix.N(); + rowIndex += SubBlockType::rows * subMatrix.N(); }); }); @@ -148,12 +139,13 @@ private: { copyValues(subMatrix, rowIndex, colIndex); - const auto numEq = std::decay_t::block_type::cols; - colIndex += numEq*subMatrix.M(); + using SubBlockType = typename std::decay_t::block_type; + + colIndex += SubBlockType::cols * subMatrix.M(); - // if we arrived at the right side of the matrix, increase the row index + // if we have arrived at the right side of the matrix, increase the row index if(colIndex == numRows) - rowIndex += numEq*subMatrix.N(); + rowIndex += SubBlockType::rows * subMatrix.N(); }); }); } @@ -165,18 +157,12 @@ private: */ static std::size_t getNumRows_(const MultiTypeBlockMatrix& A) { - // problem-agnostic version std::size_t numRows = 0; Dune::Hybrid::forEach(A[Dune::Indices::_0], [&numRows](const auto& subMatrixInFirstRow) { - // std::cout << "in matrix of first row:\n"; - // std::cout << "numRows: " << subMatrixInFirstRow.N() << ", numCols: " << subMatrixInFirstRow.M() << std::endl; - - // The number of rows of the individual submatrice's block equals the respective number of equations. - const auto numEq = std::decay_t::block_type::cols; // TODO: is this allways correct? + // the number of cols of the individual submatrice's block equals the respective number of equations. + const auto numEq = std::decay_t::block_type::cols; numRows += numEq * subMatrixInFirstRow.M(); - // std::cout << "numEq: " << numEq << std::endl; - // std::cout << "totalNumRows: " << numRows << std::endl; }); return numRows; -- GitLab From 67e13d8278875f317d6f542a087067333cbb55f8 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sat, 2 Sep 2017 21:09:20 +0200 Subject: [PATCH 03/10] [matrixconverter] Move to dumux/linear --- dumux/implicit/staggered/newtoncontroller.hh | 4 +--- dumux/{implicit/staggered => linear}/matrixconverter.hh | 0 2 files changed, 1 insertion(+), 3 deletions(-) rename dumux/{implicit/staggered => linear}/matrixconverter.hh (100%) diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh index 6e964fc768..4090447914 100644 --- a/dumux/implicit/staggered/newtoncontroller.hh +++ b/dumux/implicit/staggered/newtoncontroller.hh @@ -31,11 +31,9 @@ #include #include +#include #include "newtonconvergencewriter.hh" -#include "matrixconverter.hh" - -#include namespace Dumux { diff --git a/dumux/implicit/staggered/matrixconverter.hh b/dumux/linear/matrixconverter.hh similarity index 100% rename from dumux/implicit/staggered/matrixconverter.hh rename to dumux/linear/matrixconverter.hh -- GitLab From 96c1e24e8e9a9afffded46bde7d5e6af56f79608 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sat, 2 Sep 2017 21:11:08 +0200 Subject: [PATCH 04/10] [matrixconverter] Improve and clean-up * add includes * improve docu * use Dune::Hybrid::elementAt insteand of [] for increased flexibility --- dumux/linear/matrixconverter.hh | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/dumux/linear/matrixconverter.hh b/dumux/linear/matrixconverter.hh index d831874e40..77828b0428 100644 --- a/dumux/linear/matrixconverter.hh +++ b/dumux/linear/matrixconverter.hh @@ -24,7 +24,12 @@ #ifndef DUMUX_MATRIX_CONVERTER #define DUMUX_MATRIX_CONVERTER - +#include +#include +#include +#include +#include +#include #include namespace Dumux { @@ -43,14 +48,21 @@ class MatrixConverter public: - // copy the matrix and the vector to types the IterativeSolverBackend can handle + /*! + * \brief Converts the matrix to a type the IterativeSolverBackend can handle + * + * \param M The converted matrix + * \param A The original multitype blockmatrix + */ static auto multiTypeToBCRSMatrix(const MultiTypeBlockMatrix &A) { - // get the number of rows for the converted matrix + // get the size for the converted matrix const auto numRows = getNumRows_(A); + // create an empty BCRS matrix with 1x1 blocks auto M = BCRSMatrix(numRows, numRows, BCRSMatrix::random); + // set the occupation pattern and copy the values setOccupationPattern_(M, A); copyValues_(M, A); @@ -60,14 +72,14 @@ public: private: /*! - * \brief Sets the occupation pattern (i.e. indices) for the converted matrix + * \brief Sets the occupation pattern and indices for the converted matrix * * \param M The converted matrix - * \param A The original subMatrix of the multitype blockmatrix + * \param A The original multitype blockmatrix */ static void setOccupationPattern_(BCRSMatrix& M, const MultiTypeBlockMatrix& A) { - // set the rowsizes + // prepare the occupation pattern const auto numRows = M.N(); Dune::MatrixIndexSet occupationPattern; occupationPattern.resize(numRows, numRows); @@ -85,6 +97,7 @@ private: occupationPattern.add(startRow + row.index()*blockSizeI + i, startCol + col.index()*blockSizeJ + j); }; + // fill the pattern int rowIndex = 0; Dune::Hybrid::forEach(A, [&addIndices, &rowIndex, &occupationPattern, numRows](const auto& rowOfMultiTypeMatrix) { @@ -114,7 +127,7 @@ private: */ static void copyValues_(BCRSMatrix& M, const MultiTypeBlockMatrix& A) { - // set the rowsizes + // get number of rows const auto numRows = M.N(); // lambda function to copy the values @@ -151,14 +164,15 @@ private: } /*! - * \brief Calculates the total number of rows (== number of cols) for the converted matrix, assuming a block size of one + * \brief Calculates the total number of rows (== number of cols) for the converted matrix, assuming a block size of 1x1 * * \param A The original multitype blockmatrix */ static std::size_t getNumRows_(const MultiTypeBlockMatrix& A) { + // iterate over the first row of the multitype blockmatrix std::size_t numRows = 0; - Dune::Hybrid::forEach(A[Dune::Indices::_0], [&numRows](const auto& subMatrixInFirstRow) + Dune::Hybrid::forEach(Dune::Hybrid::elementAt(A, Dune::Indices::_0), [&numRows](const auto& subMatrixInFirstRow) { // the number of cols of the individual submatrice's block equals the respective number of equations. const auto numEq = std::decay_t::block_type::cols; -- GitLab From a4f8d42b24cfcca1ff350b725e1e0a9ffbd68b5c Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sat, 2 Sep 2017 21:11:52 +0200 Subject: [PATCH 05/10] [staggered][newtoncontroller] Clean-up * remove unused code --- dumux/implicit/staggered/newtoncontroller.hh | 85 ++------------------ 1 file changed, 5 insertions(+), 80 deletions(-) diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh index 4090447914..e297bf30ae 100644 --- a/dumux/implicit/staggered/newtoncontroller.hh +++ b/dumux/implicit/staggered/newtoncontroller.hh @@ -34,7 +34,6 @@ #include #include "newtonconvergencewriter.hh" - namespace Dumux { namespace Properties @@ -80,75 +79,6 @@ public: : ParentType(problem) {} - - auto createMatrix(JacobianMatrix &A) - { - // copy the matrix and the vector to types the IterativeSolverBackend can handle - - using MatrixBlock = typename Dune::FieldMatrix; - using SparseMatrix = typename Dune::BCRSMatrix; - - // problem-agnostic version - std::size_t numRows = 0; - Dune::Hybrid::forEach(A[Dune::Indices::_0], [&numRows](const auto& subMatrixInFirstRow) - { - std::cout << "in matrix of first row:\n"; - std::cout << "numRows: " << subMatrixInFirstRow.N() << ", numCols: " << subMatrixInFirstRow.M() << std::endl; - - // The number of rows of the individual submatrice's block equals the respective number of equations. - const auto numEq = std::decay_t::block_type::rows; // TODO: is this allways correct? - numRows += numEq * subMatrixInFirstRow.M(); - }); - - // /*const std::size_t*/ numRows = this->model_().numCellCenterDofs()*numEqCellCenter + this->model_().numFaceDofs()*numEqFace; - - - // set the rowsizes - Dune::MatrixIndexSet occupationPattern; - occupationPattern.resize(numRows, numRows); - - // lambda function to fill the occupation pattern - auto addIndices = [&occupationPattern](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol) - { - using BlockType = typename std::decay_t::block_type; - const auto blockSizeI = BlockType::rows; - const auto blockSizeJ = BlockType::cols; - for(auto row = subMatrix.begin(); row != subMatrix.end(); ++row) - for(auto col = row->begin(); col != row->end(); ++col) - for(std::size_t i = 0; i < blockSizeI; ++i) - for(std::size_t j = 0; j < blockSizeJ; ++j) - occupationPattern.add(startRow + row.index()*blockSizeI + i, startCol + col.index()*blockSizeJ + j); - - }; - - - int rowIndex = 0; - - Dune::Hybrid::forEach(A, [&rowIndex, &addIndices, numRows](const auto& rowOfMultiTypeMatrix) - { - int colIndex = 0; - - std::cout << "processing one row of multitypeblockmatrix " << std::endl; - Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [&colIndex, &addIndices, &rowIndex, numRows](const auto& subMatrix) - { - std::cout << "numRows: " << subMatrix.N() << ", numCols: " << subMatrix.M() << std::endl; - addIndices(subMatrix, rowIndex, colIndex); - - const auto numEq = std::decay_t::block_type::rows; - colIndex += numEq*subMatrix.M(); - - // if we arrived at the right side of the matrix, increase the row index - if(colIndex == numRows) - rowIndex += numEq*subMatrix.N(); - }); - }); - - auto M = SparseMatrix(numRows, numRows, SparseMatrix::random); - occupationPattern.exportIdx(M); - - return M; - } - /*! * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$. * @@ -168,28 +98,23 @@ public: SolutionVector &x, SolutionVector &b) { - // createMatrix(A); try { if (this->numSteps_ == 0) this->initialResidual_ = b.two_norm(); - // copy the matrix and the vector to types the IterativeSolverBackend can handle - using MatrixBlock = typename Dune::FieldMatrix; - using SparseMatrix = typename Dune::BCRSMatrix; - - // get the new matrix sizes - std::size_t numRows = numEqCellCenter*A[cellCenterIdx][cellCenterIdx].N() + numEqFace*A[faceIdx][cellCenterIdx].N(); - std::size_t numCols = numEqCellCenter*A[cellCenterIdx][cellCenterIdx].M() + numEqFace*A[cellCenterIdx][faceIdx].M(); - // check matrix sizes assert(A[cellCenterIdx][cellCenterIdx].N() == A[cellCenterIdx][faceIdx].N()); assert(A[faceIdx][cellCenterIdx].N() == A[faceIdx][faceIdx].N()); - assert(numRows == numCols); // create the bcrs matrix the IterativeSolver backend can handle const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); + // get the new matrix sizes + const std::size_t numRows = M.N(); + const std::size_t numCols = M.M(); + assert(numRows == numCols); + // create the vector the IterativeSolver backend can handle using VectorBlock = typename Dune::FieldVector; using BlockVector = typename Dune::BlockVector; -- GitLab From 7e56d5d79fc6c4cb15604b0bc1bd64c127775acb Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Mon, 4 Sep 2017 14:45:33 +0200 Subject: [PATCH 06/10] [matrixconverter] Add method to convert multitype blockvectors to plain blockvectors * additionally add method to copy back values to multitype blockvector * clean-up * improve docu --- dumux/linear/matrixconverter.hh | 96 +++++++++++++++++++++++++++++---- 1 file changed, 87 insertions(+), 9 deletions(-) diff --git a/dumux/linear/matrixconverter.hh b/dumux/linear/matrixconverter.hh index 77828b0428..08ac279aa3 100644 --- a/dumux/linear/matrixconverter.hh +++ b/dumux/linear/matrixconverter.hh @@ -24,22 +24,20 @@ #ifndef DUMUX_MATRIX_CONVERTER #define DUMUX_MATRIX_CONVERTER -#include -#include #include +#include +#include #include #include #include -#include namespace Dumux { /*! - *A helper classe that converts a Dune::MultiTypeBlockMatrix into a plain Dune::BCRSMatrix + * A helper classe that converts a Dune::MultiTypeBlockMatrix into a plain Dune::BCRSMatrix * TODO: allow block sizes for BCRSMatrix other than 1x1 ? * */ - template class MatrixConverter { @@ -98,10 +96,10 @@ private: }; // fill the pattern - int rowIndex = 0; + std::size_t rowIndex = 0; Dune::Hybrid::forEach(A, [&addIndices, &rowIndex, &occupationPattern, numRows](const auto& rowOfMultiTypeMatrix) { - int colIndex = 0; + std::size_t colIndex = 0; Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [&addIndices, &occupationPattern, &colIndex, &rowIndex, numRows](const auto& subMatrix) { addIndices(subMatrix, rowIndex, colIndex); @@ -144,10 +142,10 @@ private: }; - int rowIndex = 0; + std::size_t rowIndex = 0; Dune::Hybrid::forEach(A, [©Values, &rowIndex, numRows](const auto& rowOfMultiTypeMatrix) { - int colIndex = 0; + std::size_t colIndex = 0; Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [©Values, &colIndex, &rowIndex, numRows](const auto& subMatrix) { copyValues(subMatrix, rowIndex, colIndex); @@ -183,6 +181,86 @@ private: } }; + +/*! + * A helper classe that converts a Dune::MultiTypeBlockVector into a plain Dune::BlockVector and transfers back values + */ +template +class VectorConverter +{ + using VectorBlock = typename Dune::FieldVector; + using BlockVector = typename Dune::BlockVector; + +public: + + /*! + * \brief Converts a Dune::MultiTypeBlockVector to a plain 1x1 Dune::BlockVector + * + * \param A The original multitype blockvector + */ + static auto multitypeToBlockVector(const MultiTypeBlockVector& b) + { + const auto size = getSize_(b); + + BlockVector bTmp; + bTmp.resize(size); + + std::size_t startIndex = 0; + Dune::Hybrid::forEach(b, [&bTmp, &startIndex](const auto& subVector) + { + const auto numEq = std::decay_t::block_type::size(); + + for(std::size_t i = 0; i < subVector.size(); ++i) + for(std::size_t j = 0; j < numEq; ++j) + bTmp[startIndex + i*numEq + j] = subVector[i][j]; + + startIndex += numEq*subVector.size(); + }); + + return bTmp; + } + + /*! + * \brief Copys the entries of a Dune::BlockVector to a Dune::MultiTypeBlockVector + * + * \param x The multitype blockvector where the values are copied to + * \param y The regular blockvector where the values are copied from + */ + static void retrieveValues(MultiTypeBlockVector& x, const BlockVector& y) + { + std::size_t startIndex = 0; + Dune::Hybrid::forEach(x, [&y, &startIndex](auto& subVector) + { + const auto numEq = std::decay_t::block_type::size(); + + for(std::size_t i = 0; i < subVector.size(); ++i) + for(std::size_t j = 0; j < numEq; ++j) + subVector[i][j] = y[startIndex + i*numEq + j]; + + startIndex += numEq*subVector.size(); + }); + } + +private: + + /*! + * \brief Returns the size of the expanded multitype block vector + * + * \param b The multitype blockvector + */ + static std::size_t getSize_(const MultiTypeBlockVector& b) + { + std::size_t size = 0; + Dune::Hybrid::forEach(b, [&size](const auto& subVector) + { + // the size of the individual vector blocks equals the respective number of equations. + const auto numEq = std::decay_t::block_type::size(); + size += numEq * subVector.size(); + }); + + return size; + } +}; } #endif -- GitLab From a10c1f60271d4419f10365b4904990925b0a4afb Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Mon, 4 Sep 2017 14:53:34 +0200 Subject: [PATCH 07/10] [staggered][newtoncontroller] Use converter and clean-up --- dumux/implicit/staggered/newtoncontroller.hh | 24 ++++++-------------- 1 file changed, 7 insertions(+), 17 deletions(-) diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh index e297bf30ae..a6ae4dff7b 100644 --- a/dumux/implicit/staggered/newtoncontroller.hh +++ b/dumux/implicit/staggered/newtoncontroller.hh @@ -112,22 +112,17 @@ public: // get the new matrix sizes const std::size_t numRows = M.N(); - const std::size_t numCols = M.M(); - assert(numRows == numCols); + assert(numRows == M.M()); // create the vector the IterativeSolver backend can handle + const auto bTmp = VectorConverter::multitypeToBlockVector(b); + assert(bTmp.size() == numRows); + + // create a blockvector to which the linear solver writes the solution using VectorBlock = typename Dune::FieldVector; using BlockVector = typename Dune::BlockVector; - - BlockVector y, bTmp; + BlockVector y; y.resize(numRows); - bTmp.resize(numCols); - for (std::size_t i = 0; i < b[cellCenterIdx].N(); ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - bTmp[i*numEqCellCenter + j] = b[cellCenterIdx][i][j]; - for (std::size_t i = 0; i < b[faceIdx].N(); ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - bTmp[i*numEqFace + j + b[cellCenterIdx].N()*numEqCellCenter] = b[faceIdx][i][j]; // printmatrix(std::cout, M, "", ""); @@ -135,12 +130,7 @@ public: const bool converged = this->linearSolver_.solve(M, y, bTmp); // copy back the result y into x - for (std::size_t i = 0; i < x[cellCenterIdx].N(); ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - x[cellCenterIdx][i][j] = y[i*numEqCellCenter + j]; - for (std::size_t i = 0; i < x[faceIdx].N(); ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - x[faceIdx][i][j] = y[i*numEqFace + j + x[cellCenterIdx].N()*numEqCellCenter]; + VectorConverter::retrieveValues(x, y); if (!converged) DUNE_THROW(NumericalProblem, "Linear solver did not converge"); -- GitLab From d3d1e585cdd4c8ca362168bd218e5aab23e5f35e Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Mon, 4 Sep 2017 17:42:42 +0200 Subject: [PATCH 08/10] [mixeddimension][newtoncontroller] Use convertmatrix class --- dumux/mixeddimension/newtoncontroller.hh | 103 ++++------------------- 1 file changed, 15 insertions(+), 88 deletions(-) diff --git a/dumux/mixeddimension/newtoncontroller.hh b/dumux/mixeddimension/newtoncontroller.hh index e938fa248f..ecd05cb4fb 100644 --- a/dumux/mixeddimension/newtoncontroller.hh +++ b/dumux/mixeddimension/newtoncontroller.hh @@ -31,6 +31,7 @@ #include #include #include +#include namespace Dumux { @@ -346,110 +347,36 @@ public: using MatrixBlock = typename Dune::FieldMatrix; using SparseMatrix = typename Dune::BCRSMatrix; - // get the new matrix sizes - std::size_t numRows = numEqBulk*A[bulkIdx][bulkIdx].N() + numEqLowDim*A[lowDimIdx][bulkIdx].N(); - std::size_t numCols = numEqBulk*A[bulkIdx][bulkIdx].M() + numEqLowDim*A[bulkIdx][lowDimIdx].M(); - // check matrix sizes assert(A[bulkIdx][bulkIdx].N() == A[bulkIdx][lowDimIdx].N()); assert(A[lowDimIdx][bulkIdx].N() == A[lowDimIdx][lowDimIdx].N()); - assert(numRows == numCols); // create the bcrs matrix the IterativeSolver backend can handle - auto M = SparseMatrix(numRows, numCols, SparseMatrix::random); - - // set the rowsizes - // A11 and A12 - for (auto row = A[bulkIdx][bulkIdx].begin(); row != A[bulkIdx][bulkIdx].end(); ++row) - for (std::size_t i = 0; i < numEqBulk; ++i) - M.setrowsize(numEqBulk*row.index() + i, row->size()*numEqBulk); - for (auto row = A[bulkIdx][lowDimIdx].begin(); row != A[bulkIdx][lowDimIdx].end(); ++row) - for (std::size_t i = 0; i < numEqBulk; ++i) - M.setrowsize(numEqBulk*row.index() + i, M.getrowsize(numEqBulk*row.index() + i) + row->size()*numEqLowDim); - // A21 and A22 - for (auto row = A[lowDimIdx][bulkIdx].begin(); row != A[lowDimIdx][bulkIdx].end(); ++row) - for (std::size_t i = 0; i < numEqLowDim; ++i) - M.setrowsize(numEqLowDim*row.index() + i + A[bulkIdx][bulkIdx].N()*numEqBulk, row->size()*numEqBulk); - for (auto row = A[lowDimIdx][lowDimIdx].begin(); row != A[lowDimIdx][lowDimIdx].end(); ++row) - for (std::size_t i = 0; i < numEqLowDim; ++i) - M.setrowsize(numEqLowDim*row.index() + i + A[bulkIdx][bulkIdx].N()*numEqBulk, M.getrowsize(numEqLowDim*row.index() + i + A[bulkIdx][bulkIdx].N()*numEqBulk) + row->size()*numEqLowDim); - M.endrowsizes(); - - // set the indices - for (auto row = A[bulkIdx][bulkIdx].begin(); row != A[bulkIdx][bulkIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqBulk; ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - M.addindex(row.index()*numEqBulk + i, col.index()*numEqBulk + j); - - for (auto row = A[bulkIdx][lowDimIdx].begin(); row != A[bulkIdx][lowDimIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqBulk; ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - M.addindex(row.index()*numEqBulk + i, col.index()*numEqLowDim + j + A[bulkIdx][bulkIdx].M()*numEqBulk); - - for (auto row = A[lowDimIdx][bulkIdx].begin(); row != A[lowDimIdx][bulkIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqLowDim; ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - M.addindex(row.index()*numEqLowDim + i + A[bulkIdx][bulkIdx].N()*numEqBulk, col.index()*numEqBulk + j); - - for (auto row = A[lowDimIdx][lowDimIdx].begin(); row != A[lowDimIdx][lowDimIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqLowDim; ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - M.addindex(row.index()*numEqLowDim + i + A[bulkIdx][bulkIdx].N()*numEqBulk, col.index()*numEqLowDim + j + A[bulkIdx][bulkIdx].M()*numEqBulk); - M.endindices(); - - // copy values - for (auto row = A[bulkIdx][bulkIdx].begin(); row != A[bulkIdx][bulkIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqBulk; ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - M[row.index()*numEqBulk + i][col.index()*numEqBulk + j] = A[bulkIdx][bulkIdx][row.index()][col.index()][i][j]; - - for (auto row = A[bulkIdx][lowDimIdx].begin(); row != A[bulkIdx][lowDimIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqBulk; ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - M[row.index()*numEqBulk + i][col.index()*numEqLowDim + j + A[bulkIdx][bulkIdx].M()*numEqBulk] = A[bulkIdx][lowDimIdx][row.index()][col.index()][i][j]; - - for (auto row = A[lowDimIdx][bulkIdx].begin(); row != A[lowDimIdx][bulkIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqLowDim; ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - M[row.index()*numEqLowDim + i + A[bulkIdx][bulkIdx].N()*numEqBulk][col.index()*numEqBulk + j] = A[lowDimIdx][bulkIdx][row.index()][col.index()][i][j]; - - for (auto row = A[lowDimIdx][lowDimIdx].begin(); row != A[lowDimIdx][lowDimIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqLowDim; ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - M[row.index()*numEqLowDim + i + A[bulkIdx][bulkIdx].N()*numEqBulk][col.index()*numEqLowDim + j + A[bulkIdx][bulkIdx].M()*numEqBulk] = A[lowDimIdx][lowDimIdx][row.index()][col.index()][i][j]; + const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); + + // get the new matrix sizes + const std::size_t numRows = M.N(); + assert(numRows == M.M()); + + // create the vector the IterativeSolver backend can handle + const auto bTmp = VectorConverter::multitypeToBlockVector(b); + assert(bTmp.size() == numRows); // create the vector the IterativeSolver backend can handle using VectorBlock = typename Dune::FieldVector; using BlockVector = typename Dune::BlockVector; - BlockVector y, bTmp; + // create a blockvector to which the linear solver writes the solution + using VectorBlock = typename Dune::FieldVector; + using BlockVector = typename Dune::BlockVector; + BlockVector y; y.resize(numRows); - bTmp.resize(numCols); - for (std::size_t i = 0; i < b[bulkIdx].N(); ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - bTmp[i*numEqBulk + j] = b[bulkIdx][i][j]; - for (std::size_t i = 0; i < b[lowDimIdx].N(); ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - bTmp[i*numEqLowDim + j + b[bulkIdx].N()*numEqBulk] = b[lowDimIdx][i][j]; // solve bool converged = linearSolver_.solve(M, y, bTmp); // copy back the result y into x - for (std::size_t i = 0; i < x[bulkIdx].N(); ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - x[bulkIdx][i][j] = y[i*numEqBulk + j]; - for (std::size_t i = 0; i < x[lowDimIdx].N(); ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - x[lowDimIdx][i][j] = y[i*numEqLowDim + j + x[bulkIdx].N()*numEqBulk]; + VectorConverter::retrieveValues(x, y); if (!converged) DUNE_THROW(NumericalProblem, "Linear solver did not converge"); -- GitLab From 0c838b596d48a88d783a1d4f03ca1221f160af66 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Mon, 4 Sep 2017 17:48:34 +0200 Subject: [PATCH 09/10] [newtoncontroller] Use known type * get rid of decltype --- dumux/implicit/staggered/newtoncontroller.hh | 4 ++-- dumux/mixeddimension/newtoncontroller.hh | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh index a6ae4dff7b..e44e17b363 100644 --- a/dumux/implicit/staggered/newtoncontroller.hh +++ b/dumux/implicit/staggered/newtoncontroller.hh @@ -115,7 +115,7 @@ public: assert(numRows == M.M()); // create the vector the IterativeSolver backend can handle - const auto bTmp = VectorConverter::multitypeToBlockVector(b); + const auto bTmp = VectorConverter::multitypeToBlockVector(b); assert(bTmp.size() == numRows); // create a blockvector to which the linear solver writes the solution @@ -130,7 +130,7 @@ public: const bool converged = this->linearSolver_.solve(M, y, bTmp); // copy back the result y into x - VectorConverter::retrieveValues(x, y); + VectorConverter::retrieveValues(x, y); if (!converged) DUNE_THROW(NumericalProblem, "Linear solver did not converge"); diff --git a/dumux/mixeddimension/newtoncontroller.hh b/dumux/mixeddimension/newtoncontroller.hh index ecd05cb4fb..65aaba1fff 100644 --- a/dumux/mixeddimension/newtoncontroller.hh +++ b/dumux/mixeddimension/newtoncontroller.hh @@ -359,7 +359,7 @@ public: assert(numRows == M.M()); // create the vector the IterativeSolver backend can handle - const auto bTmp = VectorConverter::multitypeToBlockVector(b); + const auto bTmp = VectorConverter::multitypeToBlockVector(b); assert(bTmp.size() == numRows); // create the vector the IterativeSolver backend can handle @@ -376,7 +376,7 @@ public: bool converged = linearSolver_.solve(M, y, bTmp); // copy back the result y into x - VectorConverter::retrieveValues(x, y); + VectorConverter::retrieveValues(x, y); if (!converged) DUNE_THROW(NumericalProblem, "Linear solver did not converge"); -- GitLab From e70552db1256504b85227ce56e5c9877728d781c Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Mon, 4 Sep 2017 17:51:32 +0200 Subject: [PATCH 10/10] [matrixconverter] Rename multitypeToBlockVector -> multiTypeToBlockVector --- dumux/implicit/staggered/newtoncontroller.hh | 2 +- dumux/linear/matrixconverter.hh | 2 +- dumux/mixeddimension/newtoncontroller.hh | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh index e44e17b363..1e4cf1dd4a 100644 --- a/dumux/implicit/staggered/newtoncontroller.hh +++ b/dumux/implicit/staggered/newtoncontroller.hh @@ -115,7 +115,7 @@ public: assert(numRows == M.M()); // create the vector the IterativeSolver backend can handle - const auto bTmp = VectorConverter::multitypeToBlockVector(b); + const auto bTmp = VectorConverter::multiTypeToBlockVector(b); assert(bTmp.size() == numRows); // create a blockvector to which the linear solver writes the solution diff --git a/dumux/linear/matrixconverter.hh b/dumux/linear/matrixconverter.hh index 08ac279aa3..0615754cd0 100644 --- a/dumux/linear/matrixconverter.hh +++ b/dumux/linear/matrixconverter.hh @@ -198,7 +198,7 @@ public: * * \param A The original multitype blockvector */ - static auto multitypeToBlockVector(const MultiTypeBlockVector& b) + static auto multiTypeToBlockVector(const MultiTypeBlockVector& b) { const auto size = getSize_(b); diff --git a/dumux/mixeddimension/newtoncontroller.hh b/dumux/mixeddimension/newtoncontroller.hh index 65aaba1fff..a95a7bd1f3 100644 --- a/dumux/mixeddimension/newtoncontroller.hh +++ b/dumux/mixeddimension/newtoncontroller.hh @@ -359,7 +359,7 @@ public: assert(numRows == M.M()); // create the vector the IterativeSolver backend can handle - const auto bTmp = VectorConverter::multitypeToBlockVector(b); + const auto bTmp = VectorConverter::multiTypeToBlockVector(b); assert(bTmp.size() == numRows); // create the vector the IterativeSolver backend can handle -- GitLab