amgbackend.hh 9.07 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
#include <dune/common/parallel/mpicommunication.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/parallelhelpers.hh>
40

41
42
namespace Dumux {

43
/*!
44
45
46
 * \ingroup Linear
 * \brief A linear solver based on the ISTL AMG preconditioner
 *        and the ISTL BiCGSTAB solver.
47
 */
48
49
template <class LinearSolverTraits>
class AMGBiCGSTABBackend : public LinearSolver
50
{
51
public:
52
    /*!
53
     * \brief Construct the backend for the sequential case only
54
     *
55
     * \param paramGroup the parameter group for parameter lookup
56
     */
57
    AMGBiCGSTABBackend(const std::string& paramGroup = "")
58
    : LinearSolver(paramGroup)
59
    , isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
60
    {
61
        if (isParallel_)
62
            DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
63
64

        checkAvailabilityOfDirectSolver_();
65
66
67
68
69
70
71
72
73
    }

    /*!
     * \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
     */
74
75
    AMGBiCGSTABBackend(const typename LinearSolverTraits::GridView& gridView,
                       const typename LinearSolverTraits::DofMapper& dofMapper,
76
77
                       const std::string& paramGroup = "")
    : LinearSolver(paramGroup)
78
#if HAVE_MPI
79
    , isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
80
#endif
81
    {
82
83
84
85
#if HAVE_MPI
        if (isParallel_)
            phelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
#endif
86
        checkAvailabilityOfDirectSolver_();
87
    }
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
        return result_.converged;
105
    }
106

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

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

123
private:
124
125
126
127
128
129
130
131
132
133
    //! see https://gitlab.dune-project.org/core/dune-istl/-/issues/62
    void checkAvailabilityOfDirectSolver_()
    {
#if !HAVE_SUPERLU && !HAVE_UMFPACK
        std::cout << "\nAMGBiCGSTABBackend: No direct solver backend found. Using iterative solver as coarse grid solver.\n"
                  << "Note that dune-istl currently hard-codes a tolerance of 1e-2 for the iterative coarse grid solver.\n"
                  << "This may result in reduced accuracy or performance depending on your setup.\nConsider installing "
                  << "UMFPack (SuiteSparse) or SuperLU or apply the istl patch, see dumux/patches/README.md." << std::endl;
#endif
    }
134

135
136
137
138
139
140
141
#if HAVE_MPI
    template<class Matrix, class Vector>
    void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b)
    {
        if constexpr (LinearSolverTraits::canCommunicate)
        {
            if (isParallel_)
142
143
144
145
146
147
148
149
150
151
152
153
            {
                if (LinearSolverTraits::isNonOverlapping(phelper_->gridView()))
                {
                    using PTraits = typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>;
                    solveParallel_<PTraits>(A, x, b);
                }
                else
                {
                    using PTraits = typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>;
                    solveParallel_<PTraits>(A, x, b);
                }
            }
154
155
156
157
158
159
160
161
162
            else
                solveSequential_(A, x, b);
        }
        else
        {
            solveSequential_(A, x, b);
        }
    }

163
    template<class ParallelTraits, class Matrix, class Vector>
164
165
    void solveParallel_(Matrix& A, Vector& x, Vector& b)
    {
166
167
168
        using Comm = typename ParallelTraits::Comm;
        using LinearOperator = typename ParallelTraits::LinearOperator;
        using ScalarProduct = typename ParallelTraits::ScalarProduct;
169
170
171
172

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

        using SeqSmoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
176
        using Smoother = typename ParallelTraits::template Preconditioner<SeqSmoother>;
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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
        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_);
    }

225
#if HAVE_MPI
226
    std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> phelper_;
227
#endif
228
    Dune::InverseOperatorResult result_;
229
    bool isParallel_ = false;
230
231
};

232
} // end namespace Dumux
233

234
#endif