amgbackend.hh 8.73 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
 *                                                                           *
 *   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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
 *   (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
#include <dune/common/parallel/indexset.hh>
32
33
34

#include <dune/common/version.hh>
#if DUNE_VERSION_LT(DUNE_COMMON,2,7)
35
#include <dune/common/parallel/mpicollectivecommunication.hh>
36
37
38
39
#else
#include <dune/common/parallel/mpicommunication.hh>
#endif

40
#include <dune/grid/common/capabilities.hh>
41
42
43
44
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/solvers.hh>

45
#include <dumux/linear/solver.hh>
46
#include <dumux/linear/parallelhelpers.hh>
47

48
49
namespace Dumux {

50
/*!
51
52
53
 * \ingroup Linear
 * \brief A linear solver based on the ISTL AMG preconditioner
 *        and the ISTL BiCGSTAB solver.
54
 */
55
56
template <class LinearSolverTraits>
class AMGBiCGSTABBackend : public LinearSolver
57
{
58
public:
59
    /*!
60
     * \brief Construct the backend for the sequential case only
61
     *
62
     * \param paramGroup the parameter group for parameter lookup
63
     */
64
    AMGBiCGSTABBackend(const std::string& paramGroup = "")
65
    : LinearSolver(paramGroup)
66
    , isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
67
68
    , firstCall_(true)
    {
69
        if (isParallel_)
70
71
72
73
74
75
76
77
78
79
            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
     */
80
81
    AMGBiCGSTABBackend(const typename LinearSolverTraits::GridView& gridView,
                       const typename LinearSolverTraits::DofMapper& dofMapper,
82
83
                       const std::string& paramGroup = "")
    : LinearSolver(paramGroup)
84
85
    , phelper_(std::make_shared<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper))
    , isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
86
87
    , firstCall_(true)
    {}
88

89
90
    /*!
     * \brief Solve a linear system.
91
     *
92
93
94
95
     * \param A the matrix
     * \param x the seeked solution vector, containing the initial solution upon entry
     * \param b the right hand side vector
     */
96
97
98
    template<class Matrix, class Vector>
    bool solve(Matrix& A, Vector& x, Vector& b)
    {
99
#if HAVE_MPI
100
        solveSequentialOrParallel_(A, x, b);
101
#else
102
        solveSequential_(A, x, b);
103
#endif
104
        firstCall_ = false;
105
        return result_.converged;
106
    }
107

108
109
110
    /*!
     * \brief The name of the solver
     */
111
112
    std::string name() const
    {
113
        return "AMG-preconditioned BiCGSTAB solver";
114
115
    }

116
117
118
    /*!
     * \brief The result containing the convergence history.
     */
119
120
121
122
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }
123

124
private:
125

126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
#if HAVE_MPI
    template<class Matrix, class Vector>
    void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b)
    {
        if constexpr (LinearSolverTraits::canCommunicate)
        {
            if (isParallel_)
                solveParallel_(A, x, b);
            else
                solveSequential_(A, x, b);
        }
        else
        {
            solveSequential_(A, x, b);
        }
    }

    template<class Matrix, class Vector>
    void solveParallel_(Matrix& A, Vector& x, Vector& b)
    {
        using Traits = typename LinearSolverTraits::template Parallel<Matrix, Vector>;
        using Comm = typename Traits::Comm;
        using LinearOperator = typename Traits::LinearOperator;
        using ScalarProduct = typename Traits::ScalarProduct;

        std::shared_ptr<Comm> comm;
        std::shared_ptr<LinearOperator> linearOperator;
        std::shared_ptr<ScalarProduct> scalarProduct;
        prepareLinearAlgebraParallel<LinearSolverTraits>(A, b, comm, linearOperator, scalarProduct, *phelper_, firstCall_);

        using SeqSmoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
        using Smoother = typename Traits::template Preconditioner<SeqSmoother>;
        solveWithAmg_<Smoother>(A, x, b, linearOperator, comm, scalarProduct);
    }
#endif // HAVE_MPI

    template<class Matrix, class Vector>
    void solveSequential_(Matrix& A, Vector& x, Vector& b)
    {
        using Comm = Dune::Amg::SequentialInformation;
        using Traits = typename LinearSolverTraits::template Sequential<Matrix, Vector>;
        using LinearOperator = typename Traits::LinearOperator;
        using ScalarProduct = typename Traits::ScalarProduct;

        auto comm = std::make_shared<Comm>();
        auto linearOperator = std::make_shared<LinearOperator>(A);
        auto scalarProduct = std::make_shared<ScalarProduct>();

        using Smoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
        solveWithAmg_<Smoother>(A, x, b, linearOperator, comm, scalarProduct);
    }

    template<class Smoother, class Matrix, class Vector, class LinearOperator, class Comm, class ScalarProduct>
    void solveWithAmg_(Matrix& A, Vector& x, Vector& b,
                       std::shared_ptr<LinearOperator>& linearOperator,
                       std::shared_ptr<Comm>& comm,
                       std::shared_ptr<ScalarProduct>& scalarProduct)
    {
        using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
        using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Matrix, Dune::Amg::FirstDiagonal>>;

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

        using Amg = Dune::Amg::AMG<LinearOperator, Vector, Smoother, Comm>;
        auto amg = std::make_shared<Amg>(*linearOperator, criterion, smootherArgs, *comm);

        Dune::BiCGSTABSolver<Vector> solver(*linearOperator, *scalarProduct, *amg, this->residReduction(), this->maxIter(),
                                            comm->communicator().rank() == 0 ? this->verbosity() : 0);

        solver.apply(x, b, result_);
    }

    std::shared_ptr<ParallelISTLHelper<LinearSolverTraits>> phelper_;
207
    Dune::InverseOperatorResult result_;
208
    bool isParallel_;
209
    bool firstCall_;
210
211
};

212
213
214
215
216
// deprecation helper -> will be removed after 3.2
template<class GridView, class Traits>
using ParallelAMGBackend
= AMGBiCGSTABBackend<Traits>;

217
218
} // namespace Dumux

219
#include <dumux/common/properties.hh>
220
#include <dumux/linear/linearsolvertraits.hh>
221
222
223
224
225
226
227
228
229
230

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>
231
232
233
using AMGBackend [[deprecated("Use AMGBiCGSTABBackend instead. Will be removed after 3.2!")]]
= ParallelAMGBackend<GetPropType<TypeTag, Properties::GridView>,
                     LinearSolverTraits<GetPropType<TypeTag, Properties::GridGeometry>>>;
234
235
236

} // namespace Dumux

237
#endif // DUMUX_AMGBACKEND_HH