amgbackend.hh 12.4 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
21
 *   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
 *
22
 * \brief Provides linear solvers using the PDELab AMG backends.
23
24
25
26
 */
#ifndef DUMUX_AMGBACKEND_HH
#define DUMUX_AMGBACKEND_HH

27
28
29
30
31
32
#include <dune/common/parallel/indexset.hh>
#include <dune/common/parallel/mpicollectivecommunication.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/solvers.hh>

33
#include <dumux/linear/linearsolverproperties.hh>
34
35
#include <dumux/linear/amgproperties.hh>
#include <dumux/linear/amgparallelhelpers.hh>
36
37
#include <dumux/linear/p0fem.hh>

38
39
namespace Dumux {

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

        typedef typename Matrix::block_type MatrixBlock;
        MatrixBlock diagonal = matrix[rowIdx][rowIdx];
        diagonal.invert();

        typedef typename Vector::block_type VectorBlock;
        const VectorBlock b = rhs[rowIdx];
        diagonal.mv(b, rhs[rowIdx]);

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

70
/*!
71
 * \brief Provides a linear solver using the parallel PDELab AMG backend.
72
 */
73
template <class TypeTag>
74
75
class AMGBackend
{
76
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
77
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
78
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
79
80
    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
    enum { numEq = JacobianMatrix::block_type::rows};
81

82
83
84
85
86
87
88
89
    typedef typename GET_PROP(TypeTag, AMGPDELabBackend) PDELabBackend;
    typedef typename PDELabBackend::LinearOperator LinearOperator;
    typedef typename PDELabBackend::VType VType;
    typedef typename PDELabBackend::Comm Comm;
    typedef typename PDELabBackend::Smoother Smoother;
    typedef Dune::Amg::AMG<typename PDELabBackend::LinearOperator, VType,
                           Smoother,Comm> AMGType;
    typedef typename PDELabBackend::LinearOperator::matrix_type BCRSMat;
90

91
public:    
92
93
94
95
96
    /*!
     * \brief Construct the backend.
     * 
     * \param problem the problem at hand
     */
97
    AMGBackend(const Problem& problem)
98
    : problem_(problem), phelper_(problem_), firstCall_(true)
99
    {
100
    }
101

102
103
104
105
106
107
108
    /*!
     * \brief Solve a linear system.
     * 
     * \param A the matrix
     * \param x the seeked solution vector, containing the initial solution upon entry
     * \param b the right hand side vector
     */
109
110
111
    template<class Matrix, class Vector>
    bool solve(Matrix& A, Vector& x, Vector& b)
    {
112
113
        int maxIt = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
        int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
114
        static const double residReduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
115
        
116
117
118
#if HAVE_MPI
        Dune::SolverCategory::Category category = PDELabBackend::isNonOverlapping?
            Dune::SolverCategory::nonoverlapping : Dune::SolverCategory::overlapping;
119

120
        if(PDELabBackend::isNonOverlapping && firstCall_)
121
        {
122
            phelper_.initGhostsAndOwners();
123
124
125
126
127
128
        }

        typename PDELabBackend::Comm comm(problem_.gridView().comm(), category);
        
        if(PDELabBackend::isNonOverlapping)
        {
129
130
131
132
            // extend the matrix pattern such that it is usable for AMG
            EntityExchanger<TypeTag> exchanger(problem_);
            exchanger.getExtendedMatrix(A, phelper_);
            exchanger.sumEntries(A);
133
        }
134
135
        phelper_.createIndexSetAndProjectForAMG(A, comm);

136
137
        typename PDELabBackend::LinearOperator fop(A, comm);
        typename PDELabBackend::ScalarProduct sp(comm);
138
        int rank = comm.communicator().rank();
139
140
141
142
143
144

        // Make rhs consistent
        if(PDELabBackend::isNonOverlapping)
        {
            phelper_.makeNonOverlappingConsistent(b);
        }
145
146
147
148
149
150
151
152
153
154
155
#else
        typename PDELabBackend::Comm  comm;
        typename PDELabBackend::LinearOperator fop(A);
        typename PDELabBackend::ScalarProduct sp;
        int rank=0;
#endif
        typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments
            SmootherArgs;
        typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,
                                                                          Dune::Amg::FirstDiagonal> >
            Criterion;
156
157
        // \todo Check whether the default accumulation mode atOnceAccu is needed.
        Dune::Amg::Parameters params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu);
158
159
160
        params.setDefaultValuesIsotropic(GET_PROP_TYPE(TypeTag, GridView)::Traits::Grid::dimension);
        params.setDebugLevel(verbosity);
        Criterion criterion(params);
161
162
163
        SmootherArgs smootherArgs;
        smootherArgs.iterations = 1;
        smootherArgs.relaxationFactor = 1;
164
165

        AMGType amg(fop, criterion, smootherArgs, comm);        
166
167
        Dune::BiCGSTABSolver<typename PDELabBackend::VType> solver(fop, sp, amg, residReduction, maxIt, 
                                                                 rank==0?verbosity: 0);
168

169
170
171

