pseudo_elastic_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//====================================================================
26
28
29namespace oomph
30{
31 //=============================================================================
32 /// Functions to create instances of optimal subsidiary operators for
33 /// the PseudoElasticPreconditioner
34 //=============================================================================
35 namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper
36 {
37#ifdef OOMPH_HAS_HYPRE
38
39 /// AMG w/ GS smoothing for the augmented elastic subsidiary linear
40 /// systems
42 {
43 HyprePreconditioner* hypre_preconditioner_pt =
44 new HyprePreconditioner("Hypre for diagonal blocks in pseudo-solid");
45 hypre_preconditioner_pt->set_amg_iterations(2);
46 hypre_preconditioner_pt->amg_using_simple_smoothing();
47 if (MPI_Helpers::communicator_pt()->nproc() > 1)
48 {
49 // Jacobi in parallel
50 hypre_preconditioner_pt->amg_simple_smoother() = 0;
51 }
52 else
53 {
54 // Gauss Seidel in serial (was 3 actually...)
55 hypre_preconditioner_pt->amg_simple_smoother() = 1;
56 }
57 hypre_preconditioner_pt->hypre_method() = HyprePreconditioner::BoomerAMG;
58 hypre_preconditioner_pt->amg_damping() = 1.0;
59 hypre_preconditioner_pt->amg_coarsening() = 6;
60 return hypre_preconditioner_pt;
61 }
62
63
64 /// AMG w/ GS smoothing for the augmented elastic subsidiary linear
65 /// systems -- calls Hypre version to stay consistent with previous default
67 {
69 }
70
71#endif
72
73#ifdef OOMPH_HAS_TRILINOS
74
75 /// TrilinosML smoothing for the augmented elastic
76 /// subsidiary linear systems
78 {
80 return trilinos_prec_pt;
81 }
82
83 /// CG with diagonal preconditioner for the lagrange multiplier
84 /// subsidiary linear systems.
86 {
91 // Note: This makes CG a proper "inner iteration" for
92 // which GMRES (may) no longer converge. We should really
93 // use FGMRES or GMRESR for this. However, here the solver
94 // is so good that it'll converge very quickly anyway
95 // so there isn't much to be gained by limiting the number
96 // of iterations...
97 // prec_pt->max_iter() = 4;
98 prec_pt->solver_pt()->solver_type() = TrilinosAztecOOSolver::CG;
99 prec_pt->solver_pt()->disable_doc_time();
100 return prec_pt;
101 }
102#endif
103 } // namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper
104
105
106 /// ////////////////////////////////////////////////////////////////////////////
107 /// ////////////////////////////////////////////////////////////////////////////
108 /// ////////////////////////////////////////////////////////////////////////////
109
110
111 //=============================================================================
112 // Setup method for the PseudoElasticPreconditioner.
113 //=============================================================================
115 {
116 // clean
117 this->clean_up_memory();
118
119#ifdef PARANOID
120 // paranoid check that meshes have been set
121 if (Elastic_mesh_pt == 0)
122 {
123 std::ostringstream error_message;
124 error_message << "The elastic mesh must be set.";
125 throw OomphLibError(
126 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
127 }
129 {
130 std::ostringstream error_message;
131 error_message << "The Lagrange multiplier mesh must be set.";
132 throw OomphLibError(
133 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
134 }
135#endif
136
137 // set the meshes
138 this->set_nmesh(2);
139 this->set_mesh(0, Elastic_mesh_pt);
141
142 // Figure out number of dof types
143 unsigned n_solid_dof_types = 0;
144 unsigned n_dof_types = 0;
145
147 {
148 // get the number of solid dof types from the first element
149 n_solid_dof_types = this->ndof_types_in_mesh(0);
150
151 // get the total number of dof types
152 n_dof_types = n_solid_dof_types + this->ndof_types_in_mesh(1);
153 }
154 else
155 {
156 n_dof_types = this->ndof_types();
157 n_solid_dof_types = n_dof_types - (n_dof_types / 3);
158 }
159#ifdef PARANOID
160 if (n_dof_types % 3 != 0)
161 {
162 std::ostringstream error_message;
163 error_message << "This preconditioner requires DIM*3 types of DOF";
164 throw OomphLibError(
165 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
166 }
167#endif
168
169 // determine the dimension
170 Dim = n_dof_types / 3;
171
172#ifdef PARANOID
173 // Recast Jacobian matrix to CRDoubleMatrix
174 CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
175
176 if (cr_matrix_pt == 0)
177 {
178 std::ostringstream error_message;
179 error_message << "FSIPreconditioner only works with"
180 << " CRDoubleMatrix matrices" << std::endl;
181 throw OomphLibError(
182 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
183 }
184#endif
185
186 // Setup the dof_list scheme for block_setup.
187 // The natural DOF types order is (in 3D):
188 // 0 1 2 3 4 5 6 7 8 <- natural DOF type index
189 // x y z xc yc zc lx ly lz <- DOf type
190 //
191 // We need to group the directional displacements together:
192 // 0 1 2 3 4 5 6 7 8 <- block index
193 // x xc y yc xz zc lx ly lz <- DOF type
194 //
195 // The mapping required is:
196 // 0 2 4 1 3 5 6 7 8
197 //
198 // In general we want:
199 //
200 // dof_to_block_map[DOF type] = block type
201 Vector<unsigned> dof_list_for_block_setup(n_dof_types);
202 for (unsigned dim_i = 0; dim_i < Dim; dim_i++)
203 {
204 // bulk solid dof types
205 dof_list_for_block_setup[dim_i] = 2 * dim_i;
206
207 // constrained solid dof types
208 dof_list_for_block_setup[dim_i + Dim] = 2 * dim_i + 1;
209
210 // lagr dof types
211 dof_list_for_block_setup[dim_i + 2 * Dim] = 2 * Dim + dim_i;
212 }
213
214 // Setup the block ordering. We have the following:
215 // Block types [0, 2*Dim) are the solid blocks.
216 // Block types [2*Dim, 3*Dim) are the Lagrange multiplier blocks.
217 //
218 // For d = 0, 1, 2,
219 // Bulk solid doftypes: 2*d
220 // Constrained solid doftypes: 2*d+1
221 // Lgr doftyoes: 2*Dim + d
222 this->block_setup(dof_list_for_block_setup);
223
224 // Create dof_list for subsidiary preconditioner. This will be used later
225 // in turn_into_subsidiary_block_preconditioner(...)
226 Vector<unsigned> dof_list_for_subsidiary_prec(n_solid_dof_types);
227
228 // The dof types are currently in the order (in 3D):
229 // 0 1 2 3 4 5 6 7 8 <- DOF index
230 // x y z xc yc zc lx ly lz <- DOF type
231 //
232 // We need to group the directional displacements together:
233 // x xc y yc xz zc
234 //
235 // The mapping required is:
236 // 0 3 1 4 2 5
237 //
238 // In general we want:
239 // dof_list[subsidiary DOF type] = master DOF type
240 for (unsigned dim_i = 0; dim_i < Dim; dim_i++)
241 {
242 // bulk solid dof
243 dof_list_for_subsidiary_prec[2 * dim_i] = dim_i;
244
245 // constrained solid dof
246 dof_list_for_subsidiary_prec[2 * dim_i + 1] = dim_i + Dim;
247 }
248
249 // Get the solid blocks
250 DenseMatrix<CRDoubleMatrix*> solid_matrix_pt(
251 n_solid_dof_types, n_solid_dof_types, 0);
252
253 for (unsigned row_i = 0; row_i < n_solid_dof_types; row_i++)
254 {
255 for (unsigned col_i = 0; col_i < n_solid_dof_types; col_i++)
256 {
257 solid_matrix_pt(row_i, col_i) = new CRDoubleMatrix;
258 this->get_block(row_i, col_i, *solid_matrix_pt(row_i, col_i));
259 }
260 }
261
262 // compute the scaling (if required)
264 {
265 Scaling = CRDoubleMatrixHelpers::inf_norm(solid_matrix_pt);
266 }
267 else
268 {
269 Scaling = 1.0;
270 }
271
272 // Add the scaled identify matrix to the constrained solid blocks.
273 for (unsigned d = 0; d < Dim; d++)
274 {
275 // Where is the constrained block located?
276 unsigned block_i = 2 * d + 1;
277
278 // Data from the constrained block.
279 double* s_values = solid_matrix_pt(block_i, block_i)->value();
280 int* s_column_index = solid_matrix_pt(block_i, block_i)->column_index();
281 int* s_row_start = solid_matrix_pt(block_i, block_i)->row_start();
282 int s_nrow_local = solid_matrix_pt(block_i, block_i)->nrow_local();
283 int s_first_row = solid_matrix_pt(block_i, block_i)->first_row();
284
285 // Loop through the diagonal entries and add the scaling.
286 for (int i = 0; i < s_nrow_local; i++)
287 {
288 bool found = false;
289 for (int j = s_row_start[i]; j < s_row_start[i + 1] && !found; j++)
290 {
291 if (s_column_index[j] == i + s_first_row)
292 {
293 s_values[j] += Scaling;
294 found = true;
295 } // if
296 } // for row start
297
298 // Check if the diagonal entry is found.
299 if (!found)
300 {
301 std::ostringstream error_message;
302 error_message << "The diagonal entry for the constained block("
303 << block_i << "," << block_i << ")\n"
304 << "on local row " << i << " does not exist."
305 << std::endl;
306 throw OomphLibError(error_message.str(),
307 OOMPH_CURRENT_FUNCTION,
308 OOMPH_EXCEPTION_LOCATION);
309 }
310 } // for nrow_local
311 } // for Dim
312
313
314 // setup the solid subsidiary preconditioner
315 // //////////////////////////////// this preconditioner uses the full S
316 // matrix
318 {
321
322 // Add mesh (not actually used since this only acts as a subsidiary
323 // preconditioner...
324 // s_prec_pt->add_mesh(Elastic_mesh_pt);
325
326 Vector<Vector<unsigned>> doftype_to_doftype_map(n_solid_dof_types,
327 Vector<unsigned>(1, 0));
328
329 for (unsigned i = 0; i < n_solid_dof_types; i++)
330 {
331 doftype_to_doftype_map[i][0] = i;
332 }
333
335 this, dof_list_for_subsidiary_prec, doftype_to_doftype_map);
336
338 {
341 }
342
343 // Set the replacement blocks.
344 for (unsigned d = 0; d < Dim; d++)
345 {
346 // The dof-block located in the block-ordering.
347 unsigned block_i = 2 * d + 1;
348
349 // The dof-block located in the dof-ordering.
350 unsigned dof_block_i = Dim + d;
352 dof_block_i, dof_block_i, solid_matrix_pt(block_i, block_i));
353 }
354
355 s_prec_pt->Preconditioner::setup(matrix_pt());
356 Elastic_preconditioner_pt = s_prec_pt;
357 }
358 // otherwise it is a block based preconditioner
359 else
360 {
362
363 // set the block preconditioning method
364 switch (E_preconditioner_type)
365 {
367 {
369 }
370 break;
372 {
374 block_triangular_prec_pt =
376 block_triangular_prec_pt->upper_triangular();
377
378 s_prec_pt = block_triangular_prec_pt;
379 }
380 break;
382 {
384 block_triangular_prec_pt =
386 block_triangular_prec_pt->lower_triangular();
387
388 s_prec_pt = block_triangular_prec_pt;
389 }
390 break;
391 default:
392 {
393 std::ostringstream error_msg;
394 error_msg
395 << "There is no such block based preconditioner.\n"
396 << "Candidates are:\n"
397 << "PseudoElasticPreconditioner::Block_diagonal_preconditioner\n"
398 << "PseudoElasticPreconditioner::Block_upper_triangular_"
399 "preconditioner\n"
400 << "PseudoElasticPreconditioner::Block_lower_triangular_"
401 "preconditioner\n";
402 throw OomphLibError(
403 error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
404 }
405 break;
406 }
407
408
409 // Add mesh (not actually used since this only acts as a subsidiary
410 // preconditioner...
411 // s_prec_pt->add_mesh(Elastic_mesh_pt);
412
413 // The block to block map
414 Vector<Vector<unsigned>> doftype_to_doftype_map(Dim,
415 Vector<unsigned>(2, 0));
416 Vector<unsigned> s_prec_dof_to_block_map(Dim, 0);
417
418 unsigned tmp_index = 0;
419 for (unsigned d = 0; d < Dim; d++)
420 {
421 s_prec_dof_to_block_map[d] = d;
422 doftype_to_doftype_map[d][0] = tmp_index++;
423 doftype_to_doftype_map[d][1] = tmp_index++;
424 }
425
427 this, dof_list_for_subsidiary_prec, doftype_to_doftype_map);
428
430 {
433 }
434
435 // Set the replacement blocks.
436 for (unsigned d = 0; d < Dim; d++)
437 {
438 // The dof-block located in the block-ordering.
439 unsigned block_i = 2 * d + 1;
440
441 // Then dof-block located in the dof-ordering.
442 unsigned dof_block_i = Dim + d;
444 dof_block_i, dof_block_i, solid_matrix_pt(block_i, block_i));
445 }
446
447 s_prec_pt->set_dof_to_block_map(s_prec_dof_to_block_map);
448 s_prec_pt->Preconditioner::setup(matrix_pt());
449
450 Elastic_preconditioner_pt = s_prec_pt;
451 }
452
453 // No longer require the solid blocks
454 for (unsigned row_i = 0; row_i < n_solid_dof_types; row_i++)
455 {
456 for (unsigned col_i = 0; col_i < n_solid_dof_types; col_i++)
457 {
458 delete solid_matrix_pt(row_i, col_i);
459 solid_matrix_pt(row_i, col_i) = 0;
460 }
461 }
462
463 // next setup the lagrange multiplier preconditioners
465 for (unsigned d = 0; d < Dim; d++)
466 {
467 CRDoubleMatrix* b_pt = new CRDoubleMatrix;
468 this->get_block(2 * Dim + d, 2 * d + 1, *b_pt);
469
470 // if a non default preconditioner is specified create
471 // the preconditioners
473 {
475 (*Lagrange_multiplier_subsidiary_preconditioner_function_pt)();
476 }
477 // else use default superlu preconditioner
478 else
479 {
481 }
482
483 // and setup
485 delete b_pt;
486 b_pt = 0;
487 }
488 }
489
490 //=============================================================================
491 /// Apply the elastic subsidiary preconditioner.
492 //=============================================================================
494 const DoubleVector& r, DoubleVector& z)
495 {
496 // apply the solid preconditioner
498 }
499
500 //=============================================================================
501 /// Apply the lagrange multiplier subsidiary preconditioner.
502 //=============================================================================
504 const DoubleVector& r, DoubleVector& z)
505 {
506 // apply the lagrange multiplier preconditioner
507 for (unsigned d = 0; d < Dim; d++)
508 {
509 DoubleVector x;
510 this->get_block_vector(Dim * 2 + d, r, x);
511 DoubleVector y;
512 Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(x, y);
513 Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(y, x);
514 unsigned nrow_local = x.nrow_local();
515 double* x_pt = x.values_pt();
516 for (unsigned i = 0; i < nrow_local; i++)
517 {
518 x_pt[i] = x_pt[i] * Scaling;
519 }
520 this->return_block_vector(Dim * 2 + d, x, z);
521 }
522 }
523
524 //=============================================================================
525 /// Clears the memory.
526 //=============================================================================
528 {
529 // clean the block preconditioner base class memory
531
532 // delete the solid preconditioner
535
536 // delete the lagrange multiplier preconditioner pt
537 unsigned sz = Lagrange_multiplier_preconditioner_pt.size();
538 for (unsigned i = 0; i < sz; i++)
539 {
542 }
543 }
544
545 /// ////////////////////////////////////////////////////////////////////////////
546 /// ////////////////////////////////////////////////////////////////////////////
547 /// ////////////////////////////////////////////////////////////////////////////
548
549 //=============================================================================
550 // Setup method for the PseudoElasticPreconditioner.
551 //=============================================================================
553 {
554 // clean
555 this->clean_up_memory();
556
557#ifdef PARANOID
558 // paranoid check that meshes have been set
559 if (Elastic_mesh_pt == 0)
560 {
561 std::ostringstream error_message;
562 error_message << "The elastic mesh must be set.";
563 throw OomphLibError(
564 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
565 }
567 {
568 std::ostringstream error_message;
569 error_message << "The Lagrange multiplier mesh must be set.";
570 throw OomphLibError(
571 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
572 }
573#endif
574
575 // set the mesh
576 unsigned n_solid_dof_types = 0;
577 unsigned n_dof_types = 0;
578 this->set_mesh(0, Elastic_mesh_pt);
581 {
582 // get the number of solid dof types from the first element
583 n_solid_dof_types = this->ndof_types_in_mesh(0);
584
585 // get the total number of dof types
586 n_dof_types = n_solid_dof_types + this->ndof_types_in_mesh(1);
587 }
588 else
589 {
590 n_dof_types = this->ndof_types();
591 n_solid_dof_types = n_dof_types - (n_dof_types / 3);
592 }
593#ifdef PARANOID
594 if (n_dof_types % 3 != 0)
595 {
596 std::ostringstream error_message;
597 error_message << "This preconditioner requires DIM*3 types of DOF";
598 throw OomphLibError(
599 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
600 }
601#endif
602
603 // determine the dimension
604 Dim = n_dof_types / 3;
605
606 // Recast Jacobian matrix to CRDoubleMatrix
607 CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
608
609#ifdef PARANOID
610 if (cr_matrix_pt == 0)
611 {
612 std::ostringstream error_message;
613 error_message << "FSIPreconditioner only works with"
614 << " CRDoubleMatrix matrices" << std::endl;
615 throw OomphLibError(
616 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
617 }
618#endif
619
620 // Setup the block ordering.
621 this->block_setup();
622
623 // Create dof_list for subsidiary preconditioner. This will be used later
624 // in turn_into_subsidiary_block_preconditioner(...)
625 Vector<unsigned> dof_list_for_subsidiary_prec(n_solid_dof_types);
626 for (unsigned i = 0; i < n_solid_dof_types; i++)
627 {
628 dof_list_for_subsidiary_prec[i] = i;
629 }
630
631 // compute the scaling (if required)
633 {
634 Vector<unsigned> dof_list(n_solid_dof_types);
635 for (unsigned i = 0; i < n_solid_dof_types; i++)
636 {
637 dof_list[i] = i;
638 }
639
642 this, cr_matrix_pt, dof_list, Elastic_mesh_pt, comm_pt());
643 Scaling = helper_pt->s_inf_norm();
644 delete helper_pt;
645 helper_pt = 0;
646 }
647 else
648 {
649 Scaling = 1.0;
650 }
651
652
653 // setup the solid subsidiary preconditioner
654 // //////////////////////////////// this preconditioner uses the full S
655 // matrix
657 {
660 Vector<unsigned> dof_list(n_solid_dof_types);
661 for (unsigned i = 0; i < n_solid_dof_types; i++)
662 {
663 dof_list[i] = i;
664 }
665 s_prec_pt->turn_into_subsidiary_block_preconditioner(this, dof_list);
667 {
670 }
671 s_prec_pt->scaling() = Scaling;
672 s_prec_pt->Preconditioner::setup(matrix_pt());
673 Elastic_preconditioner_pt = s_prec_pt;
674 }
675
676 // otherwise it is a block based preconditioner
677 else
678 {
679 // create the preconditioner
682 Vector<unsigned> dof_list(n_solid_dof_types);
683 for (unsigned i = 0; i < n_solid_dof_types; i++)
684 {
685 dof_list[i] = i;
686 }
687 s_prec_pt->turn_into_subsidiary_block_preconditioner(this, dof_list);
688
689 // set the subsidiary solve method
691 {
694 }
695
696 // set the scaling
697 s_prec_pt->scaling() = Scaling;
698
699 // set the block preconditioning method
700 switch (E_preconditioner_type)
701 {
704 break;
707 break;
710 break;
711 default:
712 break;
713 }
714
715 // setup
716 s_prec_pt->Preconditioner::setup(matrix_pt());
717 Elastic_preconditioner_pt = s_prec_pt;
718 }
719
720 // next setup the lagrange multiplier preconditioners
722 for (unsigned d = 0; d < Dim; d++)
723 {
724 CRDoubleMatrix* b_pt = new CRDoubleMatrix;
725 this->get_block(2 * Dim + d, Dim + d, *b_pt);
726
727 // if a non default preconditioner is specified create
728 // the preconditioners
730 {
732 (*Lagrange_multiplier_subsidiary_preconditioner_function_pt)();
733 }
734
735 // else use default superlu preconditioner
736 else
737 {
739 }
740
741 // and setup
743 delete b_pt;
744 b_pt = 0;
745 }
746 }
747
748 //=============================================================================
749 /// Apply the elastic subsidiary preconditioner.
750 //=============================================================================
752 const DoubleVector& r, DoubleVector& z)
753 {
754 // apply the solid preconditioner
756 }
757
758
759 //=============================================================================
760 /// Apply the lagrange multiplier subsidiary preconditioner.
761 //=============================================================================
763 const DoubleVector& r, DoubleVector& z)
764 {
765 // apply the lagrange multiplier preconditioner
766 for (unsigned d = 0; d < Dim; d++)
767 {
768 DoubleVector x;
769 this->get_block_vector(Dim * 2 + d, r, x);
770 DoubleVector y;
771 Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(x, y);
772 Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(y, x);
773 unsigned nrow_local = x.nrow_local();
774 double* x_pt = x.values_pt();
775 for (unsigned i = 0; i < nrow_local; i++)
776 {
777 x_pt[i] = x_pt[i] * Scaling;
778 }
779 this->return_block_vector(Dim * 2 + d, x, z);
780 }
781 }
782
783 //=============================================================================
784 /// Clears the memory.
785 //=============================================================================
787 {
788 // clean the block preconditioner base class memory
790
791 // delete the solid preconditioner
794
795 // delete the lagrange multiplier preconditioner pt
796 unsigned sz = Lagrange_multiplier_preconditioner_pt.size();
797 for (unsigned i = 0; i < sz; i++)
798 {
801 }
802 }
803
804
805 /// ////////////////////////////////////////////////////////////////////////////
806 /// ////////////////////////////////////////////////////////////////////////////
807 /// ////////////////////////////////////////////////////////////////////////////
808
809
810 //=============================================================================
811 /// Setup the preconditioner
812 //=============================================================================
814 {
815 // clean memory
816 this->clean_up_memory();
817
818#ifdef PARANOID
819 // paranoid check that this preconditioner has an even number of DOF types
820 if (this->ndof_types() % 2 != 0)
821 {
822 std::ostringstream error_message;
823 error_message
824 << "This SUBSIDIARY preconditioner requires an even number of "
825 << "types of DOF";
826 throw OomphLibError(
827 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
828 }
829#endif
830
831 // assemble dof_to_block_map
832 unsigned ndof_types = this->ndof_types();
833 Vector<unsigned> dof_to_block_map(ndof_types, 0);
834 for (unsigned i = ndof_types / 2; i < ndof_types; i++)
835 {
836 dof_to_block_map[i] = 1;
837 }
838
839 this->block_setup(dof_to_block_map);
840
841 // get block 11
842 CRDoubleMatrix* s11_pt = new CRDoubleMatrix;
843 this->get_block(1, 1, *s11_pt);
844
845 // add the scaled identity matrix to block 11
846 double* s11_values = s11_pt->value();
847 int* s11_column_index = s11_pt->column_index();
848 int* s11_row_start = s11_pt->row_start();
849 int s11_nrow_local = s11_pt->nrow_local();
850 int s11_first_row = s11_pt->first_row();
851 for (int i = 0; i < s11_nrow_local; i++)
852 {
853 bool found = false;
854 for (int j = s11_row_start[i]; j < s11_row_start[i + 1] && !found; j++)
855 {
856 if (s11_column_index[j] == i + s11_first_row)
857 {
858 s11_values[j] += Scaling;
859 found = true;
860 }
861 }
862 }
863
864 VectorMatrix<BlockSelector> required_blocks(2, 2);
865 const bool want_block = true;
866 for (unsigned b_i = 0; b_i < 2; b_i++)
867 {
868 for (unsigned b_j = 0; b_j < 2; b_j++)
869 {
870 required_blocks[b_i][b_j].select_block(b_i, b_j, want_block);
871 }
872 }
873
874 required_blocks[1][1].set_replacement_block_pt(s11_pt);
875
876 CRDoubleMatrix s_prec_pt = this->get_concatenated_block(required_blocks);
877
878 delete s11_pt;
879 s11_pt = 0;
880
881 // setup the preconditioner
883 {
884 Preconditioner_pt = (*Subsidiary_preconditioner_function_pt)();
885 }
886 else
887 {
889 }
890 Preconditioner_pt->setup(&s_prec_pt);
891 }
892
893 //=============================================================================
894 /// Apply the preconditioner.
895 //=============================================================================
898 {
899 DoubleVector x;
901 DoubleVector y;
904 }
905
906
907 /// ////////////////////////////////////////////////////////////////////////////
908 /// ////////////////////////////////////////////////////////////////////////////
909 /// ////////////////////////////////////////////////////////////////////////////
910
911
912 //=============================================================================
913 /// clean up the memory
914 //=============================================================================
917 {
918 // number of block types
919 unsigned n_block = Diagonal_block_preconditioner_pt.size();
920
921 // delete diagonal blocks
922 for (unsigned i = 0; i < n_block; i++)
923 {
926 if (Method == 1)
927 {
928 for (unsigned j = i + 1; j < n_block; j++)
929 {
932 }
933 }
934 else if (Method == 2)
935 {
936 for (unsigned j = 0; j < i; j++)
937 {
940 }
941 }
942 }
943
944 // clean up the block preconditioner
946 }
947
948 //=============================================================================
949 /// Setup the preconditioner.
950 //=============================================================================
952 {
953 // clean the memory
954 this->clean_up_memory();
955
956 // determine the number of DOF types
957 unsigned n_dof_types = this->ndof_types();
958
959#ifdef PARANOID
960 // must be Dim*2 dof types
961 if (n_dof_types % 2 != 0)
962 {
963 std::ostringstream error_message;
964 error_message << "This preconditioner requires DIM*3 types of DOF";
965 throw OomphLibError(
966 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
967 }
968#endif
969
970 // store the dimension of the problem
971 unsigned dim = n_dof_types / 2;
972
973 // assemble the dof to block lookup scheme
974 Vector<unsigned> dof_to_block_map(n_dof_types, 0);
975 for (unsigned d = 0; d < dim; d++)
976 {
977 dof_to_block_map[d] = d;
978 dof_to_block_map[d + dim] = d;
979 }
980
981 // setup the blocks look up schemes
982 this->block_setup(dof_to_block_map);
983
984 // Storage for the diagonal block preconditioners
986
987 // storage for the off diagonal matrix vector products
988 Off_diagonal_matrix_vector_products.resize(dim, dim, 0);
989
990 // setup the subsidiary preconditioners
991 for (unsigned d = 0; d < dim; d++)
992 {
993 Vector<unsigned> dof_list(2);
994 dof_list[0] = d;
995 dof_list[1] = d + dim;
996
1000 ->turn_into_subsidiary_block_preconditioner(this, dof_list);
1002 {
1004 ->set_subsidiary_preconditioner_function(
1006 }
1008
1009 Diagonal_block_preconditioner_pt[d]->Preconditioner::setup(matrix_pt());
1010
1011 // the preconditioning method.\n
1012 // 0 - block diagonal\n
1013 // 1 - upper triangular\n
1014 // 2 - lower triangular\n
1015 // next setup the off diagonal mat vec operators if required
1016 if (Method == 1 || Method == 2)
1017 {
1018 unsigned l = d + 1;
1019 unsigned u = dim;
1020 if (Method == 2)
1021 {
1022 l = 0;
1023 u = d;
1024 }
1025 for (unsigned j = l; j < u; j++)
1026 {
1027 CRDoubleMatrix* block_matrix_pt = new CRDoubleMatrix;
1028 this->get_block(d, j, *block_matrix_pt);
1030 // Off_diagonal_matrix_vector_products(d,j)->setup(block_matrix_pt);
1032 Off_diagonal_matrix_vector_products(d, j), block_matrix_pt, j);
1033
1034 delete block_matrix_pt;
1035 block_matrix_pt = 0;
1036 }
1037 }
1038 } // setup the subsidiary preconditioner.
1039 } // preconditioner setup
1040
1041 //=============================================================================
1042 /// Apply preconditioner to r
1043 //=============================================================================
1046 {
1047 // copy r
1048 DoubleVector r(res);
1049
1050 unsigned n_block;
1051
1052 // Cache umber of block types (also the spatial DIM)
1053 n_block = this->nblock_types();
1054
1055 // loop parameters
1056 int start = n_block - 1;
1057 int end = -1;
1058 int step = -1;
1059 if (Method != 1)
1060 {
1061 start = 0;
1062 end = n_block;
1063 step = 1;
1064 }
1065
1066 // the preconditioning method.
1067 // 0 - block diagonal
1068 // 1 - upper triangular
1069 // 2 - lower triangular
1070 //
1071 // loop over the DIM
1072 //
1073 // For Method = 0 or 2 (diagonal, lower)
1074 // start = 2, end = -1, step = -1
1075 // i = 2,1,0
1076 //
1077 // For Method = 1 (upper)
1078 // start = 0, end = 3 step = 1
1079 // i = 0, 1, 2
1080 for (int i = start; i != end; i += step)
1081 {
1082 // solve
1083 Diagonal_block_preconditioner_pt[i]->preconditioner_solve(r, z);
1084
1085 // if upper or lower triangular
1086 if (Method != 0)
1087 {
1088 // substitute
1089 //
1090 for (int j = i + step; j != end; j += step)
1091 {
1092 DoubleVector x;
1093 this->get_block_vector(i, z, x);
1094 DoubleVector y;
1095 Off_diagonal_matrix_vector_products(j, i)->multiply(x, y);
1096 x.clear();
1097 this->get_block_vector(j, r, x);
1098 x -= y;
1099 this->return_block_vector(j, x, r);
1100 } // substitute
1101 } // if upper or lower
1102 } // for loop over DIM
1103 } // Block preconditioner solve
1104} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
Block diagonal preconditioner. By default SuperLU is used to solve the subsidiary systems,...
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...
unsigned ndof_types_in_mesh(const unsigned &i) const
Return the number of DOF types in mesh i. WARNING: This should only be used by the upper-most master ...
CRDoubleMatrix get_concatenated_block(const VectorMatrix< BlockSelector > &selected_block)
Returns a concatenation of the block matrices specified by the argument selected_block....
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...
bool is_master_block_preconditioner() const
Return true if this preconditioner is the master block preconditioner.
unsigned nblock_types() const
Return the number of block types.
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 ...
void return_block_ordered_preconditioner_vector(const DoubleVector &w, DoubleVector &v) const
Takes the block ordered vector, w, and reorders it in natural order. Reordered vector is returned in ...
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 setup_matrix_vector_product(MatrixVectorProduct *matvec_prod_pt, CRDoubleMatrix *block_pt, const Vector< unsigned > &block_col_indices)
Setup a matrix vector product. matvec_prod_pt is a pointer to the MatrixVectorProduct,...
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 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_ordered_preconditioner_vector(const DoubleVector &v, DoubleVector &w)
Given the naturally ordered vector, v, return the vector rearranged in block order in w....
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
void upper_triangular()
Use as an upper triangular preconditioner.
void lower_triangular()
Use as a lower triangular preconditioner.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1060
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1072
double * value()
Access to C-style value array.
Definition: matrices.h:1084
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
Definition: matrices.h:386
unsigned nrow_local() const
access function for the num of local rows on this processor.
unsigned first_row() const
access function for the first row on this processor
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
double * values_pt()
access function to the underlying values
void clear()
wipes the DoubleVector
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
Base class for general purpose block preconditioners. Deals with setting subsidiary preconditioners a...
void set_dof_to_block_map(Vector< unsigned > &dof_to_block_map)
Specify a DOF to block map.
void set_subsidiary_preconditioner_function(SubsidiaryPreconditionerFctPt sub_prec_fn)
access function to set the subsidiary preconditioner function.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
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
void amg_using_simple_smoothing()
Function to select use of 'simple' AMG smoothers as controlled by the flag AMG_simple_smoother.
Definition: hypre_solver.h:977
double & amg_damping()
Access function to AMG_damping parameter.
unsigned & amg_coarsening()
Access function to AMG_coarsening flag.
unsigned & hypre_method()
Access function to Hypre_method flag – specified via enumeration.
Definition: hypre_solver.h:945
A preconditioner for performing inner iteration preconditioner solves. The template argument SOLVER s...
static OomphCommunicator * communicator_pt()
access to global communicator. This is the oomph-lib equivalent of MPI_COMM_WORLD
Matrix-based diagonal preconditioner.
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
An OomphLibError object which should be thrown when an run-time error is encountered....
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 const OomphCommunicator * comm_pt() const
Get function for comm pointer.
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 lagrange_multiplier_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the lagrange multiplier subsidiary preconditioner.
Preconditioner * Elastic_preconditioner_pt
storage for the preconditioner for the solid system
Vector< Preconditioner * > Lagrange_multiplier_preconditioner_pt
lagrange multiplier preconditioner pt
bool Use_inf_norm_of_s_scaling
boolean indicating whether the inf-norm of S should be used as scaling. Default = true;
SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt
The solid subsidiary preconditioner function pointer.
Elastic_preconditioner_type E_preconditioner_type
An unsigned indicating which method should be used for preconditioning the solid component.
void elastic_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the elastic subsidiary preconditioner.
double Scaling
The scaling. Defaults to infinity norm of S.
unsigned Dim
the dimension of the problem
SubsidiaryPreconditionerFctPt Lagrange_multiplier_subsidiary_preconditioner_function_pt
The Lagrange multiplier subsidary preconditioner function pointer.
Mesh * Elastic_mesh_pt
Pointer to the mesh containing the solid elements.
Mesh * Lagrange_multiplier_mesh_pt
Pointer to the mesh containing the Lagrange multiplier elements.
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
///////////////////////////////////////////////////////////////////////////// ///////////////////////...
void preconditioner_solve(const DoubleVector &res, DoubleVector &z)
Apply preconditioner to r.
SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt
The SubisidaryPreconditionerFctPt.
void use_upper_triangular_approximation()
Use as an upper triangular preconditioner.
Vector< PseudoElasticPreconditionerSubsidiaryPreconditionerOld * > Diagonal_block_preconditioner_pt
Vector of SuperLU preconditioner pointers for storing the preconditioners for each diagonal block.
void use_lower_triangular_approximation()
Use as a lower triangular preconditioner.
unsigned Method
the preconditioning method. 0 - block diagonal 1 - upper triangular 2 - lower triangular
void set_subsidiary_preconditioner_function(SubsidiaryPreconditionerFctPt sub_prec_fn)
access function to set the subsidiary preconditioner function.
DenseMatrix< MatrixVectorProduct * > Off_diagonal_matrix_vector_products
Matrix of matrix vector product operators for the off diagonals.
double & scaling()
Specify the scaling. Default is 1.0 Must be set before setup(...).
///////////////////////////////////////////////////////////////////////////// ///////////////////////...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the preconditioner.
SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt
the SubisidaryPreconditionerFctPt
void set_subsidiary_preconditioner_function(SubsidiaryPreconditionerFctPt sub_prec_fn)
access function to set the subsidiary preconditioner function.
double & scaling()
Specify the scaling. Default is 1.0 Must be called before setup(...).
SubsidiaryPreconditionerFctPt Lagrange_multiplier_subsidiary_preconditioner_function_pt
The Lagrange multiplier subsidiary preconditioner function pointer.
void lagrange_multiplier_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the lagrange multiplier subsidiary preconditioner.
bool Use_inf_norm_of_s_scaling
boolean indicating whether the inf-norm of S should be used as scaling. Default = true;
double Scaling
The scaling. Defaults to infinity norm of S.
Vector< Preconditioner * > Lagrange_multiplier_preconditioner_pt
lagrange multiplier preconditioner pt
void elastic_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the elastic subsidiary preconditioner.
Mesh * Lagrange_multiplier_mesh_pt
Pointer to the mesh containing the Lagrange multiplier elements.
SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt
The solid subsidiary preconditioner function pointer.
Mesh * Elastic_mesh_pt
Pointer to the mesh containing the solid elements.
Elastic_preconditioner_type E_preconditioner_type
An unsigned indicating which method should be used for preconditioning the solid component.
Preconditioner * Elastic_preconditioner_pt
storage for the preconditioner for the solid system
unsigned Dim
the dimension of the problem
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....
An interface to the Trilinos ML class - provides a function to construct a serial ML object,...
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
void start(const unsigned &i)
(Re-)start i-th timer
Preconditioner * get_lagrange_multiplier_preconditioner()
CG with diagonal preconditioner for the lagrange multiplier subsidiary linear systems.
Preconditioner * get_elastic_preconditioner()
AMG w/ GS smoothing for the augmented elastic subsidiary linear systems – calls Hypre version to stay...
Preconditioner * get_elastic_preconditioner_trilinos_ml()
TrilinosML smoothing for the augmented elastic subsidiary linear systems.
Preconditioner * get_elastic_preconditioner_hypre()
AMG w/ GS smoothing for the augmented elastic subsidiary linear systems.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...