istlsolverfactorybackend.hh 9.92 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 <dune/common/version.hh>
30
31
32
33
34
35
36
37
38
39
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>

// preconditioners
#include <dune/istl/preconditioners.hh>
#include <dune/istl/paamg/amg.hh>

// solvers
#include <dune/istl/solvers.hh>

40
41
42
43
#include <dumux/linear/solver.hh>
#include <dumux/linear/amgtraits.hh>
#include <dumux/linear/amgparallelhelpers.hh>

44
45
#if DUNE_VERSION_NEWER_REV(DUNE_ISTL,2,7,1)

46
47
namespace Dumux {

48
49
/*!
 * \ingroup Linear
50
 * \brief A linear solver using the dune-istl solver factory
51
52
 *        allowing choosing the solver and preconditioner
 *        at runtime.
53
54
 * \note the solvers are configured via the input file
 * \note Requires Dune version 2.7.1 or newer
55
56
 */
template <class Matrix, class Vector, class GridGeometry>
57
class IstlSolverFactoryBackend : public LinearSolver
58
59
60
61
62
63
64
65
66
67
68
69
{
    using GridView = typename GridGeometry::GridView;
    using AMGTraits =  AmgTraits<Matrix, Vector, GridGeometry>;
    using Grid = typename GridView::Grid;
    using LinearOperator = typename AMGTraits::LinearOperator;
    using ScalarProduct = typename AMGTraits::ScalarProduct;
    using VType = typename AMGTraits::VType;
    using Comm = typename AMGTraits::Comm;
    using BCRSMat = typename AMGTraits::LinearOperator::matrix_type;
    using DofMapper = typename AMGTraits::DofMapper;

public:
70
71
72
    //! translation table for solver parameters
    static std::vector<std::array<std::string,2> > istlToDumuxSolverParams;

73
74
75
76
77
    /*!
     * \brief Construct the backend for the sequential case only
     *
     * \param paramGroup the parameter group for parameter lookup
     */
78
79
    IstlSolverFactoryBackend(const std::string& paramGroup = "")
    : firstCall_(true)
80
81
82
    {
        if (Dune::MPIHelper::getCollectiveCommunication().size() > 1)
            DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
83
84
85

        resetDefaultParameters();
        convertParameterTree_(paramGroup);
86
    }
87

88
89
90
91
92
93
94
    /*!
     * \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
     */
95
96
97
98
99
    IstlSolverFactoryBackend(const GridView& gridView,
                             const DofMapper& dofMapper,
                             const std::string& paramGroup = "")
    : phelper_(std::make_shared<ParallelISTLHelper<GridView, AMGTraits>>(gridView, dofMapper))
    , firstCall_(true)
100
    {
101
102
        resetDefaultParameters();
        convertParameterTree_(paramGroup);
103
    }
104

105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
    /*!
     * \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
     */
    bool solve(Matrix& A, Vector& x, Vector& b)
    {
        int rank = 0;
        std::shared_ptr<Comm> comm;
        std::shared_ptr<LinearOperator> fop;
        std::shared_ptr<ScalarProduct> sp; // not used.

#if HAVE_MPI
        if constexpr (AMGTraits::isParallel)
            prepareLinearAlgebraParallel<AMGTraits>(A, b, comm, fop, sp, *phelper_, firstCall_);
        else
            prepareLinearAlgebraSequential<AMGTraits>(A, comm, fop, sp);
#else
        prepareLinearAlgebraSequential<AMGTraits>(A, comm, fop, sp);
#endif

        if (firstCall_)
        {
            Dune::initSolverFactories<typename AMGTraits::LinearOperator>();
        }
132
        std::shared_ptr<Dune::InverseOperator<Vector, Vector>> solver;
133
        try{
134
135
136
137
138
139
140
141
142
            solver = 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;
        }
        try
        {
143
144
            solver->apply(x,b,result_);
        }catch(Dune::Exception& e){
145
            std::cerr << "Exception thrown during linear solve." << std::endl;
146
147
148
149
150
151
            std::cerr << e.what() << std::endl;
            throw e;
        }
        firstCall_ = false;
        return result_.converged;
    }
152
153
154
155
156
157
158
159
160
161
162
163

    //! reset some defaults for the solver parameters
    void resetDefaultParameters()
    {
        params_["restart"] = "10";
        params_["maxit"] = "250";
        params_["reduction"] = "1e-13";
        params_["verbose"] = "0";
        params_["preconditioner.iterations"] = "1";
        params_["preconditioner.relaxation"] = "1.0";
    }

164
165
private:

166
    void convertParameterTree_(const std::string& paramGroup="")
167
168
    {
        const auto& loggingTree = Parameters::getTree();
Timo Koch's avatar
Timo Koch committed
169
        auto matchingGroups = loggingTree.getSubGroups("LinearSolver", paramGroup);
170

171
        for (const auto& [istlKey, dumuxKey] : istlToDumuxSolverParams)
172
173
174
        {
            for (const auto fullGroup : matchingGroups)
            {
Timo Koch's avatar
Timo Koch committed
175
176
                auto istlName = fullGroup + "." + istlKey;
                auto dumuxName = fullGroup + "." + dumuxKey;
177
178
179
180
181
182
183
184
                if(loggingTree.hasKeyOrDefaultKey(dumuxName))
                {
                    if(loggingTree.hasKeyOrDefaultKey(istlName))
                    {
                        std::cerr << "Found equivalent keys " << istlName
                                  << " " << dumuxName << std::endl
                                  << "Please use only one (e.g. " << dumuxName
                                  << ")." << std::endl;
185
                        DUNE_THROW(Dune::InvalidStateException, "Ambiguous parameters used for linear solver");
186
                    }
Timo Koch's avatar
Timo Koch committed
187
                    params_[istlKey] = loggingTree.get<std::string>(dumuxName);
188
189
190
191
                    break;
                }
                else if (loggingTree.hasKey(istlName))
                {
Timo Koch's avatar
Timo Koch committed
192
                    params_[istlKey] = loggingTree.get<std::string>(istlName);
193
194
195
196
197
                    break;
                }
            }
        }

198
        // prevent throw in solve
199
200
201
202
203
204
205
206
207
208
209
        if (!params_.hasKey("type"))
            DUNE_THROW(Dune::InvalidStateException, "Solverfactory needs a specified \"type\" key to select the solver");
    }

    std::shared_ptr<ParallelISTLHelper<GridView, AMGTraits>> phelper_;
    bool firstCall_;
    Dune::InverseOperatorResult result_;
    Dune::ParameterTree params_;
};

template<class Matrix, class Vector, class Geometry>
210
211
212
213
214
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
241
242
243
244
245
246
247
248
std::vector<std::array<std::string, 2>>
IstlSolverFactoryBackend<Matrix, Vector, Geometry>::istlToDumuxSolverParams =
{
    // solver params
    {"verbose", "Verbosity"},
    {"maxit", "MaxIterations"},
    {"reduction", "ResidualReduction"},
    {"type", "Type"},
    {"restart", "Restart"}, // cycles before restarting
    {"mmax", "MaxOrthogonalizationVectors"},

    // preconditioner params
    {"preconditioner.verbosity", "PreconditionerVerbosity"},
    {"preconditioner.type", "PreconditionerType"},
    {"preconditioner.iterations", "PreconditionerIterations"},
    {"preconditioner.relaxation", "PreconditionerRelaxation"},
    {"preconditioner.n", "ILUOrder"},
    {"preconditioner.resort", "ILUResort"},
    {"preconditioner.smootherRelaxation", "AmgSmootherRelaxation"},
    {"preconditioner.smootherIterations", "AmgSmootherIterations"},
    {"preconditioner.maxLevel", "AmgMaxLevel"},
    {"preconditioner.coarsenTarget", "AmgCoarsenTarget"},
    {"preconditioner.minCoarseningRate", "MinCoarseningRate"},
    {"preconditioner.prolongationDampingFactor", "AmgProlongationDampingFactor"},
    {"preconditioner.alpha", "AmgAlpha"},
    {"preconditioner.beta", "AmgBeta"},
    {"preconditioner.additive", "AmgAdditive"},
    {"preconditioner.gamma", "AmgGamma"},
    {"preconditioner.preSteps", "AmgPreSmoothingSteps"},
    {"preconditioner.postSteps", "AmgPostSmoothingSteps"},
    {"preconditioner.criterionSymmetric", "AmgCriterionSymmetric"},
    {"preconditioner.strengthMeasure", "AmgStrengthMeasure"},
    {"preconditioner.diagonalRowIndex", "AmgDiagonalRowIndex"},
    {"preconditioner.defaultAggregationSizeMode", "DefaultAggregationSizeMode"},
    {"preconditioner.defaultAggregationDimension", "defaultAggregationDimension"},
    {"preconditioner.maxAggregateDistance", "MaxAggregateDistance"},
    {"preconditioner.minAggregateSize", "MinAggregateSize"},
    {"preconditioner.maxAggregateSize", "MaxAggregateSize"}
};
249
250

} // end namespace Dumux
251

252
#else
253
#warning "Generic dune-istl solver factory backend needs dune-istl >= 2.7.1!"
254
#endif // DUNE version check
255
#endif // header guard