        solver.apply(x, b, result_);
        firstCall_ = false;
172
        return result_.converged;
173
        }
174

175
176
177
    /*!
     * \brief The result containing the convergence history.
     */
178
179
180
181
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }
182
    
183
184
private:
    const Problem& problem_;
185
    ParallelISTLHelper<TypeTag> phelper_;
186
    Dune::InverseOperatorResult result_;
187
    bool firstCall_;
188
189
};

190
191
192
/*!
 * \brief Provides a linear solver using the sequential PDELab AMG backend.
 */
193
194
195
196
template <class TypeTag>
class SeqAMGBackend
{
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
197
198
199
200
201
202
203
204
205
    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
    enum { numEq = JacobianMatrix::block_type::rows};
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> > MType;
    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > VType;
    typedef Dune::Amg::SequentialInformation Comm;
    typedef Dune::AssembledLinearOperator<MType,VType, VType> LinearOperator;
    typedef Dune::SeqSSOR<MType,VType, VType> Smoother;
    typedef Dune::Amg::AMG<LinearOperator,VType,Smoother> AMGType;
206
207
public:

208
209
210
211
212
    /*!
     * \brief Construct the backend.
     * 
     * \param problem the problem at hand
     */
213
214
215
216
    SeqAMGBackend(const Problem& problem)
    : problem_(problem)
    {}

217
218
219
220
221
222
223
    /*!
     * \brief Solve a linear system.
     * 
     * \param A the matrix
     * \param x the seeked solution vector, containing the initial solution upon entry
     * \param b the right hand side vector
     */
224
225
226
    template<class Matrix, class Vector>
    bool solve(Matrix& A, Vector& x, Vector& b)
    {
227
228
229
        int maxIt = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
        int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
        static const double residReduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
        LinearOperator fop(A);
        typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<typename LinearOperator::matrix_type,
                                                                          Dune::Amg::FirstDiagonal> >
            Criterion;
        Criterion criterion(15,2000);
        criterion.setDefaultValuesIsotropic(2);
        typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
        SmootherArgs smootherArgs;
        smootherArgs.iterations = 1;
        smootherArgs.relaxationFactor = 1;

        AMGType amg(fop, criterion, smootherArgs, 1, 1, 1);
        Dune::BiCGSTABSolver<VType> solver(fop, amg, residReduction, maxIt, 
                                          verbosity);
        solver.apply(x, b, result_);        
245
246
247
        return result_.converged;
    }

248
249
250
    /*!
     * \brief The result containing the convergence history.
     */
251
252
253
254
255
256
257
258
259
260
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }

private:
    const Problem& problem_;
    Dune::InverseOperatorResult result_;
};

261
262
263
264
265
266
/*!
 * \brief Provides a linear solver using the sequential PDELab AMG backend.
 * 
 * The linear system is scaled beforehand, possibly improving the 
 * convergence behavior of the iterative solver.
 */
267
268
269
270
template <class TypeTag>
class ScaledSeqAMGBackend
{
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
271
272
273
274
275
276
277
278
279
    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
    enum { numEq = JacobianMatrix::block_type::rows};
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> > MType;
    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > VType;
    typedef Dune::Amg::SequentialInformation Comm;
    typedef Dune::MatrixAdapter<MType,VType, VType> LinearOperator;
    typedef Dune::SeqSSOR<MType,VType, VType> Smoother;
    typedef Dune::Amg::AMG<LinearOperator,VType,Smoother> AMGType;
280
281
public:

282
283
284
285
286
    /*!
     * \brief Construct the backend.
     * 
     * \param problem the problem at hand
     */
287
288
289
290
    ScaledSeqAMGBackend(const Problem& problem)
    : problem_(problem)
    {}

291
292
293
294
295
296
297
    /*!
     * \brief Solve a linear system.
     * 
     * \param A the matrix
     * \param x the seeked solution vector, containing the initial solution upon entry
     * \param b the right hand side vector
     */
298
299
300
301
302
    template<class Matrix, class Vector>
    bool solve(Matrix& A, Vector& x, Vector& b)
    {
        scaleLinearSystem(A, b);

303
304
305
        int maxIt = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
        int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
        static const double residReduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
306
307
308
309
310
311
312
313
314
315
316
317
318
        LinearOperator fop(A);
        typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MType,
                                                                          Dune::Amg::FirstDiagonal> >
            Criterion;
        Criterion criterion(15,2000);
        criterion.setDefaultValuesIsotropic(2);
        typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
        SmootherArgs smootherArgs;
        smootherArgs.iterations = 1;
        smootherArgs.relaxationFactor = 1;
        AMGType amg(fop, criterion, smootherArgs, 1, 1, 1);
        Dune::BiCGSTABSolver<VType> solver(fop, amg, residReduction, maxIt, 
                                          verbosity);
319
        
320
321
322
        return result_.converged;
    }

323
324
325
    /*!
     * \brief The result containing the convergence history.
     */
326
327
328
329
330
331
332
333
334
335
336
337
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }

private:
    const Problem& problem_;
    Dune::InverseOperatorResult result_;
};

} // namespace Dumux

338
#endif // DUMUX_AMGBACKEND_HH