1
//---------------------------------------------------------------------------
2
// $Id: trilinos_precondition_block.h 18547 2009-04-03 13:04:41Z kronbichler $
5
// Copyright (C) 2008 by the deal.II authors
7
// This file is subject to QPL and may not be distributed
8
// without copyright and license information. Please refer
9
// to the file deal.II/doc/license.html for the text and
10
// further information on this license.
12
//---------------------------------------------------------------------------
13
#ifndef __deal2__trilinos_precondition_block_h
14
#define __deal2__trilinos_precondition_block_h
17
#include <base/config.h>
18
#include <lac/exceptions.h>
19
#include <lac/trilinos_precondition.h>
20
#include <lac/trilinos_sparse_matrix.h>
23
#ifdef DEAL_II_USE_TRILINOS
25
#include <Teuchos_RCP.hpp>
26
#include <Thyra_LinearOperatorDecl.hpp>
27
#include <Epetra_Operator.h>
29
// some forward declarations
30
class Ifpack_Preconditioner;
33
DEAL_II_NAMESPACE_OPEN
35
/*! @addtogroup TrilinosWrappers
39
namespace TrilinosWrappers
41
// forward declarations
42
class BlockSparseMatrix;
45
* This implements an interface class of preconditioners for block
46
* matrices based on Trilinos, which are intended to be used with the
47
* solver classes SolverBlockCG and SolverBlockGMRES. This is based on
48
* the Trilinos Thyra abstract interface of linear operators. This
49
* class does not implement any method, and derived classes will need
52
* @ingroup TrilinosWrappers
53
* @ingroup Preconditioners
54
* @author Martin Kronbichler, 2008
56
class PreconditionBlockBase : public Subscriptor
62
PreconditionBlockBase();
67
~PreconditionBlockBase();
71
* Preconditioner object.
73
Teuchos::RCP<Thyra::ConstLinearOperator<double> > preconditioner;
75
friend class SolverBlockBase;
79
* Internal function that sets
80
* up an inverse matrix.
83
Thyra::ConstLinearOperator<double>
84
inverse_matrix (const SparseMatrix &M,
85
const Epetra_Operator *P,
86
const bool is_symmetric,
87
const unsigned int n_iterations,
88
const double solve_tolerance,
89
const bool output_details);
92
// Forward declarations.
95
* This class implements a black box preconditioner for saddle points
96
* systems arising from the Stokes or Navier–Stokes equations as
97
* specified by the papers <em>D. Silvester, A. Wathen, Fast iterative
98
* solution of stabilised Stokes systems part II. Using general block
99
* preconditioners, SIAM J. Numer. Anal. 31:1352–1367 (1994)</em>
100
* and <em> D. Kay, D. Loghin, A. Wathen, A preconditioner for the
101
* steady-state Navier–Stokes equations, SIAM
102
* J. Sci. Comput. 24(1):237–256 (2002)</em>, respectively.
104
* The preconditioner is based an approximation to the Schur
105
* complement of the block matrix. The Schur complement $S=B
106
* A_{\mathbf u}^{-1} B^T$ is approximated by a mass matrix $M_p$ on
107
* the pressure space in the case of the Stokes equations, and as a
108
* product $S^{-1} = L_p^{-1} F_p M_p^{-1}$ with pressure Laplace
109
* matrix $L_p$, pressure convection-diffusion operator $F_p$
110
* (corresponding to the sum of time derivative, convection and
111
* diffusion), and pressure mass matrix $M_p$ in the case of the
112
* Navier–Stokes equations.
114
* @ingroup TrilinosWrappers
115
* @ingroup Preconditioners
116
* @author Martin Kronbichler, 2008
118
class PreconditionStokes : public PreconditionBlockBase
122
* Standardized data struct to
123
* pipe additional data to the
126
struct AdditionalData
129
* Constructor. This sets the
130
* additional data to its
131
* default values. For lazy
132
* initializations, the user
133
* will just set the parameters
134
* in here and let the
135
* <tt>initialize</tt> function
136
* set up good preconditioners,
137
* where some good values have
138
* to be provided here. The
139
* default ones mean that we
140
* assume the velocity-velocity
141
* matrix to be symmetric (only
142
* true for Stokes problems!),
143
* we use a relative inner
144
* tolerance of 1e-9 for the
145
* inversion of matrices, and
146
* take in IC preconditioner
147
* for the inversion of the
149
* matrix. Moreover, there is
150
* no inner CG/GMRES solve
152
* velocity-velocity block by
153
* default. Otherwise, we do a
154
* solve with the specified
155
* solver and the number of
157
* tolerance. Note that too
158
* many iterations on Av can
159
* slow down the overall
160
* procedure considerably.
162
AdditionalData (const bool right_preconditioning = false,
163
const bool Av_is_symmetric = true,
164
const std::vector<std::vector<bool> > &Av_constant_modes = std::vector<std::vector<bool> > (1),
165
const double inner_solve_tolerance = 1e-9,
166
const bool use_ssor_on_mass_matrix = false,
167
const bool velocity_uses_higher_order_elements = false,
168
const bool pressure_uses_higher_order_elements = false,
169
const bool Av_do_solve = false,
170
const unsigned int Av_n_iters = 1,
171
const double Av_tolerance = 1e-3,
172
const bool ouput_details = false);
175
* This flag specifies whether
176
* the preconditioner should be
178
* preconditioner or right
179
* preconditioner. Note that
180
* this setting must be the
181
* same as the one used in the
184
bool right_preconditioning;
187
* This flag determines whether
188
* the velocity-velocity matrix
189
* is symmetric (and positive
190
* definite) or not. This
191
* choice will influence the
192
* AMG preconditioner that is
193
* built for that system as
194
* well as the inner solve on
195
* Av (if that flag is
198
bool Av_is_symmetric;
201
* This determines the near
202
* null space for the operator
203
* Av, the vector-valued
204
* velocity-velocity space. In
206
* performance, this vector has
207
* to be specified using the
208
* DoFTools::extract_constant_modes()
211
std::vector<std::vector<bool> > Av_constant_modes;
214
* Tolerance level that the
215
* inner solvers on the
216
* pressure mass matrix (and
217
* the pressure Laplace matrix
218
* in case of a Navier-Stokes
219
* problem) should be solve
220
* to. It is essential that
221
* this tolerance is in the
222
* same order of magnitude than
223
* the one chosen for the outer
224
* GMRES solver (or little
225
* less). Otherwise unexpected
226
* variations in the number of
227
* iterations can occur from
230
* Default value: <tt>1e-9</tt>
231
* residual relative to initial
234
double inner_solve_tolerance;
237
* Determines whether the inner
238
* solve for the pressure mass
239
* matrix should use an IC
240
* preconditioner (incomplete
241
* Cholesky factorization) or
242
* an SSOR preconditioner. The
243
* former tends to be faster in
244
* most situations, whereas the
245
* latter uses almost no
246
* additional memory besides
247
* the matrix storage.
249
bool use_ssor_on_mass_matrix;
252
* Determines whether the
253
* underlying velocity
254
* discretization uses linear
255
* elements or higher order
256
* elements, which has
257
* consequences for the set up
258
* of the internal setup of the
259
* ML AMG preconditioner.
261
bool velocity_uses_higher_order_elements;
264
* Determines whether the
265
* underlying pressure
266
* discretization uses linear
267
* elements or higher order
268
* elements, which has
269
* consequences for the set up
270
* of the internal setup of the
271
* ML AMG preconditioner.
273
bool pressure_uses_higher_order_elements;
276
* Flag to determine whether we
277
* should do an inner solve on
278
* the velocity-velocity
279
* matrix. By default, this
280
* option is not on, since it
281
* is not necessary for a good
282
* performance of the outer
283
* solver. However, it can be
284
* beneficial to perform some
285
* iterations on this block
286
* only and not on the whole
287
* block, so that the number of
288
* outer iterations is
289
* reduced. See the discussion
290
* in the @ref step_22
291
* "step-22" tutorial program.
293
unsigned int Av_do_solve;
296
* The number of iterations in
299
unsigned int Av_n_iters;
302
* The residual tolerance that
303
* is going to be used for the
305
* velocity-velocity matrix.
310
* This defines whether
311
* internal details about the
312
* solve process should be
313
* written to screen. This can
314
* be a lot of information.
320
* Lazy setup of the block
321
* preconditioner for the
323
* system. This function takes
324
* the matrices given here and
325
* firs calculates good
326
* preconditioners, i.e.,
327
* algebraic multigrid
328
* preconditioners for the
329
* velocity-velocity matrix
330
* <tt>Av</tt>, that can be
331
* specified to be different
332
* from the (0,0) block in the
333
* system matrix, IC/SSOR on
334
* the pressure mass matrix
335
* <tt>Mp_matrix</tt>, and AMG
336
* on the pressure Laplace
338
* <tt>Lp_matrix</tt>. Next,
339
* these preconditioners are
340
* used to generate a block
343
void initialize (const BlockSparseMatrix &system_matrix,
344
const SparseMatrix &Mp_matrix,
345
const AdditionalData &additional_data,
346
const SparseMatrix &Av_matrix = SparseMatrix(),
347
const SparseMatrix &Lp_matrix = SparseMatrix(),
348
const SparseMatrix &Fp_matrix = SparseMatrix());
352
* preconditioner for the
353
* (Navier-) Stokes system. In
354
* contrast to the other
355
* <tt>initialize</tt>
356
* function, this one expects
357
* the user to specify
358
* preconditioner objects for
359
* the individual matrices,
360
* that will then be used for
361
* the inner inversion of the
362
* matrices when building the
364
* preconditioner. If no
365
* preconditioner is specified,
366
* the inner solvers will use
367
* an identity preconditioner
370
void initialize (const BlockSparseMatrix &system_matrix,
371
const SparseMatrix &Mp_matrix,
372
const AdditionalData &additional_data,
373
const PreconditionBase &Mp_preconditioner = PreconditionBase(),
374
const PreconditionBase &Av_preconditioner = PreconditionBase(),
375
const SparseMatrix &Lp_matrix = SparseMatrix(),
376
const PreconditionBase &Lp_preconditioner = PreconditionBase(),
377
const SparseMatrix &Fp_matrix = SparseMatrix());
380
* This function can be used
381
* for a faster recalculation
382
* of the preconditioner
383
* construction when the system
384
* matrix entries underlying
385
* the preconditioner have
386
* changed but not the sparsity
387
* pattern (this means that the
388
* grid has remained the same
389
* and the matrix structure has
391
* reinitialized). Moreover, it
392
* is assumed that the PDE
393
* parameters have not changed
394
* drastically. This function
395
* is only needed for the lazy
397
* preconditioner. In the other
398
* case, the preconditioner can
399
* be modified outside this
400
* function, which will be
401
* recongnized in this
402
* preconditioner as well when
403
* doing a solve, without the
404
* need to restructure anything
412
DeclException1 (ExcNonMatchingMaps,
414
<< "The sparse matrix the preconditioner is based on "
415
<< "uses a map that is not compatible to the one in vector"
417
<< ". Check preconditioner and matrix setup.");
422
* Pointer to preconditioner
423
* for the mass matrix.
425
Teuchos::RCP<PreconditionSSOR> Mp_precondition_ssor;
428
* Pointer to preconditioner
429
* for the mass matrix.
431
Teuchos::RCP<PreconditionIC> Mp_precondition_ic;
434
* Pointer to preconditioner
435
* for the velocity-velocity
438
Teuchos::RCP<PreconditionAMG> Av_precondition;
441
* Pointer to preconditioner
442
* for the pressure Laplace
445
Teuchos::RCP<PreconditionAMG> Lp_precondition;
452
DEAL_II_NAMESPACE_CLOSE
454
#endif // DEAL_II_USE_TRILINOS
456
/*-------------------------- trilinos_precondition_block.h ---------------------------*/
459
/*-------------------------- trilinos_precondition_block.h ---------------------------*/