lagrange_enforced_flow_preconditioner.cc
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$
7 //LIC//
8 //LIC// $LastChangedDate$
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//=====================================================================
31 
32 namespace oomph
33 {
34 //==========================================================================
35 /// Namespace for subsidiary preconditioner creation helper functions
36 //==========================================================================
37 namespace Lagrange_Enforced_Flow_Preconditioner_Subsidiary_Operator_Helper
38 {
39  /// \short CG with diagonal preconditioner for W-block subsidiary linear
40  /// systems.
42  {
43 #ifdef OOMPH_HAS_TRILINOS
47  <TrilinosAztecOOSolver,MatrixBasedDiagPreconditioner>;
48 
49  // Note: This makes CG a proper "inner iteration" for
50  // which GMRES (may) no longer converge. We should really
51  // use FGMRES or GMRESR for this. However, here the solver
52  // is so good that it'll converge very quickly anyway
53  // so there isn't much to be gained by limiting the number
54  // of iterations...
55  prec_pt->max_iter() = 4;
56  prec_pt->solver_pt()->solver_type() = TrilinosAztecOOSolver::CG;
57  prec_pt->solver_pt()->disable_doc_time();
58  return prec_pt;
59 #else
60  std::ostringstream err_msg;
61  err_msg << "Inner CG preconditioner is unavailable.\n"
62  << "Please install Trilinos.\n";
63  throw OomphLibError(err_msg.str(),
64  OOMPH_CURRENT_FUNCTION,
65  OOMPH_EXCEPTION_LOCATION);
66 #endif
67  } // function get_w_cg_preconditioner
68 
69  /// \short Hypre Boomer AMG setting for the augmented momentum block
70  /// of a 2D Navier-Stokes problem using the simple form of the viscous
71  /// term (for serial code).
73  {
74 #ifdef OOMPH_HAS_HYPRE
75  // Create a new HyprePreconditioner
76  HyprePreconditioner* hypre_preconditioner_pt = new HyprePreconditioner;
77 
78  // Coarsening strategy
79  // 1 = classical RS with no boundary treatment (not recommended in
80  // parallel)
81  hypre_preconditioner_pt->amg_coarsening() = 1;
82 
83  // Strength of dependence = 0.25
84  hypre_preconditioner_pt->amg_strength() = 0.25;
85 
86 
87  // Set the smoothers
88  // 1 = Gauss-Seidel, sequential (very slow in parallel!)
89  hypre_preconditioner_pt->amg_simple_smoother() = 1;
90 
91  // Set smoother damping (not required, so set to -1)
92  hypre_preconditioner_pt->amg_damping() = -1;
93 
94 
95  // Set number of cycles to 1xV(2,2)
96  hypre_preconditioner_pt->set_amg_iterations(1);
97  hypre_preconditioner_pt->amg_smoother_iterations()=2;
98 
99  return hypre_preconditioner_pt;
100 #else
101  std::ostringstream err_msg;
102  err_msg << "hypre preconditioner is not available.\n"
103  << "Please install Hypre.\n";
104  throw OomphLibError(err_msg.str(),
105  OOMPH_CURRENT_FUNCTION,
106  OOMPH_EXCEPTION_LOCATION);
107 #endif
108  } // function boomer_amg_for_2D_momentum_simple_visc
109 
110  /// \short Hypre Boomer AMG setting for the augmented momentum block
111  /// of a 2D Navier-Stokes problem using the stress divergence form of the
112  /// viscous term (for serial code).
114  {
115 #ifdef OOMPH_HAS_HYPRE
116  // Create a new HyprePreconditioner
117  HyprePreconditioner* hypre_preconditioner_pt = new HyprePreconditioner;
118 
119  // Coarsening strategy
120  // 1 = classical RS with no boundary treatment (not recommended in
121  // parallel)
122  hypre_preconditioner_pt->amg_coarsening() = 1;
123 
124  // Strength of dependence = 0.668
125  hypre_preconditioner_pt->amg_strength() = 0.668;
126 
127 
128  // Set the smoothers
129  // 1 = Gauss-Seidel, sequential (very slow in parallel!)
130  hypre_preconditioner_pt->amg_simple_smoother() = 1;
131 
132  // Set smoother damping (not required, so set to -1)
133  hypre_preconditioner_pt->amg_damping() = -1;
134 
135 
136  // Set number of cycles to 1xV(2,2)
137  hypre_preconditioner_pt->set_amg_iterations(1);
138  hypre_preconditioner_pt->amg_smoother_iterations()=2;
139 
140  return hypre_preconditioner_pt;
141 #else
142  std::ostringstream err_msg;
143  err_msg << "hypre preconditioner is not available.\n"
144  << "Please install Hypre.\n";
145  throw OomphLibError(err_msg.str(),
146  OOMPH_CURRENT_FUNCTION,
147  OOMPH_EXCEPTION_LOCATION);
148 #endif
149  } // function boomer_amg_for_2D_momentum_stressdiv_visc
150 
151  /// \short Hypre Boomer AMG setting for the augmented momentum block
152  /// of a 3D Navier-Stokes problem (for serial code).
154  {
155 #ifdef OOMPH_HAS_HYPRE
156  // Create a new HyprePreconditioner
157  HyprePreconditioner* hypre_preconditioner_pt = new HyprePreconditioner;
158 
159  // Coarsening strategy
160  // 1 = classical RS with no boundary treatment (not recommended in
161  // parallel)
162  hypre_preconditioner_pt->amg_coarsening() = 1;
163 
164  // Strength of dependence = 0.668
165  hypre_preconditioner_pt->amg_strength() = 0.8;
166 
167 
168  // Set the smoothers
169  // 1 = Gauss-Seidel, sequential (very slow in parallel!)
170  hypre_preconditioner_pt->amg_simple_smoother() = 1;
171 
172  // Set smoother damping (not required, so set to -1)
173  hypre_preconditioner_pt->amg_damping() = -1;
174 
175 
176  // Set number of cycles to 1xV(2,2)
177  hypre_preconditioner_pt->set_amg_iterations(1);
178  hypre_preconditioner_pt->amg_smoother_iterations()=2;
179 
180  return hypre_preconditioner_pt;
181 #else
182  std::ostringstream err_msg;
183  err_msg << "hypre preconditioner is not available.\n"
184  << "Please install Hypre.\n";
185  throw OomphLibError(err_msg.str(),
186  OOMPH_CURRENT_FUNCTION,
187  OOMPH_EXCEPTION_LOCATION);
188 #endif
189  } // function boomer_amg_for_3D_momentum
190 
191  /// \short Hypre Boomer AMG setting for the augmented momentum block
192  /// of a 3D Navier-Stokes problem (for serial code).
194  {
195 #ifdef OOMPH_HAS_HYPRE
196  // Create a new HyprePreconditioner
197  HyprePreconditioner* hypre_preconditioner_pt = new HyprePreconditioner;
198 
199  // Coarsening strategy
200  // 1 = classical RS with no boundary treatment (not recommended in
201  // parallel)
202  hypre_preconditioner_pt->amg_coarsening() = 1;
203 
204  // Strength of dependence = 0.668
205  hypre_preconditioner_pt->amg_strength() = 0.8;
206 
207 
208  // Set the smoothers
209  // 1 = Gauss-Seidel, sequential (very slow in parallel!)
210  hypre_preconditioner_pt->amg_simple_smoother() = 1;
211 
212  // Set smoother damping (not required, so set to -1)
213  hypre_preconditioner_pt->amg_damping() = -1;
214 
215 
216  // Set number of cycles to 1xV(2,2)
217  hypre_preconditioner_pt->set_amg_iterations(2);
218  hypre_preconditioner_pt->amg_smoother_iterations()=2;
219 
220  return hypre_preconditioner_pt;
221 #else
222  std::ostringstream err_msg;
223  err_msg << "hypre preconditioner is not available.\n"
224  << "Please install Hypre.\n";
225  throw OomphLibError(err_msg.str(),
226  OOMPH_CURRENT_FUNCTION,
227  OOMPH_EXCEPTION_LOCATION);
228 #endif
229  } // function boomer_amg_for_3D_momentum
230 
231 
232 
233  /// \short Hypre Boomer AMG setting for the 2D Poisson problem
234  /// (for serial code).
236  {
237 #ifdef OOMPH_HAS_HYPRE
238  // Create a new HyprePreconditioner
239  HyprePreconditioner* hypre_preconditioner_pt = new HyprePreconditioner;
240 
241  // Coarsening strategy
242  // 1 = classical RS with no boundary treatment (not recommended in
243  // parallel)
244  hypre_preconditioner_pt->amg_coarsening() = 1;
245 
246  // Strength of dependence = 0.25
247  hypre_preconditioner_pt->amg_strength() = 0.25;
248 
249 
250  // Set the smoothers
251  // 0 = Jacobi
252  hypre_preconditioner_pt->amg_simple_smoother() = 0;
253 
254  // Set Jacobi damping = 2/3
255  hypre_preconditioner_pt->amg_damping() = 0.668;
256 
257 
258  // Set number of cycles to 1xV(2,2)
259  hypre_preconditioner_pt->set_amg_iterations(2);
260  hypre_preconditioner_pt->amg_smoother_iterations()=1;
261 
262  return hypre_preconditioner_pt;
263 #else
264  std::ostringstream err_msg;
265  err_msg << "hypre preconditioner is not available.\n"
266  << "Please install Hypre.\n";
267  throw OomphLibError(err_msg.str(),
268  OOMPH_CURRENT_FUNCTION,
269  OOMPH_EXCEPTION_LOCATION);
270 #endif
271  } // function boomer_amg_for_2D_poisson_problem
272 
273  /// \short Hypre Boomer AMG setting for the 3D Poisson problem
274  /// (for serial code).
276  {
277 #ifdef OOMPH_HAS_HYPRE
278  // Create a new HyprePreconditioner
279  HyprePreconditioner* hypre_preconditioner_pt = new HyprePreconditioner;
280 
281  // Coarsening strategy
282  // 1 = classical RS with no boundary treatment (not recommended in
283  // parallel)
284  hypre_preconditioner_pt->amg_coarsening() = 1;
285 
286  // Strength of dependence = 0.7
287  hypre_preconditioner_pt->amg_strength() = 0.7;
288 
289 
290  // Set the smoothers
291  // 0 = Jacobi
292  hypre_preconditioner_pt->amg_simple_smoother() = 0;
293 
294  // Set smoother damping = 2/3
295  hypre_preconditioner_pt->amg_damping() = 0.668;
296 
297 
298  // Set number of cycles to 2xV(1,1)
299  hypre_preconditioner_pt->set_amg_iterations(2);
300  hypre_preconditioner_pt->amg_smoother_iterations()=1;
301 
302  return hypre_preconditioner_pt;
303 #else
304  std::ostringstream err_msg;
305  err_msg << "hypre preconditioner is not available.\n"
306  << "Please install Hypre.\n";
307  throw OomphLibError(err_msg.str(),
308  OOMPH_CURRENT_FUNCTION,
309  OOMPH_EXCEPTION_LOCATION);
310 #endif
311  } // function boomer_amg_for_3D_poisson_problem
312 
313 } // namespace
314 
315 ////////////////////////////////////////////////////////////////////////////
316 ////////////////////////////////////////////////////////////////////////////
317 ////////////////////////////////////////////////////////////////////////////
318 
319 /// \short Apply the preconditioner.
320 /// r is the residual (rhs), z will contain the solution.
322  const DoubleVector& r, DoubleVector& z)
323 {
324 #ifdef PARANOID
326  {
327  std::ostringstream error_message;
328  error_message << "setup() must be called before using "
329  << "preconditioner_solve()";
330  throw OomphLibError(
331  error_message.str(),
332  OOMPH_CURRENT_FUNCTION,
333  OOMPH_EXCEPTION_LOCATION);
334  }
335  if (z.built())
336  {
337  if (z.nrow() != r.nrow())
338  {
339  std::ostringstream error_message;
340  error_message << "The vectors z and r must have the same number of "
341  << "of global rows";
342  throw OomphLibError(
343  error_message.str(),
344  OOMPH_CURRENT_FUNCTION,
345  OOMPH_EXCEPTION_LOCATION);
346  }
347  }
348 #endif
349 
350  // if z is not setup then give it the same distribution
351  if (!z.distribution_pt()->built())
352  {
353  z.build(r.distribution_pt(),0.0);
354  }
355 
356  // Working vectors.
357  DoubleVector temp_vec;
358  DoubleVector another_temp_vec;
359  DoubleVector yet_another_temp_vec;
360 
361 
362  // -----------------------------------------------------------------------
363  // Step 1 - apply approximate W block inverse to Lagrange multiplier
364  // unknowns
365  // -----------------------------------------------------------------------
366  // For each subsystem associated with each Lagrange multiplier, we loop
367  // through and:
368  // 1) extract the block entries from r
369  // 2) apply the inverse
370  // 3) return the entries to z.
371  for(unsigned l_i = 0; l_i < N_lagrange_doftypes; l_i++)
372  {
373  // The Lagrange multiplier block type.
374  const unsigned l_ii = N_fluid_doftypes + l_i;
375 
376  // Extract the block
377  this->get_block_vector(l_ii,r,temp_vec);
378 
379  // Apply the inverse.
380  const unsigned vec_nrow_local = temp_vec.nrow_local();
381  double* vec_values_pt = temp_vec.values_pt();
382  for (unsigned i = 0; i < vec_nrow_local; i++)
383  {
384  vec_values_pt[i]
385  = vec_values_pt[i]*Inv_w_diag_values[l_i][i];
386  } // for
387 
388  // Return the unknowns
389  this->return_block_vector(l_ii,temp_vec,z);
390 
391  // Clear vectors.
392  temp_vec.clear();
393  } // for
394 
395  // -----------------------------------------------------------------------
396  // Step 2 - apply the augmented Navier-Stokes matrix inverse to the
397  // velocity and pressure unknowns
398  // -----------------------------------------------------------------------
399 
400  // At this point, all vectors are cleared.
402  {
403  // Which block types corresponds to the fluid block types.
404  Vector<unsigned> fluid_block_indices(N_fluid_doftypes,0);
405  for (unsigned b = 0; b < N_fluid_doftypes; b++)
406  {
407  fluid_block_indices[b] = b;
408  }
409 
410  this->get_concatenated_block_vector(fluid_block_indices,r,temp_vec);
411 
412  // temp_vec contains the (concatenated) fluid rhs.
414  ->preconditioner_solve(temp_vec,another_temp_vec);
415 
416  temp_vec.clear();
417 
418  // Return it to the unknowns.
419  this->return_concatenated_block_vector(fluid_block_indices,
420  another_temp_vec,z);
421 
422  another_temp_vec.clear();
423  }
424  else
425  {
426  // The Navier-Stokes preconditioner is a block preconditioner.
427  // Thus is handles all of the block vector extraction and returns.
429  }
430 } // end of preconditioner_solve
431 
432 /// \short Set the meshes,
433 /// the first mesh in the vector must be the bulk mesh.
435  const Vector<Mesh*> &mesh_pt)
436 {
437  // There should be at least two meshes passed to this preconditioner.
438  const unsigned nmesh = mesh_pt.size();
439 
440 #ifdef PARANOID
441  if(nmesh < 2)
442  {
443  std::ostringstream err_msg;
444  err_msg << "There should be at least two meshes.\n";
445  throw OomphLibError(err_msg.str(),
446  OOMPH_CURRENT_FUNCTION,
447  OOMPH_EXCEPTION_LOCATION);
448  }
449 
450  // Check that all pointers are not null
451  for(unsigned mesh_i = 0; mesh_i < nmesh; mesh_i++)
452  {
453  if (mesh_pt[mesh_i]==0)
454  {
455  std::ostringstream err_msg;
456  err_msg << "The pointer mesh_pt[" << mesh_i << "] is null.\n";
457  throw OomphLibError(err_msg.str(),
458  OOMPH_CURRENT_FUNCTION,
459  OOMPH_EXCEPTION_LOCATION);
460  }
461  }
462 
463  // We assume that the first mesh is the Navier-Stokes "bulk" mesh.
464  // To check this, the elemental dimension must be the same as the
465  // nodal (spatial) dimension.
466  //
467  // We store the elemental dimension i.e. the number of local coordinates
468  // required to parametrise its geometry.
469  const unsigned elemental_dim = mesh_pt[0]->elemental_dimension();
470 
471  // The dimension of the nodes in the first element in the (supposedly)
472  // bulk mesh.
473  const unsigned nodal_dim = mesh_pt[0]->nodal_dimension();
474 
475  // Check if the first mesh is the "bulk" mesh.
476  // Here we assume only one mesh contains "bulk" elements.
477  // All subsequent meshes contain block preconditionable elements which
478  // re-classify the bulk velocity DOFs to constrained velocity DOFs.
479  if (elemental_dim != nodal_dim)
480  {
481  std::ostringstream err_msg;
482  err_msg << "In the first mesh, the elements have elemental dimension "
483  << "of " << elemental_dim << ",\n"
484  << "with a nodal dimension of " << nodal_dim << ".\n"
485  << "The first mesh does not contain 'bulk' elements.\n"
486  << "Please re-order your mesh_pt vector.\n";
487 
488  throw OomphLibError(err_msg.str(),
489  OOMPH_CURRENT_FUNCTION,
490  OOMPH_EXCEPTION_LOCATION);
491  }
492 #endif
493 
494  // Set the number of meshes
495  this->set_nmesh(nmesh);
496 
497  // Set the meshes
498  for(unsigned mesh_i = 0; mesh_i < nmesh; mesh_i++)
499  {
500  this->set_mesh(mesh_i,mesh_pt[mesh_i]);
501  }
502 
503  // We also store the meshes and number of meshes locally in this class.
504  // This may seem slightly redundant, since we always have all the meshes
505  // stored in the upper most master block preconditioner.
506  // But at the moment there is no mapping/look up scheme between
507  // master and subsidiary block preconditioners for meshes.
508  //
509  // If this is a subsidiary block preconditioner, we don't know which of
510  // the master's meshes belong to us. We need this information to set up
511  // look up lists in the function setup(...).
512  // Thus we store them local to this class.
514  My_nmesh = nmesh;
515 } // function set_meshes
516 
517 
518 //==========================================================================
519 /// Setup the Lagrange enforced flow preconditioner. This
520 /// extracts blocks corresponding to the velocity and Lagrange multiplier
521 /// unknowns, creates the matrices actually needed in the application of the
522 /// preconditioner and deletes what can be deleted... Note that
523 /// this preconditioner needs a CRDoubleMatrix.
524 //==========================================================================
526 {
527  // clean
528  this->clean_up_memory();
529 
530 #ifdef PARANOID
531  // Paranoid check that meshes have been set. In this preconditioner, we
532  // always have to set meshes.
533  if(My_nmesh == 0)
534  {
535  std::ostringstream err_msg;
536  err_msg << "There are no meshes set. Please call set_meshes(...)\n"
537  << "with at least two mesh.";
538  throw OomphLibError(err_msg.str(),
539  OOMPH_CURRENT_FUNCTION,
540  OOMPH_EXCEPTION_LOCATION);
541  }
542 #endif
543 
544  // -----------------------------------------------------------------------
545  // Step 1 - Construct the dof_to_block_map vector.
546  // -----------------------------------------------------------------------
547  // Assumption: The first mesh is always the "bulk" mesh
548  // (Navier-Stokes mesh), which contains block preconditionable
549  // Navier-Stokes elements. Thus the bulk elements classify the velocity
550  // and pressure degrees of freedom (ndof_types = 3(4) in 2(3)D).
551  // All subsequent meshes contain the constrained velocity DOF types,
552  // then the Lagrange multiplier DOF types.
553  //
554  // Thus, a general ordering of DOF types (in 3D) follows the ordering:
555  //
556  // 0 1 2 3 4 5 6 7 8 ..x x+0 x+1 x+2 x+3 x+4
557  // [u v w p] [u v w l1 l2 ...] [u v w l1 l2 ...] ...
558  //
559  // where the square brackets [] represent the DOF types in each mesh.
560  //
561  // Example:
562  // Consider the case of imposing parallel outflow (3 constrained velocity
563  // DOF types and 2 Lagrange multiplier DOF types) and tangential flow (3
564  // constrained velocity DOF types and 1 Lagrange multiplier DOF type)
565  // along two different boundaries in 3D. The resulting natural ordering of
566  // the DOF types is:
567  //
568  // [0 1 2 3] [4 5 6 7 8 ] [9 10 11 12 ]
569  // [u v w p] [up vp wp Lp1 Lp2] [ut vt wt Lt1]
570  //
571  //
572  // In our implementation, the desired block structure is:
573  // | u v w | up vp wp | ut vt wt | p | Lp1 Lp2 Lt1 |
574  //
575  // The dof_to_block_map should have the following construction:
576  //
577  // dof_to_block_map[$dof_number] = $block_number
578  //
579  // Thus, the dof_to_block_map for the example above should be:
580  //
581  // To achieve this, we use the dof_to_block_map:
582  // dof_to_block_map = [0 1 2 9 3 4 5 10 11 6 7 8 12]
583  //
584  // To generalise the construction of the dof_to_block_map vector
585  // (to work for any number of meshes), we first require some auxiliary
586  // variables to aid us in this endeavour.
587  // -----------------------------------------------------------------------
588 
589  // Set up the My_ndof_types_in_mesh vector.
590  // If this is already constructed, we reuse it instead.
591  if(My_ndof_types_in_mesh.size() == 0)
592  {
593  for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
594  {
595  My_ndof_types_in_mesh.push_back(My_mesh_pt[mesh_i]->ndof_types());
596  }
597  }
598 
599  // Get the spatial dimension of the problem.
600  unsigned spatial_dim = My_mesh_pt[0]->nodal_dimension();
601 
602  // Get the number of DOF types.
603  unsigned n_dof_types = ndof_types();
604 
605 #ifdef PARANOID
606  // We cannot check which DOF types are "correct" in the sense that there
607  // is a distinction between bulk and constrained velocity DOF tyoes.
608  // But we can at least check if the ndof_types matches the total number
609  // of DOF types from the meshes.
610  //
611  // This check is not necessary for the upper most master block
612  // preconditioner since the ndof_types() is calculated by looping
613  // through the meshes!
614  //
615  // This check is useful if this is a subsidiary block preconditioner and
616  // incorrect DOF type merging has taken place.
618  {
619  unsigned tmp_ndof_types = 0;
620  for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
621  {
622  tmp_ndof_types += My_ndof_types_in_mesh[mesh_i];
623  }
624 
625  if(tmp_ndof_types != n_dof_types)
626  {
627  std::ostringstream err_msg;
628  err_msg << "The number of DOF types are incorrect.\n"
629  << "The total DOF types from the meshes is: "
630  << tmp_ndof_types << ".\n"
631  << "The number of DOF types from "
632  << "BlockPreconditioner::ndof_types() is " << n_dof_types <<".\n";
633  throw OomphLibError(err_msg.str(),
634  OOMPH_CURRENT_FUNCTION,
635  OOMPH_EXCEPTION_LOCATION);
636  }
637  }
638 #endif
639 
640  // The number of velocity DOF types: We assume that all meshes classify
641  // at least some of the velocity DOF types (bulk/constrained), thus the
642  // total number of velocity DOF types is the spatial dimension multiplied
643  // by the number of meshes.
644  N_velocity_doftypes = My_nmesh*spatial_dim;
645 
646  // Fluid has +1 for the pressure.
648 
649  // The rest are Lagrange multiplier DOF types.
650  N_lagrange_doftypes = n_dof_types - N_fluid_doftypes;
651 
652  //////////////////////////////////////////////////////////////
653  // Now construct the DOF to block map for block_setup() //
654  //////////////////////////////////////////////////////////////
655 
656  // Now that we have
657  //
658  // (1) spatial dimension of the problem and
659  // (2) the number of DOF types in each of the meshes.
660  //
661  // We observe that the problem dimension is 3 and
662  // My_ndof_type_in_mesh = [4, 5, 4].
663  //
664  // With these information we can construct the desired block structure:
665  // | u v w | up vp wp | ut vt wt | p | Lp1 Lp2 Lt1 |
666  //
667  // The block structure is determined by the vector dof_to_block_map we
668  // give to the function block_setup(...).
669 
670  // This preconditioner permutes the DOF numbers to get the block numbers.
671  // I.e. nblock_types = ndof_types, but they are re-ordered
672  // so that we have (in order):
673  // 1) bulk velocity DOF types
674  // 2) constrained velocity DOF types
675  // 3) pressure DOF type
676  // 4) Lagrange multiplier DOF types
677  //
678  // Recall that the natural ordering of the DOF types are ordered first by
679  // their meshes, then the elemental DOF type as described by the
680  // function get_dof_numbers_for_unknowns().
681  //
682  // Also recall that (for this problem), we assume/require that every mesh
683  // (re-)classify at least some of the velocity DOF types, furthermore,
684  // the velocity DOF type classification comes before the
685  // pressure / Lagrange multiplier DOF types.
686  //
687  // Consider the same example as shown previously, repeated here for your
688  // convenience:
689  //
690  // Natural DOF type ordering:
691  // 0 1 2 3 4 5 6 7 8 9 10 11 12 <- DOF number.
692  // [u v w p] [up vp wp Lp1 Lp2] [ut vt wt Lt1] <- DOF type.
693  //
694  // Desired block structure:
695  // 0 1 2 3 4 5 6 7 8 9 10 11 12 <- block number
696  // [u v w | up vp wp | ut vt wt ] [p | Lp1 Lp2 Lt1] <- DOF type
697  //
698  // dof_to_block_map[$dof_number] = $block_number
699  //
700  // Required dof_to_block_map:
701  // 0 1 2 9 3 4 5 10 11 6 7 8 12 <- block number
702  // [u v w p up vp wp Lp1 Lp2 ut vt wt Lt1] <- DOF type
703  //
704  // Consider the 4th entry of the dof_to_block_map, this represents the
705  // pressure DOF type, from the desired block structure we see this takes
706  // the block number 9.
707  //
708  // One way to generalise the construction of the dof_to_block_map is to
709  // lump the first $spatial_dimension number of DOF types
710  // from each mesh together, then lump the remaining DOF types together.
711  //
712  // Notice that the values in the velocity DOF types of the
713  // dof_to_block_map vector increases sequentially, from 0 to 8. The values
714  // of the Lagrange multiplier DOF types follow the same pattern, from 9 to
715  // 12. We follow this construction.
716 
717  // Storage for the dof_to_block_map vector.
718  Vector<unsigned> dof_to_block_map(n_dof_types,0);
719 
720  // Index for the dof_to_block_map vector.
721  unsigned temp_index = 0;
722 
723  // Value for the velocity DOF type.
724  unsigned velocity_number = 0;
725 
726  // Value for the pressure/Lagrange multiplier DOF type.
727  unsigned pres_lgr_number = N_velocity_doftypes;
728 
729  // Loop through the number of meshes.
730  for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
731  {
732  // Fill in the velocity DOF types.
733  for (unsigned dim_i = 0; dim_i < spatial_dim; dim_i++)
734  {
735  dof_to_block_map[temp_index++] = velocity_number++;
736  } // for
737 
738  // Loop through the pressure/Lagrange multiplier DOF types.
739  unsigned ndof_type_in_mesh_i = My_ndof_types_in_mesh[mesh_i];
740  for (unsigned doftype_i = spatial_dim;
741  doftype_i < ndof_type_in_mesh_i; doftype_i++)
742  {
743  dof_to_block_map[temp_index++] = pres_lgr_number++;
744  } // for
745  } // for
746 
747  // Call the block setup
748  this->block_setup(dof_to_block_map);
749 
750 
751  // -----------------------------------------------------------------------
752  // Step 2 - Get the infinite norm of Navier-Stokes F block.
753  // -----------------------------------------------------------------------
754 
755  // Extract the velocity blocks.
758 
759  for(unsigned row_i = 0; row_i < N_velocity_doftypes; row_i++)
760  {
761  for(unsigned col_i = 0; col_i < N_velocity_doftypes; col_i++)
762  {
763  v_aug_pt(row_i,col_i) = new CRDoubleMatrix;
764  this->get_block(row_i,col_i,*v_aug_pt(row_i,col_i));
765  } // for
766  } // for
767 
768  // Now get the infinite norm.
770  {
772  } // if
773 
774 #ifdef PARANOID
775  // Warning for division by zero.
776  if(Scaling_sigma == 0.0)
777  {
778  std::ostringstream warning_stream;
779  warning_stream << "WARNING: " << std::endl
780  << "The scaling (Scaling_sigma) is "
781  << Scaling_sigma << ".\n"
782  << "Division by 0 will occur." << "\n";
783  OomphLibWarning(warning_stream.str(),
784  OOMPH_CURRENT_FUNCTION,
785  OOMPH_EXCEPTION_LOCATION);
786  }
787  if(Scaling_sigma > 0.0)
788  {
789  std::ostringstream warning_stream;
790  warning_stream << "WARNING: " << std::endl
791  << "The scaling (Scaling_sigma) is positive: "
792  << Scaling_sigma << std::endl
793  << "Performance may be degraded."
794  << std::endl;
795  OomphLibWarning(warning_stream.str(),
796  OOMPH_CURRENT_FUNCTION,
797  OOMPH_EXCEPTION_LOCATION);
798  }
799 #endif
800 
801  // -----------------------------------------------------------------------
802  // Step 3 - Augment the constrained fluid blocks.
803  // -----------------------------------------------------------------------
804  // Loop through the Lagrange multipliers and do three things:
805  // For each Lagrange block:
806  // 3.1) Extract the mass matrices and store the location of non-zero mass
807  // matrices.
808  //
809  // 3.2) a) Create inv_w (for the augmentation)
810  // b) Store the diagonal values of inv_w (for preconditioner solve)
811  //
812  // 3.3) Perform the augmentation: v_aug + m_i * inv(w_i) * m_j
813  //
814  // 3.4) Clean up memory.
815 
816  // Storage for the inv w diag values.
817  Inv_w_diag_values.clear();
818  for(unsigned l_i = 0; l_i < N_lagrange_doftypes; l_i++)
819  {
820  // Step 3.1: Location and extraction of non-zero mass matrices.
821 
822  // Storage for the current Lagrange block mass matrices.
825 
826  // Block type for the Lagrange multiplier.
827  const unsigned lgr_block_num = N_fluid_doftypes + l_i;
828 
829  // Store the mass matrix locations for the current Lagrange block.
830  Vector<unsigned> mm_locations;
831 
832  // Store the number of mass matrices.
833  unsigned n_mm = 0;
834 
835  // Go along the block columns for the current Lagrange block ROW.
836  for(unsigned col_i = 0; col_i < N_velocity_doftypes; col_i++)
837  {
838  // Get the block matrix for this block column.
839  CRDoubleMatrix* mm_temp_pt = new CRDoubleMatrix;
840  this->get_block(lgr_block_num, col_i, *mm_temp_pt);
841 
842  // Check if this is non-zero
843  if(mm_temp_pt->nnz() > 0)
844  {
845  mm_locations.push_back(col_i);
846  mm_pt.push_back(mm_temp_pt);
847  n_mm++;
848  }
849  else
850  {
851  // This is just an empty matrix. No need to keep this.
852  delete mm_temp_pt;
853  }
854  } // loop through the columns of the Lagrange row.
855 
856 #ifdef PARANOID
857  if (n_mm == 0)
858  {
859  std::ostringstream warning_stream;
860  warning_stream << "WARNING:\n"
861  << "There are no mass matrices on Lagrange block row "
862  << l_i << ".\n"
863  << "Perhaps the problem setup is incorrect." << std::endl;
864  OomphLibWarning(warning_stream.str(),
865  OOMPH_CURRENT_FUNCTION,
866  OOMPH_EXCEPTION_LOCATION);
867  }
868 #endif
869 
870  // Get the transpose of the mass matrices.
871  for(unsigned mm_i = 0; mm_i < n_mm; mm_i++)
872  {
873  // Get the block matrix for this block column.
874  CRDoubleMatrix* mm_temp_pt = new CRDoubleMatrix;
875  this->get_block(mm_locations[mm_i],lgr_block_num,*mm_temp_pt);
876 
877  if(mm_temp_pt->nnz() > 0)
878  {
879  mmt_pt.push_back(mm_temp_pt);
880  }
881  else
882  {
883  // There should be a non-zero mass matrix here, since L=(L^T)^T
884 #ifdef PARANOID
885  {
886  std::ostringstream warning_stream;
887  warning_stream << "WARNING:\n"
888  << "The mass matrix block " << mm_locations[mm_i]
889  << " in L^T block " << l_i << " is zero.\n"
890  << "Perhaps the problem setup is incorrect." << std::endl;
891  OomphLibWarning(warning_stream.str(),
892  OOMPH_CURRENT_FUNCTION,
893  OOMPH_EXCEPTION_LOCATION);
894  }
895 #endif
896  }
897  } // loop through the ROW of the Lagrange COLUMN.
898 
899  // Step 3.2: Create inv_w and store its diagonal entries.
900 
901  // Storage for inv_w
902  CRDoubleMatrix* inv_w_pt
903  = new CRDoubleMatrix(this->Block_distribution_pt[lgr_block_num]);
904 
905  // Get the number of local rows for this Lagrange block.
906  unsigned long l_i_nrow_local
907  = this->Block_distribution_pt[lgr_block_num]->nrow_local();
908 
909  // The first row, for the column offset (required in parallel).
910  unsigned l_i_first_row
911  = this->Block_distribution_pt[lgr_block_num]->first_row();
912 
913  // A vector to contain the results of mass matrices squared.
914  Vector<double> w_i_diag_values(l_i_nrow_local,0);
915 
916  // Create mm*mm^T, and component-wise add the mass matrices
917  for(unsigned m_i = 0; m_i < n_mm; m_i++)
918  {
919  // Create mm*mm^T
920  CRDoubleMatrix temp_mm_sqrd;
921  temp_mm_sqrd.build(mm_pt[m_i]->distribution_pt());
922  mm_pt[m_i]->multiply(*mmt_pt[m_i],temp_mm_sqrd);
923 
924  // Extract diagonal entries
925  Vector<double> m_diag = temp_mm_sqrd.diagonal_entries();
926 
927  // Loop through the entries, add them.
928  for(unsigned long row_i = 0; row_i < l_i_nrow_local; row_i++)
929  {
930  w_i_diag_values[row_i] += m_diag[row_i];
931  }
932  }
933 
934  // Storage for inv_w matrix vectors
935  Vector<double> invw_i_diag_values(l_i_nrow_local,0);
936  Vector<int> w_i_column_indices(l_i_nrow_local);
937  Vector<int> w_i_row_start(l_i_nrow_local+1);
938 
939  // Divide by Scaling_sigma and create the inverse of w.
940  for(unsigned long row_i = 0; row_i < l_i_nrow_local; row_i++)
941  {
942  invw_i_diag_values[row_i] = Scaling_sigma / w_i_diag_values[row_i];
943 
944  w_i_column_indices[row_i] = row_i + l_i_first_row;
945  w_i_row_start[row_i] = row_i;
946  }
947 
948  w_i_row_start[l_i_nrow_local] = l_i_nrow_local;
949 
950  // Theses are square matrices. So we can use the l_i_nrow_global for the
951  // number of columns.
952  unsigned long l_i_nrow_global
953  = this->Block_distribution_pt[lgr_block_num]->nrow();
954  inv_w_pt->build(l_i_nrow_global,
955  invw_i_diag_values,
956  w_i_column_indices,
957  w_i_row_start);
958 
959  Inv_w_diag_values.push_back(invw_i_diag_values);
960 
961 
962  // Step 3.3: Perform the augmentation: v_aug + m_i * inv(w_i) * m_j
963 
964  ////////////////////////////////////////////////////////////////////////
965  // Now we create the augmented matrix in v_aug_pt.
966  // v_aug_pt is already re-ordered
967  // Loop through the mm_locations
968  for(unsigned ii = 0; ii < n_mm; ii++)
969  {
970  // Location of the mass matrix
971  unsigned aug_i = mm_locations[ii];
972 
973  for(unsigned jj = 0; jj < n_mm; jj++)
974  {
975  // Location of the mass matrix
976  unsigned aug_j = mm_locations[jj];
977 
978  // Storage for intermediate results.
979  CRDoubleMatrix aug_block;
980  CRDoubleMatrix another_aug_block;
981 
982  // aug_block = mmt*inv_w
983  mmt_pt[ii]->multiply((*inv_w_pt),(aug_block));
984 
985  // another_aug_block = aug_block*mm = mmt*inv_w*mm
986  aug_block.multiply(*mm_pt[jj],another_aug_block);
987 
988  // Add the augmentation.
989  v_aug_pt(aug_i,aug_j)->add(another_aug_block,
990  *v_aug_pt(aug_i,aug_j));
991  } // loop jj
992  } // loop ii
993 
994  // Step 3.4: Clean up memory.
995  delete inv_w_pt; inv_w_pt = 0;
996  for(unsigned m_i = 0; m_i < n_mm; m_i++)
997  {
998  delete mm_pt[m_i]; mm_pt[m_i] = 0;
999  delete mmt_pt[m_i]; mmt_pt[m_i] = 0;
1000  }
1001  } // loop through Lagrange multipliers.
1002 
1003  // -----------------------------------------------------------------------
1004  // Step 4 - Replace all the velocity blocks
1005  // -----------------------------------------------------------------------
1006  // When we replace (DOF) blocks, the indices have to respect the DOF type
1007  // ordering. This is because only DOF type information is passed between
1008  // preconditioners. So we need to create the inverse of the
1009  // dof_to_block_map... a block_to_dof_map!
1010  //
1011  // Note: Once the DOF blocks have been replaced, further calls to
1012  // get_block at this preconditioning hierarchy level and/or lower will use
1013  // the (nearest) replaced blocks.
1014  Vector<unsigned> block_to_dof_map(n_dof_types,0);
1015  for (unsigned dof_i = 0; dof_i < n_dof_types; dof_i++)
1016  {
1017  block_to_dof_map[dof_to_block_map[dof_i]] = dof_i;
1018  }
1019 
1020  // Now do the replacement of all blocks in v_aug_pt
1021  for (unsigned block_row_i = 0;
1022  block_row_i < N_velocity_doftypes; block_row_i++)
1023  {
1024  unsigned temp_dof_row_i = block_to_dof_map[block_row_i];
1025  for (unsigned block_col_i = 0;
1026  block_col_i < N_velocity_doftypes; block_col_i++)
1027  {
1028  unsigned temp_dof_col_i = block_to_dof_map[block_col_i];
1030  temp_dof_row_i,temp_dof_col_i,
1031  v_aug_pt(block_row_i,block_col_i));
1032  }
1033  }
1034 
1035  // -----------------------------------------------------------------------
1036  // Step 5 - Set up Navier-Stokes preconditioner
1037  // If the subsidiary fluid preconditioner is a block preconditioner:
1038  // 5.1) Set up the dof_number_in_master_map
1039  // 5.2) Set up the doftype_coarsen_map
1040  // otherwise:
1041  // 5.3) Concatenate the fluid matrices into one big matrix and pass it
1042  // to the solver.
1043  // -----------------------------------------------------------------------
1044 
1045  // First we determine if we're using a block preconditioner or not.
1047  navier_stokes_block_preconditioner_pt
1048  = dynamic_cast<BlockPreconditioner<CRDoubleMatrix>* >
1051  if(navier_stokes_block_preconditioner_pt == 0)
1052  {
1053  Navier_stokes_preconditioner_is_block_preconditioner = false;
1054  }
1055 
1056  // If the Navier-Stokes preconditioner is a block preconditioner, then we
1057  // need to turn it into a subsidiary block preconditioner.
1058  if(Navier_stokes_preconditioner_is_block_preconditioner)
1059  {
1060  // Assumption: All Navier-Stokes block preconditioners take $dim
1061  // number of velocity DOF types and 1 pressure DOF type.
1062 
1063  // Step 5.1: Set up the dof_number_in_master_map
1064 
1065  // First we create the dof_number_in_master_map vector to pass to the
1066  // subsidiary preconditioner's
1067  // turn_into_subsidiary_block_preconditioner(...) function.
1068  //
1069  // This vector maps the subsidiary DOF numbers and the master DOF
1070  // numbers and has the construction:
1071  //
1072  // dof_number_in_master_map[subsidiary DOF number] = master DOF number
1073  //
1074  // Example: Using the example above, our problem has the natural
1075  // DOF type ordering:
1076  //
1077  // 0 1 2 3 4 5 6 7 8 9 10 11 12 <- DOF number in master
1078  // [u v w p] [up vp wp Lp1 Lp2] [ut vt wt Lt1] <- DOF type
1079  //
1080  // For now, we ignore the fact that this preconditioner's number of
1081  // DOF types may be more fine grain than what is assumed by the
1082  // subsidiary block preconditioner. Normally, the order of the values in
1083  // dof_number_in_master_map matters, since the indices must match the
1084  // DOF type in the subsidiary block preconditioner (see the assumption
1085  // above). For example, if the (subsidiary) LSC block preconditioner
1086  // requires the DOF type ordering:
1087  //
1088  // [0 1 2 3]
1089  // u v w p
1090  //
1091  // Then the dof_number_in_master_map vector must match the u velocity
1092  // DOF type in the subsidiary preconditioner with the u velocity in the
1093  // master preconditioner, etc...
1094  //
1095  // However, we shall see (later) that it does not matter in this
1096  // instance because the DOF type coarsening feature overrides the
1097  // ordering provided here.
1098  //
1099  // For now, we only need to ensure that the subsidiary preconditioner
1100  // knows about 10 master DOF types (9 velocity and 1 pressure), the
1101  // ordering does not matter.
1102  //
1103  // We pass to the subsidiary block preconditioner the following
1104  // dof_number_in_master_map:
1105  // [0 1 2 4 5 6 9 10 11 3] <- DOF number in master
1106  // u v w up vp wp ut vt wt p <- corresponding DOF type
1107 
1108  // Variable to keep track of the DOF number.
1109  unsigned temp_dof_number = 0;
1110 
1111  // Storage for the dof_number_in_master_map
1112  Vector<unsigned> dof_number_in_master_map;
1113  for(unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
1114  {
1115  // Store the velocity dof types.
1116  for(unsigned dim_i = 0; dim_i < spatial_dim; dim_i++)
1117  {
1118  dof_number_in_master_map.push_back(temp_dof_number + dim_i);
1119  } // for spatial_dim
1120 
1121  // Update the DOF index
1122  temp_dof_number += My_ndof_types_in_mesh[mesh_i];
1123  } // for My_nmesh
1124 
1125  // Push back the pressure DOF type
1126  dof_number_in_master_map.push_back(spatial_dim);
1127 
1128 
1129  // Step 5.2 DOF type coarsening.
1130 
1131  // Since this preconditioner works with more fine grained DOF types than
1132  // the subsidiary block preconditioner, we need to tell the subsidiary
1133  // preconditioner which DOF types to coarsen together.
1134  //
1135  // The LSC preconditioner expects 2(3) velocity dof types, and 1
1136  // pressure DOF types. Thus, we give it this list:
1137  // u [0, 3, 6]
1138  // v [1, 4, 7]
1139  // w [2, 5, 8]
1140  // p [9]
1141  //
1142  // See how the ordering of dof_number_in_master_map (constructed above)
1143  // does not matter as long as we construct the doftype_coarsen_map
1144  // correctly.
1145 
1146  // Storage for which subsidiary DOF types to coarsen
1147  Vector<Vector<unsigned> > doftype_coarsen_map;
1148  for (unsigned direction = 0; direction < spatial_dim; direction++)
1149  {
1150  Vector<unsigned> dir_doftypes_vec(My_nmesh,0);
1151  for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
1152  {
1153  dir_doftypes_vec[mesh_i] = spatial_dim*mesh_i+direction;
1154  }
1155  doftype_coarsen_map.push_back(dir_doftypes_vec);
1156  }
1157 
1158  Vector<unsigned> ns_p_vec(1,0);
1159 
1160  // This is simply the number of velocity dof types,
1161  ns_p_vec[0] = My_nmesh*spatial_dim;
1162 
1163  doftype_coarsen_map.push_back(ns_p_vec);
1164 
1165  // Turn the Navier-Stokes block preconditioner into a subsidiary block
1166  // preconditioner.
1167  navier_stokes_block_preconditioner_pt
1169  this, dof_number_in_master_map, doftype_coarsen_map);
1170 
1171  // Call the setup function
1172  navier_stokes_block_preconditioner_pt
1173  ->setup(matrix_pt());
1174  }
1175  else
1176  {
1177  // Step 5.3: This is not a block preconditioner, thus we need to
1178  // concatenate all the fluid matrices and pass them to the solver.
1179 
1180  // Select all the fluid blocks (velocity and pressure)
1181  VectorMatrix<BlockSelector> f_aug_blocks(N_fluid_doftypes,
1182  N_fluid_doftypes);
1183  for (unsigned block_i = 0; block_i < N_fluid_doftypes; block_i++)
1184  {
1185  for (unsigned block_j = 0; block_j < N_fluid_doftypes; block_j++)
1186  {
1187  f_aug_blocks[block_i][block_j].select_block(block_i,block_j,true);
1188  }
1189  }
1190 
1191  // Note: This will use the replaced blocks.
1192  CRDoubleMatrix f_aug_block
1193  = this->get_concatenated_block(f_aug_blocks);
1194 
1196  {
1198  {
1200  }
1201  }
1202  else
1203  {
1205  {
1206  std::ostringstream err_msg;
1207  err_msg << "Not using SuperLUPreconditioner for NS block,\n"
1208  << "but the Navier_stokes_preconditioner_pt is null.\n";
1209  throw OomphLibError(err_msg.str(),
1210  OOMPH_CURRENT_FUNCTION,
1211  OOMPH_EXCEPTION_LOCATION);
1212  }
1213  }
1214 
1215  // Setup the solver.
1216  Navier_stokes_preconditioner_pt->setup(&f_aug_block);
1217 
1218  } // else
1219 
1220  // Clean up memory
1221  const unsigned v_aug_nrow = v_aug_pt.nrow();
1222  const unsigned v_aug_ncol = v_aug_pt.ncol();
1223  for (unsigned v_row = 0; v_row < v_aug_nrow; v_row++)
1224  {
1225  for (unsigned v_col = 0; v_col < v_aug_ncol; v_col++)
1226  {
1227  delete v_aug_pt(v_row,v_col);
1228  v_aug_pt(v_row,v_col) = 0;
1229  }
1230  }
1231 
1233 } // func setup
1234 
1235 /// \short Function to set a new momentum matrix preconditioner
1236 /// (inexact solver)
1238  Preconditioner* new_ns_preconditioner_pt)
1239 {
1240  // Check if pointer is non-zero.
1241  if(new_ns_preconditioner_pt == 0)
1242  {
1243  std::ostringstream warning_stream;
1244  warning_stream << "WARNING: \n"
1245  << "The LSC preconditioner point is null.\n"
1246  << "Using default (SuperLU) preconditioner.\n"
1247  << std::endl;
1248  OomphLibWarning(warning_stream.str(),
1249  OOMPH_CURRENT_FUNCTION,
1250  OOMPH_EXCEPTION_LOCATION);
1251 
1254  }
1255  else
1256  {
1257  // If the default SuperLU preconditioner has been used
1258  // clean it up now...
1261  {
1264  }
1265 
1266  Navier_stokes_preconditioner_pt = new_ns_preconditioner_pt;
1268  }
1269 }// func set_navier_stokes_preconditioner
1270 
1271 //========================================================================
1272 /// \short Clears the memory.
1273 //========================================================================
1275 {
1276  // clean the block preconditioner base class memory
1278 
1279  // Delete the Navier-Stokes preconditioner pointer if we have created it.
1282  {
1285  }
1286 
1288 } // func LagrangeEnforcedFlowPreconditioner::clean_up_memory
1289 
1290 }// namespace oomph
void set_meshes(const Vector< Mesh * > &mesh_pt)
Set the meshes, the first mesh in the vector must be the bulk mesh.
unsigned & amg_simple_smoother()
Access function to AMG_simple_smoother flag.
Definition: hypre_solver.h:888
bool built() const
unsigned nrow_local() const
access function for the num of local rows on this processor.
void clear_block_preconditioner_base()
Clears all BlockPreconditioner data. Called by the destructor and the block_setup(...) methods.
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.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
void clear()
wipes the DoubleVector
void set_navier_stokes_preconditioner(Preconditioner *new_ns_preconditioner_pt=0)
Set a new Navier-Stokes matrix preconditioner (inexact solver)
double & amg_damping()
Access function to AMG_damping parameter.
Definition: hypre_solver.h:916
double * values_pt()
access function to the underlying values
Vector< LinearAlgebraDistribution * > Block_distribution_pt
The distribution for the blocks.
cstr elem_len * i
Definition: cfortran.h:607
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...
unsigned nrow() const
access function to the number of global rows.
void return_concatenated_block_vector(const Vector< unsigned > &block_vec_number, const DoubleVector &b, DoubleVector &v) const
Takes concatenated block ordered vector, b, and copies its entries to the appropriate entries in the ...
bool is_subsidiary_block_preconditioner() const
Return true if this preconditioner is a subsidiary preconditioner.
unsigned & amg_coarsening()
Access function to AMG_coarsening flag.
Definition: hypre_solver.h:910
void return_block_vector(const unsigned &n, const DoubleVector &b, DoubleVector &v) const
Takes the n-th block ordered vector, b, and copies its entries to the appropriate entries in the natu...
double Scaling_sigma
Scaling for the augmentation: Scaling_sigma*(LL^T)
An interface to the Trilinos AztecOO classes allowing it to be used as an Oomph-lib LinearSolver...
void set_amg_iterations(const unsigned &amg_iterations)
Function to set the number of times to apply BoomerAMG.
Definition: hypre_solver.h:872
Preconditioner * boomer_amg_for_2D_poisson_problem()
Hypre Boomer AMG setting for the 2D Poisson problem (for serial code).
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
Matrix-based diagonal preconditioner.
bool Navier_stokes_preconditioner_is_block_preconditioner
Flag to indicate if the preconditioner for the Navier-Stokes block is a block preconditioner or not...
void set_replacement_dof_block(const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix *replacement_dof_block_pt)
Set replacement dof-level blocks. Only dof-level blocks can be set. This is important due to how the ...
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1068
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void get_concatenated_block_vector(const Vector< unsigned > &block_vec_number, const DoubleVector &v, DoubleVector &b)
Takes the naturally ordered vector and extracts the blocks indicated by the block number (the values)...
unsigned My_nmesh
The number of meshes. This is used to create various lookup lists.
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...
unsigned nmesh() const
Return the number of meshes in Mesh_pt.
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
unsigned & amg_smoother_iterations()
Access function to AMG_smoother_iterations.
Definition: hypre_solver.h:907
A preconditioner for performing inner iteration preconditioner solves. The template argument SOLVER s...
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:507
void set_nmesh(const unsigned &n)
Specify the number of meshes required by this block preconditioner. Note: elements in different meshe...
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.
void set_mesh(const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Set the i-th mesh for this block preconditioner. Note: The method set_nmesh(...) must be called befor...
void get_block(const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const
Put block (i,j) into output_matrix. This block accounts for any coarsening of dof types and any repla...
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.
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
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...
unsigned ndof_types() const
Return the total number of DOF types.
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...
double & amg_strength()
Access function to AMG_strength.
Definition: hypre_solver.h:919
void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< MATRIX > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
Function to turn this preconditioner into a subsidiary preconditioner that operates within a bigger "...
Vector< double > diagonal_entries() const
returns a Vector of diagonal entries of this matrix. This only works with square matrices. This condition may be relaxed in the future if need be.
Definition: matrices.cc:3442
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1805
Vector< Mesh * > My_mesh_pt
Storage for the meshes. In our implementation, the first mesh must always be the Navier-Stokes (bulk)...
Class for dense matrices, storing all the values of the matrix as a pointer to a pointer with assorte...
Definition: communicator.h:50
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.
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...
CRDoubleMatrix * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*...
double inf_norm(const DenseMatrix< CRDoubleMatrix * > &matrix_pt)
Compute infinity (maximum) norm of sub blocks as if it was one matrix.
Definition: matrices.cc:3708
void get_block_vector(const unsigned &n, const DoubleVector &v, DoubleVector &b) const
Takes the naturally ordered vector, v and returns the n-th block vector, b. Here n is the block numbe...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
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...
void build(const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)
build method: vector of values, vector of column indices, vector of row starts and number of rows and...
Definition: matrices.cc:1692
CRDoubleMatrix get_concatenated_block(const VectorMatrix< BlockSelector > &selected_block)
Returns a concatenation of the block matrices specified by the argument selected_block. The VectorMatrix selected_block must be correctly sized as it is used to determine the number of sub block matrices to concatenate.