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