From d3cb7f6e6bedb8886e2a36df34dad2c1206b352c Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Thu, 12 May 2011 15:54:31 +0000
Subject: [PATCH] generalized BoxBiCGSTABILU0Solver to BoxLinearSolver, built
two derived linear solvers BoxBiCGSTABILU0Solver (default for BoxModel) and
BoxCGILU0Solver (for OnePModel). This also solves the convergence issue for
test_1p.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5815 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/1p/1pproperties.hh | 3 +
.../richards/richardspropertydefaults.hh | 2 +-
dumux/linear/boxbicgstabilu0solver.hh | 158 +-----------
dumux/linear/boxlinearsolver.hh | 241 ++++++++++++++++++
dumux/nonlinear/newtoncontroller.hh | 3 +-
test/boxmodels/1p/grids/test_1p_2d.dgf | 2 +-
6 files changed, 250 insertions(+), 159 deletions(-)
create mode 100644 dumux/linear/boxlinearsolver.hh
diff --git a/dumux/boxmodels/1p/1pproperties.hh b/dumux/boxmodels/1p/1pproperties.hh
index be1a028c8a..b934310b19 100644
--- a/dumux/boxmodels/1p/1pproperties.hh
+++ b/dumux/boxmodels/1p/1pproperties.hh
@@ -58,6 +58,9 @@ NEW_PROP_TAG(Fluid); //!< The fluid for the single-phase problems
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(UpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes
+// use the CG solver preconditioned by ILU-0
+SET_TYPE_PROP(BoxOneP, LinearSolver, Dumux::BoxCGILU0Solver<TypeTag> );
+
// \}
};
diff --git a/dumux/boxmodels/richards/richardspropertydefaults.hh b/dumux/boxmodels/richards/richardspropertydefaults.hh
index aa8f551bb0..e0f17b7554 100644
--- a/dumux/boxmodels/richards/richardspropertydefaults.hh
+++ b/dumux/boxmodels/richards/richardspropertydefaults.hh
@@ -1,4 +1,4 @@
-// $Id: richardsproperties.hh 3840 2010-07-15 10:14:15Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/linear/boxbicgstabilu0solver.hh b/dumux/linear/boxbicgstabilu0solver.hh
index 13c17cbd68..7b77c8fe43 100644
--- a/dumux/linear/boxbicgstabilu0solver.hh
+++ b/dumux/linear/boxbicgstabilu0solver.hh
@@ -26,162 +26,8 @@
#ifndef DUMUX_BICGSTAB_ILU0_SOLVER_HH
#define DUMUX_BICGSTAB_ILU0_SOLVER_HH
-#include <dumux/common/propertysystem.hh>
-#include <dumux/linear/vertexborderlistfromgrid.hh>
-#include <dumux/linear/overlappingbcrsmatrix.hh>
-#include <dumux/linear/overlappingblockvector.hh>
-#include <dumux/linear/overlappingpreconditioner.hh>
-#include <dumux/linear/overlappingscalarproduct.hh>
-#include <dumux/linear/overlappingoperator.hh>
+#warning inclusion of this file is DEPRECATED, please include boxlinearsolver.hh instead
-#include <dune/istl/solvers.hh>
-#include <dune/istl/preconditioners.hh>
-
-namespace Dumux {
-
-/*!
- * \brief Provides a linear solver for the stabilized BiCG method with
- * an ILU-0 preconditioner.
- *
- * This solver's intention is to be used in conjunction with the box
- * method, so it assumes that the vertices are the only DOFs.
- */
-template <class TypeTag>
-class BoxBiCGStabILU0Solver
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) Matrix;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) Vector;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
-
- typedef Dumux::OverlappingBCRSMatrix<Matrix> OverlappingMatrix;
- typedef typename OverlappingMatrix::Overlap Overlap;
- typedef Dumux::OverlappingBlockVector<typename Vector::block_type, Overlap> OverlappingVector;
- typedef Dune::SeqILU0<OverlappingMatrix, OverlappingVector, OverlappingVector> SeqPreconditioner;
- //typedef Dune::SeqJac<OverlappingMatrix, OverlappingVector, OverlappingVector> SeqPreconditioner;
- typedef Dumux::OverlappingPreconditioner<SeqPreconditioner, Overlap> OverlappingPreconditioner;
- typedef Dumux::OverlappingScalarProduct<OverlappingVector, Overlap> OverlappingScalarProduct;
- typedef Dumux::OverlappingOperator<OverlappingMatrix, OverlappingVector, OverlappingVector> OverlappingOperator;
- typedef Dune::BiCGSTABSolver<OverlappingVector> Solver;
-
-public:
- BoxBiCGStabILU0Solver(const Problem &problem, int overlapSize=3)
- : problem_(problem)
- , overlapSize_(overlapSize)
- {
- overlapMatrix_ = 0;
- overlapb_ = 0;
- overlapx_ = 0;
- };
-
- ~BoxBiCGStabILU0Solver()
- { cleanup_(); }
-
- /*!
- * \brief Set the structure of the linear system of equations to be solved.
- *
- * This method allocates space an does the necessary
- * communication before actually calling the solve() method. As
- * long as the structure of the linear system does not change, the
- * solve method can be called arbitrarily often.
- */
- void setStructureMatrix(const Matrix &M)
- {
- cleanup_();
- prepare_();
- };
-
- /*!
- * \brief Actually solve the linear system of equations.
- *
- * \return true if the residual reduction could be achieved, else false.
- */
- bool solve(const Matrix &M,
- Vector &x,
- const Vector &b,
- double residReduction,
- int verbosityLevel = 0)
- {
- if (!overlapMatrix_) {
- // make sure that the overlapping matrix and block vectors
- // have been created
- prepare_(M);
- };
-
- // copy the values of the non-overlapping linear system of
- // equations to the overlapping one. On ther border, we add up
- // the values of all processes (using the assignAdd() methods)
- overlapMatrix_->assignAdd(M);
- overlapb_->assignAdd(b);
- (*overlapx_) = 0.0;
-
- // create sequential and overlapping preconditioners
- //SeqPreconditioner seqPreCond(*overlapMatrix_, 1, 1.0);
- SeqPreconditioner seqPreCond(*overlapMatrix_, 1.0);
- OverlappingPreconditioner preCond(seqPreCond, overlapMatrix_->overlap());
-
- // create the scalar products and linear operators for ISTL
- OverlappingScalarProduct scalarProd(overlapMatrix_->overlap());
- OverlappingOperator opA(*overlapMatrix_);
-
- // create the actual solver
- Solver solver(opA,
- scalarProd,
- preCond,
- residReduction,
- /*maxIterations=*/250,
- verbosityLevel);
-
- // run the solver
- Dune::InverseOperatorResult result;
- solver.apply(*overlapx_, *overlapb_, result);
-
- // copy the result back to the non-overlapping vector
- overlapx_->assignTo(x);
-
- // return the result of the solver
- return result.converged;
- };
-
-private:
- void prepare_(const Matrix &M)
- {
- VertexBorderListFromGrid<GridView, VertexMapper>
- borderListCreator(problem_.gridView(), problem_.vertexMapper());
-
- // create the overlapping Jacobian matrix
- overlapMatrix_ = new OverlappingMatrix (M,
- borderListCreator.foreignBorderList(),
- borderListCreator.domesticBorderList(),
- overlapSize_);
-
- // create the overlapping vectors for the residual and the
- // solution
- overlapb_ = new OverlappingVector(overlapMatrix_->overlap());
- overlapx_ = new OverlappingVector(*overlapb_);
- };
-
- void cleanup_()
- {
- // create the overlapping Jacobian matrix and vectors
- delete overlapMatrix_;
- delete overlapb_;
- delete overlapx_;
-
- overlapMatrix_ = 0;
- overlapb_ = 0;
- overlapx_ = 0;
- };
-
- const Problem &problem_;
-
- int overlapSize_;
- OverlappingMatrix *overlapMatrix_;
- OverlappingVector *overlapb_;
- OverlappingVector *overlapx_;
-};
-
-} // namespace Dumux
+#include <dumux/linear/boxlinearsolver.hh>
#endif
diff --git a/dumux/linear/boxlinearsolver.hh b/dumux/linear/boxlinearsolver.hh
new file mode 100644
index 0000000000..ebb91fe13e
--- /dev/null
+++ b/dumux/linear/boxlinearsolver.hh
@@ -0,0 +1,241 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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 Provides a linear solver for the stabilized BiCG method with
+ * an ILU-0 preconditioner.
+ */
+#ifndef DUMUX_BOXLINEARSOLVER_HH
+#define DUMUX_BOXLINEARSOLVER_HH
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/linear/vertexborderlistfromgrid.hh>
+#include <dumux/linear/overlappingbcrsmatrix.hh>
+#include <dumux/linear/overlappingblockvector.hh>
+#include <dumux/linear/overlappingpreconditioner.hh>
+#include <dumux/linear/overlappingscalarproduct.hh>
+#include <dumux/linear/overlappingoperator.hh>
+
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+
+namespace Dumux {
+
+/*!
+ * \brief Provides a linear solver for the stabilized BiCG method with
+ * an ILU-0 preconditioner.
+ *
+ * This solver's intention is to be used in conjunction with the box
+ * method, so it assumes that the vertices are the only DOFs.
+ */
+template <class TypeTag>
+class BoxLinearSolver
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) Vector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef Dumux::OverlappingBCRSMatrix<Matrix> OverlappingMatrix;
+ typedef typename OverlappingMatrix::Overlap Overlap;
+ typedef Dumux::OverlappingBlockVector<typename Vector::block_type, Overlap> OverlappingVector;
+ typedef Dumux::OverlappingScalarProduct<OverlappingVector, Overlap> OverlappingScalarProduct;
+ typedef Dumux::OverlappingOperator<OverlappingMatrix, OverlappingVector, OverlappingVector> OverlappingOperator;
+
+public:
+ BoxLinearSolver(const Problem &problem, int overlapSize)
+ : problem_(problem)
+ , overlapSize_(overlapSize)
+ {
+ overlapMatrix_ = 0;
+ overlapb_ = 0;
+ overlapx_ = 0;
+ };
+
+ ~BoxLinearSolver()
+ { cleanup_(); }
+
+ /*!
+ * \brief Set the structure of the linear system of equations to be solved.
+ *
+ * This method allocates space an does the necessary
+ * communication before actually calling the solve() method. As
+ * long as the structure of the linear system does not change, the
+ * solve method can be called arbitrarily often.
+ */
+ void setStructureMatrix(const Matrix &M)
+ {
+ cleanup_();
+ prepare_();
+ };
+
+ /*!
+ * \brief Actually solve the linear system of equations.
+ *
+ * \return true if the residual reduction could be achieved, else false.
+ */
+ template <class SeqPreconditioner, class Solver>
+ bool solve(const Matrix &M,
+ Vector &x,
+ const Vector &b,
+ double residReduction,
+ int verbosityLevel = 0)
+ {
+ if (!overlapMatrix_) {
+ // make sure that the overlapping matrix and block vectors
+ // have been created
+ prepare_(M);
+ };
+
+ // copy the values of the non-overlapping linear system of
+ // equations to the overlapping one. On ther border, we add up
+ // the values of all processes (using the assignAdd() methods)
+ overlapMatrix_->assignAdd(M);
+ overlapb_->assignAdd(b);
+ (*overlapx_) = 0.0;
+
+ // create sequential and overlapping preconditioners
+ //SeqPreconditioner seqPreCond(*overlapMatrix_, 1, 1.0);
+ SeqPreconditioner seqPreCond(*overlapMatrix_, 1.0);
+ typedef Dumux::OverlappingPreconditioner<SeqPreconditioner, Overlap> OverlappingPreconditioner;
+ OverlappingPreconditioner preCond(seqPreCond, overlapMatrix_->overlap());
+
+ // create the scalar products and linear operators for ISTL
+ OverlappingScalarProduct scalarProd(overlapMatrix_->overlap());
+ OverlappingOperator opA(*overlapMatrix_);
+
+ // create the actual solver
+ Solver solver(opA,
+ scalarProd,
+ preCond,
+ residReduction,
+ /*maxIterations=*/250,
+ verbosityLevel);
+
+ // run the solver
+ Dune::InverseOperatorResult result;
+ solver.apply(*overlapx_, *overlapb_, result);
+
+ // copy the result back to the non-overlapping vector
+ overlapx_->assignTo(x);
+
+ // return the result of the solver
+ return result.converged;
+ };
+
+private:
+ void prepare_(const Matrix &M)
+ {
+ VertexBorderListFromGrid<GridView, VertexMapper>
+ borderListCreator(problem_.gridView(), problem_.vertexMapper());
+
+ // create the overlapping Jacobian matrix
+ overlapMatrix_ = new OverlappingMatrix (M,
+ borderListCreator.foreignBorderList(),
+ borderListCreator.domesticBorderList(),
+ overlapSize_);
+
+ // create the overlapping vectors for the residual and the
+ // solution
+ overlapb_ = new OverlappingVector(overlapMatrix_->overlap());
+ overlapx_ = new OverlappingVector(*overlapb_);
+ };
+
+ void cleanup_()
+ {
+ // create the overlapping Jacobian matrix and vectors
+ delete overlapMatrix_;
+ delete overlapb_;
+ delete overlapx_;
+
+ overlapMatrix_ = 0;
+ overlapb_ = 0;
+ overlapx_ = 0;
+ };
+
+ const Problem &problem_;
+
+ int overlapSize_;
+ OverlappingMatrix *overlapMatrix_;
+ OverlappingVector *overlapb_;
+ OverlappingVector *overlapx_;
+};
+
+template <class TypeTag>
+class BoxBiCGStabILU0Solver : public BoxLinearSolver<TypeTag>
+{
+ typedef BoxLinearSolver<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) Matrix;
+ typedef Dumux::OverlappingBCRSMatrix<Matrix> OverlappingMatrix;
+ typedef typename OverlappingMatrix::Overlap Overlap;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) Vector;
+ typedef Dumux::OverlappingBlockVector<typename Vector::block_type, Overlap> OverlappingVector;
+ typedef Dune::SeqILU0<OverlappingMatrix, OverlappingVector, OverlappingVector> SeqPreconditioner;
+ typedef Dune::BiCGSTABSolver<OverlappingVector> Solver;
+
+public:
+ template <class Problem>
+ BoxBiCGStabILU0Solver(const Problem &problem, int overlapSize = 3)
+ : ParentType(problem, overlapSize)
+ {}
+
+ bool solve(const Matrix &M,
+ Vector &x,
+ const Vector &b,
+ double residReduction,
+ int verbosityLevel = 0)
+ {
+ return ParentType::template solve<SeqPreconditioner, Solver>(M, x, b, residReduction, verbosityLevel);
+ }
+};
+
+template <class TypeTag>
+class BoxCGILU0Solver : public BoxLinearSolver<TypeTag>
+{
+ typedef BoxLinearSolver<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) Matrix;
+ typedef Dumux::OverlappingBCRSMatrix<Matrix> OverlappingMatrix;
+ typedef typename OverlappingMatrix::Overlap Overlap;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) Vector;
+ typedef Dumux::OverlappingBlockVector<typename Vector::block_type, Overlap> OverlappingVector;
+ typedef Dune::SeqILU0<OverlappingMatrix, OverlappingVector, OverlappingVector> SeqPreconditioner;
+ typedef Dune::CGSolver<OverlappingVector> Solver;
+
+public:
+ template <class Problem>
+ BoxCGILU0Solver(const Problem &problem, int overlapSize = 3)
+ : ParentType(problem, overlapSize)
+ {}
+
+ bool solve(const Matrix &M,
+ Vector &x,
+ const Vector &b,
+ double residReduction,
+ int verbosityLevel = 0)
+ {
+ return ParentType::template solve<SeqPreconditioner, Solver>(M, x, b, residReduction, verbosityLevel);
+ }
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index 8d5a8bfe1f..b8aba21e4e 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -36,7 +36,7 @@
#include <dumux/io/vtkmultiwriter.hh>
-#include <dumux/linear/boxbicgstabilu0solver.hh>
+#include <dumux/linear/boxlinearsolver.hh>
namespace Dumux
{
@@ -455,6 +455,7 @@ public:
verbosityLevel = GET_PROP_VALUE(TypeTag, PTAG(NewtonLinearSolverVerbosity));
}
+ Dune::writeMatrixToMatlab(A, "matrix.txt");
int converged = linearSolver_.solve(A, x, b, residReduction, verbosityLevel);
// make sure all processes converged
diff --git a/test/boxmodels/1p/grids/test_1p_2d.dgf b/test/boxmodels/1p/grids/test_1p_2d.dgf
index d92fa32fcc..caa0ae9193 100644
--- a/test/boxmodels/1p/grids/test_1p_2d.dgf
+++ b/test/boxmodels/1p/grids/test_1p_2d.dgf
@@ -2,7 +2,7 @@ DGF
Interval
0 0 % first corner
1 1 % second corner
-12 12 % 1 cells in x and 1 in y direction
+10 10 % cells in x and y direction
#
BOUNDARYDOMAIN
--
GitLab