solid_preconditioners.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_SOLID_PRECONDITIONERS_HEADER
27#define OOMPH_SOLID_PRECONDITIONERS_HEADER
28
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35
36// oomphlib headers
37#include "../generic/matrices.h"
38#include "../generic/assembly_handler.h"
39#include "../generic/problem.h"
40#include "../generic/block_preconditioner.h"
41#include "../generic/preconditioner.h"
42#include "../generic/SuperLU_preconditioner.h"
43#include "../generic/matrix_vector_product.h"
44
45
46namespace oomph
47{
48 //===========================================================================
49 /// The least-squares commutator (LSC; formerly BFBT)
50 /// preconditioner. It uses blocks corresponding to the displacement/position
51 /// and pressure unknowns, i.e. there are a total of 2x2 blocks,
52 /// and all displacement/position components are treated as a
53 /// single block of unknowns.
54 ///
55 /// Here are the details: An "ideal" preconditioner
56 /// would solve the saddle point system
57 /// \f[ \left( \begin{array}{cc} {\bf F} & {\bf G} \\ {\bf D} & {\bf 0} \end{array} \right) \left( \begin{array}{c} {\bf z}_u \\ {\bf z}_p \end{array} \right) = \left( \begin{array}{c} {\bf r}_u \\ {\bf r}_p \end{array} \right) \f]
58 /// where \f$ {\bf F}\f$, \f$ {\bf G} \f$, and \f$ {\bf D}\f$ are
59 /// the blocks that arise in the Jacobian of the pressure-based
60 /// equations of linear and nonlinear elasticity (with dofs in order
61 /// of displacement/position and pressure).
62 /// The use of this preconditioner would ensure the convergence
63 /// of any iterative linear solver in a single iteration but its
64 /// application is, of course, exactly as expensive as a direct solve.
65 /// The LSC/BFBT preconditioner replaces the exact Jacobian by
66 /// a block-triangular approximation
67 /// \f[ \left( \begin{array}{cc} {\bf F} & {\bf G} \\ {\bf 0} & -{\bf M}_s \end{array} \right) \left( \begin{array}{c} {\bf z}_u \\ {\bf z}_p \end{array} \right) = \left( \begin{array}{c} {\bf r}_u \\ {\bf r}_p \end{array} \right), \f]
68 /// where \f${\bf M}_s\f$ is an approximation to the pressure
69 /// Schur-complement \f$ {\bf S} = {\bf D} {\bf F}^{-1}{\bf G}. \f$
70 /// This system can be solved in two steps:
71 /// -# Solve the second row for \f$ {\bf z}_p\f$ via
72 /// \f[ {\bf z}_p = - {\bf M}_s^{-1} {\bf r}_p \f]
73 /// -# Given \f$ {\bf z}_p \f$ , solve the first row for \f$ {\bf z}_u\f$ via
74 /// \f[ {\bf z}_u = {\bf F}^{-1} \big( {\bf r}_u - {\bf G} {\bf z}_p \big) \f].
75 /// In the LSC/BFBT preconditioner, the action of the inverse pressure
76 /// Schur complement
77 /// \f[ {\bf z}_p = - {\bf M}_s^{-1} {\bf r}_p \f]
78 /// is approximated by
79 /// \f[ {\bf z}_p = - \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)^{-1} \big({\bf D} \widehat{\bf Q}^{-1}{\bf F} \widehat{\bf Q}^{-1}{\bf G}\big) \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)^{-1} {\bf r}_p, \f]
80 /// where \f$ \widehat{\bf Q} \f$ is the diagonal of the
81 /// displacement/position mass matrix. The evaluation of this expression
82 /// involves two linear solves involving the matrix
83 /// \f[ {\bf P} = \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big) \f]
84 /// which has the character of a matrix arising from the discretisation
85 /// of a Poisson problem on the pressure space. We also have
86 /// to evaluate matrix-vector products with the matrix
87 /// \f[ {\bf E}={\bf D}\widehat{\bf Q}^{-1}{\bf F}\widehat{\bf Q}^{-1}{\bf G} \f]
88 /// Details of the theory can be found in "Finite Elements and
89 /// Fast Iterative Solvers with Applications in Incompressible Fluid
90 /// Dynamics" by Howard C. Elman, David J. Silvester, and Andrew J. Wathen,
91 /// published by Oxford University Press, 2006.
92 ///
93 /// In our implementation of the preconditioner, the linear systems
94 /// can either be solved "exactly", using SuperLU (in its incarnation
95 /// as an exact preconditioner; this is the default) or by any
96 /// other Preconditioner (inexact solver) specified via the access functions
97 /// \code
98 /// PressureBasedSolidLSCPreconditioner::set_f_preconditioner(...)
99 /// \endcode
100 /// or
101 /// \code
102 /// PressureBasedSolidLSCPreconditioner::set_p_preconditioner(...)
103 /// \endcode
104 //===========================================================================
106 : public BlockPreconditioner<CRDoubleMatrix>
107 {
108 public:
109 /// Constructor - sets defaults for control flags
112 {
113 // Flag to indicate that the preconditioner has been setup
114 // previously -- if setup() is called again, data can
115 // be wiped.
117
118 // By default we use SuperLU for both p and f blocks
121
122 // resize the mesh pt
123 // note: meaningless if subsidiary preconditioner
124 this->set_nmesh(1);
125 Solid_mesh_pt = 0;
126
127 // Set default preconditioners (inexact solvers) -- they are
128 // members of this class!
131
132 // Flag to determine if mass matrix diagonal Q^{-1}
133 // is used for scaling.
135
136 // set Doc_time to false
137 Doc_time = false;
138
139 // null the off diagonal Block matrix pt
140 Bt_mat_vec_pt = 0;
141
142 // null the F matrix vector product helper
143 F_mat_vec_pt = 0;
144
145 // null the QBt matrix vector product pt
146 QBt_mat_vec_pt = 0;
147
148 // null the E matrix vector product helper
149 E_mat_vec_pt = 0;
150
151 // by default we do not form the E matrix (BQFQBt)
152 Form_BFBt_product = false;
153 }
154
155 /// Destructor
157 {
159 }
160
161 /// Broken copy constructor
164
165 /// Broken assignment operator
166 // Commented out broken assignment operator because this can lead to a
167 // conflict warning when used in the virtual inheritence hierarchy.
168 // Essentially the compiler doesn't realise that two separate
169 // implementations of the broken function are the same and so, quite
170 // rightly, it shouts.
171 /*void operator=(const PressureBasedSolidLSCPreconditioner&) = delete;*/
172
173 /// Setup the preconditioner
174 void setup();
175
176 /// Apply preconditioner to Vector r
178
179 /// specify the mesh containing the mesh containing the
180 /// block-preconditionable solid elements. The dimension of the
181 /// problem must also be specified.
183 {
185 }
186
187 /// Enable mass matrix diagonal scaling in the
188 /// Schur complement approximation
190 {
192 }
193
194 /// Enable mass matrix diagonal scaling in the
195 /// Schur complement approximation
197 {
199 }
200
201 /// Return whether the mass matrix is using diagonal
202 /// scaling or not
204 {
206 }
207
208 /// Function to set a new pressure matrix preconditioner (inexact solver)
209 void set_p_preconditioner(Preconditioner* new_p_preconditioner_pt)
210 {
211 // If the default preconditioner has been used
212 // clean it up now...
214 {
215 delete P_preconditioner_pt;
216 }
217 P_preconditioner_pt = new_p_preconditioner_pt;
219 }
220
221 /// Function to (re-)set pressure matrix preconditioner (inexact
222 /// solver) to SuperLU
224 {
226 {
229 }
230 }
231
232 /// Function to set a new momentum matrix preconditioner (inexact solver)
233 void set_f_preconditioner(Preconditioner* new_f_preconditioner_pt)
234 {
235 // If the default preconditioner has been used
236 // clean it up now...
238 {
239 delete F_preconditioner_pt;
240 }
241 F_preconditioner_pt = new_f_preconditioner_pt;
243 }
244
245 /// Function to (re-)set momentum matrix preconditioner (inexact
246 /// solver) to SuperLU
248 {
250 {
253 }
254 }
255
256
257 /// Enable documentation of time
259 {
260 Doc_time = true;
261 }
262
263 /// Disable documentation of time
265 {
266 Doc_time = false;
267 }
268
269 /// If this function is called then:
270 /// in setup(...) : BFBt is computed.
271 /// in preconditioner_solve(...) : a single matrix vector product with
272 /// BFBt is performed.
274 {
275 Form_BFBt_product = true;
276 }
277
278 /// if this function is called then:
279 /// in setup(...) : the matrices B, F are assembled and stored
280 /// (the default behaviour) .
281 /// in preconditioner_solve(...) : a sequence of matrix vector products
282 /// with B, F, and Bt is performed.
283 /// (Note: in this discussion no scaling was considered but B and Bt
284 /// are replaced with BQ and QBt with scaling)
286 {
287 Form_BFBt_product = false;
288 }
289
290 /// Helper function to delete preconditioner data.
291 void clean_up_memory();
292
293 private:
294 // oomph-lib objects
295 // -----------------
296
297 // Pointers to preconditioner (=inexact solver) objects
298 // -----------------------------------------------------
299 /// Pointer to the 'preconditioner' for the pressure matrix
301
302 /// Pointer to the 'preconditioner' for the F matrix
304
305 /// flag indicating whether the default F preconditioner is used
307
308 /// flag indicating whether the default P preconditioner is used
310
311 /// Control flag is true if the preconditioner has been setup
312 /// (used so we can wipe the data when the preconditioner is
313 /// called again)
315
316 /// Control flag is true if mass matrix diagonal scaling
317 /// is used in the Schur complement approximation
319
320 /// Helper function to assemble the diagonal of the
321 /// mass matrix from the elemental contributions defined in
322 /// PressureBasedSolidEquations<DIM>::get_mass_matrix_diagonal(...).
324
325 /// Boolean indicating whether the momentum system preconditioner
326 /// is a block preconditioner
328
329 /// Set Doc_time to true for outputting results of timings
331
332 /// MatrixVectorProduct operator for F if BFBt is not to be formed.
334
335 /// MatrixVectorProduct operator for QBt if BFBt is not to be formed.
337
338 /// MatrixVectorProduct operator for Bt;
340
341 /// MatrixVectorProduct operator for E (BFBt) if BFBt is to be formed.
343
344 /// indicates whether BFBt should be formed or the component matrices
345 /// should be retained.
346 /// If true then:
347 /// in setup(...) : BFBt is computed.
348 /// in preconditioner_solve(...) : a single matrix vector product with
349 /// BFBt is performed.
350 /// if false then:
351 /// in setup(...) : the matrices B, F are assembled and stored.
352 /// in preconditioner_solve(...) : a sequence of matrix vector products
353 /// with B, F, and Bt is performed.
354 /// (Note: in this discussion no scaling was considered but B and Bt
355 /// are replaced with BQ and QBt with scaling)
357
358 /// the pointer to the mesh of block preconditionable solid
359 /// elements.
361 };
362
363
364 /// ////////////////////////////////////////////////////////////////////////////
365 /// ////////////////////////////////////////////////////////////////////////////
366 /// ////////////////////////////////////////////////////////////////////////////
367
368
369 //============================================================================
370 /// The exact solid preconditioner. This extracts 2x2 blocks
371 /// (corresponding to the displacement/position and pressure unknowns)
372 /// and uses these to
373 /// build a single preconditioner matrix for testing purposes.
374 /// Iterative solvers should converge in a single step if this is used.
375 /// If it doesn't something is wrong in the setup of the block matrices.
376 //=============================================================================
377 template<typename MATRIX>
379 : public BlockPreconditioner<MATRIX>
380 {
381 public:
382 /// Constructor - do nothing
384
385
386 /// Destructor - do nothing
388
389
390 /// Broken copy constructor
393
394
395 /// Broken assignment operator
396 /*void operator=(const PressureBasedSolidExactPreconditioner&) = delete;*/
397
398
399 /// Setup the preconditioner
400 void setup();
401
402 /// Apply preconditioner to r
404
405 protected:
406 /// Preconditioner matrix
407 MATRIX P_matrix;
408 };
409
410} // namespace oomph
411#endif
Block Preconditioner base class. The block structure of the overall problem is determined from the Me...
const Mesh * mesh_pt(const unsigned &i) const
Access to i-th mesh (of the various meshes that contain block preconditionable elements of the same n...
void set_nmesh(const unsigned &n)
Specify the number of meshes required by this block preconditioner. Note: elements in different meshe...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
A general mesh class.
Definition: mesh.h:67
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
PressureBasedSolidExactPreconditioner()
Constructor - do nothing.
~PressureBasedSolidExactPreconditioner()
Destructor - do nothing.
PressureBasedSolidExactPreconditioner(const PressureBasedSolidExactPreconditioner &)=delete
Broken copy constructor.
void setup()
Broken assignment operator.
void preconditioner_solve(const Vector< double > &r, Vector< double > &z)
Apply preconditioner to r.
The least-squares commutator (LSC; formerly BFBT) preconditioner. It uses blocks corresponding to the...
MatrixVectorProduct * QBt_mat_vec_pt
MatrixVectorProduct operator for QBt if BFBt is not to be formed.
void enable_doc_time()
Enable documentation of time.
bool Form_BFBt_product
indicates whether BFBt should be formed or the component matrices should be retained....
CRDoubleMatrix * assemble_mass_matrix_diagonal()
Helper function to assemble the diagonal of the mass matrix from the elemental contributions defined ...
void enable_form_BFBt_product()
If this function is called then: in setup(...) : BFBt is computed. in preconditioner_solve(....
PressureBasedSolidLSCPreconditioner()
Constructor - sets defaults for control flags.
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
MatrixVectorProduct * E_mat_vec_pt
MatrixVectorProduct operator for E (BFBt) if BFBt is to be formed.
void disable_doc_time()
Disable documentation of time.
bool Doc_time
Set Doc_time to true for outputting results of timings.
PressureBasedSolidLSCPreconditioner(const PressureBasedSolidLSCPreconditioner &)=delete
Broken copy constructor.
void set_p_preconditioner(Preconditioner *new_p_preconditioner_pt)
Function to set a new pressure matrix preconditioner (inexact solver)
Mesh * Solid_mesh_pt
the pointer to the mesh of block preconditionable solid elements.
void disable_p_matrix_scaling()
Enable mass matrix diagonal scaling in the Schur complement approximation.
bool F_preconditioner_is_block_preconditioner
Boolean indicating whether the momentum system preconditioner is a block preconditioner.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to Vector r.
void setup()
Broken assignment operator.
bool is_p_matrix_using_scaling() const
Return whether the mass matrix is using diagonal scaling or not.
void set_f_preconditioner(Preconditioner *new_f_preconditioner_pt)
Function to set a new momentum matrix preconditioner (inexact solver)
MatrixVectorProduct * Bt_mat_vec_pt
MatrixVectorProduct operator for Bt;.
bool Using_default_p_preconditioner
flag indicating whether the default P preconditioner is used
MatrixVectorProduct * F_mat_vec_pt
MatrixVectorProduct operator for F if BFBt is not to be formed.
void set_solid_mesh(Mesh *mesh_pt)
specify the mesh containing the mesh containing the block-preconditionable solid elements....
void clean_up_memory()
Helper function to delete preconditioner data.
bool Using_default_f_preconditioner
flag indicating whether the default F preconditioner is used
void set_p_superlu_preconditioner()
Function to (re-)set pressure matrix preconditioner (inexact solver) to SuperLU.
void disable_form_BFBt_product()
if this function is called then: in setup(...) : the matrices B, F are assembled and stored (the defa...
void enable_p_matrix_scaling()
Enable mass matrix diagonal scaling in the Schur complement approximation.
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
bool P_matrix_using_scaling
Control flag is true if mass matrix diagonal scaling is used in the Schur complement approximation.
void set_f_superlu_preconditioner()
Function to (re-)set momentum matrix preconditioner (inexact solver) to SuperLU.
An interface to allow SuperLU to be used as an (exact) Preconditioner.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...