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