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// Version 1.0; svn revision $LastChangedRevision: 1099 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-01-07 09:05:59 +0000 (Thu, 07 Jan 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
32 
33 namespace oomph
34 {
35  //=============================================================================
36  /// \short Functions to create instances of optimal subsidiary operators for
37  /// the PseudoElasticPreconditioner
38  //=============================================================================
39  namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper
40  {
41 #ifdef OOMPH_HAS_HYPRE
42 
43  /// \short AMG w/ GS smoothing for the augmented elastic subsidiary linear
44  /// systems
46  {
47  HyprePreconditioner* hypre_preconditioner_pt =
48  new HyprePreconditioner("Hypre for diagonal blocks in pseudo-solid");
49  hypre_preconditioner_pt->set_amg_iterations(2);
50  hypre_preconditioner_pt->amg_using_simple_smoothing();
51  if (MPI_Helpers::communicator_pt()->nproc()>1)
52  {
53  // Jacobi in parallel
54  hypre_preconditioner_pt->amg_simple_smoother() = 0;
55  }
56  else
57  {
58  // Gauss Seidel in serial (was 3 actually...)
59  hypre_preconditioner_pt->amg_simple_smoother() =1;
60  }
61  hypre_preconditioner_pt->hypre_method() = HyprePreconditioner::BoomerAMG;
62  hypre_preconditioner_pt->amg_damping() = 1.0;
63  hypre_preconditioner_pt->amg_coarsening() = 6;
64  return hypre_preconditioner_pt;
65  }
66 
67 
68  /// \short AMG w/ GS smoothing for the augmented elastic subsidiary linear
69  /// systems -- calls Hypre version to stay consistent with previous default
71  {
73  }
74 
75 #endif
76 
77 #ifdef OOMPH_HAS_TRILINOS
78 
79  /// \short TrilinosML smoothing for the augmented elastic
80  /// subsidiary linear systems
82  {
83  TrilinosMLPreconditioner* trilinos_prec_pt=
85  return trilinos_prec_pt;
86  }
87 
88  /// \short CG with diagonal preconditioner for the lagrange multiplier
89  /// subsidiary linear systems.
91  {
95  <TrilinosAztecOOSolver,MatrixBasedDiagPreconditioner>;
96  // Note: This makes CG a proper "inner iteration" for
97  // which GMRES (may) no longer converge. We should really
98  // use FGMRES or GMRESR for this. However, here the solver
99  // is so good that it'll converge very quickly anyway
100  // so there isn't much to be gained by limiting the number
101  // of iterations...
102  // prec_pt->max_iter() = 4;
103  prec_pt->solver_pt()->solver_type() = TrilinosAztecOOSolver::CG;
104  prec_pt->solver_pt()->disable_doc_time();
105  return prec_pt;
106  }
107 #endif
108  }
109 
110 
111 
112  ///////////////////////////////////////////////////////////////////////////////
113  ///////////////////////////////////////////////////////////////////////////////
114  ///////////////////////////////////////////////////////////////////////////////
115 
116 
117 
118  //=============================================================================
119  // Setup method for the PseudoElasticPreconditioner.
120  //=============================================================================
122  {
123  // clean
124  this->clean_up_memory();
125 
126 #ifdef PARANOID
127  // paranoid check that meshes have been set
128  if (Elastic_mesh_pt==0)
129  {
130  std::ostringstream error_message;
131  error_message << "The elastic mesh must be set.";
132  throw OomphLibError(error_message.str(),
133  OOMPH_CURRENT_FUNCTION,
134  OOMPH_EXCEPTION_LOCATION);
135 
136  }
138  {
139  std::ostringstream error_message;
140  error_message << "The Lagrange multiplier mesh must be set.";
141  throw OomphLibError(error_message.str(),
142  OOMPH_CURRENT_FUNCTION,
143  OOMPH_EXCEPTION_LOCATION);
144  }
145 #endif
146 
147  // set the meshes
148  this->set_nmesh(2);
149  this->set_mesh(0,Elastic_mesh_pt);
151 
152  // Figure out number of dof types
153  unsigned n_solid_dof_types = 0;
154  unsigned n_dof_types = 0;
155 
156  if (this->is_master_block_preconditioner())
157  {
158  // get the number of solid dof types from the first element
159  n_solid_dof_types = this->ndof_types_in_mesh(0);
160 
161  // get the total number of dof types
162  n_dof_types = n_solid_dof_types
163  + this->ndof_types_in_mesh(1);
164  }
165  else
166  {
167  n_dof_types = this->ndof_types();
168  n_solid_dof_types = n_dof_types - (n_dof_types/3);
169  }
170 #ifdef PARANOID
171  if (n_dof_types%3 != 0)
172  {
173  std::ostringstream error_message;
174  error_message << "This preconditioner requires DIM*3 types of DOF";
175  throw OomphLibError(error_message.str(),
176  OOMPH_CURRENT_FUNCTION,
177  OOMPH_EXCEPTION_LOCATION);
178  }
179 #endif
180 
181  // determine the dimension
182  Dim = n_dof_types/3;
183 
184 #ifdef PARANOID
185  // Recast Jacobian matrix to CRDoubleMatrix
186  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
187 
188  if (cr_matrix_pt==0)
189  {
190  std::ostringstream error_message;
191  error_message << "FSIPreconditioner only works with"
192  << " CRDoubleMatrix matrices" << std::endl;
193  throw OomphLibError(error_message.str(),
194  OOMPH_CURRENT_FUNCTION,
195  OOMPH_EXCEPTION_LOCATION);
196  }
197 #endif
198 
199  // Setup the dof_list scheme for block_setup.
200  // The natural DOF types order is (in 3D):
201  // 0 1 2 3 4 5 6 7 8 <- natural DOF type index
202  // x y z xc yc zc lx ly lz <- DOf type
203  //
204  // We need to group the directional displacements together:
205  // 0 1 2 3 4 5 6 7 8 <- block index
206  // x xc y yc xz zc lx ly lz <- DOF type
207  //
208  // The mapping required is:
209  // 0 2 4 1 3 5 6 7 8
210  //
211  // In general we want:
212  //
213  // dof_to_block_map[DOF type] = block type
214  Vector<unsigned> dof_list_for_block_setup(n_dof_types);
215  for (unsigned dim_i = 0; dim_i < Dim; dim_i++)
216  {
217  // bulk solid dof types
218  dof_list_for_block_setup[dim_i] = 2*dim_i;
219 
220  // constrained solid dof types
221  dof_list_for_block_setup[dim_i + Dim] = 2*dim_i + 1;
222 
223  // lagr dof types
224  dof_list_for_block_setup[dim_i + 2*Dim] = 2*Dim + dim_i;
225  }
226 
227  // Setup the block ordering. We have the following:
228  // Block types [0, 2*Dim) are the solid blocks.
229  // Block types [2*Dim, 3*Dim) are the Lagrange multiplier blocks.
230  //
231  // For d = 0, 1, 2,
232  // Bulk solid doftypes: 2*d
233  // Constrained solid doftypes: 2*d+1
234  // Lgr doftyoes: 2*Dim + d
235  this->block_setup(dof_list_for_block_setup);
236 
237  // Create dof_list for subsidiary preconditioner. This will be used later
238  // in turn_into_subsidiary_block_preconditioner(...)
239  Vector<unsigned> dof_list_for_subsidiary_prec(n_solid_dof_types);
240 
241  // The dof types are currently in the order (in 3D):
242  // 0 1 2 3 4 5 6 7 8 <- DOF index
243  // x y z xc yc zc lx ly lz <- DOF type
244  //
245  // We need to group the directional displacements together:
246  // x xc y yc xz zc
247  //
248  // The mapping required is:
249  // 0 3 1 4 2 5
250  //
251  // In general we want:
252  // dof_list[subsidiary DOF type] = master DOF type
253  for (unsigned dim_i = 0; dim_i < Dim; dim_i++)
254  {
255  // bulk solid dof
256  dof_list_for_subsidiary_prec[2*dim_i] = dim_i;
257 
258  // constrained solid dof
259  dof_list_for_subsidiary_prec[2*dim_i+1] = dim_i + Dim;
260  }
261 
262  // Get the solid blocks
263  DenseMatrix<CRDoubleMatrix*> solid_matrix_pt(n_solid_dof_types,
264  n_solid_dof_types,0);
265 
266  for (unsigned row_i = 0; row_i < n_solid_dof_types; row_i++)
267  {
268  for (unsigned col_i = 0; col_i < n_solid_dof_types; col_i++)
269  {
270  solid_matrix_pt(row_i,col_i) = new CRDoubleMatrix;
271  this->get_block(row_i,col_i,*solid_matrix_pt(row_i,col_i));
272  }
273  }
274 
275  // compute the scaling (if required)
277  {
278  Scaling = CRDoubleMatrixHelpers::inf_norm(solid_matrix_pt);
279  }
280  else
281  {
282  Scaling = 1.0;
283  }
284 
285  // Add the scaled identify matrix to the constrained solid blocks.
286  for(unsigned d = 0; d < Dim; d++)
287  {
288  // Where is the constrained block located?
289  unsigned block_i = 2*d+1;
290 
291  // Data from the constrained block.
292  double* s_values = solid_matrix_pt(block_i,block_i)->value();
293  int* s_column_index = solid_matrix_pt(block_i,block_i)->column_index();
294  int* s_row_start = solid_matrix_pt(block_i,block_i)->row_start();
295  int s_nrow_local = solid_matrix_pt(block_i,block_i)->nrow_local();
296  int s_first_row = solid_matrix_pt(block_i,block_i)->first_row();
297 
298  // Loop through the diagonal entries and add the scaling.
299  for (int i = 0; i < s_nrow_local; i++)
300  {
301  bool found = false;
302  for (int j = s_row_start[i];
303  j < s_row_start[i+1] && !found; j++)
304  {
305  if (s_column_index[j] == i + s_first_row)
306  {
307  s_values[j] += Scaling;
308  found = true;
309  } // if
310  } // for row start
311 
312  // Check if the diagonal entry is found.
313  if(!found)
314  {
315  std::ostringstream error_message;
316  error_message << "The diagonal entry for the constained block("
317  << block_i<<","<<block_i<<")\n"
318  << "on local row " << i << " does not exist."
319  << std::endl;
320  throw OomphLibError(error_message.str(),
321  OOMPH_CURRENT_FUNCTION,
322  OOMPH_EXCEPTION_LOCATION);
323  }
324  } // for nrow_local
325  } // for Dim
326 
327 
328  // setup the solid subsidiary preconditioner ////////////////////////////////
329  // this preconditioner uses the full S matrix
331  {
334 
335  // Add mesh (not actually used since this only acts as a subsidiary
336  // preconditioner...
337  //s_prec_pt->add_mesh(Elastic_mesh_pt);
338 
339  Vector<Vector<unsigned> > doftype_to_doftype_map(n_solid_dof_types,
340  Vector<unsigned>(1,0));
341 
342  for (unsigned i = 0; i < n_solid_dof_types; i++)
343  {
344  doftype_to_doftype_map[i][0] = i;
345  }
346 
348  this,dof_list_for_subsidiary_prec,doftype_to_doftype_map);
349 
351  {
352  s_prec_pt->
353  set_subsidiary_preconditioner_function
355  }
356 
357  // Set the replacement blocks.
358  for (unsigned d = 0; d < Dim; d++)
359  {
360  // The dof-block located in the block-ordering.
361  unsigned block_i = 2*d+1;
362 
363  // The dof-block located in the dof-ordering.
364  unsigned dof_block_i = Dim + d;
365  this->set_replacement_dof_block(dof_block_i,dof_block_i,
366  solid_matrix_pt(block_i,block_i));
367  }
368 
369  s_prec_pt->Preconditioner::setup(matrix_pt());
370  Elastic_preconditioner_pt = s_prec_pt;
371  }
372  // otherwise it is a block based preconditioner
373  else
374  {
376 
377  // set the block preconditioning method
378  switch (E_preconditioner_type)
379  {
381  {
383  }
384  break;
386  {
387  BlockTriangularPreconditioner<CRDoubleMatrix>* block_triangular_prec_pt
389  block_triangular_prec_pt->upper_triangular();
390 
391  s_prec_pt = block_triangular_prec_pt;
392  }
393  break;
395  {
396  BlockTriangularPreconditioner<CRDoubleMatrix>* block_triangular_prec_pt
398  block_triangular_prec_pt->lower_triangular();
399 
400  s_prec_pt = block_triangular_prec_pt;
401  }
402  break;
403  default:
404  {
405  std::ostringstream error_msg;
406  error_msg << "There is no such block based preconditioner.\n"
407  << "Candidates are:\n"
408  << "PseudoElasticPreconditioner::Block_diagonal_preconditioner\n"
409  << "PseudoElasticPreconditioner::Block_upper_triangular_preconditioner\n"
410  << "PseudoElasticPreconditioner::Block_lower_triangular_preconditioner\n";
411  throw OomphLibError(error_msg.str(),
412  OOMPH_CURRENT_FUNCTION,
413  OOMPH_EXCEPTION_LOCATION);
414  }
415  break;
416  }
417 
418 
419  // Add mesh (not actually used since this only acts as a subsidiary
420  // preconditioner...
421  //s_prec_pt->add_mesh(Elastic_mesh_pt);
422 
423  // The block to block map
424  Vector<Vector<unsigned> > doftype_to_doftype_map(
425  Dim,Vector<unsigned>(2,0));
426  Vector<unsigned> s_prec_dof_to_block_map(Dim,0);
427 
428  unsigned tmp_index = 0;
429  for (unsigned d = 0; d < Dim; d++)
430  {
431  s_prec_dof_to_block_map[d] = d;
432  doftype_to_doftype_map[d][0] = tmp_index++;
433  doftype_to_doftype_map[d][1] = tmp_index++;
434  }
435 
437  (this,dof_list_for_subsidiary_prec,doftype_to_doftype_map);
438 
440  {
443  }
444 
445  // Set the replacement blocks.
446  for (unsigned d = 0; d < Dim; d++)
447  {
448  // The dof-block located in the block-ordering.
449  unsigned block_i = 2*d+1;
450 
451  // Then dof-block located in the dof-ordering.
452  unsigned dof_block_i = Dim + d;
453  this->set_replacement_dof_block(dof_block_i,dof_block_i,
454  solid_matrix_pt(block_i,block_i));
455  }
456 
457  s_prec_pt->set_dof_to_block_map(s_prec_dof_to_block_map);
458  s_prec_pt->Preconditioner::setup(matrix_pt());
459 
460  Elastic_preconditioner_pt = s_prec_pt;
461  }
462 
463  // No longer require the solid blocks
464  for (unsigned row_i = 0; row_i < n_solid_dof_types; row_i++)
465  {
466  for (unsigned col_i = 0; col_i < n_solid_dof_types; col_i++)
467  {
468  delete solid_matrix_pt(row_i,col_i); solid_matrix_pt(row_i,col_i) = 0;
469  }
470  }
471 
472  // next setup the lagrange multiplier preconditioners
474  for (unsigned d = 0; d < Dim; d++)
475  {
476  CRDoubleMatrix* b_pt = new CRDoubleMatrix;
477  this->get_block(2*Dim+d,2*d+1,*b_pt);
478 
479  // if a non default preconditioner is specified create
480  // the preconditioners
482  {
484  (*Lagrange_multiplier_subsidiary_preconditioner_function_pt)();
485  }
486  // else use default superlu preconditioner
487  else
488  {
490  }
491 
492  // and setup
494  delete b_pt; b_pt = 0;
495  }
496 
497  }
498 
499  //=============================================================================
500  /// \short Apply the elastic subsidiary preconditioner.
501  //=============================================================================
503  (const DoubleVector& r, DoubleVector& z)
504  {
505  // apply the solid preconditioner
506  Elastic_preconditioner_pt->preconditioner_solve(r,z);
507  }
508 
509  //=============================================================================
510  /// \short Apply the lagrange multiplier subsidiary preconditioner.
511  //=============================================================================
514  DoubleVector& z)
515  {
516  // apply the lagrange multiplier preconditioner
517  for (unsigned d = 0; d < Dim; d++)
518  {
519  DoubleVector x;
520  this->get_block_vector(Dim*2+d,r,x);
521  DoubleVector y;
522  Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(x,y);
523  Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(y,x);
524  unsigned nrow_local = x.nrow_local();
525  double* x_pt = x.values_pt();
526  for (unsigned i = 0; i < nrow_local; i++)
527  {
528  x_pt[i] = x_pt[i] * Scaling;
529  }
530  this->return_block_vector(Dim*2+d,x,z);
531  }
532  }
533 
534  //=============================================================================
535  /// \short Clears the memory.
536  //=============================================================================
538  {
539  // clean the block preconditioner base class memory
541 
542  // delete the solid preconditioner
544 
545  // delete the lagrange multiplier preconditioner pt
546  unsigned sz = Lagrange_multiplier_preconditioner_pt.size();
547  for (unsigned i = 0; i < sz; i++)
548  {
551  }
552  }
553 
554  ///////////////////////////////////////////////////////////////////////////////
555  ///////////////////////////////////////////////////////////////////////////////
556  ///////////////////////////////////////////////////////////////////////////////
557 
558  //=============================================================================
559  // Setup method for the PseudoElasticPreconditioner.
560  //=============================================================================
562  {
563  // clean
564  this->clean_up_memory();
565 
566 #ifdef PARANOID
567  // paranoid check that meshes have been set
568  if (Elastic_mesh_pt==0)
569  {
570  std::ostringstream error_message;
571  error_message << "The elastic mesh must be set.";
572  throw OomphLibError(error_message.str(),
573  OOMPH_CURRENT_FUNCTION,
574  OOMPH_EXCEPTION_LOCATION);
575 
576  }
578  {
579  std::ostringstream error_message;
580  error_message << "The Lagrange multiplier mesh must be set.";
581  throw OomphLibError(error_message.str(),
582  OOMPH_CURRENT_FUNCTION,
583  OOMPH_EXCEPTION_LOCATION);
584  }
585 #endif
586 
587  // set the mesh
588  unsigned n_solid_dof_types = 0;
589  unsigned n_dof_types = 0;
590  this->set_mesh(0,Elastic_mesh_pt);
592  if (this->is_master_block_preconditioner())
593  {
594 
595  // get the number of solid dof types from the first element
596  n_solid_dof_types = this->ndof_types_in_mesh(0);
597 
598  // get the total number of dof types
599  n_dof_types = n_solid_dof_types
600  + this->ndof_types_in_mesh(1);
601  }
602  else
603  {
604  n_dof_types = this->ndof_types();
605  n_solid_dof_types = n_dof_types - (n_dof_types/3);
606  }
607 #ifdef PARANOID
608  if (n_dof_types%3 != 0)
609  {
610  std::ostringstream error_message;
611  error_message << "This preconditioner requires DIM*3 types of DOF";
612  throw OomphLibError(error_message.str(),
613  OOMPH_CURRENT_FUNCTION,
614  OOMPH_EXCEPTION_LOCATION);
615  }
616 #endif
617 
618  // determine the dimension
619  Dim = n_dof_types/3;
620 
621  // Recast Jacobian matrix to CRDoubleMatrix
622  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
623 
624 #ifdef PARANOID
625  if (cr_matrix_pt==0)
626  {
627  std::ostringstream error_message;
628  error_message << "FSIPreconditioner only works with"
629  << " CRDoubleMatrix matrices" << std::endl;
630  throw OomphLibError(error_message.str(),
631  OOMPH_CURRENT_FUNCTION,
632  OOMPH_EXCEPTION_LOCATION);
633  }
634 #endif
635 
636  // Setup the block ordering.
637  this->block_setup();
638 
639  // Create dof_list for subsidiary preconditioner. This will be used later
640  // in turn_into_subsidiary_block_preconditioner(...)
641  Vector<unsigned> dof_list_for_subsidiary_prec(n_solid_dof_types);
642  for (unsigned i = 0; i < n_solid_dof_types; i++)
643  {
644  dof_list_for_subsidiary_prec[i] = i;
645  }
646 
647  // compute the scaling (if required)
649  {
650  Vector<unsigned> dof_list(n_solid_dof_types);
651  for (unsigned i = 0; i < n_solid_dof_types; i++)
652  {
653  dof_list[i] = i;
654  }
655 
658  cr_matrix_pt,
659  dof_list,
661  comm_pt());
662  Scaling = helper_pt->s_inf_norm();
663  delete helper_pt; helper_pt = 0;
664  }
665  else
666  {
667  Scaling = 1.0;
668  }
669 
670 
671  // setup the solid subsidiary preconditioner ////////////////////////////////
672  // this preconditioner uses the full S matrix
674  {
677  Vector<unsigned> dof_list(n_solid_dof_types);
678  for (unsigned i = 0; i < n_solid_dof_types; i++)
679  {
680  dof_list[i] = i;
681  }
682  s_prec_pt->turn_into_subsidiary_block_preconditioner(this,dof_list);
684  {
685  s_prec_pt->
686  set_subsidiary_preconditioner_function
688  }
689  s_prec_pt->scaling() = Scaling;
690  s_prec_pt->Preconditioner::setup(matrix_pt());
691  Elastic_preconditioner_pt = s_prec_pt;
692  }
693 
694  // otherwise it is a block based preconditioner
695  else
696  {
697  // create the preconditioner
699  s_prec_pt =
701  Vector<unsigned> dof_list(n_solid_dof_types);
702  for (unsigned i = 0; i < n_solid_dof_types; i++)
703  {
704  dof_list[i] = i;
705  }
706  s_prec_pt->turn_into_subsidiary_block_preconditioner(this,dof_list);
707 
708  // set the subsidiary solve method
710  {
711  s_prec_pt->
712  set_subsidiary_preconditioner_function
714  }
715 
716  // set the scaling
717  s_prec_pt->scaling() = Scaling;
718 
719  // set the block preconditioning method
720  switch (E_preconditioner_type)
721  {
724  break;
727  break;
730  break;
731  default:
732  break;
733  }
734 
735  // setup
736  s_prec_pt->Preconditioner::setup(matrix_pt());
737  Elastic_preconditioner_pt = s_prec_pt;
738 
739  }
740 
741  // next setup the lagrange multiplier preconditioners
743  for (unsigned d = 0; d < Dim; d++)
744  {
745  CRDoubleMatrix* b_pt = new CRDoubleMatrix;
746  this->get_block(2*Dim+d,Dim+d,*b_pt);
747 
748  // if a non default preconditioner is specified create
749  // the preconditioners
751  {
753  (*Lagrange_multiplier_subsidiary_preconditioner_function_pt)();
754  }
755 
756  // else use default superlu preconditioner
757  else
758  {
760  }
761 
762  // and setup
764  delete b_pt; b_pt = 0;
765  }
766 
767  }
768 
769  //=============================================================================
770  /// \short Apply the elastic subsidiary preconditioner.
771  //=============================================================================
773  (const DoubleVector& r, DoubleVector& z)
774  {
775  // apply the solid preconditioner
776  Elastic_preconditioner_pt->preconditioner_solve(r,z);
777  }
778 
779 
780  //=============================================================================
781  /// \short Apply the lagrange multiplier subsidiary preconditioner.
782  //=============================================================================
785  DoubleVector& z)
786  {
787  // apply the lagrange multiplier preconditioner
788  for (unsigned d = 0; d < Dim; d++)
789  {
790  DoubleVector x;
791  this->get_block_vector(Dim*2+d,r,x);
792  DoubleVector y;
793  Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(x,y);
794  Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(y,x);
795  unsigned nrow_local = x.nrow_local();
796  double* x_pt = x.values_pt();
797  for (unsigned i = 0; i < nrow_local; i++)
798  {
799  x_pt[i] = x_pt[i] * Scaling;
800  }
801  this->return_block_vector(Dim*2+d,x,z);
802  }
803  }
804 
805  //=============================================================================
806  /// \short Clears the memory.
807  //=============================================================================
809  {
810  // clean the block preconditioner base class memory
812 
813  // delete the solid preconditioner
816 
817  // delete the lagrange multiplier preconditioner pt
818  unsigned sz = Lagrange_multiplier_preconditioner_pt.size();
819  for (unsigned i = 0; i < sz; i++)
820  {
823  }
824  }
825 
826 
827 
828  ///////////////////////////////////////////////////////////////////////////////
829  ///////////////////////////////////////////////////////////////////////////////
830  ///////////////////////////////////////////////////////////////////////////////
831 
832 
833 
834  //=============================================================================
835  /// \short Setup the preconditioner
836  //=============================================================================
839  {
840  // clean memory
841  this->clean_up_memory();
842 
843 #ifdef PARANOID
844  // paranoid check that this preconditioner has an even number of DOF types
845  if (this->ndof_types()%2 != 0)
846  {
847  std::ostringstream error_message;
848  error_message
849  << "This SUBSIDIARY preconditioner requires an even number of "
850  << "types of DOF";
851  throw OomphLibError(
852  error_message.str(),
853  OOMPH_CURRENT_FUNCTION,
854  OOMPH_EXCEPTION_LOCATION);
855  }
856 #endif
857 
858  // assemble dof_to_block_map
859  unsigned ndof_types = this->ndof_types();
860  Vector<unsigned> dof_to_block_map(ndof_types,0);
861  for (unsigned i = ndof_types/2; i < ndof_types; i++)
862  {
863  dof_to_block_map[i] = 1;
864  }
865 
866  this->block_setup(dof_to_block_map);
867 
868  // get block 11
869  CRDoubleMatrix* s11_pt = new CRDoubleMatrix;
870  this->get_block(1,1,*s11_pt);
871 
872  // add the scaled identity matrix to block 11
873  double* s11_values = s11_pt->value();
874  int* s11_column_index = s11_pt->column_index();
875  int* s11_row_start = s11_pt->row_start();
876  int s11_nrow_local = s11_pt->nrow_local();
877  int s11_first_row = s11_pt->first_row();
878  for (int i = 0; i < s11_nrow_local; i++)
879  {
880  bool found = false;
881  for (int j = s11_row_start[i];
882  j < s11_row_start[i+1] && !found; j++)
883  {
884  if (s11_column_index[j] == i + s11_first_row)
885  {
886  s11_values[j] += Scaling;
887  found = true;
888  }
889  }
890  }
891 
892  VectorMatrix<BlockSelector> required_blocks(2,2);
893  const bool want_block = true;
894  for (unsigned b_i = 0; b_i < 2; b_i++)
895  {
896  for (unsigned b_j = 0; b_j < 2; b_j++)
897  {
898  required_blocks[b_i][b_j].select_block(b_i,b_j,want_block);
899  }
900  }
901 
902  required_blocks[1][1].set_replacement_block_pt(s11_pt);
903 
904  CRDoubleMatrix s_prec_pt = this->get_concatenated_block(required_blocks);
905 
906  delete s11_pt; s11_pt = 0;
907 
908  // setup the preconditioner
910  {
911  Preconditioner_pt = (*Subsidiary_preconditioner_function_pt)();
912  }
913  else
914  {
916  }
917  Preconditioner_pt->setup(&s_prec_pt);
918  }
919 
920  //=============================================================================
921  /// \short Apply the preconditioner.
922  //=============================================================================
925  {
926  DoubleVector x;
928  DoubleVector y;
931  }
932 
933 
934 
935  ///////////////////////////////////////////////////////////////////////////////
936  ///////////////////////////////////////////////////////////////////////////////
937  ///////////////////////////////////////////////////////////////////////////////
938 
939 
940 
941  //=============================================================================
942  /// clean up the memory
943  //=============================================================================
946  {
947  //number of block types
948  unsigned n_block = Diagonal_block_preconditioner_pt.size();
949 
950  //delete diagonal blocks
951  for (unsigned i = 0 ; i < n_block; i++)
952  {
955  if (Method == 1)
956  {
957  for (unsigned j = i+1; j < n_block; j++)
958  {
961  }
962  }
963  else if (Method == 2)
964  {
965  for (unsigned j = 0; j < i; j++)
966  {
969  }
970  }
971  }
972 
973  // clean up the block preconditioner
975  }
976 
977  //=============================================================================
978  /// \short Setup the preconditioner.
979  //=============================================================================
982  {
983  // clean the memory
984  this->clean_up_memory();
985 
986  // determine the number of DOF types
987  unsigned n_dof_types = this->ndof_types();
988 
989 #ifdef PARANOID
990  // must be Dim*2 dof types
991  if (n_dof_types%2 != 0)
992  {
993  std::ostringstream error_message;
994  error_message << "This preconditioner requires DIM*3 types of DOF";
995  throw OomphLibError(error_message.str(),
996  OOMPH_CURRENT_FUNCTION,
997  OOMPH_EXCEPTION_LOCATION);
998  }
999 #endif
1000 
1001  // store the dimension of the problem
1002  unsigned dim = n_dof_types/2;
1003 
1004  // assemble the dof to block lookup scheme
1005  Vector<unsigned> dof_to_block_map(n_dof_types,0);
1006  for (unsigned d = 0; d < dim; d++)
1007  {
1008  dof_to_block_map[d] = d;
1009  dof_to_block_map[d+dim] = d;
1010  }
1011 
1012  //setup the blocks look up schemes
1013  this->block_setup(dof_to_block_map);
1014 
1015  // Storage for the diagonal block preconditioners
1017 
1018  // storage for the off diagonal matrix vector products
1019  Off_diagonal_matrix_vector_products.resize(dim,dim,0);
1020 
1021  // setup the subsidiary preconditioners
1022  for (unsigned d = 0; d < dim; d++)
1023  {
1024  Vector<unsigned> dof_list(2);
1025  dof_list[0]=d;
1026  dof_list[1]=d+dim;
1027 
1033  {
1037  }
1038  Diagonal_block_preconditioner_pt[d]->scaling() = Scaling;
1039 
1042 
1043  // the preconditioning method.\n
1044  // 0 - block diagonal\n
1045  // 1 - upper triangular\n
1046  // 2 - lower triangular\n
1047  // next setup the off diagonal mat vec operators if required
1048  if (Method == 1 || Method == 2)
1049  {
1050  unsigned l = d+1;
1051  unsigned u = dim;
1052  if (Method==2)
1053  {
1054  l = 0;
1055  u = d;
1056  }
1057  for (unsigned j = l; j < u; j++)
1058  {
1059  CRDoubleMatrix* block_matrix_pt = new CRDoubleMatrix;
1060  this->get_block(d,j,*block_matrix_pt);
1062  = new MatrixVectorProduct();
1063 // Off_diagonal_matrix_vector_products(d,j)->setup(block_matrix_pt);
1065  block_matrix_pt,j);
1066 
1067  delete block_matrix_pt; block_matrix_pt = 0;
1068  }
1069  }
1070  } // setup the subsidiary preconditioner.
1071  } // preconditioner setup
1072 
1073  //=============================================================================
1074  /// Apply preconditioner to r
1075  //=============================================================================
1078  {
1079  // copy r
1080  DoubleVector r(res);
1081 
1082  unsigned n_block;
1083 
1084  // Cache umber of block types (also the spatial DIM)
1085  n_block = this->nblock_types();
1086 
1087  // loop parameters
1088  int start = n_block-1;
1089  int end = -1;
1090  int step = -1;
1091  if (Method != 1)
1092  {
1093  start = 0;
1094  end = n_block;
1095  step = 1;
1096  }
1097 
1098  // the preconditioning method.
1099  // 0 - block diagonal
1100  // 1 - upper triangular
1101  // 2 - lower triangular
1102  //
1103  // loop over the DIM
1104  //
1105  // For Method = 0 or 2 (diagonal, lower)
1106  // start = 2, end = -1, step = -1
1107  // i = 2,1,0
1108  //
1109  // For Method = 1 (upper)
1110  // start = 0, end = 3 step = 1
1111  // i = 0, 1, 2
1112  for (int i = start; i != end; i+=step)
1113  {
1114 
1115  // solve
1116  Diagonal_block_preconditioner_pt[i]->preconditioner_solve(r,z);
1117 
1118  // if upper or lower triangular
1119  if (Method != 0)
1120  {
1121 
1122  // substitute
1123  //
1124  for (int j = i + step; j !=end; j+=step)
1125  {
1126  DoubleVector x;
1127  this->get_block_vector(i,z,x);
1128  DoubleVector y;
1129  Off_diagonal_matrix_vector_products(j,i)->multiply(x,y);
1130  x.clear();
1131  this->get_block_vector(j,r,x);
1132  x -= y;
1133  this->return_block_vector(j,x,r);
1134  } // substitute
1135  } // if upper or lower
1136  } // for loop over DIM
1137  } // Block preconditioner solve
1138 } // namespace oomph
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1056
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
An interface to the Trilinos ML class - provides a function to construct a serial ML object...
Elastic_preconditioner_type E_preconditioner_type
An unsigned indicating which method should be used for preconditioning the solid component.
unsigned & amg_simple_smoother()
Access function to AMG_simple_smoother flag.
Definition: hypre_solver.h:888
Preconditioner * Elastic_preconditioner_pt
storage for the preconditioner for the solid system
unsigned nrow_local() const
access function for the num of local rows on this processor.
void clear_block_preconditioner_base()
Clears all BlockPreconditioner data. Called by the destructor and the block_setup(...) methods.
A helper class for PseudoElasticPreconditioner. Note that this is NOT actually a functioning precondi...
bool Use_inf_norm_of_s_scaling
boolean indicating whether the inf-norm of S should be used as scaling. Default = true; ...
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
bool Use_inf_norm_of_s_scaling
boolean indicating whether the inf-norm of S should be used as scaling. Default = true; ...
void clear()
wipes the DoubleVector
unsigned first_row() const
access function for the first row on this processor
double & amg_damping()
Access function to AMG_damping parameter.
Definition: hypre_solver.h:916
double * values_pt()
access function to the underlying values
Mesh * Elastic_mesh_pt
Pointer to the mesh containing the solid elements.
void set_subsidiary_preconditioner_function(SubsidiaryPreconditionerFctPt sub_prec_fn)
access function to set the subsidiary preconditioner function.
cstr elem_len * i
Definition: cfortran.h:607
void set_subsidiary_preconditioner_function(SubsidiaryPreconditionerFctPt sub_prec_fn)
access function to set the subsidiary preconditioner function.
Mesh * Elastic_mesh_pt
Pointer to the mesh containing the solid elements.
void elastic_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the elastic subsidiary preconditioner.
void lagrange_multiplier_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the lagrange multiplier subsidiary preconditioner.
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 ...
unsigned & amg_coarsening()
Access function to AMG_coarsening flag.
Definition: hypre_solver.h:910
SubsidiaryPreconditionerFctPt Lagrange_multiplier_subsidiary_preconditioner_function_pt
The Lagrange multiplier subsidary preconditioner function pointer.
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...
void setup()
Broken assignment operator.
static OomphCommunicator * communicator_pt()
access to the global oomph-lib communicator
SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt
The solid subsidiary preconditioner function pointer.
void lower_triangular()
Use as a lower triangular preconditioner.
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, block_pt is a pointer to the block matrix, block_col_indices is a vector indicating which block indices does the RHS vector we want to multiply the matrix by.
An interface to the Trilinos AztecOO classes allowing it to be used as an Oomph-lib LinearSolver...
bool is_master_block_preconditioner() const
Return true if this preconditioner is the master block preconditioner.
unsigned Dim
the dimension of the problem
void set_amg_iterations(const unsigned &amg_iterations)
Function to set the number of times to apply BoomerAMG.
Definition: hypre_solver.h:872
double * value()
Access to C-style value array.
Definition: matrices.h:1062
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...
Matrix-based diagonal preconditioner.
double Scaling
The scaling. Defaults to infinity norm of S.
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 ...
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 amg_using_simple_smoothing()
Function to select use of 'simple' AMG smoothers as controlled by the flag AMG_simple_smoother.
Definition: hypre_solver.h:885
void elastic_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the elastic subsidiary preconditioner.
Vector< Preconditioner * > Lagrange_multiplier_preconditioner_pt
lagrange multiplier preconditioner pt
Preconditioner * get_elastic_preconditioner_hypre()
AMG w/ GS smoothing for the augmented elastic subsidiary linear systems.
Mesh * Lagrange_multiplier_mesh_pt
Pointer to the mesh containing the Lagrange multiplier elements.
Preconditioner * get_elastic_preconditioner()
AMG w/ GS smoothing for the augmented elastic subsidiary linear systems – calls Hypre version to stay...
void start(const unsigned &i)
(Re-)start i-th timer
SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt
The solid subsidiary preconditioner function pointer.
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
A preconditioner for performing inner iteration preconditioner solves. The template argument SOLVER s...
double & scaling()
Specify the scaling. Default is 1.0 Must be set before setup(...).
void preconditioner_solve(const DoubleVector &res, DoubleVector &z)
Apply preconditioner to r.
Vector< Preconditioner * > Lagrange_multiplier_preconditioner_pt
lagrange multiplier preconditioner pt
void set_nmesh(const unsigned &n)
Specify the number of meshes required by this block preconditioner. Note: elements in different meshe...
DenseMatrix< MatrixVectorProduct * > Off_diagonal_matrix_vector_products
Matrix of matrix vector product operators for the off diagonals.
unsigned & hypre_method()
Access function to Hypre_method flag – specified via enumeration.
Definition: hypre_solver.h:862
SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt
the SubisidaryPreconditionerFctPt
void set_mesh(const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Set the i-th mesh for this block preconditioner. Note: The method set_nmesh(...) must be called befor...
void get_block(const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const
Put block (i,j) into output_matrix. This block accounts for any coarsening of dof types and any repla...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the preconditioner.
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.
unsigned ndof_types() const
Return the total number of DOF types.
An interface to allow SuperLU to be used as an (exact) Preconditioner.
void use_lower_triangular_approximation()
Use as a lower triangular preconditioner.
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
double & scaling()
Specify the scaling. Default is 1.0 Must be called before setup(...).
Block diagonal preconditioner. By default SuperLU is used to solve the subsidiary systems...
void upper_triangular()
Use as an upper triangular preconditioner.
void use_upper_triangular_approximation()
Use as an upper triangular preconditioner.
unsigned Dim
the dimension of the problem
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 "...
Preconditioner * get_elastic_preconditioner_trilinos_ml()
TrilinosML smoothing for the augmented elastic subsidiary linear systems.
double Scaling
The scaling. Defaults to infinity norm of S.
SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt
The SubisidaryPreconditionerFctPt.
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
Class for dense matrices, storing all the values of the matrix as a pointer to a pointer with assorte...
Definition: communicator.h:50
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
void set_dof_to_block_map(Vector< unsigned > &dof_to_block_map)
Specify a DOF to block map.
General purpose block triangular preconditioner By default this is Upper triangular. By default SuperLUPreconditioner (or SuperLUDistPreconditioner) is used to solve the subsidiary systems, but other preconditioners can be used by setting them using passing a pointer to a function of type SubsidiaryPreconditionerFctPt to the method subsidiary_preconditioner_function_pt().
CRDoubleMatrix * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*...
double inf_norm(const DenseMatrix< CRDoubleMatrix * > &matrix_pt)
Compute infinity (maximum) norm of sub blocks as if it was one matrix.
Definition: matrices.cc:3708
Elastic_preconditioner_type E_preconditioner_type
An unsigned indicating which method should be used for preconditioning the solid component.
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1050
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...
unsigned nblock_types() const
Return the number of block types.
Mesh * Lagrange_multiplier_mesh_pt
Pointer to the mesh containing the Lagrange multiplier elements.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
SubsidiaryPreconditionerFctPt Lagrange_multiplier_subsidiary_preconditioner_function_pt
The Lagrange multiplier subsidiary preconditioner function pointer.
Preconditioner * get_lagrange_multiplier_preconditioner()
CG with diagonal preconditioner for the lagrange multiplier subsidiary linear systems.
CRDoubleMatrix get_concatenated_block(const VectorMatrix< BlockSelector > &selected_block)
Returns a concatenation of the block matrices specified by the argument selected_block. The VectorMatrix selected_block must be correctly sized as it is used to determine the number of sub block matrices to concatenate.
Vector< PseudoElasticPreconditionerSubsidiaryPreconditionerOld * > Diagonal_block_preconditioner_pt
Vector of SuperLU preconditioner pointers for storing the preconditioners for each diagonal block...
void use_block_diagonal_approximation()
use as a block diagonal preconditioner