general_purpose_space_time_subsidiary_block_preconditioner.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// 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//====================================================================
26// Header file
28
29#ifdef OOMPH_HAS_HYPRE
30#include "../../generic/hypre_solver.h"
31#endif
32
33/// /////////////////////////////////////////////////////////////////////
34/// /////////////////////////////////////////////////////////////////////
35/// /////////////////////////////////////////////////////////////////////
36
37namespace oomph
38{
39#ifdef OOMPH_HAS_HYPRE
40 //======start_of_HypreSubsidiaryPreconditionerHelper_namespace==============
41 /// Helper method for the block diagonal F block preconditioner to allow
42 /// hypre to be used as a subsidiary block preconditioner
43 //==========================================================================
44 namespace HypreSubsidiaryPreconditionerHelper
45 {
46 /// Assign the Hypre preconditioner pointer
48 {
49 return new HyprePreconditioner;
50 }
51 } // namespace HypreSubsidiaryPreconditionerHelper
52#endif
53
54 /// /////////////////////////////////////////////////////////////////////
55 /// /////////////////////////////////////////////////////////////////////
56 /// /////////////////////////////////////////////////////////////////////
57
58 //============================================================================
59 /// Setup for the block triangular preconditioner
60 //============================================================================
62 {
63 // For debugging...
64 bool doc_block_matrices = false;
65
66 // Clean up any previously created data
68
69 // If we're meant to build silently
70 if (this->Silent_preconditioner_setup == true)
71 {
72 // Store the output stream pointer
74
75 // Now set the oomph_info stream pointer to the null stream to
76 // disable all possible output
78 }
79
80#ifdef PARANOID
81 // If the upcast was unsuccessful
82 if (dynamic_cast<CRDoubleMatrix*>(this->matrix_pt()) == 0)
83 {
84 // Allocate space for an error message
85 std::ostringstream error_message;
86
87 // Create the error message
88 error_message << "NavierStokesSchurComplementPreconditioner only works "
89 << "with CRDoubleMatrix matrices" << std::endl;
90
91 // Throw the error message
92 throw OomphLibError(
93 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
94 }
95#endif
96
97 // Start the timer for the block setup
98 double t_block_setup_start = TimingHelpers::timer();
99
100 // Set up block look up schemes (done automatically in the
101 // BlockPreconditioner base class, based on the information
102 // provided in the block-preconditionable elements in the problem)
103
104 // This preconditioner has two types of block:
105 // Type 0: Velocity - Corresponding to DOF type 0 & 1
106 // Type 1: Pressure - Corresponding to DOF type 2
107 unsigned n_dof_types = 0;
108
109 // If this has been set up as a subsidiary preconditioner
111 {
112 // Get the number of dof types (will be told by the master preconditioner)
113 n_dof_types = this->ndof_types();
114
115 // DRAIG: Delete
116 oomph_info << "Number of dof types: " << n_dof_types << std::endl;
117
118 // Make sure there's only 3 dof types for now!
119 if (n_dof_types != 3)
120 {
121 // Allocate space for an error message
122 std::ostringstream error_message;
123
124 // Create the error message
125 error_message << "Should only be 3 dof types! You have " << n_dof_types
126 << " dof types!" << std::endl;
127
128 // Throw an error
129 throw OomphLibError(error_message.str(),
130 OOMPH_CURRENT_FUNCTION,
131 OOMPH_EXCEPTION_LOCATION);
132 }
133 }
134 else
135 {
136 // Throw an error
137 throw OomphLibError("Can only be used as a subsidiary preconditioner!",
138 OOMPH_CURRENT_FUNCTION,
139 OOMPH_EXCEPTION_LOCATION);
140 } // if (this->is_subsidiary_block_preconditioner())
141
142 // Allocate storage for the dof-to-block mapping
143 Vector<unsigned> dof_to_block_map(n_dof_types);
144
145 // The last dof type (the pressure) should have it's own block type
146 dof_to_block_map[n_dof_types - 1] = 1;
147
148 // Call the generic block setup function
149 this->block_setup(dof_to_block_map);
150
151 // Stop the timer
152 double t_block_setup_finish = TimingHelpers::timer();
153
154 // Calculate the time difference
155 double block_setup_time = t_block_setup_finish - t_block_setup_start;
156
157 // Document the time taken
158 oomph_info << "Time for block_setup(...) [sec]: " << block_setup_time
159 << std::endl;
160
161 // Get the number of block types
162 unsigned n_block_type = this->nblock_types();
163
164 // How many block types are there? Should only be 2...
165 oomph_info << "Number of block types: " << n_block_type << std::endl;
166
167 // Allocate storage for the momentum block
168 CRDoubleMatrix* f_pt = new CRDoubleMatrix;
169
170 // Allocate storage for the divergence block
171 CRDoubleMatrix* d_pt = new CRDoubleMatrix;
172
173 // Allocate storage for the gradient block
174 CRDoubleMatrix* g_pt = new CRDoubleMatrix;
175
176 // Allocate storage for the pressure Poisson matrix
177 CRDoubleMatrix* p_matrix_pt = new CRDoubleMatrix;
178
179 // Get the momentum block F
180 this->get_block(0, 0, *f_pt);
181
182 // Get the divergence block, D
183 this->get_block(1, 0, *d_pt);
184
185 // Get gradient matrix, G=D^T
186 this->get_block(0, 1, *g_pt);
187
188 // Calculate the pressure Poisson-like matrix
189 d_pt->multiply(*g_pt, *p_matrix_pt);
190
191 // If the user wants the block matrices themselves
192 if (doc_block_matrices)
193 {
194 // Filename suffix (defaults should be an empty string)
195 std::string suffix = "";
196
197 // Create space for all of the blocks
198 DenseMatrix<CRDoubleMatrix*> block_matrix_pt(
199 n_block_type, n_block_type, 0);
200
201 // Loop over the block rows
202 for (unsigned i = 0; i < n_block_type; i++)
203 {
204 // Loop over the block columns
205 for (unsigned j = 0; j < n_block_type; j++)
206 {
207 // Allocate space for the matrix
208 block_matrix_pt(i, j) = new CRDoubleMatrix;
209
210 // Get the block
211 this->get_block(i, j, *block_matrix_pt(i, j));
212 }
213 } // for (unsigned i=0;i<n_block_type;i++)
214
215 // Space for the full matrix
216 CRDoubleMatrix* j_pt = new CRDoubleMatrix;
217
218 // Concatenate the sub-blocks to get the full matrix
219 CRDoubleMatrixHelpers::concatenate(block_matrix_pt, *j_pt);
220
221 // Output the matrix in MATLAB form
222 j_pt->sparse_indexed_output_with_offset("j_" + suffix + ".csv");
223
224 // Tell the user
225 oomph_info << "Done output of j_" + suffix + ".csv" << std::endl;
226
227 // Loop over the block rows
228 for (unsigned i = 0; i < n_block_type; i++)
229 {
230 // Loop over the block columns
231 for (unsigned j = 0; j < n_block_type; j++)
232 {
233 // Clear up this block
234 delete block_matrix_pt(i, j);
235
236 // Make it a null pointer
237 block_matrix_pt(i, j) = 0;
238 }
239 } // for (unsigned i=0;i<n_block_type;i++)
240
241 // Now clear up j_pt
242 delete j_pt;
243
244 // Make it a null pointer
245 j_pt = 0;
246
247 // Output the matrix in MATLAB form
248 f_pt->sparse_indexed_output_with_offset("f_matrix" + suffix + ".dat");
249
250 // Tell the user
251 oomph_info << "Done output of f_matrix" + suffix + ".dat" << std::endl;
252
253 // Output the matrix in MATLAB form
254 d_pt->sparse_indexed_output_with_offset("d_matrix" + suffix + ".dat");
255
256 // Tell the user
257 oomph_info << "Done output of d_matrix" + suffix + ".dat" << std::endl;
258
259 // Output the matrix in MATLAB form
260 g_pt->sparse_indexed_output_with_offset("g_matrix" + suffix + ".dat");
261
262 // Tell the user
263 oomph_info << "Done output of g_matrix" + suffix + ".dat" << std::endl;
264
265 // Output the matrix in MATLAB form
266 p_matrix_pt->sparse_indexed_output_with_offset("p_matrix" + suffix +
267 ".dat");
268
269 // Tell the user
270 oomph_info << "Done output of p_matrix" + suffix + ".dat" << std::endl;
271
272 // Exit
273 exit(0);
274 }
275
276 // Kill divergence matrix because we don't need it any more
277 delete d_pt;
278
279 // Make it a null pointer
280 d_pt = 0;
281
282 // Make a new matrix-vector product
284
285 // Make a new matrix-vector product
287
288 // Set the matrix-vector product up
290
291 // Set the matrix-vector product up
293
294 // Do we need to doc. the memory statistics?
296 {
297 // Initialise the memory usage variable
299
300 // How many rows are there in the momentum block?
301 unsigned n_row_f = f_pt->nrow();
302
303 // How many nonzeros are there in the momentum block?
304 unsigned n_nnz_f = f_pt->nnz();
305
306 // Update the memory usage variable. NOTE: This calculation is done by
307 // taking into account the storage scheme for a compressed row matrix
308 Memory_usage_in_bytes += ((2 * ((n_row_f + 1) * sizeof(int))) +
309 (n_nnz_f * (sizeof(int) + sizeof(double))));
310
311 // How many rows are there in the gradient block?
312 unsigned n_row_g = g_pt->nrow();
313
314 // How many nonzeros are there in the gradient block?
315 unsigned n_nnz_g = g_pt->nnz();
316
317 // Update the memory usage variable
318 Memory_usage_in_bytes += ((2 * ((n_row_g + 1) * sizeof(int))) +
319 (n_nnz_g * (sizeof(int) + sizeof(double))));
320 }
321
322 // Kill the gradient matrix because we don't need it any more
323 delete g_pt;
324
325 // Make it a null pointer
326 g_pt = 0;
327
328#ifdef OOMPH_HAS_HYPRE
329 // Assign the preconditioner for the P block
331
332 // Set parameters for use as preconditioner on Poisson-type problem
335
336 // Indicate that we're not using the default preconditioner
338
339 // Assign the preconditioner for the P block
341
342 // Indicate that we're not using the default preconditioner
344#endif
345
346 /*
347 // Assign the preconditioner for the F block
348 F_preconditioner_pt=new BlockDiagonalPreconditioner<CRDoubleMatrix>;
349
350 // Number of velocity dof types
351 unsigned n_velocity_dof_types=2;
352
353 // Create space for the dof-to-block map
354 Vector<unsigned> dof_map(n_velocity_dof_types,0);
355
356 // Loop over the velocity dofs
357 for (unsigned i=0;i<n_velocity_dof_types;i++)
358 {
359 // Seperate the velocity dofs out (by giving each component its own block)
360 dof_map[i]=i;
361 }
362
363 // Turn it into a subsidiary preconditioner
364 dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(F_preconditioner_pt)->
365 turn_into_subsidiary_block_preconditioner(this,dof_map);
366
367 // Pass a pointer to the global system matrix
368 F_preconditioner_pt->setup(this->matrix_pt());
369 */
370
371 // If the P preconditioner has not been set
372 if (P_preconditioner_pt == 0)
373 {
374 // Just use SuperLU
376
377 // Indicate that we're using the default preconditioner
379 } // if (P_preconditioner_pt==0)
380
381 // Compute the LU decomposition of P
382 P_preconditioner_pt->setup(p_matrix_pt);
383
384 // If we need to document the memory statistics
386 {
387 // If we're using SuperLU as the preconditioner
389 {
390 // Add on the memory needed to store and calculate the LU factors of P
393 ->get_total_memory_needed_for_superlu();
394 }
395 // Using Hypre which just really means we just store the matrix P
396 else
397 {
398 // How many rows are there in the momentum block?
399 unsigned n_row_p = p_matrix_pt->nrow();
400
401 // How many nonzeros are there in the momentum block?
402 unsigned n_nnz_p = p_matrix_pt->nnz();
403
404 // Update the memory usage variable. NOTE: This calculation is done by
405 // taking into account the storage scheme for a compressed row matrix
406 Memory_usage_in_bytes += ((2 * ((n_row_p + 1) * sizeof(int))) +
407 (n_nnz_p * (sizeof(int) + sizeof(double))));
408 }
409 } // if (Compute_memory_statistics)
410
411 // Now delete the actual matrix P, we don't need it anymore
412 delete p_matrix_pt;
413
414 // Make it a null pointer
415 p_matrix_pt = 0;
416
417 // If the F preconditioner has not been setup
418 if (F_preconditioner_pt == 0)
419 {
420 // Just use SuperLU
422
423 // Indicate that we're using the default preconditioner
425 } // if (F_preconditioner_pt==0)
426
427 // Compute the LU decomposition of F
429
430 // If we need to document the memory statistics
432 {
433 // If we're using SuperLU as the preconditioner
435 {
436 // Add on the memory needed to store and calculate the LU factors of P
439 ->get_total_memory_needed_for_superlu();
440 }
441 // Using Hypre which just really means we store the matrix F (again, as
442 // it is duplicated inside Hypre)
443 else
444 {
445 // How many rows are there in the momentum block?
446 unsigned n_row_f = f_pt->nrow();
447
448 // How many nonzeros are there in the momentum block?
449 unsigned n_nnz_f = f_pt->nnz();
450
451 // Update the memory usage variable. NOTE: This calculation is done by
452 // taking into account the storage scheme for a compressed row matrix
453 Memory_usage_in_bytes += ((2 * ((n_row_f + 1) * sizeof(int))) +
454 (n_nnz_f * (sizeof(int) + sizeof(double))));
455 }
456 } // if (Compute_memory_statistics)
457
458 // Now delete the matrix, we don't need it anymore
459 delete f_pt;
460
461 // Make it a null pointer
462 f_pt = 0;
463
464 // Remember that the preconditioner has been setup so the stored
465 // information can be wiped when we come here next...
467
468 // If we're meant to build silently, reassign the oomph stream pointer
469 if (this->Silent_preconditioner_setup == true)
470 {
471 // Store the output stream pointer
473
474 // Reset our own stream pointer
475 this->Stream_pt = 0;
476 }
477 } // End of setup
478
479 //=============================================================================
480 /// Preconditioner solve for the block triangular preconditioner
481 //=============================================================================
483 const DoubleVector& r, DoubleVector& z)
484 {
485 // If z is not set up then give it the same distribution
486 if (!z.distribution_pt()->built())
487 {
488 // Build z with the same distribution as r
489 z.build(r.distribution_pt(), 0.0);
490 }
491
492 // -----------------------------------------------------------------------
493 // Step 1 - apply approximate Schur inverse to pressure unknowns (block 1)
494 // -----------------------------------------------------------------------
495 // First working vectors
496 DoubleVector temp_vec;
497
498 // Second working vector
499 DoubleVector another_temp_vec;
500
501 // Copy pressure values from residual vector to temp_vec:
502 // Loop over all entries in the global vector (this one
503 // includes velocity and pressure dofs in some random fashion)
504 this->get_block_vector(1, r, temp_vec);
505
506 // Compute the vector P^{-1}r_p
507 P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
508
509 // Now clear the vector storing r_p
510 temp_vec.clear();
511
512 // Now compute G(P^{-1}r_p)
513 G_mat_vec_pt->multiply(another_temp_vec, temp_vec);
514
515 // Clear the vector storing P^{-1}r_p
516 another_temp_vec.clear();
517
518 // Now update to get F(GP^{-1}r_p)
519 F_mat_vec_pt->multiply(temp_vec, another_temp_vec);
520
521 // Clear the vector storing GP^{-1}r_p
522 temp_vec.clear();
523
524 // Now compute D(FGP^{-1}r_p)=G^{T}(FGP^{-1}r_p)
525 G_mat_vec_pt->multiply_transpose(another_temp_vec, temp_vec);
526
527 // Clear the vector storing FGP^{-1}r_p
528 another_temp_vec.clear();
529
530 // Finally, compute the full approximation to -z_p;
531 // -z_p=P^{-1}(DFG)P^{-1}r_p NOTE: We'll move the sign over in the next step
532 P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
533
534 // Rebuild temp_vec with another_temp_vec's distribution with all
535 // entries initialised to zero
536 temp_vec.build(another_temp_vec.distribution_pt(), 0.0);
537
538 // Now fix the sign and store the result in temp_vec
539 temp_vec -= another_temp_vec;
540
541 // Now copy temp_vec (i.e. z_p) back into the global vector z
542 return_block_vector(1, temp_vec, z);
543
544 // ------------------------------------------------------------
545 // Step 2 - apply preconditioner to velocity unknowns (block 0)
546 // ------------------------------------------------------------
547 // Clear the vector storing z_p=-P^{-1}(DFG)P^{-1}r_p
548 temp_vec.clear();
549
550 // Recall that another_temp_vec (computed above) contains the negative
551 // of the solution of the Schur complement system, -z_p. Multiply by G
552 // and store the result in temp_vec (vector resizes itself).
553 G_mat_vec_pt->multiply(another_temp_vec, temp_vec);
554
555 // Now clear the vector storing P^{-1}(DFG)P^{-1}r_p
556 another_temp_vec.clear();
557
558 // Get the residual vector entries associated with the velocities, r_u
559 get_block_vector(0, r, another_temp_vec);
560
561 // Update the result to get r_u-Gz_p
562 another_temp_vec += temp_vec;
563
564 // Compute the solution z_u which satisfies z_u=F^{-1}(r_u-Gz_p)
565 F_preconditioner_pt->preconditioner_solve(another_temp_vec, temp_vec);
566
567 // Now store the result in the global solution vector
568 return_block_vector(0, temp_vec, z);
569 } // End of preconditioner_solve
570
571
572 //============================================================================
573 /// Setup for the GMRES block preconditioner
574 //============================================================================
576 {
577 // Clean up any previously created data
579
580 // If we're meant to build silently
581 if (this->Silent_preconditioner_setup == true)
582 {
583 // Store the output stream pointer
585
586 // Now set the oomph_info stream pointer to the null stream to
587 // disable all possible output
589 }
590
591 // Start the timer for the block setup
592 double t_block_setup_start = TimingHelpers::timer();
593
594 // Set up block look up schemes (done automatically in the
595 // BlockPreconditioner base class, based on the information
596 // provided in the block-preconditionable elements in the problem)
597
598 // This preconditioner has two types of block:
599 // Type 0: Velocity - Corresponding to DOF type 0 & 1
600 // Type 1: Pressure - Corresponding to DOF type 2
601 unsigned n_dof_types = 0;
602
603 // If this has been set up as a subsidiary preconditioner
605 {
606 // Get the number of dof types (will be told by the master preconditioner)
607 n_dof_types = this->ndof_types();
608
609 // Output the number of dof types
610 oomph_info << "Number of dof types: " << n_dof_types << std::endl;
611
612 // Make sure there's only 3 dof types for now!
613 if (n_dof_types != 3)
614 {
615 // Allocate space for an error message
616 std::ostringstream error_message;
617
618 // Create the error message
619 error_message << "Should only be 3 dof types! You have " << n_dof_types
620 << " dof types!" << std::endl;
621
622 // Throw an error
623 throw OomphLibError(error_message.str(),
624 OOMPH_CURRENT_FUNCTION,
625 OOMPH_EXCEPTION_LOCATION);
626 }
627 }
628 // If it's being used as a master preconditioner
629 else
630 {
631 // Throw an error
632 throw OomphLibError("Currently only used as a subsidiary preconditioner!",
633 OOMPH_CURRENT_FUNCTION,
634 OOMPH_EXCEPTION_LOCATION);
635 } // if (this->is_subsidiary_block_preconditioner())
636
637 // Aggregrate everything together since GMRES itself isn't going to take
638 // advantage of the block structure...
639 Vector<unsigned> dof_to_block_map(n_dof_types, 0);
640
641 // Call the generic block setup function
642 this->block_setup(dof_to_block_map);
643
644 // Stop the timer
645 double t_block_setup_finish = TimingHelpers::timer();
646
647 // Calculate the time difference
648 double block_setup_time = t_block_setup_finish - t_block_setup_start;
649
650 // Document the time taken
651 oomph_info << "Time for block_setup(...) [sec]: " << block_setup_time
652 << std::endl;
653
654 // Get the number of block types
655 unsigned n_block_type = this->nblock_types();
656
657 // How many block types are there? Should only be 2...
658 oomph_info << "Number of block types: " << n_block_type << std::endl;
659
660 // Space for the concatenated block matrix
662
663 // Get the (one and only) block to be used by this block preconditioner
664 this->get_block(0, 0, *Matrix_pt);
665
666 // Create a new instance of the NS subsidiary preconditioner
669
670 // Do we need to doc. the memory statistics?
672 {
673 // How many rows are there in this block?
674 unsigned n_row = Matrix_pt->nrow();
675
676 // How many nonzeros are there in this block?
677 unsigned n_nnz = Matrix_pt->nnz();
678
679 // Compute the memory usage. The only memory overhead here (for the
680 // setup phase) is storing the system matrix
681 Memory_usage_in_bytes = ((2 * ((n_row + 1) * sizeof(int))) +
682 (n_nnz * (sizeof(int) + sizeof(double))));
683
684 // We might as well document the memory statistics of the Navier-Stokes
685 // subsidiary block preconditioner while we're at it...
687 }
688
689 // Create a map to declare which of the 3 dof types in the GMRES
690 // subsidiary preconditioner correspond to the dof types in the NS
691 // subsidiary block preconditioner
692 Vector<unsigned> dof_map(n_dof_types);
693
694 // Loop over the dof types associated with the GMRES subsidiary block
695 // preconditioner
696 for (unsigned i = 0; i < n_dof_types; i++)
697 {
698 // There is, essentially, an identity mapping between the GMRES subsidiary
699 // block preconditioner and the NSSP w.r.t. the dof types
700 dof_map[i] = i;
701 }
702
703 // Now turn it into an "proper" subsidiary preconditioner
706
707 // Setup: The NS subsidiary preconditioner is itself a block preconditioner
708 // so we need to pass it a pointer to full-size matrix!
710
711 // Remember that the preconditioner has been setup so the stored
712 // information can be wiped when we come here next...
714
715 // If we're meant to build silently, reassign the oomph stream pointer
716 if (this->Silent_preconditioner_setup == true)
717 {
718 // Store the output stream pointer
720
721 // Reset our own stream pointer
722 this->Stream_pt = 0;
723 }
724 } // End of setup
725
726 //=============================================================================
727 /// Preconditioner solve for the GMRES block preconditioner
728 //=============================================================================
730 DoubleVector& solution)
731 {
732 // Get the number of block types
733 unsigned n_block_types = this->nblock_types();
734
735 // If there is more than one block type (there shouldn't be!)
736 if (n_block_types != 1)
737 {
738 // Throw an error
739 throw OomphLibError("Can only deal with one dof type!",
740 OOMPH_CURRENT_FUNCTION,
741 OOMPH_EXCEPTION_LOCATION);
742 }
743
744 // Allocate space for the solution blocks
745 Vector<DoubleVector> block_z(n_block_types);
746
747 // Allocate space for the residual blocks
748 Vector<DoubleVector> block_r(n_block_types);
749
750 // Get the reordered RHS vector
751 this->get_block_vectors(rhs, block_r);
752
753 // Get the reordered solution vector
754 this->get_block_vectors(solution, block_z);
755
756 // Initialise the entries in the solution vector to zero
757 block_z[0].initialise(0.0);
758
759 // Get number of dofs
760 unsigned n_dof = block_r[0].nrow();
761
762 // Time solver
763 double t_start = TimingHelpers::timer();
764
765 // Relative residual
766 double resid;
767
768 // Iteration counter
769 unsigned iter = 1;
770
771 // Initialise vectors
772 Vector<double> s(n_dof + 1, 0);
773 Vector<double> cs(n_dof + 1);
774 Vector<double> sn(n_dof + 1);
775 DoubleVector w(block_r[0].distribution_pt(), 0.0);
776
777 // Storage for the time taken for all preconditioner applications
778 double t_prec_application_total = 0.0;
779
780 // Storage for the RHS vector
781 DoubleVector r(block_r[0].distribution_pt(), 0.0);
782
783 // If we're using LHS preconditioning then solve Mr=b-Jx for r (assumes x=0)
785 {
786 // Initialise timer
787 double t_prec_application_start = TimingHelpers::timer();
788
789 // This is pretty inefficient but there is no other way to do this with
790 // block subsidiary preconditioners as far as I can tell because they
791 // demand to have the full r and z vectors...
792 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
793 DoubleVector r_updated(rhs.distribution_pt(), 0.0);
794
795 // Reorder the entries in block_r[0] and put them back into the
796 // global-sized vector. The subsidiary preconditioner should only
797 // operate on the dofs that this preconditioner operates on so
798 // don't worry about all the other entries
799 this->return_block_vector(0, block_r[0], r_updated);
800
801 // Solve on the global-sized vectors
803 r_updated, block_z_with_size_of_full_z);
804
805 // Now grab the part of the solution associated with this preconditioner
806 this->get_block_vector(0, block_z_with_size_of_full_z, block_z[0]);
807
808 // Doc time for setup of preconditioner
809 double t_prec_application_end = TimingHelpers::timer();
810
811 // Update the preconditioner application time total
812 t_prec_application_total +=
813 (t_prec_application_end - t_prec_application_start);
814 }
815 // If we're using RHS preconditioning, do nothing to the actual RHS vector
816 else
817 {
818 // Copy the entries into r
819 r = block_r[0];
820 }
821
822 // Calculate the initial residual norm (assumes x=0 initially)
823 double normb = r.norm();
824
825 // Set beta (the initial residual)
826 double beta = normb;
827
828 // Compute initial relative residual
829 if (normb == 0.0)
830 {
831 // To avoid dividing by zero, set normb to 1
832 normb = 1;
833 }
834
835 // Now calculate the relative residual norm
836 resid = beta / normb;
837
838 // If required, document convergence history to screen or file (if
839 // stream is open)
841 {
842 // If a file has not been provided to output to
843 if (!Output_file_stream.is_open())
844 {
845 // Output the initial relative residual norm to screen
846 oomph_info << 0 << " " << resid << std::endl;
847 }
848 // If there is a file provided to output to
849 else
850 {
851 // Output the initial relative residual norm to file
852 Output_file_stream << 0 << " " << resid << std::endl;
853 }
854 } // if (Doc_convergence_history)
855
856 // If GMRES converges immediately
857 if (resid <= Tolerance)
858 {
859 // If the user wishes GMRES to output information about the solve to
860 // screen
861 if (Doc_time)
862 {
863 // Tell the user
864 oomph_info << "GMRES block preconditioner converged immediately. "
865 << "Normalised residual norm: " << resid << std::endl;
866 }
867
868 // Now reorder the values in block_z[0] and put them back into the
869 // global solution vector
870 this->return_block_vector(0, block_z[0], solution);
871
872 // Doc time for solver
873 double t_end = TimingHelpers::timer();
874 Solution_time = t_end - t_start;
875
876 // If the user wishes GMRES to output information about the solve to
877 // screen
878 if (Doc_time)
879 {
880 // Tell the user
881 oomph_info << "Time for all preconditioner applications [sec]: "
882 << t_prec_application_total
883 << "\nTotal time for solve with GMRES block preconditioner "
884 << "[sec]: " << Solution_time << std::endl;
885 }
886
887 // We're finished; end here!
888 return;
889 } // if (resid<=Tolerance)
890
891 // Initialise vector of orthogonal basis vectors, v
892 Vector<DoubleVector> v(n_dof + 1);
893
894 // Initialise the upper Hessenberg matrix H.
895 // NOTE: For implementation purposes the upper hessenberg matrix indices
896 // are swapped so the matrix is effectively transposed
897 Vector<Vector<double>> H(n_dof + 1);
898
899 // Loop as many times as we're allowed
900 while (iter <= Max_iter)
901 {
902 // Copy the entries of the residual vector, r, into v[0]
903 v[0] = r;
904
905 // Now update the zeroth basis vector, v[0], to have values r/beta
906 v[0] /= beta;
907
908 // Storage for ...?
909 s[0] = beta;
910
911 // Inner iteration counter for restarted version (we don't actually use
912 // a restart in this implementation of GMRES...)
913 unsigned iter_restart;
914
915 // Perform iterations
916 for (iter_restart = 0; iter_restart < n_dof && iter <= Max_iter;
917 iter_restart++, iter++)
918 {
919 // Resize next column of the upper Hessenberg matrix
920 H[iter_restart].resize(iter_restart + 2);
921
922 // Solve for w inside the scope of these braces
923 {
924 // Temporary vector
925 DoubleVector temp(block_r[0].distribution_pt(), 0.0);
926
927 // If we're using LHS preconditioning solve Mw=Jv[i] for w
929 {
930 // Calculate temp=Jv[i]
931 Matrix_pt->multiply(v[iter_restart], temp);
932
933 // Get the current time
934 double t_prec_application_start = TimingHelpers::timer();
935
936 // This is pretty inefficient but there is no other way to do this
937 // with block sub preconditioners as far as I can tell because they
938 // demand to have the full r and z vectors...
939 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(),
940 0.0);
941 DoubleVector r_updated(rhs.distribution_pt(), 0.0);
942
943 // Put the values in temp into the global-sized vector
944 this->return_block_vector(0, temp, r_updated);
945
946 // Apply the NSSP to the global-sized vectors (this will extract
947 // the appropriate sub-vectors and solve on the reduced system)
949 r_updated, block_z_with_size_of_full_z);
950
951 // Now grab the part of the solution associated with this
952 // preconditioner
953 this->get_block_vector(0, block_z_with_size_of_full_z, w);
954
955 // End timer
956 double t_prec_application_end = TimingHelpers::timer();
957
958 // Update the preconditioner application time total
959 t_prec_application_total +=
960 (t_prec_application_end - t_prec_application_start);
961 }
962 // If we're using RHS preconditioning solve
963 else
964 {
965 // Initialise timer
966 double t_prec_application_start = TimingHelpers::timer();
967
968 // This is pretty inefficient but there is no other way to do this
969 // with block sub preconditioners as far as I can tell because they
970 // demand to have the full r and z vectors...
971 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt());
972 DoubleVector r_updated(rhs.distribution_pt());
973
974 // Put the values in v[iter_restart] back into the global-sized
975 // vector
976 this->return_block_vector(0, v[iter_restart], r_updated);
977
978 // Solve on the global-sized vectors
980 r_updated, block_z_with_size_of_full_z);
981
982 // Now grab the part of the solution associated with this
983 // preconditioner
984 this->get_block_vector(0, block_z_with_size_of_full_z, temp);
985
986 // Doc time for setup of preconditioner
987 double t_prec_application_end = TimingHelpers::timer();
988
989 // Update the preconditioner application time total
990 t_prec_application_total +=
991 (t_prec_application_end - t_prec_application_start);
992
993 // Calculate w=Jv_m where v_m=M^{-1}v
994 Matrix_pt->multiply(temp, w);
995 }
996 } // End of solve for w
997
998 // Get a pointer to the values in w
999 double* w_pt = w.values_pt();
1000
1001 // Loop over the columns of the (transposed) Hessenberg matrix
1002 for (unsigned k = 0; k <= iter_restart; k++)
1003 {
1004 // Initialise the entry in the upper Hessenberg matrix
1005 H[iter_restart][k] = 0.0;
1006
1007 // Get the pointer to the values of the k-th basis vector
1008 double* vk_pt = v[k].values_pt();
1009
1010 // Update H with the dot product of w and v[k]
1011 H[iter_restart][k] += w.dot(v[k]);
1012
1013 // Loop over the entries in w and v[k] again
1014 for (unsigned i = 0; i < n_dof; i++)
1015 {
1016 // Now update the values in w
1017 w_pt[i] -= H[iter_restart][k] * vk_pt[i];
1018 }
1019 } // for (unsigned k=0;k<=iter_restart;k++)
1020
1021 // Calculate the (iter_restart+1,iter_restart)-th entry in the upper
1022 // Hessenberg matrix (remember, H is actually the transpose of the
1023 // proper upper Hessenberg matrix)
1024 H[iter_restart][iter_restart + 1] = w.norm();
1025
1026 // Build the (iter_restart+1)-th basis vector by copying w
1027 v[iter_restart + 1] = w;
1028
1029 // Now rescale the basis vector
1030 v[iter_restart + 1] /= H[iter_restart][iter_restart + 1];
1031
1032 // Loop over the columns of the transposed upper Hessenberg matrix
1033 for (unsigned k = 0; k < iter_restart; k++)
1034 {
1035 // Apply a Givens rotation to entries of H
1037 H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
1038 }
1039
1040 // Now generate the cs and sn entries for a plane rotation
1041 generate_plane_rotation(H[iter_restart][iter_restart],
1042 H[iter_restart][iter_restart + 1],
1043 cs[iter_restart],
1044 sn[iter_restart]);
1045
1046 // Use the newly generated entries in cs and sn to apply a Givens
1047 // rotation to those same entries
1048 apply_plane_rotation(H[iter_restart][iter_restart],
1049 H[iter_restart][iter_restart + 1],
1050 cs[iter_restart],
1051 sn[iter_restart]);
1052
1053 // Now apply the Givens rotation to the entries in s
1054 apply_plane_rotation(s[iter_restart],
1055 s[iter_restart + 1],
1056 cs[iter_restart],
1057 sn[iter_restart]);
1058
1059 // Compute the current residual
1060 beta = std::fabs(s[iter_restart + 1]);
1061
1062 // Compute the relative residual norm
1063 resid = beta / normb;
1064
1065 // If required, document the convergence history to screen or file
1067 {
1068 // If an output file has not been provided
1069 if (!Output_file_stream.is_open())
1070 {
1071 // Output the convergence history to screen
1072 oomph_info << iter << " " << resid << std::endl;
1073 }
1074 // An output file has been provided so output to that
1075 else
1076 {
1077 // Output the convergence history to file
1078 Output_file_stream << iter << " " << resid << std::endl;
1079 }
1080 } // if (Doc_convergence_history)
1081
1082 // If the required tolerance was met
1083 if (resid < Tolerance)
1084 {
1085 // Storage for the global-sized solution vector. Strictly speaking
1086 // we could actually just use the vector, solution but this can be
1087 // done after we know we've implemented everything correctly...
1088 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
1089
1090 // Take the current solution, reorder the entries appropriately
1091 // and stick them into the global sized vector
1092 this->return_block_vector(0, block_z[0], block_z_with_size_of_full_z);
1093
1094 // This will update block_z[0] and won't touch the global-size vector
1095 // block_x_with_size_of_full_x. We're in charge of reordering the
1096 // entries and putting it in solution
1097 update(
1098 iter_restart, H, s, v, block_z_with_size_of_full_z, block_z[0]);
1099
1100 // Now go into block_z[0], reorder the entries in the appropriate
1101 // manner and put them into the vector, solution
1102 this->return_block_vector(0, block_z[0], solution);
1103
1104 // If the user wishes GMRES to document the convergence
1105 if (Doc_time)
1106 {
1107 // Document the convergence to screen
1109 << "\nGMRES block preconditioner converged (1). Normalised "
1110 << "residual norm: " << resid
1111 << "\nNumber of iterations to convergence: " << iter << "\n"
1112 << std::endl;
1113 }
1114
1115 // Get the current time
1116 double t_end = TimingHelpers::timer();
1117
1118 // Calculate the time difference, i.e. the time taken for the whole
1119 // solve
1120 Solution_time = t_end - t_start;
1121
1122 // Storage the iteration count
1123 Iterations = iter;
1124
1125 // If the user wishes GMRES to document the convergence
1126 if (Doc_time)
1127 {
1128 // Document the convergence to screen
1130 << "Time for all preconditioner applications [sec]: "
1131 << t_prec_application_total
1132 << "\nTotal time for solve with GMRES block preconditioner "
1133 << "[sec]: " << Solution_time << std::endl;
1134 }
1135
1136 // We're done now; finish here
1137 return;
1138 } // if (resid<Tolerance)
1139 } // for (iter_restart=0;iter_restart<n_dof&&iter<=Max_iter;...
1140
1141 // Update
1142 if (iter_restart > 0)
1143 {
1144 // Storage for the global-sized solution vector
1145 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
1146
1147 // Take the current solution, reorder the entries appropriately
1148 // and stick them into the global sized vector
1149 this->return_block_vector(0, block_z[0], block_z_with_size_of_full_z);
1150
1151 // Update the (reordered block) solution vector
1152 update(
1153 (iter_restart - 1), H, s, v, block_z_with_size_of_full_z, block_z[0]);
1154
1155 // Now go into block_z[0], reorder the entries in the appropriate manner
1156 // and put them into the vector, solution
1157 this->return_block_vector(0, block_z[0], solution);
1158 }
1159
1160 // Solve Mr=(b-Jx) for r
1161 {
1162 // Temporary vector to store b-Jx
1163 DoubleVector temp(block_r[0].distribution_pt(), 0.0);
1164
1165 // Calculate temp=b-Jx
1166 Matrix_pt->residual(block_z[0], block_r[0], temp);
1167
1168 // If we're using LHS preconditioning
1170 {
1171 // Initialise timer
1172 double t_prec_application_start = TimingHelpers::timer();
1173
1174 // This is pretty inefficient but there is no other way to do this
1175 // with block sub preconditioners as far as I can tell because they
1176 // demand to have the full r and z vectors...
1177 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt());
1178 DoubleVector r_updated(rhs.distribution_pt());
1179
1180 // Reorder the entries of temp and put them into the global-sized
1181 // vector
1182 this->return_block_vector(0, temp, r_updated);
1183
1184 // Solve on the global-sized vectors
1186 r_updated, block_z_with_size_of_full_z);
1187
1188 // Now grab the part of the solution associated with this
1189 // preconditioner
1190 this->get_block_vector(0, block_z_with_size_of_full_z, r);
1191
1192 // Doc time for setup of preconditioner
1193 double t_prec_application_end = TimingHelpers::timer();
1194
1195 // Update the preconditioner application time total
1196 t_prec_application_total +=
1197 (t_prec_application_end - t_prec_application_start);
1198 }
1199 } // End of solve Mr=(b-Jx) for r
1200
1201 // Compute the current residual norm
1202 beta = r.norm();
1203
1204 // if relative residual within tolerance
1205 resid = beta / normb;
1206 if (resid < Tolerance)
1207 {
1208 // Now store the result in the global solution vector
1209 this->return_block_vector(0, block_z[0], solution);
1210
1211 if (Doc_time)
1212 {
1214 << "\nGMRES block preconditioner converged (2). Normalised "
1215 << "residual norm: " << resid
1216 << "\nNumber of iterations to convergence: " << iter << "\n"
1217 << std::endl;
1218 }
1219
1220 // Doc time for solver
1221 double t_end = TimingHelpers::timer();
1222 Solution_time = t_end - t_start;
1223
1224 Iterations = iter;
1225
1226 if (Doc_time)
1227 {
1229 << "Time for all preconditioner applications [sec]: "
1230 << t_prec_application_total
1231 << "\nTotal time for solve with GMRES block preconditioner "
1232 << "[sec]: " << Solution_time << std::endl;
1233 }
1234 return;
1235 }
1236 }
1237
1238
1239 // otherwise GMRES failed convergence
1240 oomph_info << "\nGMRES block preconditioner did not converge to required "
1241 << "tolerance! \nReturning with normalised residual norm: "
1242 << resid << "\nafter " << Max_iter << " iterations.\n"
1243 << std::endl;
1244
1246 {
1247 std::ostringstream error_message_stream;
1248 error_message_stream << "Solver failed to converge and you requested an "
1249 << "error on convergence failures.";
1250 throw OomphLibError(error_message_stream.str(),
1251 OOMPH_EXCEPTION_LOCATION,
1252 OOMPH_CURRENT_FUNCTION);
1253 }
1254 } // End of solve_helper
1255} // End of namespace oomph
static char t char * s
Definition: cfortran.h:568
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...
void get_block(const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const
Put block (i,j) into output_matrix. This block accounts for any coarsening of dof types and any repla...
void get_block_vector(const unsigned &n, const DoubleVector &v, DoubleVector &b) const
Takes the naturally ordered vector, v and returns the n-th block vector, b. Here n is the block numbe...
unsigned nblock_types() const
Return the number of block types.
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 get_block_vectors(const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
Takes the naturally ordered vector and rearranges it into a vector of sub vectors corresponding to th...
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*,...
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
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
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1096
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1002
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
Definition: matrices.h:386
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
virtual void residual(const DoubleVector &x, const DoubleVector &b, DoubleVector &residual_)
Find the residual, i.e. r=b-Ax the residual.
Definition: matrices.h:326
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.
double norm() const
compute the 2 norm of this vector
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
double * values_pt()
access function to the underlying values
void clear()
wipes the DoubleVector
double Memory_usage_in_bytes
Storage for the memory usage of the solver if the flag above is set to true (in bytes)
void apply_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Apply plane rotation. This is done using the update:
virtual void clean_up_memory()
Clean up the memory (empty for now...)
SpaceTimeNavierStokesSubsidiaryPreconditioner * Navier_stokes_subsidiary_preconditioner_pt
Pointer to the preconditioner for the block matrix.
bool Compute_memory_statistics
Flag to indicate whether or not to record the memory statistics this preconditioner.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
void update(const unsigned &k, const Vector< Vector< double > > &H, const Vector< double > &s, const Vector< DoubleVector > &v, const DoubleVector &block_x_with_size_of_full_x, DoubleVector &x)
Helper function to update the result vector using the result, x=x_0+V_m*y.
void generate_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Generate a plane rotation. This is done by finding the values of (i....
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: hypre_solver.h:826
double Tolerance
Convergence tolerance.
bool Throw_error_after_max_iter
Should we throw an error instead of just returning when we hit the max iterations?
unsigned Max_iter
Maximum number of iterations.
double Solution_time
linear solver solution time
std::ofstream Output_file_stream
Output file stream for convergence history.
bool Doc_convergence_history
Flag indicating if the convergence history is to be documented.
bool built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:77
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.
std::ostream *& stream_pt()
Access function for the stream pointer.
An OomphLibError object which should be thrown when an run-time error is encountered....
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
std::ostream * Stream_pt
Pointer to the output stream – defaults to std::cout.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
bool Silent_preconditioner_setup
Boolean to indicate whether or not the build should be done silently.
General purpose block triangular preconditioner. By default this is Upper triangular....
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
double Memory_usage_in_bytes
Storage for the memory usage of the solver if the flag above is set to true (in bytes)
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
bool Compute_memory_statistics
Flag to indicate whether or not to record the memory statistics this preconditioner.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
bool Using_default_f_preconditioner
Flag indicating whether the default F preconditioner is used.
bool Using_default_p_preconditioner
Flag indicating whether the default P preconditioner is used.
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
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
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:167
void concatenate(const DenseMatrix< CRDoubleMatrix * > &matrix_pt, CRDoubleMatrix &result_matrix)
Concatenate CRDoubleMatrix matrices. The in matrices are concatenated such that the block structure o...
Definition: matrices.cc:4349
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
Preconditioner * set_hypre_preconditioner()
Assign the Hypre preconditioner pointer.
void set_defaults_for_2D_poisson_problem(HyprePreconditioner *hypre_preconditioner_pt)
Set default parameters for use as preconditioner in 2D Poisson-type problem.
Definition: hypre_solver.cc:45
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Nullstream oomph_nullstream
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...