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