istlsolverfactorybackend.hh 11.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   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 3 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
 * \ingroup Linear
 * \brief Provides a generic linear solver based on the ISTL that chooses the
23
 *        solver and preconditioner at runtime
24
25
 */

26
27
#ifndef DUMUX_LINEAR_ISTL_SOLVERFACTORYBACKEND_HH
#define DUMUX_LINEAR_ISTL_SOLVERFACTORYBACKEND_HH
28

29
#include <memory>
30

31
#include <dune/common/version.hh>
32
33
34
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>

35
36
#include "linearsolverparameters.hh"

37
38
39
// preconditioners
#include <dune/istl/preconditioners.hh>
#include <dune/istl/paamg/amg.hh>
40
#include "preconditioners.hh"
41
42
43

// solvers
#include <dune/istl/solvers.hh>
44
#include <dune/istl/solverfactory.hh>
45

46
#include <dumux/common/typetraits/matrix.hh>
47
#include <dumux/linear/solver.hh>
48
#include <dumux/linear/parallelhelpers.hh>
49
#include <dumux/linear/istlsolverregistry.hh>
50
51
52

namespace Dumux {

53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
// initSolverFactoriesForMultiTypeBlockMatrix differs in different compilation units,
// so we have it in an anonymous namespace
namespace {

/*!
 * \brief Initializes the direct solvers, preconditioners and iterative solvers in
 *        the factories with the corresponding Matrix and Vector types.
 * \note  We currently consider only direct solvers and preconditioners provided by
 *        Dumux which, unlike the ones implemented in Dune, also support MultiTypeBlockMatrices.
 * \note This function could be removed once Dune::initSolverFactories supports MultiTypeBlockMatrices.
 * \tparam LinearOperator the linear operator type
 */
template<class LinearOperator>
int initSolverFactoriesForMultiTypeBlockMatrix()
{
    using M  = typename LinearOperator::matrix_type;
    using X  = typename LinearOperator::range_type;
    using Y  = typename LinearOperator::domain_type;
    using TL = Dune::TypeList<M,X,Y>;
    auto& dsfac = Dune::DirectSolverFactory<M,X,Y>::instance();
    Dune::addRegistryToFactory<TL>(dsfac, Dumux::MultiTypeBlockMatrixDirectSolverTag{});
74
#if DUNE_VERSION_GT(DUNE_ISTL,2,7)
75
    auto& pfac = Dune::PreconditionerFactory<LinearOperator,X,Y>::instance();
76
77
78
#else
    auto& pfac = Dune::PreconditionerFactory<M,X,Y>::instance();
#endif
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    Dune::addRegistryToFactory<TL>(pfac, Dumux::MultiTypeBlockMatrixPreconditionerTag{});
    using TLS = Dune::TypeList<X,Y>;
    auto& isfac = Dune::IterativeSolverFactory<X,Y>::instance();
    return Dune::addRegistryToFactory<TLS>(isfac, Dune::IterativeSolverTag{});
}
} // end namespace

/*!
 * \brief Initialize the solver factories for regular matrices or MultiTypeBlockMatrices
 * \tparam Matrix the matrix
 * \tparam LinearOperator the linear operator
 *
 * \note This function could be removed once Dune::initSolverFactories supports MultiTypeBlockMatrices.
 */
template<class Matrix, class LinearOperator>
void initSolverFactories()
{
    if constexpr (isMultiTypeBlockMatrix<Matrix>::value)
        initSolverFactoriesForMultiTypeBlockMatrix<LinearOperator>();
    else
99
#if DUNE_VERSION_GT(DUNE_ISTL,2,7)
100
        Dune::initSolverFactories<LinearOperator>();
101
102
103
104
105
106
107
#else
    {
        using X  = typename LinearOperator::range_type;
        using Y  = typename LinearOperator::domain_type;
        Dune::initSolverFactories<Matrix, X, Y>();
    }
#endif
108
109
}

110
111
/*!
 * \ingroup Linear
112
 * \brief A linear solver using the dune-istl solver factory
Timo Koch's avatar
Timo Koch committed
113
 *        to choose the solver and preconditioner at runtime.
114
 * \note the solvers are configured via the input file
115
 * \note requires Dune version 2.7.1 or newer
116
 */
117
template <class LinearSolverTraits>
118
class IstlSolverFactoryBackend : public LinearSolver
119
120
{
public:
121

122
123
124
125
126
    /*!
     * \brief Construct the backend for the sequential case only
     *
     * \param paramGroup the parameter group for parameter lookup
     */
127
    IstlSolverFactoryBackend(const std::string& paramGroup = "")
128
    : paramGroup_(paramGroup)
129
    , isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
130
    {
131
        if (isParallel_)
132
            DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
133

134
        reset();
135
    }
136

137
138
139
    /*!
     * \brief Construct the backend for parallel or sequential runs
     *
Timo Koch's avatar
Timo Koch committed
140
     * \param gridView the grid view for parallel communication via the grid
141
142
143
     * \param dofMapper an index mapper for dof entities
     * \param paramGroup the parameter group for parameter lookup
     */
144
145
    IstlSolverFactoryBackend(const typename LinearSolverTraits::GridView& gridView,
                             const typename LinearSolverTraits::DofMapper& dofMapper,
146
                             const std::string& paramGroup = "")
147
    : paramGroup_(paramGroup)
148
#if HAVE_MPI
149
150
    , parallelHelper_(std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper))
    , isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
151
#endif
152
    {
153
        reset();
154
    }
155

156
157
158
159
160
161
162
    /*!
     * \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
     */
163
    template<class Matrix, class Vector>
164
165
166
    bool solve(Matrix& A, Vector& x, Vector& b)
    {
#if HAVE_MPI
167
        solveSequentialOrParallel_(A, x, b);
168
#else
169
        solveSequential_(A, x, b);
170
171
172
173
#endif
        firstCall_ = false;
        return result_.converged;
    }
174

175
176
177
    //! reset the linear solver factory
    void reset()
    {
178
        firstCall_ = true;
179
        params_ = LinearSolverParameters<LinearSolverTraits>::createParameterTree(paramGroup_);
180
181
182
183
        checkMandatoryParameters_();
        name_ = params_.get<std::string>("preconditioner.type") + "-preconditioned " + params_.get<std::string>("type");
        if (params_.get<int>("verbose", 0) > 0)
            std::cout << "Initialized linear solver of type: " << name_ << std::endl;
184
185
    }

186
187
188
189
190
191
192
193
194
195
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }

    const std::string& name() const
    {
        return name_;
    }

