lagrange_enforced_flow_preconditioner.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// Version 1.0; svn revision $LastChangedRevision: 1277 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-12-30 06:46:52 +0000 (Fri, 30 Dec 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//=====================================================================
30 #ifndef OOMPH_LAGRANGE_ENFORCED_FLOW_PRECONDITIONERS_HEADER
31 #define OOMPH_LAGRANGE_ENFORCED_FLOW_PRECONDITIONERS_HEADER
32 
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36 #include <oomph-lib-config.h>
37 #endif
38 
39 // oomphlib headers
40 #include "../generic/matrices.h"
41 #include "../generic/assembly_handler.h"
42 #include "../generic/problem.h"
43 #include "../generic/block_preconditioner.h"
44 #include "../generic/preconditioner.h"
45 #include "../generic/SuperLU_preconditioner.h"
46 #include "../generic/matrix_vector_product.h"
47 #include "../generic/general_purpose_preconditioners.h"
48 #include "../generic/general_purpose_block_preconditioners.h"
49 #ifdef OOMPH_HAS_HYPRE
50 #include "../generic/hypre_solver.h"
51 #endif
52 #ifdef OOMPH_HAS_TRILINOS
53 #include "../generic/trilinos_solver.h"
54 #endif
55 
56 namespace oomph
57 {
58 //==========================================================================
59 /// Namespace for subsidiary preconditioner creation helper functions
60 //==========================================================================
61 namespace Lagrange_Enforced_Flow_Preconditioner_Subsidiary_Operator_Helper
62 {
63  /// \short CG with diagonal preconditioner for W-block subsidiary linear
64  /// systems.
65  extern Preconditioner* get_w_cg_preconditioner();
66 
67  /// \short Hypre Boomer AMG setting for the augmented momentum block
68  /// of a 2D Navier-Stokes problem using the simple form of the viscous
69  /// term (for serial code).
70  extern Preconditioner* boomer_amg_for_2D_momentum_simple_visc();
71 
72  /// \short Hypre Boomer AMG setting for the augmented momentum block
73  /// of a 2D Navier-Stokes problem using the stress divergence form of the
74  /// viscous term (for serial code).
75  extern Preconditioner* boomer_amg_for_2D_momentum_stressdiv_visc();
76 
77  /// \short Hypre Boomer AMG setting for the augmented momentum block
78  /// of a 3D Navier-Stokes problem (for serial code).
79  extern Preconditioner* boomer_amg_for_3D_momentum();
80 
81  /// \short Hypre Boomer AMG setting for the augmented momentum block
82  /// of a 3D Navier-Stokes problem (for serial code).
83  extern Preconditioner* boomer_amg2v22_for_3D_momentum();
84 
85  /// \short Hypre Boomer AMG setting for the 2D Poisson problem
86  /// (for serial code).
87  extern Preconditioner* boomer_amg_for_2D_poisson_problem();
88 
89  /// \short Hypre Boomer AMG setting for the 3D Poisson problem
90  /// (for serial code).
91  extern Preconditioner* boomer_amg_for_3D_poisson_problem();
92 
93 } // namespace
94 
95 
96 //==========================================================================
97 /// \short The preconditioner for the Lagrange multiplier constrained
98 /// Navier-Stokes equations. The velocity components are constrained by
99 /// Lagrange multiplier, which are applied via OOMPH-LIB's FACE elements.
100 ///
101 ///
102 /// The linearised Jacobian takes the block form:
103 ///
104 /// | F_ns | L^T |
105 /// |------------|
106 /// | L | 0 |
107 ///
108 /// where L correspond to the constrained block,
109 /// F_ns is the Navier-Stokes block with the following block structure
110 ///
111 /// | F | G^T |
112 /// |----------|
113 /// | D | 0 |
114 ///
115 /// Here F is the momentum block, G the discrete gradient operator,
116 /// and D the discrete divergence operator. For unstabilised elements,
117 /// we have D = G^T and in much of the literature the divergence matrix is
118 /// denoted by B.
119 ///
120 /// The Lagrange enforced flow preconditioner takes the form:
121 /// | F_aug | |
122 /// |-------|----|
123 /// | | Wd |
124 ///
125 /// where F_aug = F_ns + L^T*inv(Wd)*L is an augmented Navier-Stokes block
126 /// and Wd=(1/Scaling_sigma)*diag(LL^T).
127 ///
128 /// In our implementation of the preconditioner, the linear systems
129 /// associated with the (1,1) block can either be solved "exactly",
130 /// using SuperLU (in its incarnation as an exact preconditioner;
131 /// this is the default) or by any other Preconditioner (inexact solver)
132 /// specified via the access functions
133 ///
134 /// LagrangeEnforcedFlowPreconditioner::set_navier_stokes_preconditioner(...)
135 ///
136 /// Access to the elements is provided via meshes. However, a Vector of
137 /// meshes is taken, each mesh contains a different type of block
138 /// preconditionable element. This allows the (re-)classification of the
139 /// constrained velocity DOF types.
140 ///
141 /// The first mesh in the Vector Mesh_pt must be the 'bulk' mesh.
142 /// The rest are assumed to contain FACEELMENTS.
143 ///
144 /// Thus, the most general block structure (in 3D) is:
145 ///
146 /// 0 1 2 3 4 5 6 7 8 ..x x+0 x+1 x+2 x+3 x+4
147 /// [u v w p] [u v w l1 l2 ...] [u v w l1 l2 ...] ...
148 /// Bulk Surface 1 Surface 2 ...
149 ///
150 /// where the DOF types in [] are the DOF types associated with each mesh.
151 ///
152 /// For example, consider a unit cube domain [0,1]^3 with parallel outflow
153 /// imposed (in mesh 0) and tangential flow imposed (in mesh 1), then there
154 /// are 13 DOF types and our implementation respects the following
155 /// (natural) DOF type order:
156 ///
157 /// bulk mesh 0 mesh 1
158 /// [0 1 2 3] [4 5 6 7 8 ] [9 10 11 12 ]
159 /// [u v w p] [up vp wp Lp1 Lp2] [ut vt wt Lt1]
160 ///
161 /// Via the appropriate mapping, the block_setup(...) function will
162 /// re-order the DOF types into the following block types:
163 ///
164 /// 0 1 2 3 4 5 6 7 8 9 10 11 12 <- Block type
165 /// 0 4 9 1 5 10 2 6 11 3 7 8 12 <- DOF type
166 /// [u up ut v vp vt w wp wt ] [p] [Lp1 Lp2 Lt1]
167 ///
168 //==========================================================================
170 : public BlockPreconditioner<CRDoubleMatrix>
171 {
172 public:
173 
174  /// \short This preconditioner includes the option to use subsidiary
175  /// operators other than SuperLUPreconditioner for this problem.
176  /// This is the typedef of a function that should return an instance
177  /// of a subsidiary preconditioning operator. This preconditioner is
178  /// responsible for the destruction of the subsidiary preconditioners.
179  typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
180 
181  /// Constructor - initialise variables.
184  {
185  // The null pointer.
187 
188  // By default, the linear systems associated with the diagonal blocks
189  // are solved "exactly" using SuperLU (in its incarnation as an exact
190  // preconditioner. This is not a block preconditioner.
192 
193  // Flag to indicate to use SuperLU or not.
195 
196  // Empty vector of meshes and set the number of meshes to zero.
197  My_mesh_pt.resize(0,0);
198  My_nmesh = 0;
199 
200  // The number of DOF types within the meshes.
201  My_ndof_types_in_mesh.resize(0,0);
202 
203  // Initialise other variables.
205  Scaling_sigma = 0.0;
207  N_fluid_doftypes = 0;
210  } // constructor
211 
212  /// \short Destructor
214  {
215  this->clean_up_memory();
216  }
217 
218  /// \short Broken copy constructor
221  {
222  BrokenCopy::broken_copy("LagrangeEnforcedFlowPreconditioner");
223  }
224 
225  /// \short Broken assignment operator
227  {
228  BrokenCopy::broken_assign(" LagrangeEnforcedFlowPreconditioner");
229  }
230 
231  /// \short Setup method for the LagrangeEnforcedFlowPreconditioner.
232  void setup();
233 
234  /// \short Apply the preconditioner.
235  /// r is the residual (rhs), z will contain the solution.
237 
238  /// \short Set the meshes,
239  /// the first mesh in the vector must be the bulk mesh.
240  void set_meshes(const Vector<Mesh*> &mesh_pt);
241 
242  /// \short Set flag to use the infinite norm of the Navier-Stokes F matrix
243  /// as the scaling sigma. This is the default behaviour. Note: the norm of
244  /// the NS F matrix positive, however, we actually use the negative of
245  /// the norm. This is because the underlying Navier-Stokes Jacobian is
246  /// multiplied by -1. Ask Andrew/Matthias for more detail.
248  {
250  }
251 
252  /// \short Access function to set the scaling sigma.
253  /// Note: this also sets the flag to use the infinite norm of
254  /// the Navier-Stokes F matrix as the scaling sigma to false.
255  /// Warning is given if trying to set scaling sigma to be equal to
256  /// or greater than zero.
257  void set_scaling_sigma(const double& scaling_sigma)
258  {
259  // Check if scaling sigma is zero or positive.
260 #ifdef PARANOID
261  if(scaling_sigma == 0.0)
262  {
263  std::ostringstream warning_stream;
264  warning_stream << "WARNING: \n"
265  << "Setting scaling_sigma = 0.0 may cause values.\n";
266  OomphLibWarning(warning_stream.str(),
267  OOMPH_CURRENT_FUNCTION,
268  OOMPH_EXCEPTION_LOCATION);
269  }
270  if(scaling_sigma > 0.0)
271  {
272  std::ostringstream warning_stream;
273  warning_stream << "WARNING: " << std::endl
274  << "The scaling (scaling_sigma) is positive: "
275  << Scaling_sigma << "\n"
276  << "Performance may be degraded.\n";
277  OomphLibWarning(warning_stream.str(),
278  OOMPH_CURRENT_FUNCTION,
279  OOMPH_EXCEPTION_LOCATION);
280  }
281 #endif
282 
285  }
286 
287  /// \short Read (const) function to get the scaling sigma.
288  double scaling_sigma() const
289  {
290  return Scaling_sigma;
291  }
292 
293  /// \short Set a new Navier-Stokes matrix preconditioner
294  /// (inexact solver)
296  Preconditioner* new_ns_preconditioner_pt = 0);
297 
298  ///\short Set Navier-Stokes matrix preconditioner (inexact
299  /// solver) to SuperLU
301  {
303  {
307  }
308  }
309 
310  /// \short Clears the memory.
311  void clean_up_memory();
312 
313 private:
314 
315  /// \short Control flag is true if the preconditioner has been setup
316  /// (used so we can wipe the data when the preconditioner is
317  /// called again)
319 
320  /// \short Scaling for the augmentation: Scaling_sigma*(LL^T)
322 
323  /// \short Flag to indicate if we want to use the infinite norm of the
324  /// Navier-Stokes momentum block for the scaling sigma.
326 
327  /// \short Inverse W values
329 
330  /// \short Pointer to the 'preconditioner' for the Navier-Stokes block
332 
333  /// \short Flag to indicate if the preconditioner for the Navier-Stokes
334  /// block is a block preconditioner or not.
336 
337  /// \short Flag to indicate whether the default NS preconditioner is used
339 
340  /// \short Storage for the meshes. In our implementation, the first mesh
341  /// must always be the Navier-Stokes (bulk) mesh, followed by surface
342  /// meshes.
344 
345  /// \short The number of DOF types in each mesh. This is used create
346  /// various lookup lists.
348 
349  /// \short The number of meshes. This is used to create various lookup
350  /// lists.
351  unsigned My_nmesh;
352 
353  /// \short The number of Lagrange multiplier DOF types.
355 
356  /// \short The number of fluid DOF types (including pressure).
358 
359  /// \short The number of velocity DOF types.
361 
362 }; // end of LagrangeEnforcedFlowPreconditioner class
363 
364 }
365 #endif
void set_meshes(const Vector< Mesh * > &mesh_pt)
Set the meshes, the first mesh in the vector must be the bulk mesh.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Vector< unsigned > My_ndof_types_in_mesh
The number of DOF types in each mesh. This is used create various lookup lists.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the preconditioner. r is the residual (rhs), z will contain the solution.
void set_navier_stokes_preconditioner(Preconditioner *new_ns_preconditioner_pt=0)
Set a new Navier-Stokes matrix preconditioner (inexact solver)
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...
double Scaling_sigma
Scaling for the augmentation: Scaling_sigma*(LL^T)
void set_scaling_sigma(const double &scaling_sigma)
Access function to set the scaling sigma. Note: this also sets the flag to use the infinite norm of t...
Preconditioner * boomer_amg_for_2D_poisson_problem()
Hypre Boomer AMG setting for the 2D Poisson problem (for serial code).
bool Navier_stokes_preconditioner_is_block_preconditioner
Flag to indicate if the preconditioner for the Navier-Stokes block is a block preconditioner or not...
unsigned My_nmesh
The number of meshes. This is used to create various lookup lists.
void set_superlu_for_navier_stokes_preconditioner()
Set Navier-Stokes matrix preconditioner (inexact solver) to SuperLU.
void setup()
Setup method for the LagrangeEnforcedFlowPreconditioner.
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...
Preconditioner * boomer_amg_for_3D_poisson_problem()
Hypre Boomer AMG setting for the 3D Poisson problem (for serial code).
Preconditioner * get_w_cg_preconditioner()
CG with diagonal preconditioner for W-block subsidiary linear systems.
Preconditioner * Navier_stokes_preconditioner_pt
Pointer to the 'preconditioner' for the Navier-Stokes block.
unsigned N_lagrange_doftypes
The number of Lagrange multiplier DOF types.
bool Using_superlu_ns_preconditioner
Flag to indicate whether the default NS preconditioner is used.
Vector< Vector< double > > Inv_w_diag_values
Inverse W values.
bool Use_norm_f_for_scaling_sigma
Flag to indicate if we want to use the infinite norm of the Navier-Stokes momentum block for the scal...
double scaling_sigma() const
Read (const) function to get the scaling sigma.
An interface to allow SuperLU to be used as an (exact) Preconditioner.
unsigned N_fluid_doftypes
The number of fluid DOF types (including pressure).
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
Vector< Mesh * > My_mesh_pt
Storage for the meshes. In our implementation, the first mesh must always be the Navier-Stokes (bulk)...
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
unsigned N_velocity_doftypes
The number of velocity DOF types.
void operator=(const LagrangeEnforcedFlowPreconditioner &)
Broken assignment operator.
Preconditioner * boomer_amg_for_2D_momentum_simple_visc()
Hypre Boomer AMG setting for the augmented momentum block of a 2D Navier-Stokes problem using the sim...
void use_norm_f_for_scaling_sigma()
Set flag to use the infinite norm of the Navier-Stokes F matrix as the scaling sigma. This is the default behaviour. Note: the norm of the NS F matrix positive, however, we actually use the negative of the norm. This is because the underlying Navier-Stokes Jacobian is multiplied by -1. Ask Andrew/Matthias for more detail.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
LagrangeEnforcedFlowPreconditioner()
Constructor - initialise variables.
The preconditioner for the Lagrange multiplier constrained Navier-Stokes equations. The velocity components are constrained by Lagrange multiplier, which are applied via OOMPH-LIB's FACE elements.
Preconditioner * boomer_amg2v22_for_3D_momentum()
Hypre Boomer AMG setting for the augmented momentum block of a 3D Navier-Stokes problem (for serial c...
Preconditioner * boomer_amg_for_3D_momentum()
Hypre Boomer AMG setting for the augmented momentum block of a 3D Navier-Stokes problem (for serial c...
Preconditioner * boomer_amg_for_2D_momentum_stressdiv_visc()
Hypre Boomer AMG setting for the augmented momentum block of a 2D Navier-Stokes problem using the str...