amgbackend.hh 8.79 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
3
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
6
7
8
9
10
11
12
 *                                                                           *
 *   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          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
18
19
20
 *   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
Thomas Fetzer's avatar
Thomas Fetzer committed
21
 * \ingroup Linear
22
 * \brief Provides a parallel linear solver based on the ISTL AMG preconditioner
23
 *        and the ISTL BiCGSTAB solver.
24
 */
25
26
#ifndef DUMUX_PARALLEL_AMGBACKEND_HH
#define DUMUX_PARALLEL_AMGBACKEND_HH
27

28
29
30
#include <memory>

#include <dune/common/exceptions.hh>
31
32
#include <dune/common/parallel/indexset.hh>
#include <dune/common/parallel/mpicollectivecommunication.hh>
33
#include <dune/grid/common/capabilities.hh>
34
35
36
37
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/solvers.hh>

38
#include <dumux/linear/solver.hh>
39
#include <dumux/linear/amgparallelhelpers.hh>
40

41
42
namespace Dumux {

43
/*!
44
 * \ingroup Linear
45
 * \brief Scale the linear system by the inverse of
46
 * its (block-)diagonal entries.
47
 *
48
49
50
 * \param matrix the matrix to scale
 * \param rhs the right hand side vector to scale
 */
51
52
53
54
55
56
template <class Matrix, class Vector>
void scaleLinearSystem(Matrix& matrix, Vector& rhs)
{
    typename Matrix::RowIterator row = matrix.begin();
    for(; row != matrix.end(); ++row)
    {
57
        using size_type = typename Matrix::size_type;
58
59
        size_type rowIdx = row.index();

60
        using MatrixBlock = typename Matrix::block_type;
61
62
63
        MatrixBlock diagonal = matrix[rowIdx][rowIdx];
        diagonal.invert();

64
        using VectorBlock = typename Vector::block_type;
65
66
67
68
69
70
71
72
73
        const VectorBlock b = rhs[rowIdx];
        diagonal.mv(b, rhs[rowIdx]);

        typename Matrix::ColIterator col = row->begin();
        for (; col != row->end(); ++col)
            col->leftmultiply(diagonal);
    }
}

74
/*!
75
76
77
 * \ingroup Linear
 * \brief A linear solver based on the ISTL AMG preconditioner
 *        and the ISTL BiCGSTAB solver.
78
 */
79
80
template <class GridView, class AmgTraits>
class ParallelAMGBackend : public LinearSolver
81
{
82
    using Grid = typename GridView::Grid;
83
84
85
86
87
88
89
90
    using LinearOperator = typename AmgTraits::LinearOperator;
    using ScalarProduct = typename AmgTraits::ScalarProduct;
    using VType = typename AmgTraits::VType;
    using Comm = typename AmgTraits::Comm;
    using Smoother = typename AmgTraits::Smoother;
    using AMGType = Dune::Amg::AMG<typename AmgTraits::LinearOperator, VType, Smoother,Comm>;
    using BCRSMat = typename AmgTraits::LinearOperator::matrix_type;
    using DofMapper = typename AmgTraits::DofMapper;
91
public:
92
    /*!
93
     * \brief Construct the backend for the sequential case only
94
     *
95
     * \param paramGroup the parameter group for parameter lookup
96
     */
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    ParallelAMGBackend(const std::string& paramGroup = "")
    : LinearSolver(paramGroup)
    , firstCall_(true)
    {
        if (Dune::MPIHelper::getCollectiveCommunication().size() > 1)
            DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
    }

    /*!
     * \brief Construct the backend for parallel or sequential runs
     *
     * \param gridView the grid view on which we are performing the multi-grid
     * \param dofMapper an index mapper for dof entities
     * \param paramGroup the parameter group for parameter lookup
     */
    ParallelAMGBackend(const GridView& gridView,
                       const DofMapper& dofMapper,
                       const std::string& paramGroup = "")
    : LinearSolver(paramGroup)
    , phelper_(std::make_shared<ParallelISTLHelper<GridView, AmgTraits>>(gridView, dofMapper))
117
118
    , firstCall_(true)
    {}
119

120
121
    /*!
     * \brief Solve a linear system.
122
     *
123
124
125
126
     * \param A the matrix
     * \param x the seeked solution vector, containing the initial solution upon entry
     * \param b the right hand side vector
     */
127
128
129
    template<class Matrix, class Vector>
    bool solve(Matrix& A, Vector& x, Vector& b)
    {
130
        int rank = 0;
131
132
133
        std::shared_ptr<Comm> comm;
        std::shared_ptr<LinearOperator> fop;
        std::shared_ptr<ScalarProduct> sp;
134
        static const bool isParallel = AmgTraits::isParallel;
135
136
        prepareLinearAlgebra_<Matrix, Vector, isParallel>(A, b, rank, comm, fop, sp);

137
138
139
        using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
        using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat, Dune::Amg::FirstDiagonal>>;

140
        //! \todo Check whether the default accumulation mode atOnceAccu is needed.
141
        //! \todo make parameters changeable at runtime from input file / parameter tree
142
        Dune::Amg::Parameters params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu);
143
        params.setDefaultValuesIsotropic(Grid::dimension);
144
        params.setDebugLevel(this->verbosity());
145
        Criterion criterion(params);
146
147
148
        SmootherArgs smootherArgs;
        smootherArgs.iterations = 1;
        smootherArgs.relaxationFactor = 1;
149

150
        AMGType amg(*fop, criterion, smootherArgs, *comm);
151
152
        Dune::BiCGSTABSolver<VType> solver(*fop, *sp, amg, this->residReduction(), this->maxIter(),
                                           rank == 0 ? this->verbosity() : 0);
153
154
155

        solver.apply(x, b, result_);
        firstCall_ = false;
156
        return result_.converged;
157
    }
158

159
160
161
    /*!
     * \brief The name of the solver
     */
162
163
164
165
166
    std::string name() const
    {
        return "AMG preconditioned BiCGSTAB solver";
    }

167
168
169
    /*!
     * \brief The result containing the convergence history.
     */
170
171
172
173
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }
174