196
197
private:

198
199
    void checkMandatoryParameters_()
    {
200
        if (!params_.hasKey("type"))
201
202
203
204
205
206
            DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Type\" parameter to select the solver");

        if (!params_.hasKey("preconditioner.type"))
            DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Preconditioner.Type\" parameter to select the preconditioner");
    }

207
208
209
210
#if HAVE_MPI
    template<class Matrix, class Vector>
    void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b)
    {
211
212
213
214
        // Dune::MultiTypeBlockMatrix does not provide a ColIterator which is needed by Dune::NonoverlappingSchwarzOperator.
        // We therefore can only solve these types of systems sequentially.
        // TODO: This can be adapted once the situation in Dune ISTL changes.
        if constexpr (LinearSolverTraits::canCommunicate && !isMultiTypeBlockMatrix<Matrix>::value)
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
        {
            if (isParallel_)
            {
                if (LinearSolverTraits::isNonOverlapping(parallelHelper_->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);
                }
            }
            else
                solveSequential_(A, x, b);
        }
        else
        {
            solveSequential_(A, x, b);
        }
    }

    template<class ParallelTraits, class Matrix, class Vector>
    void solveParallel_(Matrix& A, Vector& x, Vector& b)
    {
241
#if DUNE_VERSION_GT_REV(DUNE_ISTL,2,7,0)
242
243
244
245
246
247
        using Comm = typename ParallelTraits::Comm;
        using LinearOperator = typename ParallelTraits::LinearOperator;
        using ScalarProduct = typename ParallelTraits::ScalarProduct;

        if (firstCall_)
        {
248
            initSolverFactories<Matrix, LinearOperator>();
249
250
251
252
253
254
255
256
257
258
259
260
261
            parallelHelper_->initGhostsAndOwners();
        }

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

        // construct solver
        auto solver = getSolverFromFactory_(linearOperator);

        // solve linear system
        solver->apply(x, b, result_);
262
#else
263
        DUNE_THROW(Dune::NotImplemented, "Parallel solvers only available for dune-istl > 2.7.0");
264
#endif
265
266
267
268
269
270
271
272
273
274
275
276
    }
#endif // HAVE_MPI

    template<class Matrix, class Vector>
    void solveSequential_(Matrix& A, Vector& x, Vector& b)
    {
        // construct linear operator
        using Traits = typename LinearSolverTraits::template Sequential<Matrix, Vector>;
        using LinearOperator = typename Traits::LinearOperator;
        auto linearOperator = std::make_shared<LinearOperator>(A);

        if (firstCall_)
277
            initSolverFactories<Matrix, LinearOperator>();
278
279
280
281
282
283
284
285
286

        // construct solver
        auto solver = getSolverFromFactory_(linearOperator);

        // solve linear system
        solver->apply(x, b, result_);
    }

    template<class LinearOperator>
287
288
289
290
291
292
293
294
295
    auto getSolverFromFactory_(std::shared_ptr<LinearOperator>& fop)
    {
        try { return Dune::getSolverFromFactory(fop, params_); }
        catch(Dune::Exception& e)
        {
            std::cerr << "Could not create solver with factory" << std::endl;
            std::cerr << e.what() << std::endl;
            throw e;
        }
296
297
    }

298
    const std::string paramGroup_;
299
#if HAVE_MPI
300
    std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> parallelHelper_;
301
302
#endif
    bool isParallel_ = false;
303
    bool firstCall_;
304

305
306
    Dune::InverseOperatorResult result_;
    Dune::ParameterTree params_;
307
    std::string name_;
308
309
310
};

} // end namespace Dumux
311
312

#endif // header guard