solid_preconditioners.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: 1162 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-04-18 13:27:54 +0100 (Mon, 18 Apr 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 #include "solid_preconditioners.h"
31 
32 
33 namespace oomph
34 {
35 
36 
37 
38 
39 //===========================================================================
40 /// Setup the least-squares commutator solid preconditioner. This
41 /// extracts blocks corresponding to the position/displacement and pressure
42 /// unknowns, creates the matrices actually needed in the application of the
43 /// preconditioner and deletes what can be deleted... Note that
44 /// this preconditioner needs a CRDoubleMatrix.
45 //============================================================================
48  {
49 
50  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51  // NOTE: In the interest of minimising memory usage, several containers
52  // are recycled, therefore their content/meaning changes
53  // throughout this function. The code is carefully annotated
54  // but you'll have to read it line by line!
55  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 
57  // make sure any old data is deleted
59 
60 #ifdef PARANOID
61  // paranoid check that the solid mesh pt has been set
62  if (Solid_mesh_pt == 0)
63  {
64  std::ostringstream error_message;
65  error_message << "The solid mesh pointer must be set.\n"
66  << "Use method set_solid_mesh(...)";
67  throw OomphLibError(error_message.str(),
68  OOMPH_CURRENT_FUNCTION,
69  OOMPH_EXCEPTION_LOCATION);
70  }
71 #endif
72 
73  // set the mesh
74  this->set_mesh(0,Solid_mesh_pt);
75 
76  // Get blocks
77  // ----------
78 
79  // Set up block look up schemes (done automatically in the
80  // BlockPreconditioner base class, based on the information
81  // provided in the block-preconditionable elements in the problem)
82 
83  // this preconditioner has two types of block:
84  // type 0: displacement/positions - corresponding to DOFs 0 to n-2
85  // type 1: pressure - corresponding to DOF n-1
86  double t_block_start = TimingHelpers::timer();
87  unsigned ndof_types = 0;
89  {
90  ndof_types = this->ndof_types();
91  }
92  else
93  {
94  ndof_types = this->ndof_types_in_mesh(0);
95  }
96  Vector<unsigned> dof_to_block_map(ndof_types);
97  dof_to_block_map[ndof_types-1]=1;
98  this->block_setup(dof_to_block_map);
99  double t_block_finish = TimingHelpers::timer();
100  double block_setup_time = t_block_finish - t_block_start;
101  if(Doc_time)
102  {
103  oomph_info << "Time for block_setup(...) [sec]: "
104  << block_setup_time << "\n";
105  }
106 
107  // determine whether the F preconditioner is a block preconditioner (and
108  // therefore a subsidiary preconditioner)
109  BlockPreconditioner<CRDoubleMatrix>* F_block_preconditioner_pt =
112  if (F_block_preconditioner_pt == 0)
113  {
114  F_preconditioner_is_block_preconditioner = false;
115  }
116 
117  // Get B (the divergence block)
118  double t_get_B_start = TimingHelpers::timer();
119  CRDoubleMatrix* b_pt = new CRDoubleMatrix;
120  this->get_block(1,0,*b_pt);
121  double t_get_B_finish = TimingHelpers::timer();
122  if(Doc_time)
123  {
124  double get_B_time = t_get_B_finish - t_get_B_start;
125  oomph_info << "Time to get B [sec]: "
126  << get_B_time << "\n";
127  }
128 
129  // the pointer for f
130  CRDoubleMatrix* f_pt = new CRDoubleMatrix;
131 
132  // the pointer for the p matrix
133  CRDoubleMatrix* p_matrix_pt = new CRDoubleMatrix;
134 
135  // the pointer for bt
136  CRDoubleMatrix* bt_pt = new CRDoubleMatrix;
137 
138  // if BFBt is to be formed
139  if (Form_BFBt_product)
140  {
141 
142  // If using scaling then replace B with Bt
143  CRDoubleMatrix* ivmm_pt = 0;
145  {
146 
147  // Assemble the ivmm diagonal
148  double ivmm_assembly_start_t = TimingHelpers::timer();
149  ivmm_pt = assemble_mass_matrix_diagonal();
150  double ivmm_assembly_finish_t = TimingHelpers::timer();
151  if (Doc_time)
152  {
153  double
154  ivmm_assembly_time = ivmm_assembly_finish_t - ivmm_assembly_start_t;
155  oomph_info << "Time to assemble inverse mass matrix [sec]: "
156  << ivmm_assembly_time << "\n";
157  }
158 
159  // assemble BQ (stored in B)
160  double t_BQ_start = TimingHelpers::timer();
161  CRDoubleMatrix* temp_matrix_pt = new CRDoubleMatrix;
162  b_pt->multiply(*ivmm_pt,*temp_matrix_pt);
163  delete b_pt; b_pt = 0;
164  b_pt = temp_matrix_pt;
165  double t_BQ_finish = TimingHelpers::timer();
166  if(Doc_time)
167  {
168  double t_BQ_time = t_BQ_finish - t_BQ_start;
169  oomph_info << "Time to generate BQ [sec]: "
170  << t_BQ_time << std::endl;
171  }
172  }
173 
174  // Get Bt
175  double t_get_Bt_start = TimingHelpers::timer();
176  this->get_block(0,1,*bt_pt);
177  double t_get_Bt_finish = TimingHelpers::timer();
178  if(Doc_time)
179  {
180  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
181  oomph_info << "Time to get Bt [sec]: "
182  << t_get_Bt_time << std::endl;
183  }
184 
185  // now form the P matrix by multiplying B (which if using scaling will be
186  // BQ) with Bt
187  double t_P_start = TimingHelpers::timer();
188  b_pt->multiply(*bt_pt,*p_matrix_pt);
189  double t_P_finish = TimingHelpers::timer();
190  if(Doc_time)
191  {
192  double t_P_time = t_P_finish - t_P_start;
193  oomph_info << "Time to generate P matrix [sec]: "
194  << t_P_time << std::endl;
195  }
196 
197  // Multiply auxiliary matrix by diagonal of mass matrix if
198  // required
200  {
201  CRDoubleMatrix* temp_matrix_pt = new CRDoubleMatrix;
202  double t_QBt_start = TimingHelpers::timer();
203  ivmm_pt->multiply(*bt_pt,*temp_matrix_pt);
204  delete bt_pt; bt_pt = 0;
205  bt_pt = temp_matrix_pt;
206  double t_QBt_finish = TimingHelpers::timer();
207 
208  // Output times
209  if(Doc_time)
210  {
211  double t_QBt_time = t_QBt_finish - t_QBt_start;
212  oomph_info << "Time to generate QBt [sec]: "
213  << t_QBt_time << std::endl;
214  }
215  }
216 
217  // Clean up memory
218  delete ivmm_pt;
219 
220  // get block 0 0
221  double t_get_F_start = TimingHelpers::timer();
222  this->get_block(0,0,*f_pt);
223  double t_get_F_finish = TimingHelpers::timer();
224  if(Doc_time)
225  {
226  double t_get_F_time = t_get_F_finish - t_get_F_start;
227  oomph_info << "Time to get F [sec]: "
228  << t_get_F_time << std::endl;
229  }
230 
231  // Auxiliary matrix for intermediate results
232  double t_aux_matrix_start = TimingHelpers::timer();
233  CRDoubleMatrix* aux_matrix_pt = new CRDoubleMatrix;
234  f_pt->multiply(*bt_pt, *aux_matrix_pt);
235  double t_aux_matrix_finish = TimingHelpers::timer();
236  if(Doc_time)
237  {
238  double t_aux_time = t_aux_matrix_finish - t_aux_matrix_start;
239  oomph_info << "Time to generate FQBt [sec]: "
240  << t_aux_time << std::endl;
241  }
242 
243  // can now delete Bt (or QBt if using scaling)
244  delete bt_pt;
245  bt_pt = 0;
246 
247  // and if F requires a block preconditioner then we can delete F
248  if (F_preconditioner_is_block_preconditioner)
249  {
250  delete f_pt;
251  }
252 
253  // now form BFBt
254  double t_E_matrix_start = TimingHelpers::timer();
255  CRDoubleMatrix* e_matrix_pt = new CRDoubleMatrix;
256  b_pt->multiply(*aux_matrix_pt,*e_matrix_pt);
257  delete aux_matrix_pt;
258  delete b_pt;
259  double t_E_matrix_finish = TimingHelpers::timer();
260  if(Doc_time)
261  {
262  double t_E_time = t_E_matrix_finish - t_E_matrix_start;
263  oomph_info << "Time to generate E (B*(F*Bt)) [sec]: "
264  << t_E_time << std::endl;
265  }
266  double t_E_matvec_start = TimingHelpers::timer();
268  //E_mat_vec_pt->setup(e_matrix_pt);
269  this->setup_matrix_vector_product(E_mat_vec_pt,e_matrix_pt,1);
270  double t_E_matvec_finish = TimingHelpers::timer();
271  delete e_matrix_pt;
272  if(Doc_time)
273  {
274  double t_E_time = t_E_matvec_finish - t_E_matvec_start;
275  oomph_info << "Time to build E (BFBt) matrix vector operator E [sec]: "
276  << t_E_time << std::endl;
277  }
278 
279  // and rebuild Bt
280  t_get_Bt_start = TimingHelpers::timer();
281  bt_pt = new CRDoubleMatrix;
282  this->get_block(0,1,*bt_pt);
283  t_get_Bt_finish = TimingHelpers::timer();
284  if(Doc_time)
285  {
286  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
287  oomph_info << "Time to get Bt [sec]: "
288  << t_get_Bt_time << std::endl;
289  }
290  }
291 
292 
293  /////////////////////////////////////////////////////////////////////////////
294 
295 
296  // otherwise we are not forming BFBt
297  else
298  {
299 
300  // get the inverse mass matrix (Q)
301  CRDoubleMatrix* ivmm_pt = 0;
303  {
304  double ivmm_assembly_start_t = TimingHelpers::timer();
305  ivmm_pt = assemble_mass_matrix_diagonal();
306  double ivmm_assembly_finish_t = TimingHelpers::timer();
307  if (Doc_time)
308  {
309  double
310  ivmm_assembly_time = ivmm_assembly_finish_t - ivmm_assembly_start_t;
311  oomph_info << "Time to assemble Q (inverse diagonal "
312  << "mass matrix) [sec]: "
313  << ivmm_assembly_time << "\n";
314  }
315  }
316 
317  // Get Bt
318  double t_get_Bt_start = TimingHelpers::timer();
319  this->get_block(0,1,*bt_pt);
320  double t_get_Bt_finish = TimingHelpers::timer();
321  if(Doc_time)
322  {
323  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
324  oomph_info << "Time to get Bt [sec]: "
325  << t_get_Bt_time << std::endl;
326  }
327 
328  // next QBt
330  {
331  double t_QBt_matrix_start = TimingHelpers::timer();
332  CRDoubleMatrix* qbt_pt = new CRDoubleMatrix;
333  ivmm_pt->multiply(*bt_pt, *qbt_pt);
334  delete bt_pt; bt_pt = 0;
335  bt_pt = qbt_pt;
336  double t_QBt_matrix_finish = TimingHelpers::timer();
337  if(Doc_time)
338  {
339  double t_QBt_time = t_QBt_matrix_finish - t_QBt_matrix_start;
340  oomph_info << "Time to generate QBt [sec]: "
341  << t_QBt_time << std::endl;
342  }
343  delete ivmm_pt;
344  }
345 
346  // form P
347  double t_p_matrix_start = TimingHelpers::timer();
348  b_pt->multiply(*bt_pt, *p_matrix_pt);
349  double t_p_matrix_finish = TimingHelpers::timer();
350  if(Doc_time)
351  {
352  double t_p_time = t_p_matrix_finish - t_p_matrix_start;
353  oomph_info << "Time to generate P [sec]: "
354  << t_p_time << std::endl;
355  }
356  delete b_pt; b_pt = 0;
357 
358  // build the matvec operator for QBt
359  double t_QBt_MV_start = TimingHelpers::timer();
361  //QBt_mat_vec_pt->setup(bt_pt);
363  double t_QBt_MV_finish = TimingHelpers::timer();
364  if(Doc_time)
365  {
366  double t_p_time = t_QBt_MV_finish - t_QBt_MV_start;
367  oomph_info << "Time to build QBt matrix vector operator [sec]: "
368  << t_p_time << std::endl;
369  }
370  delete bt_pt; bt_pt = 0;
371 
372  // get F
373  double t_get_F_start = TimingHelpers::timer();
374  this->get_block(0,0,*f_pt);
375  double t_get_F_finish = TimingHelpers::timer();
376  if(Doc_time)
377  {
378  double t_get_F_time = t_get_F_finish - t_get_F_start;
379  oomph_info << "Time to get F [sec]: "
380  << t_get_F_time << std::endl;
381  }
382 
383  // form the matrix vector product helper
384  double t_F_MV_start = TimingHelpers::timer();
386  //F_mat_vec_pt->setup(f_pt);
388  double t_F_MV_finish = TimingHelpers::timer();
389  if(Doc_time)
390  {
391  double t_F_MV_time = t_F_MV_finish - t_F_MV_start;
392  oomph_info << "Time to build F Matrix Vector Operator [sec]: "
393  << t_F_MV_time << std::endl;
394  }
395 
396  // if F is a block preconditioner then we can delete the F matrix
397  if (F_preconditioner_is_block_preconditioner)
398  {
399  delete f_pt; f_pt = 0;
400  }
401 
402  // and rebuild Bt
403  t_get_Bt_start = TimingHelpers::timer();
404  bt_pt = new CRDoubleMatrix;
405  this->get_block(0,1,*bt_pt);
406  t_get_Bt_finish = TimingHelpers::timer();
407  if(Doc_time)
408  {
409  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
410  oomph_info << "Time to get Bt [sec]: "
411  << t_get_Bt_time << std::endl;
412  }
413  }
414 
415 
416  /////////////////////////////////////////////////////////////////////////////
417 
418 
419  // form the matrix vector operator for Bt
420  double t_Bt_MV_start = TimingHelpers::timer();
422  //Bt_mat_vec_pt->setup(bt_pt);
424  double t_Bt_MV_finish = TimingHelpers::timer();
425  if(Doc_time)
426  {
427  double t_Bt_MV_time = t_Bt_MV_finish - t_Bt_MV_start;
428  oomph_info << "Time to build Bt Matrix Vector Operator [sec]: "
429  << t_Bt_MV_time << std::endl;
430  }
431  delete bt_pt; bt_pt = 0;
432 
433  // if the P preconditioner has not been setup
434  if (P_preconditioner_pt == 0)
435  {
438  }
439 
440  // std::stringstream junk;
441  // junk << "p_matrix" << comm_pt()->my_rank()
442  // << ".dat";
443  // p_matrix_pt->sparse_indexed_output_with_offset(junk.str());
444  // oomph_info << "Done output of " << junk.str() << std::endl;
445 
446  // Setup the preconditioner for the Pressure matrix
447  double t_p_prec_start = TimingHelpers::timer();
448  P_preconditioner_pt->setup(p_matrix_pt);
449  delete p_matrix_pt; p_matrix_pt = 0;
450  p_matrix_pt=0;
451  double t_p_prec_finish = TimingHelpers::timer();
452  if(Doc_time)
453  {
454  double t_p_prec_time = t_p_prec_finish - t_p_prec_start;
455  oomph_info << "P sub-preconditioner setup time [sec]: "
456  << t_p_prec_time << "\n";
457  }
458 
459  // Set up solver for solution of system with momentum matrix
460  // ----------------------------------------------------------
461 
462  // if the F preconditioner has not been setup
463  if (F_preconditioner_pt == 0)
464  {
467  }
468 
469  // if F is a block preconditioner
470  double t_f_prec_start = TimingHelpers::timer();
471  if (F_preconditioner_is_block_preconditioner)
472  {
473  unsigned ndof_types = this->ndof_types();
474  ndof_types--;
475  Vector<unsigned> dof_map(ndof_types);
476  for (unsigned i = 0; i < ndof_types; i++)
477  {
478  dof_map[i] = i;
479  }
480  F_block_preconditioner_pt->
482  F_block_preconditioner_pt->setup(matrix_pt());
483  }
484  // otherwise F is not a block preconditioner
485  else
486  {
487  F_preconditioner_pt->setup(f_pt);
488  delete f_pt; f_pt = 0;
489  }
490  double t_f_prec_finish = TimingHelpers::timer();
491  if(Doc_time)
492  {
493  double t_f_prec_time = t_f_prec_finish - t_f_prec_start;
494  oomph_info << "F sub-preconditioner setup time [sec]: "
495  << t_f_prec_time << "\n";
496  }
497 
498  // Remember that the preconditioner has been setup so
499  // the stored information can be wiped when we
500  // come here next...
502  }
503 
504 
505 
506 //=======================================================================
507  /// Apply preconditioner to r.
508 //=======================================================================
511  {
512 
513 #ifdef PARANOID
515  {
516  std::ostringstream error_message;
517  error_message << "setup must be called before using preconditioner_solve";
518  throw OomphLibError(
519  error_message.str(),
520  OOMPH_CURRENT_FUNCTION,
521  OOMPH_EXCEPTION_LOCATION);
522  }
523  if (z.built())
524  {
525  if (z.nrow() != r.nrow())
526  {
527  std::ostringstream error_message;
528  error_message << "The vectors z and r must have the same number of "
529  << "of global rows";
530  throw OomphLibError(
531  error_message.str(),
532  OOMPH_CURRENT_FUNCTION,
533  OOMPH_EXCEPTION_LOCATION);
534  }
535  }
536 #endif
537 
538  double t_start_overall = TimingHelpers::timer();
539  double t_start = TimingHelpers::timer();
540  double t_end=0;
541 
542  // if z is not setup then give it the same distribution
543  if (!z.built())
544  {
545  z.build(r.distribution_pt(),0.0);
546  }
547 
548  // Step 1 - apply approximate Schur inverse to pressure unknowns (block 1)
549  // -----------------------------------------------------------------------
550 
551  // Working vectors
552  DoubleVector temp_vec;
553  DoubleVector another_temp_vec;
554 
555  // Copy pressure values from residual vector to temp_vec:
556  // Loop over all entries in the global vector (this one
557  // includes displacement/position and pressure dofs in some random fashion)
558  this->get_block_vector(1,r,temp_vec);
559 
560 
561  if(Doc_time)
562  {
563  t_end=TimingHelpers::timer();
564  oomph_info << "LSC prec solve: Time for get block vector: "
565  << t_end-t_start << std::endl;
566  t_start=TimingHelpers::timer();
567  }
568 
569  // NOTE: The vector temp_vec now contains the vector r_p.
570 
571  // Solve first pressure Poisson system
572 #ifdef PARANOID
573  // check a solver has been set
574  if (P_preconditioner_pt==0)
575  {
576  std::ostringstream error_message;
577  error_message << "P_preconditioner_pt has not been set.";
578  throw OomphLibError(
579  error_message.str(),
580  OOMPH_CURRENT_FUNCTION,
581  OOMPH_EXCEPTION_LOCATION);
582  }
583 #endif
584 
585  // use some Preconditioner's preconditioner_solve function
586  P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
587 
588 
589  if(Doc_time)
590  {
591  t_end=TimingHelpers::timer();
592  oomph_info << "LSC prec solve: First P solve [nrow="
593  << P_preconditioner_pt->nrow() << "]: "
594  << t_end-t_start << std::endl;
595  t_start=TimingHelpers::timer();
596  }
597 
598 
599  // NOTE: The vector another_temp_vec now contains the vector P^{-1} r_p
600 
601  // Multiply another_temp_vec by matrix E and stick the result into temp_vec
602  temp_vec.clear();
603  if (Form_BFBt_product)
604  {
605  E_mat_vec_pt->multiply(another_temp_vec, temp_vec);
606  }
607  else
608  {
609  QBt_mat_vec_pt->multiply(another_temp_vec, temp_vec);
610  another_temp_vec.clear();
611  F_mat_vec_pt->multiply(temp_vec,another_temp_vec);
612  temp_vec.clear();
613  QBt_mat_vec_pt->multiply_transpose(another_temp_vec, temp_vec);
614  }
615 
616 
617  if(Doc_time)
618  {
619  t_end=TimingHelpers::timer();
620  oomph_info << "LSC prec solve: E matrix vector product: "
621  << t_end-t_start << std::endl;
622  t_start=TimingHelpers::timer();
623  }
624 
625  // NOTE: The vector temp_vec now contains E P^{-1} r_p
626 
627  // Solve second pressure Poisson system using preconditioner_solve
628  another_temp_vec.clear();
629  P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
630 
631 
632  if(Doc_time)
633  {
634  t_end=TimingHelpers::timer();
635  oomph_info << "LSC prec solve: Second P solve [nrow="
636  << P_preconditioner_pt->nrow() << "]: "
637  << t_end-t_start << std::endl;
638  t_start=TimingHelpers::timer();
639  }
640 
641 
642  // NOTE: The vector another_temp_vec now contains z_p = P^{-1} E P^{-1} r_p
643  // as required (apart from the sign which we'll fix in the
644  // next step.
645 
646  // Now copy another_temp_vec (i.e. z_p) back into the global vector z.
647  // Loop over all entries in the global results vector z:
648  temp_vec.build(another_temp_vec.distribution_pt(),0.0);
649  temp_vec -= another_temp_vec;
650  return_block_vector(1,temp_vec,z);
651 
652  // Step 2 - apply preconditioner to displacement/positon unknowns (block 0)
653  // ------------------------------------------------------------------------
654 
655  // Recall that another_temp_vec (computed above) contains the
656  // negative of the solution of the Schur complement systen, -z_p.
657  // Multiply by G (stored in Block_matrix_pt(0,1) and store
658  // result in temp_vec (vector resizes itself).
659  temp_vec.clear();
660  Bt_mat_vec_pt->multiply(another_temp_vec, temp_vec);
661 
662 
663  if(Doc_time)
664  {
665  t_end=TimingHelpers::timer();
666  oomph_info << "LSC prec solve: G matrix vector product: "
667  << t_end-t_start << std::endl;
668  t_start=TimingHelpers::timer();
669  }
670 
671  // NOTE: temp_vec now contains -G z_p
672 
673  // The vector another_temp_vec is no longer needed -- re-use it to store
674  // displacement/position quantities:
675  another_temp_vec.clear();
676 
677  // Loop over all enries in the global vector and find the
678  // entries associated with the displacement/position:
679  get_block_vector(0,r,another_temp_vec);
680  another_temp_vec += temp_vec;
681 
682  // NOTE: The vector another_temp_vec now contains r_u - G z_p
683 
684  // Solve momentum system
685 #ifdef PARANOID
686  // check a solver has been set
687  if (F_preconditioner_pt==0)
688  {
689  std::ostringstream error_message;
690  error_message << "F_preconditioner_pt has not been set.";
691  throw OomphLibError(
692  error_message.str(),
693  OOMPH_CURRENT_FUNCTION,
694  OOMPH_EXCEPTION_LOCATION);
695  }
696 #endif
697 
698  // use some Preconditioner's preconditioner solve
699  // and return
701  {
702  return_block_vector(0,another_temp_vec,z);
704  }
705  else
706  {
707  F_preconditioner_pt->preconditioner_solve(another_temp_vec, temp_vec);
708  return_block_vector(0,temp_vec,z);
709  }
710 
711  if(Doc_time)
712  {
713  t_end=TimingHelpers::timer();
714  oomph_info << "LSC prec solve: F solve [nrow="
715  << P_preconditioner_pt->nrow() << "]: "
716  << t_end-t_start << std::endl;
717  oomph_info << "LSC prec solve: Overall "
718  << t_end-t_start_overall << std::endl;
719  }
720 
721  }
722 
723 
724 //========================================================================
725 /// Helper function to assemble the diagonal of the
726 /// mass matrix from the elemental contributions defined in
727 /// SolidElementWithDiagonalMassMatrix::get_mass_matrix_diagonal(...).
728 //========================================================================
731  {
732  // determine the rows required by this processor
733  unsigned first_row = this->block_distribution_pt(0)->first_row();
734  unsigned nrow_local = this->block_distribution_pt(0)->nrow_local();
735  unsigned nrow = this->block_distribution_pt(0)->nrow();
736 
737  // create storage for the diagonals
738  double* m_values = new double[nrow_local];
739  for (unsigned i = 0; i < nrow_local; i++)
740  {
741  m_values[i] = 0;
742  }
743 
744  // if the problem is distributed
745  bool distributed = false;
746 #ifdef OOMPH_HAS_MPI
747  if (any_mesh_distributed() ||
749  {
750  distributed = true;
751  }
752 #endif
753 
754  // next we get the diagonal mass matrix data
755  if (distributed)
756  {
757 #ifdef OOMPH_HAS_MPI
758  // the number of processors
759  unsigned nproc = this->comm_pt()->nproc();
760 
761  // and my rank
762  unsigned my_rank = this->comm_pt()->my_rank();
763 
764  // determine the rows for which we have lookup rows
765  // if the problem is NOT distributed then we only classify global equation
766  // on this processor to avoid duplication (as every processor holds
767  // every element)
768  unsigned first_lookup_row = 0;
769  unsigned last_lookup_row = 0;
770  first_lookup_row = this->master_distribution_pt()->first_row();
771  last_lookup_row = first_lookup_row +
772  this->master_distribution_pt()->nrow_local() - 1;
773 
774  // find number of local elements
775  unsigned n_el = Solid_mesh_pt->nelement();
776 
777  // the diagonal mass matrix contributions that have been
778  // classified and should be sent to another processor
779  Vector<double>* classified_contributions_send
780  = new Vector<double>[nproc];
781 
782  // the corresponding block indices
783  Vector<unsigned>* classified_indices_send
784  = new Vector<unsigned>[nproc];
785 
786  // the maitrix contributions that cannot be classified by this processor
787  // and therefore must be sent to another for classification
788  Vector<double>* unclassified_contributions_send
789  = new Vector<double>[nproc];
790 
791  // the corresponding global indices that require classification
792  Vector<unsigned>* unclassified_indices_send
793  = new Vector<unsigned>[nproc];
794 
795  // get the master distribution pt
797  this->master_distribution_pt();
798 
799  // get the displacement/position distribution pt
800  const LinearAlgebraDistribution* displ_dist_pt
801  = this->block_distribution_pt(0);
802 
803  // get the contribution for each element
804  for (unsigned e = 0; e < n_el; e++)
805  {
806 
807  // check that the element is not halo d
808  if (!Solid_mesh_pt->element_pt(e)->is_halo())
809  {
810 
811  // find number of degrees of freedom in the element
812  // (this is slightly too big because it includes the
813  // pressure dofs but this doesn't matter)
814  unsigned el_dof = Solid_mesh_pt->element_pt(e)->ndof();
815 
816  // allocate local storage for the element's contribution to the
817  // mass matrix diagonal
818  Vector<double> el_vmm_diagonal(el_dof);
819 
821  cast_el_pt=dynamic_cast<
823  ( Solid_mesh_pt->element_pt(e) );
824 
825 
826  if (cast_el_pt==0)
827  {
828 #ifdef PARANOID
829  std::ostringstream error_message;
830  error_message
831  << "Failed cast to "
832  << "SolidElementWithDiagonalMassMatrix*\n"
833  << "Element is of type: "
834  << typeid(*(Solid_mesh_pt->element_pt(e))).name() << "\n"
835  << typeid(Solid_mesh_pt->element_pt(e)).name() << std::endl;
837  error_message.str(),
838  "PressureBasedSolidLSCPreconditioner::assemble_mass_matrix_diagonal()",
839  OOMPH_EXCEPTION_LOCATION);
840 #endif
841  }
842  else
843  {
844  cast_el_pt->get_mass_matrix_diagonal(el_vmm_diagonal);
845  }
846 
847  // get the contribution for each dof
848  for (unsigned i = 0; i < el_dof; i++)
849  {
850 
851  //Get the equation number
852  unsigned eqn_number = Solid_mesh_pt
853  ->element_pt(e)->eqn_number(i);
854 
855  // if I have lookup information on this processor
856  if (eqn_number >= first_lookup_row &&
857  eqn_number <= last_lookup_row)
858  {
859 
860  // bypass non displacement/position DOFs
861  if ( this->block_number(eqn_number)==0 )
862  {
863 
864  // get the index in the block
865  unsigned index = this->index_in_block(eqn_number);
866 
867  // determine which processor requires the block index
868  for (unsigned p = 0; p < nproc; p++)
869  {
870  if (index >= displ_dist_pt->first_row(p) &&
871  (index <
872  (displ_dist_pt->first_row(p)
873  +displ_dist_pt->nrow_local(p))))
874  {
875 
876  // if it is required by this processor then add the
877  // contribution
878  if (p == my_rank)
879  {
880  m_values[index-first_row] += el_vmm_diagonal[i];
881  }
882  // other wise store it for communication
883  else
884  {
885  classified_contributions_send[p]
886  .push_back(el_vmm_diagonal[i]);
887  classified_indices_send[p].push_back(index);
888  }
889  }
890  }
891  }
892  }
893 
894  // if we do not have the lookup information on this processor
895  // then we send the mass matrix contribution to a processor
896  // which we know does have the lookup information
897  // the assumption: the processor for which the master block
898  // preconditioner distribution will definitely hold the lookup
899  // data for eqn_number (although others may)
900  else if (any_mesh_distributed())
901  {
902 
903  // determine which processor requires the block index
904  unsigned p = 0;
905  while (!(eqn_number >= master_distribution_pt->first_row(p) &&
906  (eqn_number < (master_distribution_pt->first_row(p)
907  +master_distribution_pt->nrow_local(p)))))
908  {
909  p++;
910  }
911 
912  // store the data
913  unclassified_contributions_send[p]
914  .push_back(el_vmm_diagonal[i]);
915  unclassified_indices_send[p].push_back(eqn_number);
916  }
917  }
918  }
919  }
920 
921  // next the unclassified contributions are communicated to
922  // processors that can classify them
923 
924  // first determine how many unclassified rows are to be sent to
925  // each processor
926  unsigned* n_unclassified_send = new unsigned[nproc];
927  for (unsigned p = 0; p < nproc; p++)
928  {
929  if (p == my_rank)
930  {
931  n_unclassified_send[p] = 0;
932  }
933  else
934  {
935  n_unclassified_send[p]
936  = unclassified_contributions_send[p].size();
937  }
938  }
939 
940  // then all-to-all number of unclassified to be sent / recv
941  unsigned* n_unclassified_recv = new unsigned[nproc];
942  MPI_Alltoall(n_unclassified_send,1,MPI_UNSIGNED,
943  n_unclassified_recv,1,MPI_UNSIGNED,
944  this->comm_pt()->mpi_comm());
945 
946  // the base displacement for the sends
947  MPI_Aint base_displacement;
948  MPI_Address(m_values,&base_displacement);
949 
950  // allocate storage for the data to be recieved
951  // and post the sends and recvs
952  Vector<double*> unclassified_contributions_recv(nproc);
953  Vector<unsigned*> unclassified_indices_recv(nproc);
954  Vector<MPI_Request> unclassified_recv_requests;
955  Vector<MPI_Request> unclassified_send_requests;
956  Vector<unsigned> unclassified_recv_proc;
957  for (unsigned p = 0; p < nproc; p++)
958  {
959  if (p != my_rank)
960  {
961  // recv
962  if (n_unclassified_recv[p] > 0)
963  {
964  unclassified_contributions_recv[p]
965  = new double[n_unclassified_recv[p]];
966  unclassified_indices_recv[p] = new
967  unsigned[n_unclassified_recv[p]];
968 
969  // data for the struct data type
970  MPI_Datatype recv_types[2];
971  MPI_Aint recv_displacements[2];
972  int recv_sz[2];
973 
974  // contributions
975  MPI_Type_contiguous(n_unclassified_recv[p],MPI_DOUBLE,
976  &recv_types[0]);
977  MPI_Type_commit(&recv_types[0]);
978  MPI_Address(unclassified_contributions_recv[p],
979  &recv_displacements[0]);
980  recv_displacements[0] -= base_displacement;
981  recv_sz[0] = 1;
982 
983  // indices
984  MPI_Type_contiguous(n_unclassified_recv[p],MPI_UNSIGNED,
985  &recv_types[1]);
986  MPI_Type_commit(&recv_types[1]);
987  MPI_Address(unclassified_indices_recv[p],
988  &recv_displacements[1]);
989  recv_displacements[1] -= base_displacement;
990  recv_sz[1] = 1;
991 
992  // build the final recv type
993  MPI_Datatype final_recv_type;
994  MPI_Type_struct(2,recv_sz,recv_displacements,recv_types,
995  &final_recv_type);
996  MPI_Type_commit(&final_recv_type);
997 
998  // and recv
999  MPI_Request req;
1000  MPI_Irecv(m_values,1,final_recv_type,p,0,
1001  comm_pt()->mpi_comm(),&req);
1002  unclassified_recv_requests.push_back(req);
1003  unclassified_recv_proc.push_back(p);
1004  MPI_Type_free(&recv_types[0]);
1005  MPI_Type_free(&recv_types[1]);
1006  MPI_Type_free(&final_recv_type);
1007  }
1008 
1009  // send
1010  if (n_unclassified_send[p] > 0)
1011  {
1012  // data for the struct data type
1013  MPI_Datatype send_types[2];
1014  MPI_Aint send_displacements[2];
1015  int send_sz[2];
1016 
1017  // contributions
1018  MPI_Type_contiguous(n_unclassified_send[p],MPI_DOUBLE,
1019  &send_types[0]);
1020  MPI_Type_commit(&send_types[0]);
1021  MPI_Address(&unclassified_contributions_send[p][0],
1022  &send_displacements[0]);
1023  send_displacements[0] -= base_displacement;
1024  send_sz[0] = 1;
1025 
1026  // indices
1027  MPI_Type_contiguous(n_unclassified_send[p],MPI_UNSIGNED,
1028  &send_types[1]);
1029  MPI_Type_commit(&send_types[1]);
1030  MPI_Address(&unclassified_indices_send[p][0],
1031  &send_displacements[1]);
1032  send_displacements[1] -= base_displacement;
1033  send_sz[1] = 1;
1034 
1035  // build the final send type
1036  MPI_Datatype final_send_type;
1037  MPI_Type_struct(2,send_sz,send_displacements,send_types,
1038  &final_send_type);
1039  MPI_Type_commit(&final_send_type);
1040 
1041  // and send
1042  MPI_Request req;
1043  MPI_Isend(m_values,1,final_send_type,p,0,
1044  comm_pt()->mpi_comm(),&req);
1045  unclassified_send_requests.push_back(req);
1046  MPI_Type_free(&send_types[0]);
1047  MPI_Type_free(&send_types[1]);
1048  MPI_Type_free(&final_send_type);
1049  }
1050  }
1051  }
1052 
1053  // next classify the data as it is received
1054  unsigned n_unclassified_recv_req = unclassified_recv_requests.size();
1055  while (n_unclassified_recv_req > 0)
1056  {
1057  // get the processor number and remove the completed request
1058  // for the vector of requests
1059  int req_num;
1060  MPI_Waitany(n_unclassified_recv_req,&unclassified_recv_requests[0],
1061  &req_num,MPI_STATUS_IGNORE);
1062  unsigned p = unclassified_recv_proc[req_num];
1063  unclassified_recv_requests.erase(unclassified_recv_requests.begin()
1064  +req_num);
1065  unclassified_recv_proc.erase(unclassified_recv_proc.begin()+req_num);
1066  n_unclassified_recv_req--;
1067 
1068  // next classify the dofs
1069  // and store them for sending to other processors if required
1070  unsigned n_recv = n_unclassified_recv[p];
1071  for (unsigned i = 0; i < n_recv; i++)
1072  {
1073  unsigned eqn_number = unclassified_indices_recv[p][i];
1074  // bypass non displacement/position DOFs
1075  if ( this->block_number(eqn_number)==0 )
1076  {
1077 
1078  // get the index in the block
1079  unsigned index = this->index_in_block(eqn_number);
1080 
1081  // determine which processor requires the block index
1082  for (unsigned pp = 0; pp < nproc; pp++)
1083  {
1084 
1085 
1086  if (index >= displ_dist_pt->first_row(pp) &&
1087  (index <
1088  (displ_dist_pt->first_row(pp)
1089  +displ_dist_pt->nrow_local(pp))))
1090  {
1091 
1092  // if it is required by this processor then add the
1093  // contribution
1094  if (pp == my_rank)
1095  {
1096  m_values[index-first_row]
1097  += unclassified_contributions_recv[p][i];
1098  }
1099  // other wise store it for communication
1100  else
1101  {
1102  double v = unclassified_contributions_recv[p][i];
1103  classified_contributions_send[pp].push_back(v);
1104  classified_indices_send[pp].push_back(index);
1105  }
1106  }
1107  }
1108  }
1109  }
1110 
1111  // clean up
1112  delete[] unclassified_contributions_recv[p];
1113  delete[] unclassified_indices_recv[p];
1114  }
1115  delete[] n_unclassified_recv;
1116 
1117  // now all indices have been classified
1118 
1119  // next the classified contributions are communicated to
1120  // processors that require them
1121 
1122  // first determine how many classified rows are to be sent to
1123  // each processor
1124  unsigned* n_classified_send = new unsigned[nproc];
1125  for (unsigned p = 0; p < nproc; p++)
1126  {
1127  if (p == my_rank)
1128  {
1129  n_classified_send[p] = 0;
1130  }
1131  else
1132  {
1133  n_classified_send[p]
1134  = classified_contributions_send[p].size();
1135  }
1136  }
1137 
1138  // then all-to-all com number of classified to be sent / recv
1139  unsigned* n_classified_recv = new unsigned[nproc];
1140  MPI_Alltoall(n_classified_send,1,MPI_UNSIGNED,
1141  n_classified_recv,1,MPI_UNSIGNED,
1142  this->comm_pt()->mpi_comm());
1143 
1144  // allocate storage for the data to be recieved
1145  // and post the sends and recvs
1146  Vector<double*> classified_contributions_recv(nproc);
1147  Vector<unsigned*> classified_indices_recv(nproc);
1148  Vector<MPI_Request> classified_recv_requests;
1149  Vector<MPI_Request> classified_send_requests;
1150  Vector<unsigned> classified_recv_proc;
1151  for (unsigned p = 0; p < nproc; p++)
1152  {
1153  if (p != my_rank)
1154  {
1155  // recv
1156  if (n_classified_recv[p] > 0)
1157  {
1158  classified_contributions_recv[p]
1159  = new double[n_classified_recv[p]];
1160  classified_indices_recv[p] = new unsigned[n_classified_recv[p]];
1161 
1162  // data for the struct data type
1163  MPI_Datatype recv_types[2];
1164  MPI_Aint recv_displacements[2];
1165  int recv_sz[2];
1166 
1167  // contributions
1168  MPI_Type_contiguous(n_classified_recv[p],MPI_DOUBLE,
1169  &recv_types[0]);
1170  MPI_Type_commit(&recv_types[0]);
1171  MPI_Address(classified_contributions_recv[p],
1172  &recv_displacements[0]);
1173  recv_displacements[0] -= base_displacement;
1174  recv_sz[0] = 1;
1175 
1176  // indices
1177  MPI_Type_contiguous(n_classified_recv[p],MPI_UNSIGNED,
1178  &recv_types[1]);
1179  MPI_Type_commit(&recv_types[1]);
1180  MPI_Address(classified_indices_recv[p],
1181  &recv_displacements[1]);
1182  recv_displacements[1] -= base_displacement;
1183  recv_sz[1] = 1;
1184 
1185  // build the final recv type
1186  MPI_Datatype final_recv_type;
1187  MPI_Type_struct(2,recv_sz,recv_displacements,recv_types,
1188  &final_recv_type);
1189  MPI_Type_commit(&final_recv_type);
1190 
1191  // and recv
1192  MPI_Request req;
1193  MPI_Irecv(m_values,1,final_recv_type,p,0,
1194  comm_pt()->mpi_comm(),&req);
1195  classified_recv_requests.push_back(req);
1196  classified_recv_proc.push_back(p);
1197  MPI_Type_free(&recv_types[0]);
1198  MPI_Type_free(&recv_types[1]);
1199  MPI_Type_free(&final_recv_type);
1200  }
1201 
1202  // send
1203  if (n_classified_send[p] > 0)
1204  {
1205  // data for the struct data type
1206  MPI_Datatype send_types[2];
1207  MPI_Aint send_displacements[2];
1208  int send_sz[2];
1209 
1210  // contributions
1211  MPI_Type_contiguous(n_classified_send[p],MPI_DOUBLE,
1212  &send_types[0]);
1213  MPI_Type_commit(&send_types[0]);
1214  MPI_Address(&classified_contributions_send[p][0],
1215  &send_displacements[0]);
1216  send_displacements[0] -= base_displacement;
1217  send_sz[0] = 1;
1218 
1219  // indices
1220  MPI_Type_contiguous(n_classified_send[p],MPI_UNSIGNED,
1221  &send_types[1]);
1222  MPI_Type_commit(&send_types[1]);
1223  MPI_Address(&classified_indices_send[p][0],
1224  &send_displacements[1]);
1225  send_displacements[1] -= base_displacement;
1226  send_sz[1] = 1;
1227 
1228  // build the final send type
1229  MPI_Datatype final_send_type;
1230  MPI_Type_struct(2,send_sz,send_displacements,send_types,
1231  &final_send_type);
1232  MPI_Type_commit(&final_send_type);
1233 
1234  // and send
1235  MPI_Request req;
1236  MPI_Isend(m_values,1,final_send_type,p,0,
1237  comm_pt()->mpi_comm(),&req);
1238  classified_send_requests.push_back(req);
1239  MPI_Type_free(&send_types[0]);
1240  MPI_Type_free(&send_types[1]);
1241  MPI_Type_free(&final_send_type);
1242  }
1243  }
1244  }
1245 
1246  // next classify the data as it is received
1247  unsigned n_classified_recv_req = classified_recv_requests.size();
1248  while (n_classified_recv_req > 0)
1249  {
1250  // get the processor number and remove the completed request
1251  // for the vector of requests
1252  int req_num;
1253  MPI_Waitany(n_classified_recv_req,&classified_recv_requests[0],
1254  &req_num,MPI_STATUS_IGNORE);
1255  unsigned p = classified_recv_proc[req_num];
1256  classified_recv_requests.erase(classified_recv_requests.begin()
1257  +req_num);
1258  classified_recv_proc.erase(classified_recv_proc.begin()+req_num);
1259  n_classified_recv_req--;
1260 
1261  // next classify the dofs
1262  // and store them for sending to other processors if required
1263  unsigned n_recv = n_classified_recv[p];
1264  for (unsigned i = 0; i < n_recv; i++)
1265  {
1266  m_values[classified_indices_recv[p][i]-first_row]
1267  += classified_contributions_recv[p][i];
1268  }
1269 
1270  // clean up
1271  delete[] classified_contributions_recv[p];
1272  delete[] classified_indices_recv[p];
1273  }
1274 
1275  // wait for the unclassified sends to complete
1276  unsigned n_unclassified_send_req = unclassified_send_requests.size();
1277  if (n_unclassified_send_req > 0)
1278  {
1279  MPI_Waitall(n_unclassified_send_req,&unclassified_send_requests[0],
1280  MPI_STATUS_IGNORE);
1281  }
1282  delete[] unclassified_contributions_send;
1283  delete[] unclassified_indices_send;
1284  delete[] n_unclassified_send;
1285 
1286  // wait for the classified sends to complete
1287  unsigned n_classified_send_req = classified_send_requests.size();
1288  if (n_classified_send_req > 0)
1289  {
1290  MPI_Waitall(n_classified_send_req,&classified_send_requests[0],
1291  MPI_STATUS_IGNORE);
1292  }
1293  delete[] classified_indices_send;
1294  delete[] classified_contributions_send;
1295  delete[] n_classified_recv;
1296  delete[] n_classified_send;
1297 #endif
1298  }
1299 
1300  // or if the problem is not distributed
1301  else
1302  {
1303 
1304  // find number of elements
1305  unsigned n_el = Solid_mesh_pt->nelement();
1306 
1307  // get the contribution for each element
1308  for (unsigned e = 0; e < n_el; e++)
1309  {
1310 
1311  // find number of degrees of freedom in the element
1312  // (this is slightly too big because it includes the
1313  // pressure dofs but this doesn't matter)
1314  unsigned el_dof = Solid_mesh_pt->element_pt(e)->ndof();
1315 
1316  // allocate local storage for the element's contribution to the
1317  // mass matrix diagonal
1318  Vector<double> el_vmm_diagonal(el_dof);
1319 
1321  cast_el_pt=dynamic_cast<
1323  ( Solid_mesh_pt->element_pt(e) );
1324 
1325  if (cast_el_pt==0)
1326  {
1327 #ifdef PARANOID
1328 // #pragma clang diagnostic push
1329 // #pragma clang diagnostic ignored "-Wpotentially-evaluated-expression"
1330  std::ostringstream error_message;
1331  error_message
1332  << "Failed cast to "
1333  << "SolidElementWithDiagonalMassMatrix*\n"
1334  << "Element is of type: "
1335  << typeid(*(Solid_mesh_pt->element_pt(e))).name() << "\n"
1336  << typeid(Solid_mesh_pt->element_pt(e)).name() << std::endl;
1338  error_message.str(),
1339  "PressureBasedSolidLSCPreconditioner::assemble_mass_matrix_diagonal()",
1340  OOMPH_EXCEPTION_LOCATION);
1341 //#pragma clang diagnostic pop
1342 #endif
1343  }
1344  else
1345  {
1346  cast_el_pt->get_mass_matrix_diagonal(el_vmm_diagonal);
1347  }
1348 
1349  // get the contribution for each dof
1350  for (unsigned i = 0; i < el_dof; i++)
1351  {
1352  //Get the equation number
1353  unsigned eqn_number = Solid_mesh_pt
1354  ->element_pt(e)->eqn_number(i);
1355 
1356  // bypass non displacement/position DOFs
1357  if ( this->block_number(eqn_number)==0 )
1358  {
1359 
1360  // get the index in the block
1361  unsigned index = this->index_in_block(eqn_number);
1362 
1363  // if it is required on this processor
1364  if (index >= first_row &&
1365  index < first_row + nrow_local)
1366  {
1367  m_values[index-first_row] += el_vmm_diagonal[i];
1368  }
1369  }
1370  }
1371  }
1372  }
1373 
1374  // create column index and row start
1375  int* m_column_index = new int[nrow_local];
1376  int* m_row_start = new int[nrow_local+1];
1377  for (unsigned i = 0; i < nrow_local; i++)
1378  {
1379  m_values[i] = 1 / m_values[i];
1380  m_column_index[i] = first_row + i;
1381  m_row_start[i] = i;
1382  }
1383  m_row_start[nrow_local] = nrow_local;
1384 
1385  // build the matrix
1386  CRDoubleMatrix* m_pt = new CRDoubleMatrix(this->block_distribution_pt(0));
1387  m_pt->build_without_copy(nrow,nrow_local,m_values,m_column_index,
1388  m_row_start);
1389 
1390  // return the matrix;
1391  return m_pt;
1392  }
1393 
1394 //=========================================================================
1395 /// Helper function to delete preconditioner data.
1396 //=========================================================================
1398  {
1400  {
1401  // delete matvecs
1402  delete Bt_mat_vec_pt; Bt_mat_vec_pt = 0;
1403 
1404  if (!Form_BFBt_product)
1405  {
1406  delete F_mat_vec_pt; F_mat_vec_pt = 0;
1407  delete QBt_mat_vec_pt; QBt_mat_vec_pt = 0;
1408  }
1409  else
1410  {
1411  delete E_mat_vec_pt; E_mat_vec_pt = 0;
1412  }
1413 
1414  // delete stuff from displacement/position solve
1416  {
1418  }
1419 
1420  // delete stuff from Schur complement approx
1422  {
1424  }
1425  }
1426  }
1427 }
bool Form_BFBt_product
indicates whether BFBt should be formed or the component matrices should be retained. If true then: in setup(...) : BFBt is computed. in preconditioner_solve(...) : a single matrix vector product with BFBt is performed. if false then: in setup(...) : the matrices B, F are assembled and stored. in preconditioner_solve(...) : a sequence of matrix vector products with B, F, and Bt is performed. (Note: in this discussion no scaling was considered but B and Bt are replaced with BQ and QBt with scaling)
bool built() const
unsigned nrow_local() const
access function for the num of local rows on this processor.
virtual void get_mass_matrix_diagonal(Vector< double > &mass_diag)=0
Get the diagonal of whatever represents the mass matrix in the specific preconditionable element...
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero...
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
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 clean_up_memory()
Helper function to delete preconditioner data.
void clear()
wipes the DoubleVector
bool any_mesh_distributed() const
Check if any of the meshes are distributed. This is equivalent to problem.distributed() and is used a...
MatrixVectorProduct * Bt_mat_vec_pt
MatrixVectorProduct operator for Bt;.
unsigned first_row() const
access function for the first row on this processor
cstr elem_len * i
Definition: cfortran.h:607
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. If no MPI then Nrow is returned...
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
bool is_subsidiary_block_preconditioner() const
Return true if this preconditioner is a 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 ...
bool distributed() const
distribution is serial or distributed
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...
OomphInfo oomph_info
bool P_matrix_using_scaling
Control flag is true if mass matrix diagonal scaling is used in the Schur complement approximation...
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.
e
Definition: cfortran.h:575
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to Vector r.
MatrixVectorProduct * F_mat_vec_pt
MatrixVectorProduct operator for F if BFBt is not to be formed.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
MatrixVectorProduct * QBt_mat_vec_pt
MatrixVectorProduct operator for QBt if BFBt is not to be formed.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool Using_default_f_preconditioner
flag indicating whether the default F preconditioner is used
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
void multiply(const DoubleVector &x, DoubleVector &y) const
Apply the operator to the vector x and return the result in the vector y.
MatrixVectorProduct * E_mat_vec_pt
MatrixVectorProduct operator for E (BFBt) if BFBt is to be formed.
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
const LinearAlgebraDistribution * block_distribution_pt(const unsigned &b) const
Access function to the block distributions (const version).
unsigned nrow() const
access function to the number of global rows.
bool F_preconditioner_is_block_preconditioner
Boolean indicating whether the momentum system preconditioner is a block preconditioner.
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...
double timer()
returns the time in seconds after some point in past
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.
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
unsigned ndof_types() const
Return the total number of DOF types.
An interface to allow SuperLU to be used as an (exact) Preconditioner.
void setup()
Broken assignment operator.
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< CRDoubleMatrix > *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 multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1805
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
void multiply_transpose(const DoubleVector &x, DoubleVector &y) const
Apply the transpose of the operator to the vector x and return the result in the vector y...
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
bool Using_default_p_preconditioner
flag indicating whether the default P preconditioner is used
CRDoubleMatrix * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*...
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...
Mesh * Solid_mesh_pt
the pointer to the mesh of block preconditionable solid elements.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
int block_number(const unsigned &i_dof) const
Return the block number corresponding to a global index i_dof.
bool is_halo() const
Is this element a halo?
Definition: elements.h:1141
const LinearAlgebraDistribution * master_distribution_pt() const
Access function to the distribution of the master preconditioner. If this preconditioner does not hav...
int index_in_block(const unsigned &i_dof) const
Given a global dof number, returns the index in the block it belongs to. This is the overall index...
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
bool Doc_time
Set Doc_time to true for outputting results of timings.
void build_without_copy(const unsigned &ncol, const unsigned &nnz, double *value, int *column_index, int *row_start)
keeps the existing distribution and just matrix that is stored without copying the matrix data ...
Definition: matrices.cc:1731