175
private:
176
177
178
179
180
181
182
183
184
185
186
187
188
189

    /*!
     * \brief Prepare the linear algebra member variables.
     *
     * At compile time, correct constructor calls have to be chosen,
     * depending on whether the setting is parallel or sequential.
     * Since several template parameters are present, this cannot be solved
     * by a full function template specialization. Instead, class template
     * specialization has to be used.
     * This adapts example 4 from http://www.gotw.ca/publications/mill17.htm.

     * The function is called from the solve function. The call is
     * forwarded to the corresponding function of a class template.
     *
190
191
     * \tparam Matrix the matrix type
     * \tparam Vector the vector type
192
193
194
195
     * \tparam isParallel decides if the setting is parallel or sequential
     */
    template<class Matrix, class Vector, bool isParallel>
    void prepareLinearAlgebra_(Matrix& A, Vector& b, int& rank,
196
197
198
                               std::shared_ptr<Comm>& comm,
                               std::shared_ptr<LinearOperator>& fop,
                               std::shared_ptr<ScalarProduct>& sp)
199
    {
200
        LinearAlgebraPreparator<GridView, AmgTraits, isParallel>
201
          ::prepareLinearAlgebra(A, b, rank, comm, fop, sp,
202
                                 *phelper_, firstCall_);
203
204
    }

205
    std::shared_ptr<ParallelISTLHelper<GridView, AmgTraits>> phelper_;
206
    Dune::InverseOperatorResult result_;
207
    bool firstCall_;
208
209
210
211
};

} // namespace Dumux

212
213
214
215
216
217
218
219
220
221
222
223
#include <dumux/common/properties.hh>
#include <dumux/linear/amgtraits.hh>

namespace Dumux {

/*!
 * \ingroup Linear
 * \brief A linear solver based on the ISTL AMG preconditioner
 *        and the ISTL BiCGSTAB solver.
 * \note This is an adaptor using a TypeTag
 */
template<class TypeTag>
224
225
226
using AMGBackend = ParallelAMGBackend<GetPropType<TypeTag, Properties::GridView>, AmgTraits<GetPropType<TypeTag, Properties::JacobianMatrix>,
                                                                                           GetPropType<TypeTag, Properties::SolutionVector>,
                                                                                           GetPropType<TypeTag, Properties::FVGridGeometry>>>;
227
228
229

} // namespace Dumux

230
#endif // DUMUX_AMGBACKEND_HH