navier_stokes_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: 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//====================================================================
31 
32 namespace oomph
33 {
34 
35 
36 ///////////////////////////////////////////////////////////////////////
37 ///////////////////////////////////////////////////////////////////////
38 ///////////////////////////////////////////////////////////////////////
39 
40 
41 //======start_of_namespace============================================
42 /// Namespace for exact solution for pressure advection diffusion
43 /// problem
44 //====================================================================
45  namespace PressureAdvectionDiffusionValidation
46  {
47 
48 
49  /// Flag for solution
50  unsigned Flag=0;
51 
52  /// Peclet number -- overwrite with actual Reynolds number
53  double Peclet=0.0;
54 
55  /// Wind
57  {
58  if (Flag==0)
59  {
60  wind[0]=sin(6.0*x[1]);
61  wind[1]=cos(6.0*x[0]);
62  }
63  else
64  {
65  wind[0]=1.0;
66  wind[1]=0.0;
67  }
68  }
69 
70  /// Exact solution as a Vector
72  {
73  u.resize(3);
74  wind_function(x,u);
75  if (Flag==0)
76  {
77  u[2]=x[0]*x[0]*pow(1.0-x[0],2.0)*x[1]*x[1]*pow(1.0-x[1],2.0);
78  }
79  else
80  {
81  u[2]=0.1E1-Peclet*x[0]*(0.1E1-0.5*x[0]);
82  }
83  }
84 
85  /// Exact solution as a scalar
86  void get_exact_u(const Vector<double>& x, double& u)
87  {
88  if (Flag==0)
89  {
90  u=x[0]*x[0]*pow(1.0-x[0],2.0)*x[1]*x[1]*pow(1.0-x[1],2.0);
91  }
92  else
93  {
94  u=0.1E1-Peclet*x[0]*(0.1E1-0.5*x[0]);
95  }
96  }
97 
98  /// Source function required to make the solution above an exact solution
99  double source_function(const Vector<double>& x_vect)
100  {
101 
102  double x[2];
103  x[0]=x_vect[0];
104  x[1]=x_vect[1];
105 
106 
107  double source=0.0;
108 
109  if (Flag==0)
110  {
111  source=
112  Peclet*(sin(0.6E1*x[1])*(2.0*x[0]*pow(1.0-x[0],2.0)*x[1]*x[1]*pow(
113  1.0-x[1],2.0)-2.0*x[0]*x[0]*(1.0-x[0])*x[1]*x[1]*pow(1.0-x[1],2.0))+cos(0.6E1*x
114  [0])*(2.0*x[0]*x[0]*pow(1.0-x[0],2.0)*x[1]*pow(1.0-x[1],2.0)-2.0*x[0]*x[0]*pow(
115  1.0-x[0],2.0)*x[1]*x[1]*(1.0-x[1])))-2.0*pow(1.0-x[0],2.0)*x[1]*x[1]*pow(1.0-x
116  [1],2.0)+8.0*x[0]*(1.0-x[0])*x[1]*x[1]*pow(1.0-x[1],2.0)-2.0*x[0]*x[0]*x[1]*x
117  [1]*pow(1.0-x[1],2.0)-2.0*x[0]*x[0]*pow(1.0-x[0],2.0)*pow(1.0-x[1],2.0)+8.0*x
118  [0]*x[0]*pow(1.0-x[0],2.0)*x[1]*(1.0-x[1])-2.0*x[0]*x[0]*pow(1.0-x[0],2.0)*x[1]
119  *x[1];
120  }
121  else
122  {
123  source=Peclet*(-0.1E1*Peclet*(0.1E1-0.5*x[0])+0.5*Peclet*x[0])-0.1E1*Peclet
124 ;
125 
126  }
127 
128  return source;
129  }
130 
131 
132  } // end of namespace
133 
134 
135 ////////////////////////////////////////////////////////////////
136 ////////////////////////////////////////////////////////////////
137 ////////////////////////////////////////////////////////////////
138 
139 
140 //===========================================================================
141 /// Setup the least-squares commutator Navier Stokes preconditioner. This
142 /// extracts blocks corresponding to the velocity and pressure unknowns,
143 /// creates the matrices actually needed in the application of the
144 /// preconditioner and deletes what can be deleted... Note that
145 /// this preconditioner needs a CRDoubleMatrix.
146 //============================================================================
149  {
150  // For debugging...
151  bool doc_block_matrices=false;
152 
153  // For output timing results - to be removed soon. Ray
154  bool raytime_flag = false;
155 
156  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
157  // NOTE: In the interest of minimising memory usage, several containers
158  // are recycled, therefore their content/meaning changes
159  // throughout this function. The code is carefully annotated
160  // but you'll have to read it line by line!
161  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
162  double t_clean_up_memory_start = TimingHelpers::timer();
163  // make sure any old data is deleted
164  clean_up_memory();
165  double t_clean_up_memory_end = TimingHelpers::timer();
166  double clean_up_memory_time = t_clean_up_memory_end
167  - t_clean_up_memory_start;
168  if(raytime_flag)
169  {
170  oomph_info << "LSC: clean_up_memory_time: "
171  << clean_up_memory_time << std::endl;
172  }
173 
174 
175 #ifdef PARANOID
176  // paranoid check that the navier stokes mesh pt has been set
177  if (Navier_stokes_mesh_pt == 0)
178  {
179  std::ostringstream error_message;
180  error_message << "The navier stokes elements mesh pointer must be set.\n"
181  << "Use method set_navier_stokes_mesh(...)";
182  throw OomphLibError(error_message.str(),
183  OOMPH_CURRENT_FUNCTION,
184  OOMPH_EXCEPTION_LOCATION);
185  }
186 #endif
187 
188 
189  // Set the mesh
190  this->set_nmesh(1);
193 
194  // Get blocks
195  // ----------
196 
197  // In comes the current Jacobian. Recast it to a CR double matrix;
198  // shout if that can't be done.
199  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
200 
201 
202 #ifdef PARANOID
203  if (cr_matrix_pt==0)
204  {
205  std::ostringstream error_message;
206  error_message
207  << "NavierStokesSchurComplementPreconditioner only works with "
208  << "CRDoubleMatrix matrices" << std::endl;
209  throw OomphLibError(error_message.str(),
210  OOMPH_CURRENT_FUNCTION,
211  OOMPH_EXCEPTION_LOCATION);
212  }
213 #endif
214 
215 
216  if (doc_block_matrices)
217  {
218  std::stringstream junk;
219  junk << "j_matrix" << comm_pt()->my_rank()
220  << ".dat";
221  oomph_info << "About to output " << junk.str() << std::endl;
222  cr_matrix_pt->sparse_indexed_output_with_offset(junk.str());
223  oomph_info << "Done output of " << junk.str() << std::endl;
224  }
225 
226 
227  // Set up block look up schemes (done automatically in the
228  // BlockPreconditioner base class, based on the information
229  // provided in the block-preconditionable elements in the problem)
230 
231  // this preconditioner has two types of block:
232  // type 0: velocity - corresponding to DOFs 0 to n-2
233  // type 1: pressure - corresponding to DOF n-1
234  double t_block_setup_start = TimingHelpers::timer();
235  unsigned ndof_types = 0;
236 
238  {
239  ndof_types = this->ndof_types();
240  }
241  else
242  {
243  // This is the upper-most master block preconditioner, the Navier-Stokes
244  // mesh is in position 0
245  ndof_types = this->ndof_types_in_mesh(0);
246  }
247 
248  Vector<unsigned> dof_to_block_map(ndof_types);
249  dof_to_block_map[ndof_types-1]=1;
250 
251  this->block_setup(dof_to_block_map);
252 
253  double t_block_setup_finish = TimingHelpers::timer();
254  double block_setup_time = t_block_setup_finish - t_block_setup_start;
255  if(Doc_time)
256  {
257  oomph_info << "Time for block_setup(...) [sec]: "
258  << block_setup_time << "\n";
259  }
260 
261  if(raytime_flag)
262  {
263  oomph_info << "LSC: block_setup: " << block_setup_time << std::endl;
264  }
265 
266 
267  // determine whether the F preconditioner is a block preconditioner (and
268  // therefore a subsidiary preconditioner)
269  BlockPreconditioner<CRDoubleMatrix>* F_block_preconditioner_pt =
272  if (F_block_preconditioner_pt == 0)
273  {
274  F_preconditioner_is_block_preconditioner = false;
275  }
276 
277  // Get B (the divergence block)
278  double t_get_B_start = TimingHelpers::timer();
279  CRDoubleMatrix* b_pt = new CRDoubleMatrix;
280  this->get_block(1,0,*b_pt);
281 
282  double t_get_B_finish = TimingHelpers::timer();
283  double get_B_time = t_get_B_finish - t_get_B_start;
284  if(Doc_time)
285  {
286  oomph_info << "Time to get B [sec]: "
287  << get_B_time << "\n";
288  }
289 
290  if(raytime_flag)
291  {
292  oomph_info << "LSC: get block B get_B_time: " << get_B_time << std::endl;
293  }
294 
295  if (doc_block_matrices)
296  {
297  std::stringstream junk;
298  junk << "b_matrix" << comm_pt()->my_rank()
299  << ".dat";
300  b_pt->sparse_indexed_output_with_offset(junk.str());
301  oomph_info << "Done output of " << junk.str() << std::endl;
302  }
303 
304 
305  // get the inverse velocity and pressure mass matrices
306  CRDoubleMatrix* inv_v_mass_pt = 0;
307  CRDoubleMatrix* inv_p_mass_pt = 0;
308 
309  double ivmm_assembly_start_t = TimingHelpers::timer();
310  if (Use_LSC)
311  {
312  // We only need the velocity mass matrix
314  inv_v_mass_pt,
315  false);
316  }
317  else
318  {
319  // We need both mass matrices
321  inv_v_mass_pt,
322  true);
323  }
324 
325  double ivmm_assembly_finish_t = TimingHelpers::timer();
326 
327  double
328  ivmm_assembly_time = ivmm_assembly_finish_t - ivmm_assembly_start_t;
329  if (Doc_time)
330  {
331 
332  oomph_info << "Time to assemble inverse diagonal velocity and pressure"
333  << "mass matrices) [sec]: "
334  << ivmm_assembly_time << "\n";
335  }
336  if(raytime_flag)
337  {
338  oomph_info << "LSC: ivmm_assembly_time: "
339  << ivmm_assembly_time << std::endl;
340  }
341 
342 
343  if (doc_block_matrices)
344  {
345  std::stringstream junk;
346  junk << "inv_v_mass_matrix"
347  << comm_pt()->my_rank()
348  << ".dat";
349  inv_v_mass_pt->sparse_indexed_output_with_offset(junk.str());
350  oomph_info << "Done output of " << junk.str() << std::endl;
351  }
352 
353 
354  // Get gradient matrix Bt
355  CRDoubleMatrix* bt_pt = new CRDoubleMatrix;
356  double t_get_Bt_start = TimingHelpers::timer();
357  this->get_block(0,1,*bt_pt);
358  double t_get_Bt_finish = TimingHelpers::timer();
359 
360  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
361  if(Doc_time)
362  {
363  oomph_info << "Time to get Bt [sec]: "
364  << t_get_Bt_time << std::endl;
365  }
366  if(raytime_flag)
367  {
368  oomph_info << "LSC: get block Bt: "
369  << t_get_Bt_time << std::endl;
370  }
371 
372  if (doc_block_matrices)
373  {
374  std::stringstream junk;
375  junk << "bt_matrix" << comm_pt()->my_rank()
376  << ".dat";
377  bt_pt->sparse_indexed_output_with_offset(junk.str());
378  oomph_info << "Done output of " << junk.str() << std::endl;
379  }
380 
381 
382  // Build pressure Poisson matrix
383  CRDoubleMatrix* p_matrix_pt = new CRDoubleMatrix;
384 
385  // Multiply inverse velocity mass matrix by gradient matrix B^T
386  double t_QBt_matrix_start = TimingHelpers::timer();
387  CRDoubleMatrix* qbt_pt = new CRDoubleMatrix;
388  inv_v_mass_pt->multiply(*bt_pt, *qbt_pt);
389  delete bt_pt; bt_pt = 0;
390 
391  // Store product in bt_pt
392  bt_pt = qbt_pt;
393  double t_QBt_matrix_finish = TimingHelpers::timer();
394 
395  double t_QBt_time = t_QBt_matrix_finish - t_QBt_matrix_start;
396  if(Doc_time)
397  {
398  oomph_info << "Time to generate QBt [sec]: "
399  << t_QBt_time << std::endl;
400  }
401  delete inv_v_mass_pt; inv_v_mass_pt = 0;
402  if(raytime_flag)
403  {
404  oomph_info << "LSC: t_QBt_time (matrix multiplicaton): "
405  << t_QBt_time << std::endl;
406  }
407 
408  // Multiply B from left by divergence matrix B and store result in
409  // pressure Poisson matrix.
410  double t_p_matrix_start = TimingHelpers::timer();
411  b_pt->multiply(*bt_pt, *p_matrix_pt);
412  double t_p_matrix_finish = TimingHelpers::timer();
413 
414  double t_p_time = t_p_matrix_finish - t_p_matrix_start;
415  if(Doc_time)
416  {
417  oomph_info << "Time to generate P [sec]: "
418  << t_p_time << std::endl;
419  }
420  // Kill divergence matrix because we don't need it any more
421  delete b_pt; b_pt = 0;
422 
423  if(raytime_flag)
424  {
425  oomph_info << "LSC: t_p_time (matrix multiplication): "
426  << t_p_time << std::endl;
427  }
428 
429 
430  // Build the matvec operator for QBt
431  double t_QBt_MV_start = TimingHelpers::timer();
434  double t_QBt_MV_finish = TimingHelpers::timer();
435 
436  double t_p_time2 = t_QBt_MV_finish - t_QBt_MV_start;
437  if(Doc_time)
438  {
439  oomph_info << "Time to build QBt matrix vector operator [sec]: "
440  << t_p_time2 << std::endl;
441  }
442 
443  // Kill gradient matrix B^T (it's been overwritten anyway and
444  // needs to be recomputed afresh below)
445  delete bt_pt; bt_pt = 0;
446 
447  if(raytime_flag)
448  {
449  oomph_info << "LSC: QBt (setup MV product): " << t_p_time2 << std::endl;
450  }
451 
452  // Do we need the Fp stuff?
453  if (!Use_LSC)
454  {
455  // Get pressure advection diffusion matrix Fp and store in
456  // a "big" matrix (same size as the problem's Jacobian)
457  double t_get_Fp_start = TimingHelpers::timer();
458  CRDoubleMatrix full_fp_matrix;
460 
461  // Now extract the pressure pressure block
462  CRDoubleMatrix* fp_matrix_pt = new CRDoubleMatrix;
463  this->get_block_other_matrix(1,1,&full_fp_matrix,*fp_matrix_pt);
464  double t_get_Fp_finish = TimingHelpers::timer();
465  if(Doc_time)
466  {
467  double t_get_Fp_time = t_get_Fp_finish - t_get_Fp_start;
468  oomph_info << "Time to get Fp [sec]: "
469  << t_get_Fp_time << std::endl;
470  }
471 
472  // Build vector product of pressure advection diffusion matrix with
473  // inverse pressure mass matrix
474  CRDoubleMatrix* fp_qp_inv_pt = new CRDoubleMatrix;
475  fp_matrix_pt->multiply(*inv_p_mass_pt, *fp_qp_inv_pt);
476 
477  // Build the matvec operator for E = F_p Q_p^{-1}
478  double t_Fp_Qp_inv_MV_start = TimingHelpers::timer();
480  this->setup_matrix_vector_product(E_mat_vec_pt,fp_qp_inv_pt,1);
481  double t_Fp_Qp_inv_MV_finish = TimingHelpers::timer();
482  if(Doc_time)
483  {
484  double t_p_time = t_Fp_Qp_inv_MV_finish - t_Fp_Qp_inv_MV_start;
485  oomph_info << "Time to build Fp Qp^{-1} matrix vector operator [sec]: "
486  << t_p_time << std::endl;
487  }
488  // Kill pressure advection diffusion and inverse pressure mass matrices
489  delete inv_p_mass_pt; inv_p_mass_pt = 0;
490  delete fp_qp_inv_pt; fp_qp_inv_pt = 0;
491  }
492 
493 
494  // Get momentum block F
495  CRDoubleMatrix* f_pt = new CRDoubleMatrix;
496  double t_get_F_start = TimingHelpers::timer();
497  this->get_block(0,0,*f_pt);
498  double t_get_F_finish = TimingHelpers::timer();
499 
500  double t_get_F_time = t_get_F_finish - t_get_F_start;
501  if(Doc_time)
502  {
503  oomph_info << "Time to get F [sec]: "
504  << t_get_F_time << std::endl;
505  }
506  if(raytime_flag)
507  {
508  oomph_info << "LSC: get_block t_get_F_time: "
509  << t_get_F_time << std::endl;
510  }
511 
512  // form the matrix vector product helper
513  double t_F_MV_start = TimingHelpers::timer();
516  double t_F_MV_finish = TimingHelpers::timer();
517 
518  double t_F_MV_time = t_F_MV_finish - t_F_MV_start;
519  if(Doc_time)
520  {
521  oomph_info << "Time to build F Matrix Vector Operator [sec]: "
522  << t_F_MV_time << std::endl;
523  }
524  if(raytime_flag)
525  {
526  oomph_info << "LSC: MV product setup t_F_MV_time: "
527  << t_F_MV_time << std::endl;
528  }
529 
530 
531  // if F is a block preconditioner then we can delete the F matrix
532  if (F_preconditioner_is_block_preconditioner)
533  {
534  delete f_pt; f_pt = 0;
535  }
536 
537  // Rebuild Bt (remember that we temporarily overwrote
538  // it by its product with the inverse velocity mass matrix)
539  t_get_Bt_start = TimingHelpers::timer();
540  bt_pt = new CRDoubleMatrix;
541  this->get_block(0,1,*bt_pt);
542  t_get_Bt_finish = TimingHelpers::timer();
543  double t_get_Bt_time2 = t_get_Bt_finish - t_get_Bt_start;
544  if(Doc_time)
545  {
546 
547  oomph_info << "Time to get Bt [sec]: "
548  << t_get_Bt_time2 << std::endl;
549  }
550  if(raytime_flag)
551  {
552  oomph_info << "LSC: get_block t_get_Bt_time2: "
553  << t_get_Bt_time2 << std::endl;
554  }
555 
556 
557  // form the matrix vector operator for Bt
558  double t_Bt_MV_start = TimingHelpers::timer();
561 
562 // if(Doc_time)
563 // {
564 // oomph_info << "Time to build Bt Matrix Vector Operator [sec]: "
565 // << t_Bt_MV_time << std::endl;
566 // }
567 
568  delete bt_pt; bt_pt = 0;
569 
570  double t_Bt_MV_finish = TimingHelpers::timer();
571 
572  double t_Bt_MV_time = t_Bt_MV_finish - t_Bt_MV_start;
573  if(raytime_flag)
574  {
575  oomph_info << "LSC: MV product setup t_Bt_MV_time: "
576  << t_Bt_MV_time << std::endl;
577  }
578 
579  // if the P preconditioner has not been setup
580  if (P_preconditioner_pt == 0)
581  {
584  }
585 
586  // Setup the preconditioner for the Pressure matrix
587  double t_p_prec_start = TimingHelpers::timer();
588 
589  if (doc_block_matrices)
590  {
591  std::stringstream junk;
592  junk << "p_matrix" << comm_pt()->my_rank()
593  << ".dat";
594  p_matrix_pt->sparse_indexed_output_with_offset(junk.str());
595  oomph_info << "Done output of " << junk.str() << std::endl;
596  }
597 
598  P_preconditioner_pt->setup(p_matrix_pt);
599  delete p_matrix_pt; p_matrix_pt = 0;
600  double t_p_prec_finish = TimingHelpers::timer();
601 
602  double t_p_prec_time = t_p_prec_finish - t_p_prec_start;
603  if(Doc_time)
604  {
605  oomph_info << "P sub-preconditioner setup time [sec]: "
606  << t_p_prec_time << "\n";
607  }
608  if(raytime_flag)
609  {
610  oomph_info << "LSC: p_prec setup time: " << t_p_prec_time << std::endl;
611  }
612 
613 
614  // Set up solver for solution of system with momentum matrix
615  // ----------------------------------------------------------
616 
617  // if the F preconditioner has not been setup
618  if (F_preconditioner_pt == 0)
619  {
622  }
623 
624  // if F is a block preconditioner
625  double t_f_prec_start = TimingHelpers::timer();
626  if (F_preconditioner_is_block_preconditioner)
627  {
628  unsigned nvelocity_dof_types
630 
631  Vector<unsigned> dof_map(nvelocity_dof_types);
632  for (unsigned i = 0; i < nvelocity_dof_types; i++)
633  {
634  dof_map[i] = i;
635  }
636 
637  F_block_preconditioner_pt->
639 
640  F_block_preconditioner_pt->setup(matrix_pt());
641  }
642  // otherwise F is not a block preconditioner
643  else
644  {
645  F_preconditioner_pt->setup(f_pt);
646  delete f_pt; f_pt = 0;
647  }
648  double t_f_prec_finish = TimingHelpers::timer();
649  double t_f_prec_time = t_f_prec_finish - t_f_prec_start;
650  if(Doc_time)
651  {
652 
653  oomph_info << "F sub-preconditioner setup time [sec]: "
654  << t_f_prec_time << "\n";
655  }
656  if(raytime_flag)
657  {
658  oomph_info << "LSC: f_prec setup time: " << t_f_prec_time << std::endl;
659  }
660 
661  // Remember that the preconditioner has been setup so
662  // the stored information can be wiped when we
663  // come here next...
665  }
666 
667 
668 
669 //=======================================================================
670  /// Apply preconditioner to r.
671 //=======================================================================
674  {
675 #ifdef PARANOID
677  {
678  std::ostringstream error_message;
679  error_message << "setup must be called before using preconditioner_solve";
680  throw OomphLibError(
681  error_message.str(),
682  OOMPH_CURRENT_FUNCTION,
683  OOMPH_EXCEPTION_LOCATION);
684  }
685  if (z.built())
686  {
687  if (z.nrow() != r.nrow())
688  {
689  std::ostringstream error_message;
690  error_message << "The vectors z and r must have the same number of "
691  << "of global rows";
692  throw OomphLibError(
693  error_message.str(),
694  OOMPH_CURRENT_FUNCTION,
695  OOMPH_EXCEPTION_LOCATION);
696  }
697  }
698 #endif
699 
700  // if z is not setup then give it the same distribution
701  if (!z.distribution_pt()->built())
702  {
703  z.build(r.distribution_pt(),0.0);
704  }
705 
706  // Step 1 - apply approximate Schur inverse to pressure unknowns (block 1)
707  // -----------------------------------------------------------------------
708 
709  // Working vectors
710  DoubleVector temp_vec;
711  DoubleVector another_temp_vec;
712  DoubleVector yet_another_temp_vec;
713 
714  // Copy pressure values from residual vector to temp_vec:
715  // Loop over all entries in the global vector (this one
716  // includes velocity and pressure dofs in some random fashion)
717  this->get_block_vector(1,r,temp_vec);
718 
719  // NOTE: The vector temp_vec now contains the vector r_p.
720 
721  // LSC version
722  if (Use_LSC)
723  {
724  // Solve first pressure Poisson system
725 #ifdef PARANOID
726  // check a solver has been set
727  if (P_preconditioner_pt==0)
728  {
729  std::ostringstream error_message;
730  error_message << "P_preconditioner_pt has not been set.";
731  throw OomphLibError(
732  error_message.str(),
733  OOMPH_CURRENT_FUNCTION,
734  OOMPH_EXCEPTION_LOCATION);
735  }
736 #endif
737 
738  // use some Preconditioner's preconditioner_solve function
739  P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
740 
741  // NOTE: The vector another_temp_vec now contains the vector P^{-1} r_p
742 
743  // Multiply another_temp_vec by matrix E and stick the result into temp_vec
744  temp_vec.clear();
745  QBt_mat_vec_pt->multiply(another_temp_vec, temp_vec);
746  another_temp_vec.clear();
747  F_mat_vec_pt->multiply(temp_vec,another_temp_vec);
748  temp_vec.clear();
749  QBt_mat_vec_pt->multiply_transpose(another_temp_vec, temp_vec);
750 
751 
752  // NOTE: The vector temp_vec now contains E P^{-1} r_p
753 
754  // Solve second pressure Poisson system using preconditioner_solve
755  another_temp_vec.clear();
756  P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
757 
758  // NOTE: The vector another_temp_vec now contains z_p = P^{-1} E P^{-1} r_p
759  // as required (apart from the sign which we'll fix in the
760  // next step.
761  }
762  // Fp version
763  else
764  {
765 
766  // Multiply temp_vec by matrix E and stick the result into
767  // yet_another_temp_vec
768  E_mat_vec_pt->multiply(temp_vec,yet_another_temp_vec);
769 
770  // NOTE: The vector yet_another_temp_vec now contains Fp Qp^{-1} r_p
771 
772  // Solve pressure Poisson system
773 #ifdef PARANOID
774  // check a solver has been set
775  if (P_preconditioner_pt==0)
776  {
777  std::ostringstream error_message;
778  error_message << "P_preconditioner_pt has not been set.";
779  throw OomphLibError(
780  error_message.str(),
781  OOMPH_CURRENT_FUNCTION,
782  OOMPH_EXCEPTION_LOCATION);
783  }
784 #endif
785 
786  // Solve second pressure Poisson system using preconditioner_solve
787  another_temp_vec.clear();
788  P_preconditioner_pt->preconditioner_solve(yet_another_temp_vec,
789  another_temp_vec);
790 
791  // NOTE: The vector another_temp_vec now contains
792  // z_p = P^{-1} Fp Qp^{-1} r_p
793  // as required (apart from the sign which we'll fix in the
794  // next step.
795 
796  }
797 
798  // Now copy another_temp_vec (i.e. z_p) back into the global vector z.
799  // Loop over all entries in the global results vector z:
800  temp_vec.build(another_temp_vec.distribution_pt(),0.0);
801  temp_vec -= another_temp_vec;
802  return_block_vector(1,temp_vec,z);
803 
804 
805  // Step 2 - apply preconditioner to velocity unknowns (block 0)
806  // ------------------------------------------------------------
807 
808  // Recall that another_temp_vec (computed above) contains the
809  // negative of the solution of the Schur complement systen, -z_p.
810  // Multiply by G (stored in Block_matrix_pt(0,1) and store
811  // result in temp_vec (vector resizes itself).
812  temp_vec.clear();
813  Bt_mat_vec_pt->multiply(another_temp_vec, temp_vec);
814 
815  // NOTE: temp_vec now contains -G z_p
816 
817  // The vector another_temp_vec is no longer needed -- re-use it to store
818  // velocity quantities:
819  another_temp_vec.clear();
820 
821  // Loop over all enries in the global vector and find the
822  // entries associated with the velocities:
823  get_block_vector(0,r,another_temp_vec);
824  another_temp_vec += temp_vec;
825 
826  // NOTE: The vector another_temp_vec now contains r_u - G z_p
827 
828  // Solve momentum system
829 #ifdef PARANOID
830  // check a solver has been set
831  if (F_preconditioner_pt==0)
832  {
833  std::ostringstream error_message;
834  error_message << "F_preconditioner_pt has not been set.";
835  throw OomphLibError(
836  error_message.str(),
837  OOMPH_CURRENT_FUNCTION,
838  OOMPH_EXCEPTION_LOCATION);
839  }
840 #endif
841 
842  // use some Preconditioner's preconditioner solve
843  // and return
845  {
846  return_block_vector(0,another_temp_vec,z);
848  }
849  else
850  {
851  F_preconditioner_pt->preconditioner_solve(another_temp_vec, temp_vec);
852  return_block_vector(0,temp_vec,z);
853  }
854  }
855 
856 
857 //========================================================================
858 /// Helper function to assemble the inverse diagonals of the pressure and
859 /// velocity mass matrix from the elemental contributions defined in
860 /// NavierStokesElementWithDiagonalMassMatrices::
861 /// get_pressure_and_velocity_mass_matrix_diagonal(...)
862 /// If do_both=true, both are computed, otherwise only the velocity
863 /// mass matrix (the LSC version of the preconditioner only needs
864 /// that one)
865 //========================================================================
868  CRDoubleMatrix*& inv_p_mass_pt,
869  CRDoubleMatrix*& inv_v_mass_pt,
870  const bool& do_both)
871  {
872 
873  // determine the velocity rows required by this processor
874  unsigned v_first_row = this->block_distribution_pt(0)->first_row();
875  unsigned v_nrow_local = this->block_distribution_pt(0)->nrow_local();
876  unsigned v_nrow = this->block_distribution_pt(0)->nrow();
877 
878  // create storage for the diagonals
879  double* v_values = new double[v_nrow_local];
880  for (unsigned i = 0; i < v_nrow_local; i++)
881  {
882  v_values[i] = 0.0;
883  }
884 
885  // Equivalent information for pressure mass matrix (only needed for
886  // Fp version)
887  unsigned p_first_row=0;
888  unsigned p_nrow_local=0;
889  unsigned p_nrow=0;
890  double* p_values = 0;
891  if (!Use_LSC)
892  {
893  // determine the pressure rows required by this processor
894  p_first_row = this->block_distribution_pt(1)->first_row();
895  p_nrow_local = this->block_distribution_pt(1)->nrow_local();
896  p_nrow = this->block_distribution_pt(1)->nrow();
897 
898  // create storage for the diagonals
899  p_values = new double[p_nrow_local];
900  for (unsigned i = 0; i < p_nrow_local; i++)
901  {
902  p_values[i] = 0.0;
903  }
904  }
905 
906  // if the problem is distributed
907  bool distributed = false;
908 #ifdef OOMPH_HAS_MPI
909  if (problem_pt()->distributed() ||
911  {
912  distributed = true;
913  }
914 #endif
915 
916  // next we get the diagonal velocity mass matrix data
917  if (distributed)
918  {
919 
920 #ifdef OOMPH_HAS_MPI
921 
922  // the number of processors
923  unsigned nproc = comm_pt()->nproc();
924 
925  // and my rank
926  unsigned my_rank = comm_pt()->my_rank();
927 
928  // determine the rows for which we have lookup rows
929 
930  // if the problem is NOT distributed then we only classify global equations
931  // on this processor to avoid duplication (as every processor holds
932  // every element)
933  unsigned first_lookup_row = 0;
934  unsigned last_lookup_row = 0;
935  first_lookup_row = this->master_distribution_pt()->first_row();
936  last_lookup_row = first_lookup_row +
937  this->master_distribution_pt()->nrow_local() - 1;
938 
939  // find number of local elements
940  unsigned n_el = Navier_stokes_mesh_pt->nelement();
941 
942  // get the master distribution pt
944  this->master_distribution_pt();
945 
946  // Do the two blocks (0: veloc; 1: press)
947  unsigned max_block=0;
948  if (!Use_LSC) max_block=1;
949  for (unsigned block_index=0;block_index<=max_block;block_index++)
950  {
951 
952  // Local working variables: Default to velocity
953  unsigned v_or_p_first_row=v_first_row;
954  double* v_or_p_values=v_values;
955  // Switch to pressure
956  if (block_index==1)
957  {
958  v_or_p_first_row=p_first_row;
959  v_or_p_values=p_values;
960  }
961 
962 
963  // the diagonal mass matrix contributions that have been
964  // classified and should be sent to another processor
965  Vector<double>* classified_contributions_send
966  = new Vector<double>[nproc];
967 
968  // the corresponding block indices
969  Vector<unsigned>* classified_indices_send
970  = new Vector<unsigned>[nproc];
971 
972  // the matrix contributions that cannot be classified by this processor
973  // and therefore must be sent to another for classification
974  Vector<double>* unclassified_contributions_send
975  = new Vector<double>[nproc];
976 
977  // the corresponding global indices that require classification
978  Vector<unsigned>* unclassified_indices_send
979  = new Vector<unsigned>[nproc];
980 
981  // get the velocity or pressure distribution pt
982  const LinearAlgebraDistribution* velocity_or_press_dist_pt
983  = this->block_distribution_pt(block_index);
984 
985  // get the contribution for each element
986  for (unsigned e = 0; e < n_el; e++)
987  {
988 
989  // Get element
991 
992  // check that the element is not halo
993  if (!el_pt->is_halo())
994  {
995 
996  // find number of degrees of freedom in the element
997  // (this is slightly too big because it includes the
998  // pressure dofs but this doesn't matter)
999  unsigned el_dof = el_pt->ndof();
1000 
1001  // Allocate local storage for the element's contribution to the
1002  // mass matrix diagonal
1003  Vector<double> el_vmm_diagonal(el_dof,0.0);
1004  Vector<double> el_pmm_diagonal(el_dof,0.0);
1005 
1006  unsigned which_one=2;
1007  if (block_index==1) which_one=1;
1008 
1010  cast_el_pt=dynamic_cast<NavierStokesElementWithDiagonalMassMatrices*>
1011  (el_pt);
1012  if (cast_el_pt!=0)
1013  {
1015  el_pmm_diagonal,el_vmm_diagonal,which_one);
1016  }
1017 
1018  // get the contribution for each dof
1019  for (unsigned i = 0; i < el_dof; i++)
1020  {
1021 
1022  //Get the equation number
1023  unsigned eqn_number = el_pt->eqn_number(i);
1024 
1025  // if I have lookup information on this processor
1026  if ((eqn_number >= first_lookup_row) &&
1027  (eqn_number <= last_lookup_row) )
1028  {
1029 
1030  // Only use the dofs that we're dealing with here
1031  if ( this->block_number(eqn_number)==int(block_index) )
1032  {
1033 
1034  // get the index in the block
1035  unsigned index = this->index_in_block(eqn_number);
1036 
1037  // determine which processor requires the block index
1038  for (unsigned p = 0; p < nproc; p++)
1039  {
1040  if ( (index >= velocity_or_press_dist_pt->first_row(p)) &&
1041  (index < (velocity_or_press_dist_pt->first_row(p)
1042  +velocity_or_press_dist_pt->nrow_local(p)) ) )
1043  {
1044 
1045  // if it is required by this processor then add the
1046  // contribution
1047  if (p == my_rank)
1048  {
1049  if (block_index==0)
1050  {
1051  v_or_p_values[index-v_or_p_first_row]
1052  += el_vmm_diagonal[i];
1053  }
1054  else if (block_index==1)
1055  {
1056  v_or_p_values[index-v_or_p_first_row]
1057  += el_pmm_diagonal[i];
1058  }
1059  }
1060  // otherwise store it for communication
1061  else
1062  {
1063  if (block_index==0)
1064  {
1065  classified_contributions_send[p]
1066  .push_back(el_vmm_diagonal[i]);
1067  classified_indices_send[p].push_back(index);
1068  }
1069  else if (block_index==1)
1070  {
1071  classified_contributions_send[p]
1072  .push_back(el_pmm_diagonal[i]);
1073  classified_indices_send[p].push_back(index);
1074  }
1075  }
1076  }
1077  }
1078  }
1079  }
1080  // if we do not have the lookup information on this processor
1081  // then we send the mass matrix contribution to a processor
1082  // which we know to have the lookup information
1083  // the assumption: the processor for which the master block
1084  // preconditioner distribution will definitely hold the lookup
1085  // data for eqn_number (although others may)
1086  else if (problem_pt()->distributed())
1087  {
1088 
1089  // determine which processor requires the block index
1090  unsigned p = 0;
1091  while (!(eqn_number >=master_distribution_pt->first_row(p) &&
1092  (eqn_number<(master_distribution_pt->first_row(p)
1093  +master_distribution_pt->nrow_local(p)))))
1094  {
1095  p++;
1096  }
1097 
1098  // store the data
1099  if (block_index==0)
1100  {
1101  unclassified_contributions_send[p]
1102  .push_back(el_vmm_diagonal[i]);
1103  unclassified_indices_send[p].push_back(eqn_number);
1104  }
1105  else if (block_index==1)
1106  {
1107  unclassified_contributions_send[p]
1108  .push_back(el_pmm_diagonal[i]);
1109  unclassified_indices_send[p].push_back(eqn_number);
1110  }
1111 
1112  }
1113  }
1114  }
1115  }
1116 
1117  //next the unclassified contributions are communicated to
1118  //processors that can classify them
1119 
1120  //first determine how many unclassified rows are to be sent to
1121  //each processor
1122  unsigned* n_unclassified_send = new unsigned[nproc];
1123  for (unsigned p = 0; p < nproc; p++)
1124  {
1125  if (p == my_rank)
1126  {
1127  n_unclassified_send[p] = 0;
1128  }
1129  else
1130  {
1131  n_unclassified_send[p]
1132  = unclassified_contributions_send[p].size();
1133  }
1134  }
1135 
1136  //then all-to-all com number of unclassified to be sent / recv
1137  unsigned* n_unclassified_recv = new unsigned[nproc];
1138  MPI_Alltoall(n_unclassified_send,1,MPI_UNSIGNED,
1139  n_unclassified_recv,1,MPI_UNSIGNED,
1140  comm_pt()->mpi_comm());
1141 
1142  //the base displacement for the sends
1143  MPI_Aint base_displacement;
1144  MPI_Address(v_or_p_values,&base_displacement);
1145 
1146  //allocate storage for the data to be received
1147  //and post the sends and recvs
1148  Vector<double*> unclassified_contributions_recv(nproc);
1149  Vector<unsigned*> unclassified_indices_recv(nproc);
1150  Vector<MPI_Request> unclassified_recv_requests;
1151  Vector<MPI_Request> unclassified_send_requests;
1152  Vector<unsigned> unclassified_recv_proc;
1153  for (unsigned p = 0; p < nproc; p++)
1154  {
1155  if (p != my_rank)
1156  {
1157  //recv
1158  if (n_unclassified_recv[p] > 0)
1159  {
1160  unclassified_contributions_recv[p]
1161  = new double[n_unclassified_recv[p]];
1162  unclassified_indices_recv[p] = new
1163  unsigned[n_unclassified_recv[p]];
1164 
1165  //data for the struct data type
1166  MPI_Datatype recv_types[2];
1167  MPI_Aint recv_displacements[2];
1168  int recv_sz[2];
1169 
1170  //contributions
1171  MPI_Type_contiguous(n_unclassified_recv[p],MPI_DOUBLE,
1172  &recv_types[0]);
1173  MPI_Type_commit(&recv_types[0]);
1174  MPI_Address(unclassified_contributions_recv[p],
1175  &recv_displacements[0]);
1176  recv_displacements[0] -= base_displacement;
1177  recv_sz[0] = 1;
1178 
1179  //indices
1180  MPI_Type_contiguous(n_unclassified_recv[p],MPI_UNSIGNED,
1181  &recv_types[1]);
1182  MPI_Type_commit(&recv_types[1]);
1183  MPI_Address(unclassified_indices_recv[p],
1184  &recv_displacements[1]);
1185  recv_displacements[1] -= base_displacement;
1186  recv_sz[1] = 1;
1187 
1188  //build the final recv type
1189  MPI_Datatype final_recv_type;
1190  MPI_Type_struct(2,recv_sz,recv_displacements,recv_types,
1191  &final_recv_type);
1192  MPI_Type_commit(&final_recv_type);
1193 
1194  //and recv
1195  MPI_Request req;
1196  MPI_Irecv(v_or_p_values,1,final_recv_type,p,0,
1197  comm_pt()->mpi_comm(),&req);
1198  unclassified_recv_requests.push_back(req);
1199  unclassified_recv_proc.push_back(p);
1200  MPI_Type_free(&recv_types[0]);
1201  MPI_Type_free(&recv_types[1]);
1202  MPI_Type_free(&final_recv_type);
1203  }
1204 
1205  //send
1206  if (n_unclassified_send[p] > 0)
1207  {
1208  //data for the struct data type
1209  MPI_Datatype send_types[2];
1210  MPI_Aint send_displacements[2];
1211  int send_sz[2];
1212 
1213  //contributions
1214  MPI_Type_contiguous(n_unclassified_send[p],MPI_DOUBLE,
1215  &send_types[0]);
1216  MPI_Type_commit(&send_types[0]);
1217  MPI_Address(&unclassified_contributions_send[p][0],
1218  &send_displacements[0]);
1219  send_displacements[0] -= base_displacement;
1220  send_sz[0] = 1;
1221 
1222  //indices
1223  MPI_Type_contiguous(n_unclassified_send[p],MPI_UNSIGNED,
1224  &send_types[1]);
1225  MPI_Type_commit(&send_types[1]);
1226  MPI_Address(&unclassified_indices_send[p][0],
1227  &send_displacements[1]);
1228  send_displacements[1] -= base_displacement;
1229  send_sz[1] = 1;
1230 
1231  //build the final send type
1232  MPI_Datatype final_send_type;
1233  MPI_Type_struct(2,send_sz,send_displacements,send_types,
1234  &final_send_type);
1235  MPI_Type_commit(&final_send_type);
1236 
1237  //and send
1238  MPI_Request req;
1239  MPI_Isend(v_or_p_values,1,final_send_type,p,0,
1240  comm_pt()->mpi_comm(),&req);
1241  unclassified_send_requests.push_back(req);
1242  MPI_Type_free(&send_types[0]);
1243  MPI_Type_free(&send_types[1]);
1244  MPI_Type_free(&final_send_type);
1245  }
1246  }
1247  }
1248 
1249  //next classify the data as it is received
1250  unsigned n_unclassified_recv_req = unclassified_recv_requests.size();
1251  while (n_unclassified_recv_req > 0)
1252  {
1253  //get the processor number and remove the completed request
1254  //for the vector of requests
1255  int req_num;
1256  MPI_Waitany(n_unclassified_recv_req,&unclassified_recv_requests[0],
1257  &req_num,MPI_STATUS_IGNORE);
1258  unsigned p = unclassified_recv_proc[req_num];
1259  unclassified_recv_requests.erase(unclassified_recv_requests.begin()
1260  +req_num);
1261  unclassified_recv_proc.erase(unclassified_recv_proc.begin()+req_num);
1262  n_unclassified_recv_req--;
1263 
1264  //next classify the dofs
1265  //and store them for sending to other processors if required
1266  unsigned n_recv = n_unclassified_recv[p];
1267  for (unsigned i = 0; i < n_recv; i++)
1268  {
1269  unsigned eqn_number = unclassified_indices_recv[p][i];
1270  //Only deal with our block unknowns
1271  if ( this->block_number(eqn_number)==int(block_index) )
1272  {
1273 
1274  //get the index in the block
1275  unsigned index = this->index_in_block(eqn_number);
1276 
1277  //determine which processor requires the block index
1278  for (unsigned pp = 0; pp < nproc; pp++)
1279  {
1280 
1281 
1282  if ( (index >= velocity_or_press_dist_pt->first_row(pp)) &&
1283  (index < (velocity_or_press_dist_pt->first_row(pp)
1284  +velocity_or_press_dist_pt->nrow_local(pp)) ) )
1285  {
1286 
1287  //if it is required by this processor then add the
1288  //contribution
1289  if (pp == my_rank)
1290  {
1291  v_or_p_values[index-v_or_p_first_row]
1292  += unclassified_contributions_recv[p][i];
1293  }
1294  //otherwise store it for communication
1295  else
1296  {
1297  double v = unclassified_contributions_recv[p][i];
1298  classified_contributions_send[pp].push_back(v);
1299  classified_indices_send[pp].push_back(index);
1300  }
1301  }
1302  }
1303  }
1304  }
1305 
1306  //clean up
1307  delete[] unclassified_contributions_recv[p];
1308  delete[] unclassified_indices_recv[p];
1309  }
1310  delete[] n_unclassified_recv;
1311 
1312  //now all indices have been classified
1313 
1314  //next the classified contributions are communicated to
1315  //processors that require them
1316 
1317  //first determine how many classified rows are to be sent to
1318  //each processor
1319  unsigned* n_classified_send = new unsigned[nproc];
1320  for (unsigned p = 0; p < nproc; p++)
1321  {
1322  if (p == my_rank)
1323  {
1324  n_classified_send[p] = 0;
1325  }
1326  else
1327  {
1328  n_classified_send[p]
1329  = classified_contributions_send[p].size();
1330  }
1331  }
1332 
1333  //then all-to-all number of classified to be sent / recv
1334  unsigned* n_classified_recv = new unsigned[nproc];
1335  MPI_Alltoall(n_classified_send,1,MPI_UNSIGNED,
1336  n_classified_recv,1,MPI_UNSIGNED,
1337  comm_pt()->mpi_comm());
1338 
1339  //allocate storage for the data to be received
1340  //and post the sends and recvs
1341  Vector<double*> classified_contributions_recv(nproc);
1342  Vector<unsigned*> classified_indices_recv(nproc);
1343  Vector<MPI_Request> classified_recv_requests;
1344  Vector<MPI_Request> classified_send_requests;
1345  Vector<unsigned> classified_recv_proc;
1346  for (unsigned p = 0; p < nproc; p++)
1347  {
1348  if (p != my_rank)
1349  {
1350  //recv
1351  if (n_classified_recv[p] > 0)
1352  {
1353  classified_contributions_recv[p]
1354  = new double[n_classified_recv[p]];
1355  classified_indices_recv[p] = new unsigned[n_classified_recv[p]];
1356 
1357  //data for the struct data type
1358  MPI_Datatype recv_types[2];
1359  MPI_Aint recv_displacements[2];
1360  int recv_sz[2];
1361 
1362  //contributions
1363  MPI_Type_contiguous(n_classified_recv[p],MPI_DOUBLE,
1364  &recv_types[0]);
1365  MPI_Type_commit(&recv_types[0]);
1366  MPI_Address(classified_contributions_recv[p],
1367  &recv_displacements[0]);
1368  recv_displacements[0] -= base_displacement;
1369  recv_sz[0] = 1;
1370 
1371  //indices
1372  MPI_Type_contiguous(n_classified_recv[p],MPI_UNSIGNED,
1373  &recv_types[1]);
1374  MPI_Type_commit(&recv_types[1]);
1375  MPI_Address(classified_indices_recv[p],
1376  &recv_displacements[1]);
1377  recv_displacements[1] -= base_displacement;
1378  recv_sz[1] = 1;
1379 
1380  //build the final recv type
1381  MPI_Datatype final_recv_type;
1382  MPI_Type_struct(2,recv_sz,recv_displacements,recv_types,
1383  &final_recv_type);
1384  MPI_Type_commit(&final_recv_type);
1385 
1386  //and recv
1387  MPI_Request req;
1388  MPI_Irecv(v_or_p_values,1,final_recv_type,p,0,
1389  comm_pt()->mpi_comm(),&req);
1390  classified_recv_requests.push_back(req);
1391  classified_recv_proc.push_back(p);
1392  MPI_Type_free(&recv_types[0]);
1393  MPI_Type_free(&recv_types[1]);
1394  MPI_Type_free(&final_recv_type);
1395  }
1396 
1397  //send
1398  if (n_classified_send[p] > 0)
1399  {
1400  //data for the struct data type
1401  MPI_Datatype send_types[2];
1402  MPI_Aint send_displacements[2];
1403  int send_sz[2];
1404 
1405  //contributions
1406  MPI_Type_contiguous(n_classified_send[p],MPI_DOUBLE,
1407  &send_types[0]);
1408  MPI_Type_commit(&send_types[0]);
1409  MPI_Address(&classified_contributions_send[p][0],
1410  &send_displacements[0]);
1411  send_displacements[0] -= base_displacement;
1412  send_sz[0] = 1;
1413 
1414  //indices
1415  MPI_Type_contiguous(n_classified_send[p],MPI_UNSIGNED,
1416  &send_types[1]);
1417  MPI_Type_commit(&send_types[1]);
1418  MPI_Address(&classified_indices_send[p][0],
1419  &send_displacements[1]);
1420  send_displacements[1] -= base_displacement;
1421  send_sz[1] = 1;
1422 
1423  //build the final send type
1424  MPI_Datatype final_send_type;
1425  MPI_Type_struct(2,send_sz,send_displacements,send_types,
1426  &final_send_type);
1427  MPI_Type_commit(&final_send_type);
1428 
1429  //and send
1430  MPI_Request req;
1431  MPI_Isend(v_or_p_values,1,final_send_type,p,0,
1432  comm_pt()->mpi_comm(),&req);
1433  classified_send_requests.push_back(req);
1434  MPI_Type_free(&send_types[0]);
1435  MPI_Type_free(&send_types[1]);
1436  MPI_Type_free(&final_send_type);
1437  }
1438  }
1439  }
1440 
1441  //next classify the data as it is received
1442  unsigned n_classified_recv_req = classified_recv_requests.size();
1443  while (n_classified_recv_req > 0)
1444  {
1445  //get the processor number and remove the completed request
1446  //for the vector of requests
1447  int req_num;
1448  MPI_Waitany(n_classified_recv_req,&classified_recv_requests[0],
1449  &req_num,MPI_STATUS_IGNORE);
1450  unsigned p = classified_recv_proc[req_num];
1451  classified_recv_requests.erase(classified_recv_requests.begin()
1452  +req_num);
1453  classified_recv_proc.erase(classified_recv_proc.begin()+req_num);
1454  n_classified_recv_req--;
1455 
1456  //next classify the dofs
1457  //and store them for sending to other processors if required
1458  unsigned n_recv = n_classified_recv[p];
1459  for (unsigned i = 0; i < n_recv; i++)
1460  {
1461  v_or_p_values[classified_indices_recv[p][i]-v_or_p_first_row]
1462  += classified_contributions_recv[p][i];
1463  }
1464 
1465  //clean up
1466  delete[] classified_contributions_recv[p];
1467  delete[] classified_indices_recv[p];
1468  }
1469 
1470  //wait for the unclassified sends to complete
1471  unsigned n_unclassified_send_req = unclassified_send_requests.size();
1472  if (n_unclassified_send_req > 0)
1473  {
1474  MPI_Waitall(n_unclassified_send_req,&unclassified_send_requests[0],
1475  MPI_STATUS_IGNORE);
1476  }
1477  delete[] unclassified_contributions_send;
1478  delete[] unclassified_indices_send;
1479  delete[] n_unclassified_send;
1480 
1481  //wait for the classified sends to complete
1482  unsigned n_classified_send_req = classified_send_requests.size();
1483  if (n_classified_send_req > 0)
1484  {
1485  MPI_Waitall(n_classified_send_req,&classified_send_requests[0],
1486  MPI_STATUS_IGNORE);
1487  }
1488  delete[] classified_indices_send;
1489  delete[] classified_contributions_send;
1490  delete[] n_classified_recv;
1491  delete[] n_classified_send;
1492 
1493  // Copy the values back where they belong
1494  if (block_index==0)
1495  {
1496  v_values=v_or_p_values;
1497  }
1498  else if (block_index==1)
1499  {
1500  p_values=v_or_p_values;
1501  }
1502 
1503  }
1504 
1505 #endif
1506 
1507  }
1508  // or if the problem is not distributed
1509  else
1510  {
1511 
1512  // find number of elements
1513  unsigned n_el = Navier_stokes_mesh_pt->nelement();
1514 
1515  // Fp needs pressure and velocity mass matrices
1516  unsigned which_one=0;
1517  if (Use_LSC) which_one=2;
1518 
1519  // get the contribution for each element
1520  for (unsigned e = 0; e < n_el; e++)
1521  {
1522 
1523  // Get element
1525 
1526  // find number of degrees of freedom in the element
1527  // (this is slightly too big because it includes the
1528  // pressure dofs but this doesn't matter)
1529  unsigned el_dof = el_pt->ndof();
1530 
1531  // allocate local storage for the element's contribution to the
1532  // pressure and velocity mass matrix diagonal
1533  Vector<double> el_vmm_diagonal(el_dof,0.0);
1534  Vector<double> el_pmm_diagonal(el_dof,0.0);
1535 
1537  cast_el_pt=dynamic_cast<NavierStokesElementWithDiagonalMassMatrices*>(
1538  el_pt);
1539  if (cast_el_pt!=0)
1540  {
1542  el_pmm_diagonal,el_vmm_diagonal,which_one);
1543  }
1544  else
1545  {
1546 #ifdef PARANOID
1547  std::ostringstream error_message;
1548  error_message
1549  << "Navier-Stokes mesh contains element that is not of type \n"
1550  << "NavierStokesElementWithDiagonalMassMatrices. \n"
1551  << "The element is in fact of type "
1552  << typeid(*el_pt).name()
1553  << "\nWe'll assume that it does not make a used contribution \n"
1554  << "to the inverse diagonal mass matrix used in the preconditioner\n"
1555  << "and (to avoid divisions by zero) fill in dummy unit entries.\n"
1556  << "[This case currently arises with flux control problems\n"
1557  << "where for simplicity the NetFluxControlElement has been added \n"
1558  << "to the Navier Stokes mesh -- this should probably be changed at\n"
1559  << "some point -- if you get this warning in any other context\n"
1560  << "you should check your code very carefully]\n";
1562  error_message.str(),
1563  "NavierStokesSchurComplementPreconditioner::assemble_inv_press_and_veloc_mass_matrix_diagonal()",
1564  OOMPH_EXCEPTION_LOCATION);
1565 #endif
1566 
1567  // Fill in dummy entries to stop division by zero below
1568  for (unsigned j=0;j<el_dof;j++)
1569  {
1570  el_vmm_diagonal[j]=1.0;
1571  el_pmm_diagonal[j]=1.0;
1572  }
1573  }
1574 
1575  // Get the contribution for each dof
1576  for (unsigned i = 0; i < el_dof; i++)
1577  {
1578  //Get the equation number
1579  unsigned eqn_number = el_pt->eqn_number(i);
1580 
1581  // Get the velocity dofs
1582  if (this->block_number(eqn_number)==0)
1583  {
1584  // get the index in the block
1585  unsigned index = this->index_in_block(eqn_number);
1586 
1587  // if it is required on this processor
1588  if ((index >= v_first_row) &&
1589  (index < (v_first_row + v_nrow_local) ) )
1590  {
1591  v_values[index-v_first_row] += el_vmm_diagonal[i];
1592  }
1593  }
1594  // Get the pressure dofs
1595  else if (this->block_number(eqn_number)==1)
1596  {
1597  if (!Use_LSC)
1598  {
1599  // get the index in the block
1600  unsigned index = this->index_in_block(eqn_number);
1601 
1602  // if it is required on this processor
1603  if ((index >= p_first_row)&&
1604  (index < (p_first_row + p_nrow_local)) )
1605  {
1606  p_values[index-p_first_row] += el_pmm_diagonal[i];
1607  }
1608  }
1609  }
1610  }
1611  }
1612  }
1613 
1614  // Create column index and row start for velocity mass matrix
1615  int* v_column_index = new int[v_nrow_local];
1616  int* v_row_start = new int[v_nrow_local+1];
1617  for (unsigned i = 0; i < v_nrow_local; i++)
1618  {
1619 #ifdef PARANOID
1620  if (v_values[i]==0.0)
1621  {
1622  std::ostringstream error_message;
1623  error_message << "Zero entry in diagonal of velocity mass matrix\n"
1624  << "Index: " << i << std::endl;
1625  throw OomphLibError(
1626  error_message.str(),
1627  OOMPH_CURRENT_FUNCTION,
1628  OOMPH_EXCEPTION_LOCATION);
1629  }
1630 #endif
1631  v_values[i] = 1.0/v_values[i];
1632  v_column_index[i] = v_first_row + i;
1633  v_row_start[i] = i;
1634  }
1635  v_row_start[v_nrow_local] = v_nrow_local;
1636 
1637  // Build the velocity mass matrix
1638  inv_v_mass_pt = new CRDoubleMatrix(this->block_distribution_pt(0));
1639  inv_v_mass_pt->build_without_copy(v_nrow,v_nrow_local,
1640  v_values,v_column_index,
1641  v_row_start);
1642 
1643  // Create pressure mass matrix
1644  if (!Use_LSC)
1645  {
1646  // Create column index and row start for pressure mass matrix
1647  int* p_column_index = new int[p_nrow_local];
1648  int* p_row_start = new int[p_nrow_local+1];
1649  for (unsigned i = 0; i < p_nrow_local; i++)
1650  {
1651 
1652 #ifdef PARANOID
1653  if (p_values[i]==0.0)
1654  {
1655  std::ostringstream error_message;
1656  error_message << "Zero entry in diagonal of pressure mass matrix\n"
1657  << "Index: " << i << std::endl;
1658  throw OomphLibError(
1659  error_message.str(),
1660  OOMPH_CURRENT_FUNCTION,
1661  OOMPH_EXCEPTION_LOCATION);
1662  }
1663 #endif
1664  p_values[i] = 1.0/p_values[i];
1665 
1666  p_column_index[i] = p_first_row + i;
1667  p_row_start[i] = i;
1668  }
1669  p_row_start[p_nrow_local] = p_nrow_local;
1670 
1671  // Build the pressure mass matrix
1672  inv_p_mass_pt = new CRDoubleMatrix(this->block_distribution_pt(1));
1673  inv_p_mass_pt->build_without_copy(p_nrow,p_nrow_local,
1674  p_values,p_column_index,
1675  p_row_start);
1676 
1677  }
1678 
1679  }
1680 
1681 //=========================================================================
1682 /// Helper function to delete preconditioner data.
1683 //=========================================================================
1685  {
1687  {
1688  // delete matvecs
1689  delete Bt_mat_vec_pt;
1690  Bt_mat_vec_pt = 0;
1691 
1692  delete F_mat_vec_pt;
1693  F_mat_vec_pt = 0;
1694 
1695  delete QBt_mat_vec_pt;
1696  QBt_mat_vec_pt = 0;
1697 
1698  delete E_mat_vec_pt;
1699  E_mat_vec_pt = 0;
1700 
1701  // delete stuff from velocity solve
1703  {
1704  delete F_preconditioner_pt;
1705  F_preconditioner_pt = 0;
1706  }
1707 
1708  // delete stuff from Schur complement approx
1710  {
1711  delete P_preconditioner_pt;
1712  P_preconditioner_pt = 0;
1713  }
1714  }
1715  }
1716 
1717 
1718 }// end oomph namespace
A Generalised Element class.
Definition: elements.h:76
bool built() const
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to Vector r.
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero...
MatrixVectorProduct * Bt_mat_vec_pt
MatrixVectorProduct operator for Bt.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
void clear()
wipes the DoubleVector
MatrixVectorProduct * QBt_mat_vec_pt
MatrixVectorProduct operator for Qv^{-1} Bt.
virtual void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)=0
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
void clean_up_memory()
Helper function to delete preconditioner data.
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
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
void sparse_indexed_output_with_offset(std::string filename)
Indexed output function to print a matrix to a file as i,j,a(i,j) for a(i,j)!=0 only. Specify filename. This uses acual global row numbers.
Definition: matrices.h:1021
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 Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
bool 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...
double Peclet
Peclet number – overwrite with actual Reynolds number.
OomphInfo oomph_info
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
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
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
bool Using_default_p_preconditioner
flag indicating whether the default P preconditioner is used
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
bool Doc_time
Set Doc_time to true for outputting results of timings.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void get_pressure_advection_diffusion_matrix(CRDoubleMatrix &fp_matrix)
Get the pressure advection diffusion matrix.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
MatrixVectorProduct * F_mat_vec_pt
MatrixVectorProduct operator for F.
void multiply(const DoubleVector &x, DoubleVector &y) const
Apply the operator to the vector x and return the result in the vector y.
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).
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
void get_block_other_matrix(const unsigned &i, const unsigned &j, CRDoubleMatrix *in_matrix_pt, CRDoubleMatrix &output_matrix)
Get a block from a different matrix using the blocking scheme that has already been set up...
unsigned nrow() const
access function to the number of global rows.
void set_nmesh(const unsigned &n)
Specify the number of meshes required by this block preconditioner. Note: elements in different meshe...
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.
unsigned ndof_types() const
Return the total number of DOF types.
An interface to allow SuperLU to be used as an (exact) Preconditioner.
bool Allow_multiple_element_type_in_navier_stokes_mesh
Flag to indicate if there are multiple element types in the Navier-Stokes mesh.
MatrixVectorProduct * E_mat_vec_pt
MatrixVectorProduct operator for E = Fp Qp^{-1} (only for Fp variant)
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 "...
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
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...
bool Use_LSC
Boolean to indicate use of LSC (true) or Fp (false) variant.
bool F_preconditioner_is_block_preconditioner
Boolean indicating whether the momentum system preconditioner is a block preconditioner.
Mesh * Navier_stokes_mesh_pt
the pointer to the mesh of block preconditionable Navier Stokes elements.
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
void assemble_inv_press_and_veloc_mass_matrix_diagonal(CRDoubleMatrix *&inv_p_mass_pt, CRDoubleMatrix *&inv_v_mass_pt, const bool &do_both)
Helper function to assemble the inverse diagonals of the pressure and velocity mass matrices from the...
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...
bool Using_default_f_preconditioner
flag indicating whether the default F preconditioner is used
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
double source_function(const Vector< double > &x_vect)
Source function required to make the solution above an exact solution.
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...
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
void wind_function(const Vector< double > &x, Vector< double > &wind)
Wind.