From 2e76b22626f8a8c95bacce6d3605d96fc8738f9c Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Fri, 17 Feb 2017 07:56:44 +0100
Subject: [PATCH] [staggeredGrid][linear] Remove directSolverBackend class
* Convert matrix (if required) in newtoncontroller, as done in
mixeddimension
---
dumux/implicit/staggered/newtoncontroller.hh | 212 +++++++++++++++----
dumux/implicit/staggered/propertydefaults.hh | 6 +-
dumux/linear/directsolverbackend.hh | 200 -----------------
3 files changed, 173 insertions(+), 245 deletions(-)
delete mode 100644 dumux/linear/directsolverbackend.hh
diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh
index 3425808ecb..2d5fdc22d2 100644
--- a/dumux/implicit/staggered/newtoncontroller.hh
+++ b/dumux/implicit/staggered/newtoncontroller.hh
@@ -30,9 +30,21 @@
#include "properties.hh"
#include <dumux/nonlinear/newtoncontroller.hh>
+#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
#include "newtonconvergencewriter.hh"
namespace Dumux {
+
+namespace Properties
+{
+ SET_PROP(StaggeredModel, LinearSolverBlockSize)
+ {
+ // LinearSolverAcceptsMultiTypeMatrix<T>::value
+ // TODO: make somehow dependend? or only relevant for direct solvers?
+ public:
+ static constexpr auto value = 1;
+ };
+}
/*!
* \ingroup PNMModel
* \brief A PNM specific controller for the newton solver.
@@ -55,77 +67,195 @@ class StaggeredNewtonController : public NewtonController<TypeTag>
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
+ enum {
+ numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter),
+ numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace)
+ };
+
public:
StaggeredNewtonController(const Problem &problem)
: ParentType(problem)
{}
- /*!
+ /*!
* \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
*
* Throws Dumux::NumericalProblem if the linear solver didn't
* converge.
*
+ * If the linear solver doesn't accept multitype matrices we copy the matrix
+ * into a 1x1 block BCRS matrix for solving.
+ *
* \param A The matrix of the linear system of equations
* \param x The vector which solves the linear system
* \param b The right hand side of the linear system
*/
- void newtonSolveLinear(JacobianMatrix &A,
- SolutionVector &x,
- SolutionVector &b)
+ template<typename T = TypeTag>
+ typename std::enable_if<!LinearSolverAcceptsMultiTypeMatrix<T>::value, void>::type
+ newtonSolveLinear(JacobianMatrix &A,
+ SolutionVector &x,
+ SolutionVector &b)
{
- try {
+ try
+ {
if (this->numSteps_ == 0)
- {
- Scalar norm2 = b.two_norm2();
- if (this->gridView_().comm().size() > 1)
- norm2 = this->gridView_().comm().sum(norm2);
+ this->initialResidual_ = b.two_norm();
- this->initialResidual_ = std::sqrt(norm2);
- }
+ // copy the matrix and the vector to types the IterativeSolverBackend can handle
+ using MatrixBlock = typename Dune::FieldMatrix<Scalar, 1, 1>;
+ using SparseMatrix = typename Dune::BCRSMatrix<MatrixBlock>;
- int converged = this->linearSolver_.solve(A, x, b);
+ // 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();
- // make sure all processes converged
- int convergedRemote = converged;
- if (this->gridView_().comm().size() > 1)
- convergedRemote = this->gridView_().comm().min(converged);
+ // check matrix sizes
+ assert(A[cellCenterIdx][cellCenterIdx].N() == A[cellCenterIdx][faceIdx].N());
+ assert(A[faceIdx][cellCenterIdx].N() == A[faceIdx][faceIdx].N());
+ assert(numRows == numCols);
- if (!converged) {
- DUNE_THROW(NumericalProblem,
- "Linear solver did not converge");
- }
- else if (!convergedRemote) {
- DUNE_THROW(NumericalProblem,
- "Linear solver did not converge on a remote process");
- }
- }
- catch (Dune::MatrixBlockError e) {
- // make sure all processes converged
- int converged = 0;
- if (this->gridView_().comm().size() > 1)
- converged = this->gridView_().comm().min(converged);
+ // 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];
+
+ // create the vector the IterativeSolver backend can handle
+ using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
+ using BlockVector = typename Dune::BlockVector<VectorBlock>;
+
+ BlockVector y, bTmp;
+ 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];
+
+ // solve
+ 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];
+
+ if (!converged)
+ DUNE_THROW(NumericalProblem, "Linear solver did not converge");
+ }
+ catch (const Dune::Exception &e)
+ {
Dumux::NumericalProblem p;
- std::string msg;
- std::ostringstream ms(msg);
-// ms << e.what() << "M=" << A[e.r][e.c];
-// p.message(ms.str());
-// throw p;
+ p.message(e.what());
+ throw p;
}
- catch (const Dune::Exception &e) {
- // make sure all processes converged
- int converged = 0;
- if (this->gridView_().comm().size() > 1)
- converged = this->gridView_().comm().min(converged);
+ }
+ /*!
+ * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
+ *
+ * Throws Dumux::NumericalProblem if the linear solver didn't
+ * converge.
+ *
+ * \param A The matrix of the linear system of equations
+ * \param x The vector which solves the linear system
+ * \param b The right hand side of the linear system
+ */
+ template<typename T = TypeTag>
+ typename std::enable_if<LinearSolverAcceptsMultiTypeMatrix<T>::value, void>::type
+ newtonSolveLinear(JacobianMatrix &A,
+ SolutionVector &x,
+ SolutionVector &b)
+ {
+ try
+ {
+ if (this->numSteps_ == 0)
+ this->initialResidual_ = b.two_norm();
+
+ bool converged = this->linearSolver_.solve(A, x, b);
+
+ if (!converged)
+ DUNE_THROW(NumericalProblem, "Linear solver did not converge");
+ }
+ catch (const Dune::Exception &e)
+ {
Dumux::NumericalProblem p;
p.message(e.what());
throw p;
}
}
- /*!
+ /*!
* \brief Update the current solution with a delta vector.
*
* The error estimates required for the newtonConverged() and
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
index 57dca34ef3..8b6caf0818 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/dumux/implicit/staggered/propertydefaults.hh
@@ -47,9 +47,7 @@
#include <dune/istl/multitypeblockvector.hh>
#include <dune/istl/multitypeblockmatrix.hh>
-#include <dumux/linear/directsolverbackend.hh>
-
-
+#include <dumux/linear/seqsolverbackend.hh>
#include "assembler.hh"
#include "localresidual.hh"
@@ -213,7 +211,7 @@ SET_PROP(StaggeredModel, DofTypeIndices)
//! set default solver
// SET_TYPE_PROP(StaggeredModel, LinearSolver, Dumux::GSBiCGSTABBackend<TypeTag>);
-SET_TYPE_PROP(StaggeredModel, LinearSolver, Dumux::StaggeredGridUMFPackBackend<TypeTag>);
+SET_TYPE_PROP(StaggeredModel, LinearSolver, Dumux::UMFPackBackend<TypeTag>);
//! set the block level to 2, suitable for e.g. the Dune::MultiTypeBlockMatrix
SET_INT_PROP(StaggeredModel, LinearSolverPreconditionerBlockLevel, 2);
diff --git a/dumux/linear/directsolverbackend.hh b/dumux/linear/directsolverbackend.hh
deleted file mode 100644
index b5b2f90b7a..0000000000
--- a/dumux/linear/directsolverbackend.hh
+++ /dev/null
@@ -1,200 +0,0 @@
-// -*- 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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Dumux multidimension direct solver backend
- */
-#ifndef DUMUX_STAGGEREDGRID_DIRECT_SOLVER_BACKEND_HH
-#define DUMUX_STAGGEREDGRID_DIRECT_SOLVER_BACKEND_HH
-
-//#include <dune/istl/superlu.hh>
-#include <dune/istl/umfpack.hh>
-
-#include <dumux/common/parameters.hh>
-#include <dumux/common/basicproperties.hh>
-#include <dumux/linear/linearsolverproperties.hh>
-
-namespace Dumux {
-
-#if HAVE_UMFPACK
-/*! \brief UMFPackBackend for MultiTypeMatrices
- * Copies the coupled matrix into a single BCRSMatrix.
- * Very slow!! Only wise to use for verification purposes.
- */
-template <class TypeTag>
-class StaggeredGridUMFPackBackend
-{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-
- enum {
- numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter),
- numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace)
- };
-
- typename GET_PROP(TypeTag, DofTypeIndices)::CellCenterIdx cellCenterIdx;
- typename GET_PROP(TypeTag, DofTypeIndices)::FaceIdx faceIdx;
-
-public:
-
- StaggeredGridUMFPackBackend(const Problem& problem) {}
-
- template<class Matrix, class Vector>
- bool solve(const Matrix& A, Vector& x, const Vector& b)
- {
- // copy the matrix and the vector to types the UMFPackBackend can handle
- using MatrixBlock = typename Dune::FieldMatrix<Scalar, 1, 1>;
- using SparseMatrix = typename Dune::BCRSMatrix<MatrixBlock>;
-
- // 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 UMFPack 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];
-
- // create the vector the UMFPack backend can handle
- using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
- using BlockVector = typename Dune::BlockVector<VectorBlock>;
-
- BlockVector y, bTmp;
- 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];
-
- int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
- Dune::UMFPack<SparseMatrix> solver(M, verbosity > 0);
-
- solver.apply(y, bTmp, result_);
-
- std::size_t size = y.size();
- for (std::size_t i = 0; i < size; i++)
- {
- if (std::isnan(y[i][0]) || std::isinf(y[i][0]))
- {
- result_.converged = false;
- break;
- }
- }
-
- // 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];
-
-// printmatrix(std::cout,M, "umfpack", "");
-
- return result_.converged;
- }
-
- const Dune::InverseOperatorResult& result() const
- {
- return result_;
- }
-
-private:
- Dune::InverseOperatorResult result_;
-};
-
-#endif // HAVE_UMFPACK
-
-} // end namespace Dumux
-
-#endif
--
GitLab