hypre_solver.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#include "hypre_solver.h"
27
28// For problem->get_jacobian(...)
29#include "problem.h"
30
31
32namespace oomph
33{
34 /// /////////////////////////////////////////////////////////////////////
35 /// /////////////////////////////////////////////////////////////////////
36 /// /////////////////////////////////////////////////////////////////////
37
38 //==================================================================
39 /// Default settings for various uses of the HYPRE solver
40 //==================================================================
41 namespace Hypre_default_settings
42 {
43 /// Set default parameters for use as preconditioner in
44 /// 2D Poisson-type problem.
46 HyprePreconditioner* hypre_preconditioner_pt)
47 {
48 // Set iterations to 1
49 hypre_preconditioner_pt->set_amg_iterations(1);
50
51 // Use simple smoother
52 hypre_preconditioner_pt->amg_using_simple_smoothing();
53
54 // Smoother types:
55 // 0=Jacobi
56 // 1=Gauss-Seidel
57 hypre_preconditioner_pt->amg_simple_smoother() = 1;
58
59 // AMG preconditioner
60 hypre_preconditioner_pt->hypre_method() = HyprePreconditioner::BoomerAMG;
61
62 // Choose strength parameter for amg
63 hypre_preconditioner_pt->amg_strength() = 0.25;
64
65 // Coarsening type
66 hypre_preconditioner_pt->amg_coarsening() = 0;
67 }
68
69
70 /// Set default parameters for use as preconditioner in
71 /// 3D Poisson-type problem.
73 HyprePreconditioner* hypre_preconditioner_pt)
74 {
75 // Set default settings as for 2D Poisson
76 set_defaults_for_2D_poisson_problem(hypre_preconditioner_pt);
77
78 // Change strength parameter for amg
79 hypre_preconditioner_pt->amg_strength() = 0.7;
80 }
81
82
83 /// Set default parameters for use as preconditioner in
84 /// for momentum block in Navier-Stokes problem
86 HyprePreconditioner* hypre_preconditioner_pt)
87 {
88 // Set default settings as for 2D Poisson
89 set_defaults_for_2D_poisson_problem(hypre_preconditioner_pt);
90
91 // Change smoother type:
92 // 0=Jacobi
93 // 1=Gauss-Seidel
94 hypre_preconditioner_pt->amg_simple_smoother() = 0;
95
96 // Set smoother damping
97 hypre_preconditioner_pt->amg_damping() = 0.5;
98
99 // Change strength parameter for amg
100 hypre_preconditioner_pt->amg_strength() = 0.75;
101 }
102
103 } // namespace Hypre_default_settings
104
105
106 /// /////////////////////////////////////////////////////////////////////
107 /// /////////////////////////////////////////////////////////////////////
108 // functions for HypreHelpers namespace
109 /// /////////////////////////////////////////////////////////////////////
110 /// /////////////////////////////////////////////////////////////////////
111
112 namespace HypreHelpers
113 {
114 //========================================================================
115 /// Default for AMG strength (0.25 recommended for 2D problems;
116 /// larger (0.5-0.75, say) for 3D
117 //========================================================================
118 double AMG_strength = 0.25;
119
120 //========================================================================
121 /// Default AMG coarsening strategy. Coarsening types include:
122 /// 0 = CLJP (parallel coarsening using independent sets)
123 /// 1 = classical RS with no boundary treatment (not recommended
124 /// in parallel)
125 /// 3 = modified RS with 3rd pass to add C points on the boundaries
126 /// 6 = Falgout (uses 1 then CLJP using interior coarse points as
127 /// first independent set)
128 /// 8 = PMIS (parallel coarsening using independent sets - lower
129 /// complexities than 0, maybe also slower convergence)
130 /// 10= HMIS (one pass RS on each processor then PMIS on interior
131 /// coarse points as first independent set)
132 /// 11= One pass RS on each processor (not recommended)
133 //========================================================================
134 unsigned AMG_coarsening = 6;
135
136
137 //========================================================================
138 /// AMG interpolation truncation factor
139 //========================================================================
140 double AMG_truncation = 0.0;
141
142 //========================================================================
143 /// Helper function to check the Hypre error flag, return the message
144 /// associated with any error, and reset the global error flag to zero.
145 /// This function also returns the error value.
146 //========================================================================
147 int check_HYPRE_error_flag(std::ostringstream& message)
148 {
149 // get the Hypre error flag
150 int err = HYPRE_GetError();
151
152 // Tell us all about it...
153 if (err)
154 {
155 oomph_info << "Hypre error flag=" << err << std::endl;
156 char* error_message = new char[128];
157 HYPRE_DescribeError(err, error_message);
158 message << "WARNING: " << std::endl
159 << "HYPRE error message: " << error_message << std::endl;
160 delete[] error_message;
161 }
162
163 // reset Hypre's global error flag
165
166 return err;
167 }
168
169
170 //========================================================================
171 /// Helper function to create a HYPRE_IJVector and HYPRE_ParVector.
172 /// + If no MPI then serial vectors are created
173 /// + If MPI and serial input vector then distributed hypre vectors are
174 /// created
175 /// + If MPI and distributed input vector the distributed output vectors
176 /// are created.
177 //========================================================================
178 void create_HYPRE_Vector(const DoubleVector& oomph_vec,
179 const LinearAlgebraDistribution* dist_pt,
180 HYPRE_IJVector& hypre_ij_vector,
181 HYPRE_ParVector& hypre_par_vector)
182 {
183 // the lower and upper row of the vector on this processor
184 unsigned lower = dist_pt->first_row();
185 unsigned upper = dist_pt->first_row() + dist_pt->nrow_local() - 1;
186
187 // number of local rows
188 unsigned nrow_local = dist_pt->nrow_local();
189
190 // initialize Hypre vector
191#ifdef OOMPH_HAS_MPI
192 HYPRE_IJVectorCreate(
193 oomph_vec.distribution_pt()->communicator_pt()->mpi_comm(),
194 lower,
195 upper,
196 &hypre_ij_vector);
197#else
198 HYPRE_IJVectorCreate(MPI_COMM_WORLD, lower, upper, &hypre_ij_vector);
199#endif
200 HYPRE_IJVectorSetObjectType(hypre_ij_vector, HYPRE_PARCSR);
201 HYPRE_IJVectorInitialize(hypre_ij_vector);
202
203 // set up array containing indices
204 int* indices = new int[nrow_local];
205 double* values = new double[nrow_local];
206 const unsigned hypre_first_row = dist_pt->first_row();
207 unsigned j = 0;
208 if (!oomph_vec.distributed() && dist_pt->distributed())
209 {
210 j = hypre_first_row;
211 }
212 const double* o_pt = oomph_vec.values_pt();
213 for (unsigned i = 0; i < nrow_local; i++)
214 {
215 indices[i] = hypre_first_row + i;
216 values[i] = o_pt[j + i];
217 }
218
219 // insert values
220 HYPRE_IJVectorSetValues(hypre_ij_vector, nrow_local, indices, values);
221
222 // assemble vectors
223 HYPRE_IJVectorAssemble(hypre_ij_vector);
224 HYPRE_IJVectorGetObject(hypre_ij_vector, (void**)&hypre_par_vector);
225
226 // clean up
227 delete[] indices;
228 delete[] values;
229 }
230
231
232 //========================================================================
233 /// Helper function to create a HYPRE_IJVector and HYPRE_ParVector.
234 /// + If no MPI then serial vectors are created
235 /// + If MPI and serial input vector then distributed hypre vectors are
236 /// created
237 /// + If MPI and distributed input vector the distributed output vectors
238 /// are created.
239 //========================================================================
241 HYPRE_IJVector& hypre_ij_vector,
242 HYPRE_ParVector& hypre_par_vector)
243 {
244 // the lower and upper row of the vector on this processor
245 unsigned lower = dist_pt->first_row();
246 unsigned upper = dist_pt->first_row() + dist_pt->nrow_local() - 1;
247
248 // initialize Hypre vector
249#ifdef OOMPH_HAS_MPI
250 HYPRE_IJVectorCreate(
251 dist_pt->communicator_pt()->mpi_comm(), lower, upper, &hypre_ij_vector);
252#else
253 HYPRE_IJVectorCreate(MPI_COMM_WORLD, lower, upper, &hypre_ij_vector);
254#endif
255 HYPRE_IJVectorSetObjectType(hypre_ij_vector, HYPRE_PARCSR);
256 HYPRE_IJVectorInitialize(hypre_ij_vector);
257
258 // assemble vectors
259 HYPRE_IJVectorAssemble(hypre_ij_vector);
260 HYPRE_IJVectorGetObject(hypre_ij_vector, (void**)&hypre_par_vector);
261 }
262
263
264 //========================================================================
265 /// Helper function to create a serial HYPRE_IJMatrix and
266 /// HYPRE_ParCSRMatrix from a CRDoubleMatrix
267 /// NOTE: dist_pt is rebuilt to match the distribution of the hypre solver
268 /// which is not necassarily the same as the oomph lib matrix
269 //========================================================================
271 HYPRE_IJMatrix& hypre_ij_matrix,
272 HYPRE_ParCSRMatrix& hypre_par_matrix,
274 {
275#ifdef PARANOID
276 // check that the matrix is built
277 if (!oomph_matrix->built())
278 {
279 std::ostringstream error_message;
280 error_message << "The matrix has not been built";
281 throw OomphLibError(error_message.str(),
282 OOMPH_CURRENT_FUNCTION,
283 OOMPH_EXCEPTION_LOCATION);
284 }
285 // check the matrix is square
286 if (oomph_matrix->nrow() != oomph_matrix->ncol())
287 {
288 std::ostringstream error_message;
289 error_message << "create_HYPRE_Matrix require a square matrix. "
290 << "Matrix is " << oomph_matrix->nrow() << " by "
291 << oomph_matrix->ncol() << std::endl;
292 throw OomphLibError(error_message.str(),
293 OOMPH_CURRENT_FUNCTION,
294 OOMPH_EXCEPTION_LOCATION);
295 }
296#endif
297
298#ifdef OOMPH_HAS_MPI
299 // Trap the case when we have compiled with MPI,
300 // but are running in serial
302 {
303 std::ostringstream error_stream;
304 error_stream
305 << "Oomph-lib has been compiled with MPI support and "
306 << "you are using HYPRE.\n"
307 << "For this combination of flags, MPI must be initialised.\n"
308 << "Call MPI_Helpers::init() in the "
309 << "main() function of your driver code\n";
310 throw OomphLibError(
311 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
312 }
313#endif
314
315 // find number of rows/columns
316 const unsigned nrow = int(oomph_matrix->nrow());
317
318 // get pointers to the matrix
319
320 // column indices of matrix
321 const int* matrix_cols = oomph_matrix->column_index();
322
323 // entries of matrix
324 const double* matrix_vals = oomph_matrix->value();
325
326 // row starts
327 const int* matrix_row_start = oomph_matrix->row_start();
328
329 // build the distribution
330 if (oomph_matrix->distribution_pt()->distributed())
331 {
332 dist_pt->build(oomph_matrix->distribution_pt());
333 }
334 else
335 {
336 bool distributed = true;
337 if (oomph_matrix->distribution_pt()->communicator_pt()->nproc() == 1)
338 {
339 distributed = false;
340 }
341 dist_pt->build(oomph_matrix->distribution_pt()->communicator_pt(),
342 nrow,
343 distributed);
344 }
345
346 // initialize hypre matrix
347 unsigned lower = dist_pt->first_row();
348 unsigned upper = lower + dist_pt->nrow_local() - 1;
349
350#ifdef OOMPH_HAS_MPI
351 HYPRE_IJMatrixCreate(dist_pt->communicator_pt()->mpi_comm(),
352 lower,
353 upper,
354 lower,
355 upper,
356 &hypre_ij_matrix);
357#else
358 HYPRE_IJMatrixCreate(
359 MPI_COMM_WORLD, lower, upper, lower, upper, &hypre_ij_matrix);
360#endif
361 HYPRE_IJMatrixSetObjectType(hypre_ij_matrix, HYPRE_PARCSR);
362 HYPRE_IJMatrixInitialize(hypre_ij_matrix);
363
364 // set up a row map
365 // and first row / nrow_local
366 const unsigned hypre_nrow_local = dist_pt->nrow_local();
367 const unsigned hypre_first_row = dist_pt->first_row();
368 int* ncols_per_row = new int[hypre_nrow_local];
369 int* row_map = new int[hypre_nrow_local];
370 for (unsigned i = 0; i < hypre_nrow_local; i++)
371 {
372 unsigned j = i;
373 if (!oomph_matrix->distributed() && dist_pt->distributed())
374 {
375 j += hypre_first_row;
376 }
377 ncols_per_row[i] = matrix_row_start[j + 1] - matrix_row_start[j];
378 row_map[i] = hypre_first_row + i;
379 }
380
381 // put values in HYPRE matrix
382 int local_start = 0;
383 if (!oomph_matrix->distributed() && dist_pt->distributed())
384 {
385 local_start += matrix_row_start[hypre_first_row];
386 }
387
388
389 HYPRE_IJMatrixSetValues(hypre_ij_matrix,
390 hypre_nrow_local,
391 ncols_per_row,
392 row_map,
393 matrix_cols + local_start,
394 matrix_vals + local_start);
395
396 // assemble matrix
397 HYPRE_IJMatrixAssemble(hypre_ij_matrix); // hierher leak?
398 HYPRE_IJMatrixGetObject(hypre_ij_matrix, (void**)&hypre_par_matrix);
399
400 // tidy up memory
401 delete[] ncols_per_row;
402 delete[] row_map;
403 }
404
405 //====================================================================
406 /// Helper function to set Euclid options using a command line
407 /// like array.
408 //=====================================================================
409 void euclid_settings_helper(const bool& use_block_jacobi,
410 const bool& use_row_scaling,
411 const bool& use_ilut,
412 const int& level,
413 const double& drop_tol,
414 const int& print_level,
415 HYPRE_Solver& euclid_object)
416 {
417 // Easier to use C-arrays rather than std::strings because Euclid takes
418 // char** as argument and because C++ doesn't provide decent number to
419 // std::string conversion functions.
420
421 int n_args = 0;
422 const char* args[22];
423
424 // first argument is empty string
425 args[n_args++] = "";
426
427 // switch on/off Block Jacobi ILU
428 args[n_args++] = "-bj";
429 if (use_block_jacobi)
430 {
431 args[n_args++] = "1";
432 }
433 else
434 {
435 args[n_args++] = "0";
436 }
437
438 // switch on/off row scaling
439 args[n_args++] = "-rowScale";
440 if (use_row_scaling)
441 {
442 args[n_args++] = "1";
443 }
444 else
445 {
446 args[n_args++] = "0";
447 }
448
449 // set level for ILU(k)
450 args[n_args++] = "-level";
451 char level_value[10];
452 sprintf(level_value, "%d", level);
453 args[n_args++] = level_value;
454
455 // // set drop tol for ILU(k) factorization
456 // args[n_args++] = "-sparseA";
457 // char droptol[20];
458 // sprintf(droptol,"%f",drop_tol);
459 // args[n_args++] = droptol;
460
461 // // set ILUT factorization if required
462 // if (use_ilut)
463 // {
464 // args[n_args++] = "-ilut";
465 // args[n_args++] = droptol;
466 // }
467
468 // set printing of Euclid data
469 if (print_level == 0)
470 {
471 args[n_args++] = "-eu_stats";
472 args[n_args++] = "0";
473 args[n_args++] = "-eu_mem";
474 args[n_args++] = "0";
475 }
476 if (print_level == 1)
477 {
478 args[n_args++] = "-eu_stats";
479 args[n_args++] = "1";
480 args[n_args++] = "-eu_mem";
481 args[n_args++] = "0";
482 }
483 if (print_level == 2)
484 {
485 args[n_args++] = "-eu_stats";
486 args[n_args++] = "1";
487 args[n_args++] = "-eu_mem";
488 args[n_args++] = "1";
489 }
490
491 // set next entry in array to null
492 args[n_args] = 0;
493
494 // Send the parameters
495 HYPRE_EuclidSetParams(euclid_object, n_args, const_cast<char**>(args));
496 }
497
498 } // namespace HypreHelpers
499
500
501 /// ////////////////////////////////////////////////////////////////////
502 /// ////////////////////////////////////////////////////////////////////
503 // functions for HypreInterface class
504 /// ////////////////////////////////////////////////////////////////////
505 /// ////////////////////////////////////////////////////////////////////
506
507
508 //=============================================================================
509 /// Helper function which creates a Hypre matrix from a CRDoubleMatrix
510 /// If OOMPH-LIB has been set up for MPI use, the Hypre matrix is
511 /// distributed over the available processors.
512 //=============================================================================
514 {
515 // reset Hypre's global error flag
517
518 // issue warning if the matrix is small compared to the number of processors
519 if (unsigned(2 * matrix_pt->distribution_pt()->communicator_pt()->nproc()) >
520 matrix_pt->nrow())
521 {
523 << "Warning: HYPRE based solvers may fail if 2*number of processors "
524 << "is greater than the number of unknowns!" << std::endl;
525 }
526
527 // store the distribution
528 // generate the Hypre matrix
531
532 // Output error messages if required
534 {
535 std::ostringstream message;
536 int err = HypreHelpers::check_HYPRE_error_flag(message);
537 if (err)
538 {
539 OomphLibWarning(message.str(),
540 "HypreSolver::hypre_matrix_setup()",
541 OOMPH_EXCEPTION_LOCATION);
542 }
543 }
544
545 // delete CRDoubleMatrix if required
547 {
548 matrix_pt->clear();
549 }
550 }
551
552
553 //=============================================================================
554 /// Sets up the solver data required for use in an oomph-lib
555 /// LinearSolver or Preconditioner, once the Hypre matrix has been
556 /// generated using hypre_matrix_setup(...).
557 //=============================================================================
559 {
560 // Store time
561 double t_start = TimingHelpers::timer();
562 double t_end = 0;
563
564
565 // reset Hypre's global error flag
567
568 // create dummy Hypre vectors which are required for setup
569 HYPRE_IJVector dummy_sol_ij;
570 HYPRE_ParVector dummy_sol_par;
571 HYPRE_IJVector dummy_rhs_ij;
572 HYPRE_ParVector dummy_rhs_par;
574 Hypre_distribution_pt, dummy_sol_ij, dummy_sol_par);
576 Hypre_distribution_pt, dummy_rhs_ij, dummy_rhs_par);
577
578 // Set up internal preconditioner for CG, GMRES or BiCGSTAB
579 // --------------------------------------------------------
580 if ((Hypre_method >= CG) && (Hypre_method <= BiCGStab))
581 {
582 // AMG preconditioner
584 {
585 // set up BoomerAMG
586 HYPRE_BoomerAMGCreate(&Preconditioner);
587 HYPRE_BoomerAMGSetPrintLevel(Preconditioner, AMG_print_level);
588 HYPRE_BoomerAMGSetMaxLevels(Preconditioner, AMG_max_levels);
589 HYPRE_BoomerAMGSetMaxIter(Preconditioner, 1);
590 HYPRE_BoomerAMGSetTol(Preconditioner, 0.0);
591 HYPRE_BoomerAMGSetCoarsenType(Preconditioner, AMG_coarsening);
592 HYPRE_BoomerAMGSetStrongThreshold(Preconditioner, AMG_strength);
593 HYPRE_BoomerAMGSetMaxRowSum(Preconditioner, AMG_max_row_sum);
594 HYPRE_BoomerAMGSetTruncFactor(Preconditioner, AMG_truncation);
595
597 {
598 HYPRE_BoomerAMGSetRelaxType(Preconditioner, AMG_simple_smoother);
599 HYPRE_BoomerAMGSetNumSweeps(Preconditioner, AMG_smoother_iterations);
600
601 // This one gives a memory leak
602 // double * relaxweight = new double[AMG_max_levels];
603
604 // This is how they do it in a hypre demo code
605 double* relaxweight = hypre_CTAlloc(double, AMG_max_levels);
606
607 for (unsigned i = 0; i < AMG_max_levels; i++)
608 {
609 relaxweight[i] = AMG_damping;
610 }
611 HYPRE_BoomerAMGSetRelaxWeight(Preconditioner, relaxweight);
612 }
613 else
614 {
615 HYPRE_BoomerAMGSetSmoothType(Preconditioner, AMG_complex_smoother);
616 HYPRE_BoomerAMGSetSmoothNumLevels(Preconditioner, AMG_max_levels);
617 HYPRE_BoomerAMGSetSmoothNumSweeps(Preconditioner,
619
620 // If we are using Euclid then set up additional Euclid only options
621 if (AMG_complex_smoother == 9)
622 {
631 }
632 }
633
635 }
636
637 // Euclid preconditioner
639 {
640#ifdef OOMPH_HAS_MPI
641 HYPRE_EuclidCreate(Hypre_distribution_pt->communicator_pt()->mpi_comm(),
643#else
644 HYPRE_EuclidCreate(MPI_COMM_WORLD, &Preconditioner);
645#endif
646
647 // Set parameters
655
657 }
658
659 // ParaSails preconditioner
661 {
662#ifdef OOMPH_HAS_MPI
663 HYPRE_ParaSailsCreate(
666#else
667 HYPRE_ParaSailsCreate(MPI_COMM_WORLD, &Preconditioner);
668#endif
669 HYPRE_ParaSailsSetSym(Preconditioner, ParaSails_symmetry);
670 HYPRE_ParaSailsSetParams(
672 HYPRE_ParaSailsSetFilter(Preconditioner, ParaSails_filter);
674 }
675
676 // check error flag
678 {
679 std::ostringstream message;
680 int err = HypreHelpers::check_HYPRE_error_flag(message);
681 if (err)
682 {
683 OomphLibWarning(message.str(),
684 "HypreSolver::hypre_setup()",
685 OOMPH_EXCEPTION_LOCATION);
686 }
687 }
688 } // end of setting up internal preconditioner
689
690
691 // set up solver
692 // -------------
693 t_start = TimingHelpers::timer();
694
695 // AMG solver
696 if (Hypre_method == BoomerAMG)
697 {
698 if (Output_info)
699 {
700 oomph_info << "Setting up BoomerAMG, ";
701 }
702
703 // set up BoomerAMG
704 HYPRE_BoomerAMGCreate(&Solver);
705 HYPRE_BoomerAMGSetPrintLevel(Solver, AMG_print_level);
706 HYPRE_BoomerAMGSetMaxLevels(Solver, AMG_max_levels);
707 HYPRE_BoomerAMGSetMaxIter(Solver, Max_iter);
708 HYPRE_BoomerAMGSetTol(Solver, Tolerance);
709 HYPRE_BoomerAMGSetCoarsenType(Solver, AMG_coarsening);
710 HYPRE_BoomerAMGSetStrongThreshold(Solver, AMG_strength);
711 HYPRE_BoomerAMGSetMaxRowSum(Solver, AMG_max_row_sum);
712 HYPRE_BoomerAMGSetTruncFactor(Solver, AMG_truncation);
713
715 {
716 HYPRE_BoomerAMGSetRelaxType(Solver, AMG_simple_smoother);
717 HYPRE_BoomerAMGSetNumSweeps(Solver, AMG_smoother_iterations);
718
719 // This one gives a memory leak
720 // double * relaxweight = new double[AMG_max_levels];
721
722 // This is how they do it in a hypre demo code
723 double* relaxweight = hypre_CTAlloc(double, AMG_max_levels);
724
725 for (unsigned i = 0; i < AMG_max_levels; i++)
726 {
727 relaxweight[i] = AMG_damping;
728 }
729 HYPRE_BoomerAMGSetRelaxWeight(Solver, relaxweight);
730 }
731 else
732 {
733 HYPRE_BoomerAMGSetSmoothType(Solver, AMG_complex_smoother);
734 HYPRE_BoomerAMGSetSmoothNumLevels(Solver, AMG_max_levels);
735 HYPRE_BoomerAMGSetSmoothNumSweeps(Solver, AMG_smoother_iterations);
736
737 /* Other settings
738 * 6 & Schwarz smoothers & HYPRE_BoomerAMGSetDomainType,
739 * HYPRE_BoomerAMGSetOverlap, \\
740 * & & HYPRE_BoomerAMGSetVariant, HYPRE_BoomerAMGSetSchwarzRlxWeight
741 * \\
742 * 7 & Pilut & HYPRE_BoomerAMGSetDropTol, HYPRE_BoomerAMGSetMaxNzPerRow
743 * \\
744 * 8 & ParaSails & HYPRE_BoomerAMGSetSym, HYPRE_BoomerAMGSetLevel, \\
745 * & & HYPRE_BoomerAMGSetFilter, HYPRE_BoomerAMGSetThreshold \\
746 * 9 & Euclid & HYPRE_BoomerAMGSetEuclidFile \\
747 */
748
749 // If we are using Euclid then set up additional Euclid only options
750 if (AMG_complex_smoother == 9)
751 {
760 }
761
762 // Add any others here as required...
763 }
764
765 // MemoryUsage::doc_memory_usage("before amg setup [solver]");
766 // MemoryUsage::insert_comment_to_continous_top("BEFORE AMG SETUP
767 // [SOLVER]");
768
769 HYPRE_BoomerAMGSetup(Solver, Matrix_par, dummy_rhs_par, dummy_sol_par);
770
771 // MemoryUsage::doc_memory_usage("after amg setup [solver]");
772 // MemoryUsage::insert_comment_to_continous_top("AFTER AMG SETUP
773 // [SOLVER]");
774
776 }
777
778 // Euclid solver
779 else if (Hypre_method == Euclid)
780 {
781 if (Output_info)
782 {
783 oomph_info << "Setting up Euclid, ";
784 }
785#ifdef OOMPH_HAS_MPI
786 HYPRE_EuclidCreate(Hypre_distribution_pt->communicator_pt()->mpi_comm(),
787 &Solver);
788#else
789 HYPRE_EuclidCreate(MPI_COMM_WORLD, &Solver);
790#endif
791
792 // Set parameters
800
801 HYPRE_EuclidSetup(Solver, Matrix_par, dummy_rhs_par, dummy_sol_par);
803 }
804
805 // ParaSails preconditioner
806 else if (Hypre_method == ParaSails)
807 {
808 if (Output_info)
809 {
810 oomph_info << "Setting up ParaSails, ";
811 }
812#ifdef OOMPH_HAS_MPI
813 HYPRE_ParaSailsCreate(
815#else
816 HYPRE_ParaSailsCreate(MPI_COMM_WORLD, &Solver);
817#endif
818 HYPRE_ParaSailsSetSym(Solver, ParaSails_symmetry);
819 HYPRE_ParaSailsSetParams(Solver, ParaSails_thresh, ParaSails_nlevel);
820 HYPRE_ParaSailsSetFilter(Solver, ParaSails_filter);
821
822 HYPRE_ParaSailsSetup(Solver, Matrix_par, dummy_rhs_par, dummy_sol_par);
824 }
825
826 // CG solver
827 else if (Hypre_method == CG)
828 {
829 if (Output_info)
830 {
831 oomph_info << "Setting up CG, ";
832 }
833
834#ifdef OOMPH_HAS_MPI
835 HYPRE_ParCSRPCGCreate(
837#else
838 HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &Solver);
839#endif
840 HYPRE_PCGSetTol(Solver, Tolerance);
841 HYPRE_PCGSetLogging(Solver, 0);
842 HYPRE_PCGSetPrintLevel(Solver, Krylov_print_level);
843 HYPRE_PCGSetMaxIter(Solver, Max_iter);
844
845 // set preconditioner
847 {
848 if (Output_info)
849 {
850 oomph_info << " with BoomerAMG preconditioner, ";
851 }
852
853 HYPRE_PCGSetPrecond(Solver,
854 (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSolve,
855 (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSetup,
857 }
858 else if (Internal_preconditioner == Euclid) // Euclid
859 {
860 if (Output_info)
861 {
862 oomph_info << " with Euclid ILU preconditioner, ";
863 }
864
865 HYPRE_PCGSetPrecond(Solver,
866 (HYPRE_PtrToSolverFcn)HYPRE_EuclidSolve,
867 (HYPRE_PtrToSolverFcn)HYPRE_EuclidSetup,
869 }
870 else if (Internal_preconditioner == ParaSails) // ParaSails
871 {
872 if (Output_info)
873 {
874 oomph_info << " with ParaSails approximate inverse preconditioner, ";
875 }
876
877 HYPRE_PCGSetPrecond(Solver,
878 (HYPRE_PtrToSolverFcn)HYPRE_ParaSailsSolve,
879 (HYPRE_PtrToSolverFcn)HYPRE_ParaSailsSetup,
881 }
882 else
883 {
884 if (Output_info)
885 {
886 oomph_info << " with no preconditioner";
887 }
888 }
889
890
891 HYPRE_PCGSetup(Solver,
892 (HYPRE_Matrix)Matrix_par,
893 (HYPRE_Vector)dummy_rhs_par,
894 (HYPRE_Vector)dummy_sol_par);
895
896
898 }
899
900 // GMRES solver
901 else if (Hypre_method == GMRES)
902 {
903 if (Output_info)
904 {
905 oomph_info << "Setting up GMRES";
906 }
907
908#ifdef OOMPH_HAS_MPI
909 HYPRE_ParCSRGMRESCreate(
911#else
912 HYPRE_ParCSRGMRESCreate(MPI_COMM_WORLD, &Solver);
913#endif
914 HYPRE_GMRESSetTol(Solver, Tolerance);
915 HYPRE_GMRESSetKDim(Solver, Max_iter);
916 HYPRE_GMRESSetLogging(Solver, 0);
917 HYPRE_GMRESSetPrintLevel(Solver, Krylov_print_level);
918 HYPRE_GMRESSetMaxIter(Solver, Max_iter);
919
920 // set preconditioner
922 {
923 if (Output_info)
924 {
925 oomph_info << " with BoomerAMG preconditioner, ";
926 }
927
928 HYPRE_GMRESSetPrecond(Solver,
929 (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSolve,
930 (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSetup,
932 }
933 else if (Internal_preconditioner == Euclid) // Euclid
934 {
935 if (Output_info)
936 {
937 oomph_info << " with Euclid ILU preconditioner, ";
938 }
939
940 HYPRE_GMRESSetPrecond(Solver,
941 (HYPRE_PtrToSolverFcn)HYPRE_EuclidSolve,
942 (HYPRE_PtrToSolverFcn)HYPRE_EuclidSetup,
944 }
945 else if (Internal_preconditioner == ParaSails) // ParaSails
946 {
947 if (Output_info)
948 {
949 oomph_info << " with ParaSails approximate inverse preconditioner, ";
950 }
951
952 HYPRE_GMRESSetPrecond(Solver,
953 (HYPRE_PtrToSolverFcn)HYPRE_ParaSailsSolve,
954 (HYPRE_PtrToSolverFcn)HYPRE_ParaSailsSetup,
956 }
957 else
958 {
959 if (Output_info)
960 {
961 oomph_info << " with no preconditioner";
962 }
963 }
964
965 HYPRE_GMRESSetup(Solver,
966 (HYPRE_Matrix)Matrix_par,
967 (HYPRE_Vector)dummy_rhs_par,
968 (HYPRE_Vector)dummy_sol_par);
969
971 }
972
973 // BiCGStab solver
974 else if (Hypre_method == BiCGStab)
975 {
976 if (Output_info)
977 {
978 oomph_info << "Setting up BiCGStab";
979 }
980#ifdef OOMPH_HAS_MPI
981 HYPRE_ParCSRBiCGSTABCreate(
983#else
984 HYPRE_ParCSRBiCGSTABCreate(MPI_COMM_WORLD, &Solver);
985#endif
986 HYPRE_BiCGSTABSetTol(Solver, Tolerance);
987 HYPRE_BiCGSTABSetLogging(Solver, 0);
988 HYPRE_BiCGSTABSetPrintLevel(Solver, Krylov_print_level);
989 HYPRE_BiCGSTABSetMaxIter(Solver, Max_iter);
990
991 // set preconditioner
993 {
994 if (Output_info)
995 {
996 oomph_info << " with BoomerAMG preconditioner, ";
997 }
998
999 HYPRE_BiCGSTABSetPrecond(Solver,
1000 (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSolve,
1001 (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSetup,
1003 }
1004 else if (Internal_preconditioner == Euclid) // Euclid
1005 {
1006 if (Output_info)
1007 {
1008 oomph_info << " with Euclid ILU preconditioner, ";
1009 }
1010
1011 HYPRE_BiCGSTABSetPrecond(Solver,
1012 (HYPRE_PtrToSolverFcn)HYPRE_EuclidSolve,
1013 (HYPRE_PtrToSolverFcn)HYPRE_EuclidSetup,
1015 }
1016 else if (Internal_preconditioner == ParaSails) // ParaSails
1017 {
1018 if (Output_info)
1019 {
1020 oomph_info << " with ParaSails approximate inverse preconditioner, ";
1021 }
1022
1023 HYPRE_BiCGSTABSetPrecond(Solver,
1024 (HYPRE_PtrToSolverFcn)HYPRE_ParaSailsSolve,
1025 (HYPRE_PtrToSolverFcn)HYPRE_ParaSailsSetup,
1027 }
1028 else
1029 {
1030 if (Output_info)
1031 {
1032 oomph_info << " with no preconditioner, ";
1033 }
1034 }
1035
1036 HYPRE_BiCGSTABSetup(Solver,
1037 (HYPRE_Matrix)Matrix_par,
1038 (HYPRE_Vector)dummy_rhs_par,
1039 (HYPRE_Vector)dummy_sol_par);
1040
1042 }
1043
1044 // no solver exists for this value of Solver flag
1045 else
1046 {
1047 std::ostringstream error_message;
1048 error_message << "Solver has been set to an invalid value. "
1049 << "current value=" << Solver;
1050 throw OomphLibError(
1051 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1052 }
1053
1054 t_end = TimingHelpers::timer();
1055 double solver_setup_time = t_end - t_start;
1056
1057 // destroy dummy hypre vectors
1058 HYPRE_IJVectorDestroy(dummy_sol_ij);
1059 HYPRE_IJVectorDestroy(dummy_rhs_ij);
1060
1061 // check error flag
1063 {
1064 std::ostringstream message;
1065 int err = HypreHelpers::check_HYPRE_error_flag(message);
1066 if (err)
1067 {
1068 OomphLibWarning(message.str(),
1069 "HypreSolver::hypre_solver_setup()",
1070 OOMPH_EXCEPTION_LOCATION);
1071 }
1072 }
1073
1074 // output times
1075 if (Output_info)
1076 {
1077 oomph_info << "time for setup [s] : " << solver_setup_time << std::endl;
1078 }
1079 }
1080
1081
1082 //===================================================================
1083 /// Helper function performs a solve if solver data has been set
1084 /// up using hypre_solver_setup(...).
1085 //====================================================================
1087 DoubleVector& solution)
1088 {
1089 // Record time
1090 double t_start = TimingHelpers::timer();
1091
1092 // Set up hypre vectors
1093 // --------------------
1094
1095 // Hypre vector for rhs values
1096 HYPRE_IJVector rhs_ij;
1097 HYPRE_ParVector rhs_par;
1098
1099 // Hypre vector for solution
1100 HYPRE_IJVector solution_ij;
1101 HYPRE_ParVector solution_par;
1102
1103 // set up rhs_values and vec_indices
1105 rhs, Hypre_distribution_pt, rhs_ij, rhs_par);
1106
1108 Hypre_distribution_pt, solution_ij, solution_par);
1109
1110 // check error flag
1112 {
1113 std::ostringstream message;
1114 int err = HypreHelpers::check_HYPRE_error_flag(message);
1115 if (err)
1116 {
1117 OomphLibWarning(message.str(),
1118 "HypreSolver::hypre_solve()",
1119 OOMPH_EXCEPTION_LOCATION);
1120 }
1121 }
1122
1123 // solve
1124 // -----
1125
1126 // for solver stats
1127 int iterations = 0;
1128 double norm = 0;
1129
1130 // Get the norm of rhs
1131 const double rhs_norm = rhs.norm();
1132 bool do_solving = false;
1133 if (rhs_norm > 0.0)
1134 {
1135 do_solving = true;
1136 }
1137
1138#ifdef OOMPH_HAS_MPI
1139 // We need to check whether any processor requires to solve, if that
1140 // is the case then do the solving
1142 {
1143 if (MPI_Helpers::communicator_pt()->nproc() > 1)
1144 {
1145 unsigned this_processor_do_solving = 0;
1146 unsigned all_processors_do_solving = 0;
1147 if (do_solving)
1148 {
1149 this_processor_do_solving = 1;
1150 }
1151 // Get the communicator
1153 // Communicate with all procesoors
1154 MPI_Allreduce(&this_processor_do_solving,
1155 &all_processors_do_solving,
1156 1,
1157 MPI_UNSIGNED,
1158 MPI_SUM,
1159 comm_pt->mpi_comm());
1160 if (all_processors_do_solving > 0)
1161 {
1162 do_solving = true;
1163 }
1164 }
1165 }
1166#endif
1167
1168 if (do_solving)
1169 {
1171 {
1172 HYPRE_BoomerAMGSolve(Solver, Matrix_par, rhs_par, solution_par);
1173 HYPRE_BoomerAMGGetNumIterations(Solver, &iterations);
1174 HYPRE_BoomerAMGGetFinalRelativeResidualNorm(Solver, &norm);
1175 }
1176 else if (Existing_solver == CG)
1177 {
1178 HYPRE_PCGSolve(Solver,
1179 (HYPRE_Matrix)Matrix_par,
1180 (HYPRE_Vector)rhs_par,
1181 (HYPRE_Vector)solution_par);
1182 HYPRE_PCGGetNumIterations(Solver, &iterations);
1183 HYPRE_PCGGetFinalRelativeResidualNorm(Solver, &norm);
1184 }
1185 else if (Existing_solver == GMRES)
1186 {
1187 HYPRE_GMRESSolve(Solver,
1188 (HYPRE_Matrix)Matrix_par,
1189 (HYPRE_Vector)rhs_par,
1190 (HYPRE_Vector)solution_par);
1191 HYPRE_GMRESGetNumIterations(Solver, &iterations);
1192 HYPRE_GMRESGetFinalRelativeResidualNorm(Solver, &norm);
1193 }
1194 else if (Existing_solver == BiCGStab)
1195 {
1196 HYPRE_BiCGSTABSolve(Solver,
1197 (HYPRE_Matrix)Matrix_par,
1198 (HYPRE_Vector)rhs_par,
1199 (HYPRE_Vector)solution_par);
1200 HYPRE_BiCGSTABGetNumIterations(Solver, &iterations);
1201 HYPRE_BiCGSTABGetFinalRelativeResidualNorm(Solver, &norm);
1202 }
1203 else if (Existing_solver == Euclid)
1204 {
1205 HYPRE_EuclidSolve(Solver, Matrix_par, rhs_par, solution_par);
1206 }
1207 else if (Existing_solver == ParaSails)
1208 {
1209 HYPRE_ParaSailsSolve(Solver, Matrix_par, rhs_par, solution_par);
1210 }
1211
1212 // output any error message
1214 {
1215 std::ostringstream message;
1216 int err = HypreHelpers::check_HYPRE_error_flag(message);
1217 if (err)
1218 {
1219 OomphLibWarning(message.str(),
1220 "HypreSolver::hypre_solve()",
1221 OOMPH_EXCEPTION_LOCATION);
1222 }
1223 }
1224
1225 } // if (do_solving)
1226
1227 // Copy result to solution
1228 unsigned nrow_local = Hypre_distribution_pt->nrow_local();
1229 unsigned first_row = Hypre_distribution_pt->first_row();
1230 int* indices = new int[nrow_local];
1231 for (unsigned i = 0; i < nrow_local; i++)
1232 {
1233 indices[i] = first_row + i;
1234 }
1235 LinearAlgebraDistribution* soln_dist_pt;
1236 if (solution.built())
1237 {
1238 soln_dist_pt = new LinearAlgebraDistribution(solution.distribution_pt());
1239 }
1240 else
1241 {
1242 soln_dist_pt = new LinearAlgebraDistribution(rhs.distribution_pt());
1243 }
1244 solution.build(Hypre_distribution_pt, 0.0);
1245 HYPRE_IJVectorGetValues(
1246 solution_ij, nrow_local, indices, solution.values_pt());
1247 solution.redistribute(soln_dist_pt);
1248 delete[] indices;
1249 delete soln_dist_pt;
1250
1251 // output any error message
1253 {
1254 std::ostringstream message;
1255 int err = HypreHelpers::check_HYPRE_error_flag(message);
1256 if (err)
1257 {
1258 OomphLibWarning(message.str(),
1259 "HypreSolver::hypre_solve()",
1260 OOMPH_EXCEPTION_LOCATION);
1261 }
1262 }
1263
1264 // deallocation
1265 HYPRE_IJVectorDestroy(solution_ij);
1266 HYPRE_IJVectorDestroy(rhs_ij);
1267
1268 // Record time
1269 double solve_time = 0;
1270 if (Output_info)
1271 {
1272 double t_end = TimingHelpers::timer();
1273 solve_time = t_end - t_start;
1274 }
1275
1276 // output timings and info
1277 if (Output_info)
1278 {
1279 oomph_info << "Time for HYPRE solve [s] : " << solve_time << std::endl;
1280 }
1281
1282 // for iterative solvers output iterations and final norm
1283 if ((Hypre_method >= CG) && (Hypre_method <= BoomerAMG))
1284 {
1285 if (iterations > 1)
1286 {
1287 if (Output_info)
1288 oomph_info << "Number of iterations : " << iterations
1289 << std::endl;
1290 if (Output_info)
1291 oomph_info << "Final Relative Residual Norm : " << norm << std::endl;
1292 }
1293 }
1294 }
1295
1296
1297 //===================================================================
1298 /// hypre_clean_up_memory() deletes any existing Hypre solver and
1299 /// Hypre matrix
1300 //====================================================================
1302 {
1303 // is there an existing solver
1304 if (Existing_solver != None)
1305 {
1306 // delete matrix
1307 HYPRE_IJMatrixDestroy(Matrix_ij);
1308
1309 // delete solver
1311 {
1312 HYPRE_BoomerAMGDestroy(Solver);
1313 }
1314 else if (Existing_solver == CG)
1315 {
1316 HYPRE_ParCSRPCGDestroy(Solver);
1317 }
1318 else if (Existing_solver == GMRES)
1319 {
1320 HYPRE_ParCSRGMRESDestroy(Solver);
1321 }
1322 else if (Existing_solver == BiCGStab)
1323 {
1324 HYPRE_ParCSRBiCGSTABDestroy(Solver);
1325 }
1326 else if (Existing_solver == Euclid)
1327 {
1328 HYPRE_EuclidDestroy(Solver);
1329 }
1330 else if (Existing_solver == ParaSails)
1331 {
1332 HYPRE_ParaSailsDestroy(Solver);
1333 }
1335
1336 // delete preconditioner
1338 {
1339 HYPRE_BoomerAMGDestroy(Preconditioner);
1340 }
1341 else if (Existing_preconditioner == Euclid)
1342 {
1343 HYPRE_EuclidDestroy(Preconditioner);
1344 }
1346 {
1347 HYPRE_ParaSailsDestroy(Preconditioner);
1348 }
1350
1351 // check error flag
1353 {
1354 std::ostringstream message;
1355 int err = HypreHelpers::check_HYPRE_error_flag(message);
1356 if (err)
1357 {
1358 OomphLibWarning(message.str(),
1359 "HypreSolver::clean_up_memory()",
1360 OOMPH_EXCEPTION_LOCATION);
1361 }
1362 }
1363 }
1364 }
1365
1366
1367 /// ////////////////////////////////////////////////////////////////////
1368 /// ////////////////////////////////////////////////////////////////////
1369 // functions for HypreSolver class
1370 /// ////////////////////////////////////////////////////////////////////
1371 /// ////////////////////////////////////////////////////////////////////
1372
1373
1374 //===================================================================
1375 /// Problem-based solve function to generate the Jacobian matrix and
1376 /// residual vector and use HypreInterface::hypre_solver_setup(...)
1377 /// and HypreInterface::hypre_solve(...) to solve the linear system.
1378 /// This function will delete any existing data.
1379 /// Note: the returned solution vector is NOT distributed, i.e. all
1380 /// processors hold all values because this is what the Newton solver
1381 /// requires.
1382 //====================================================================
1383 void HypreSolver::solve(Problem* const& problem_pt, DoubleVector& solution)
1384 {
1385 double t_start = TimingHelpers::timer();
1386
1387 // Set Output_time flag for HypreInterface
1389
1390 // Delete any existing solver data
1392
1393 // Set flag to allow deletion of the oomphlib Jacobian matrix
1394 // (we're in control)
1395 Delete_input_data = true;
1396
1397 // Get oomph-lib Jacobian matrix and residual vector
1398 DoubleVector residual;
1399 CRDoubleMatrix* matrix = new CRDoubleMatrix;
1400 problem_pt->get_jacobian(residual, *matrix);
1401
1402 // Output times
1403 if (Doc_time)
1404 {
1405 oomph_info << "Time to generate Jacobian and residual [s] : ";
1406 double t_end = TimingHelpers::timer();
1407 oomph_info << t_end - t_start << std::endl;
1408 }
1409
1410 // generate hypre matrix
1411 hypre_matrix_setup(matrix);
1412
1413 // call hypre_solver_setup to generate linear solver data
1415
1416 // perform hypre_solve
1417 hypre_solve(residual, solution);
1418
1419 // delete solver data if required
1420 if (!Enable_resolve)
1421 {
1423 }
1424 }
1425
1426
1427 //===================================================================
1428 /// Uses HypreInterface::hypre_solve(...) to solve the linear system
1429 /// for a CRDoubleMatrix or a DistributedCRDoubleMatrix.
1430 /// In the latter case, the rhs needs to be distributed too,
1431 /// i.e. the length of the rhs vector must be equal to the number of
1432 /// rows stored locally.
1433 /// Note: the returned solution vector is never distributed, i.e. all
1434 /// processors hold all values
1435 //====================================================================
1436 void HypreSolver::solve(DoubleMatrixBase* const& matrix_pt,
1437 const DoubleVector& rhs,
1438 DoubleVector& solution)
1439 {
1440#ifdef PARANOID
1441 // check the matrix is square
1442 if (matrix_pt->nrow() != matrix_pt->ncol())
1443 {
1444 std::ostringstream error_message;
1445 error_message << "HypreSolver require a square matrix. "
1446 << "Matrix is " << matrix_pt->nrow() << " by "
1447 << matrix_pt->ncol() << std::endl;
1448 throw OomphLibError(
1449 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1450 }
1451#endif
1452
1453 // Set Output_time flag for HypreInterface
1455
1456 // Clean up existing solver data
1458
1459 // Set flag to decide if oomphlib matrix can be deleted
1460 // (Recall that Delete_matrix defaults to false).
1462
1463 // Try cast to a CRDoubleMatrix
1464 CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1465
1466 // If cast successful set things up for a serial solve
1467 if (cr_matrix_pt)
1468 {
1469 // rebuild the distribution
1470 this->build_distribution(cr_matrix_pt->distribution_pt());
1471
1472#ifdef PARANOID
1473 // check that rhs has the same distribution as the matrix (now stored
1474 // in Distribution_pt)
1475 if (*this->distribution_pt() != *rhs.distribution_pt())
1476 {
1477 std::ostringstream error_message;
1478 error_message << "The distribution of the rhs vector and the matrix "
1479 << " must be the same" << std::endl;
1480 throw OomphLibError(error_message.str(),
1481 OOMPH_CURRENT_FUNCTION,
1482 OOMPH_EXCEPTION_LOCATION);
1483 }
1484 // if the solution is setup make sure it has the same distribution as
1485 // the matrix as well
1486 if (solution.built())
1487 {
1488 if (*this->distribution_pt() != *solution.distribution_pt())
1489 {
1490 std::ostringstream error_message;
1491 error_message << "The distribution of the solution vector is setup "
1492 << "there it must be the same as the matrix."
1493 << std::endl;
1494 throw OomphLibError(error_message.str(),
1495 OOMPH_CURRENT_FUNCTION,
1496 OOMPH_EXCEPTION_LOCATION);
1497 }
1498 }
1499#endif
1500
1501 hypre_matrix_setup(cr_matrix_pt);
1502 }
1503 else
1504 {
1505#ifdef PARANOID
1506 std::ostringstream error_message;
1507 error_message << "HypreSolver only work with "
1508 << "CRDoubleMatrix matrices" << std::endl;
1509 throw OomphLibError(
1510 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1511#endif
1512 }
1513
1514 // call hypre_setup to generate Hypre matrix and linear solver data
1516
1517 // perform hypre_solve
1518 hypre_solve(rhs, solution);
1519
1520 // delete solver data if required
1521 if (!Enable_resolve)
1522 {
1524 }
1525 }
1526
1527
1528 //===================================================================
1529 /// Resolve performs a linear solve using current solver data (if
1530 /// such data exists).
1531 //====================================================================
1533 {
1534#ifdef PARANOID
1535 // check solver data exists
1536 if (existing_solver() == None)
1537 {
1538 std::ostringstream error_message;
1539 error_message << "resolve(...) requires that solver data has been "
1540 << "set up by a previous call to solve(...) after "
1541 << "a call to enable_resolve()" << std::endl;
1542 throw OomphLibError(
1543 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1544 }
1545 // check that rhs has the same distribution as the matrix (now stored
1546 // in Distribution_pt)
1547 if (*this->distribution_pt() != *rhs.distribution_pt())
1548 {
1549 std::ostringstream error_message;
1550 error_message << "The distribution of the rhs vector and the matrix "
1551 << " must be the same" << std::endl;
1552 throw OomphLibError(
1553 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1554 }
1555 // if the solution is setup make sure it has the same distribution as
1556 // the matrix as well
1557 if (solution.built())
1558 {
1559 if (*this->distribution_pt() != *solution.distribution_pt())
1560 {
1561 std::ostringstream error_message;
1562 error_message << "The distribution of the solution vector is setup "
1563 << "there it must be the same as the matrix."
1564 << std::endl;
1565 throw OomphLibError(error_message.str(),
1566 OOMPH_CURRENT_FUNCTION,
1567 OOMPH_EXCEPTION_LOCATION);
1568 }
1569 }
1570#endif
1571
1572 // Set Output_info flag for HypreInterface
1574
1575 // solve
1576 hypre_solve(rhs, solution);
1577
1578 // Note: never delete solver data as the preconditioner is typically
1579 // called repeatedly.
1580 }
1581
1582
1583 //===================================================================
1584 /// clean_up_memory() deletes any existing solver data
1585 //====================================================================
1587 {
1589 }
1590
1591
1592 /// ////////////////////////////////////////////////////////////////////
1593 /// ////////////////////////////////////////////////////////////////////
1594 // functions for HyprePreconditioner class
1595 /// ////////////////////////////////////////////////////////////////////
1596 /// ////////////////////////////////////////////////////////////////////
1597
1598
1599 //=============================================================================
1600 /// Static double that accumulates the preconditioner
1601 /// solve time of all instantiations of this class. Reset
1602 /// it manually, e.g. after every Newton solve, using
1603 /// reset_cumulative_solve_times().
1604 //=============================================================================
1606
1607 //=============================================================================
1608 /// map of static doubles that accumulates the preconditioner
1609 /// solve time of all instantiations of this class, labeled by
1610 /// context string. Reset
1611 /// it manually, e.g. after every Newton solve, using
1612 /// reset_cumulative_solve_times().
1613 //=============================================================================
1614 std::map<std::string, double>
1616
1617 //=============================================================================
1618 /// Static unsigned that accumulates the number of preconditioner
1619 /// solves of all instantiations of this class. Reset
1620 /// it manually, e.g. after every Newton solve, using
1621 /// reset_cumulative_solve_times().
1622 //=============================================================================
1624
1625 //=============================================================================
1626 /// Static unsigned that accumulates the number of preconditioner
1627 /// solves of all instantiations of this class, labeled by
1628 /// context string. Reset
1629 /// it manually, e.g. after every Newton solve, using
1630 /// reset_cumulative_solve_times().
1631 //=============================================================================
1632 std::map<std::string, unsigned>
1634
1635 //=============================================================================
1636 /// Static unsigned that stores nrow for the most recent
1637 /// instantiations of this class, labeled by
1638 /// context string. Reset
1639 /// it manually, e.g. after every Newton solve, using
1640 /// reset_cumulative_solve_times().
1641 //=============================================================================
1642 std::map<std::string, unsigned> HyprePreconditioner::Context_based_nrow;
1643
1644 //=============================================================================
1645 /// Report cumulative solve times of all instantiations of this
1646 /// class
1647 //=============================================================================
1649 {
1650 oomph_info << "\n\n=====================================================\n";
1651 oomph_info << "Cumulative HyprePreconditioner solve time "
1653 << " for " << Cumulative_npreconditioner_solve << " solves";
1655 {
1656 oomph_info << " ( "
1659 << " per solve )";
1660 }
1661 oomph_info << std::endl << std::endl;
1663 {
1664 oomph_info << "Breakdown by context: " << std::endl;
1665 for (std::map<std::string, double>::iterator it =
1668 it++)
1669 {
1671 << (*it).first << " " << (*it).second << " for "
1673 << " solves";
1675 {
1677 << " ( "
1678 << (*it).second /
1679 double(
1681 << " per solve; "
1682 << (*it).second /
1683 double(
1685 .first]) /
1686 double(Context_based_nrow[(*it).first])
1687 << " per solve per dof )";
1688 }
1689 oomph_info << std::endl;
1690 }
1691 }
1692 oomph_info << "\n=====================================================\n";
1693 oomph_info << std::endl;
1694 }
1695
1696 //=============================================================================
1697 /// Reset cumulative solve times
1698 //=============================================================================
1700 {
1705 Context_based_nrow.clear();
1706 }
1707
1708
1709 //=============================================================================
1710 /// An interface to allow HypreSolver to be used as a Preconditioner
1711 /// for the oomph-lib IterativeLinearSolver class.
1712 /// Matrix has to be of type CRDoubleMatrix or DistributedCRDoubleMatrix.
1713 //=============================================================================
1715 {
1716 // Set Output_info flag for HypreInterface
1718
1719#ifdef PARANOID
1720 // check the matrix is square
1721 if (matrix_pt()->nrow() != matrix_pt()->ncol())
1722 {
1723 std::ostringstream error_message;
1724 error_message << "HyprePreconditioner require a square matrix. "
1725 << "Matrix is " << matrix_pt()->nrow() << " by "
1726 << matrix_pt()->ncol() << std::endl;
1727 throw OomphLibError(
1728 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1729 }
1730#endif
1731
1732 // clean up any previous solver data
1734
1735 // set flag to decide if oomphlib matrix can be deleted
1736 // (Recall that Delete_matrix defaults to false).
1738
1739 // Try casting to a CRDoubleMatrix
1740 CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
1741
1742 // If cast successful set things up for a serial solve
1743 if (cr_matrix_pt)
1744 {
1745 this->build_distribution(cr_matrix_pt->distribution_pt());
1746 hypre_matrix_setup(cr_matrix_pt);
1747 }
1748 else
1749 {
1750#ifdef PARANOID
1751 std::ostringstream error_message;
1752 error_message << "HyprePreconditioner only work with "
1753 << "CRDoubleMatrix matrices" << std::endl;
1754 throw OomphLibError(
1755 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1756#endif
1757 }
1758
1759 if (Context_string != "")
1760 {
1761 oomph_info << "Setup of HyprePreconditioner in context \" "
1762 << Context_string << "\": nrow, nrow_local, nnz "
1763 << cr_matrix_pt->nrow() << " " << cr_matrix_pt->nrow_local()
1764 << " " << cr_matrix_pt->nnz() << std::endl;
1765 }
1766 Context_based_nrow[Context_string] = cr_matrix_pt->nrow();
1767
1768 // call hypre_solver_setup
1770 }
1771
1772 //===================================================================
1773 /// Preconditioner_solve uses a hypre solver to precondition vector r
1774 //====================================================================
1776 DoubleVector& z)
1777 {
1778 // Store time
1779 double t_start = TimingHelpers::timer();
1780
1781#ifdef PARANOID
1782 // check solver data exists
1783 if (existing_solver() == None)
1784 {
1785 std::ostringstream error_message;
1786 error_message << "preconditioner_solve(...) requires that data has "
1787 << "been set up using the function setup(...)" << std::endl;
1788 throw OomphLibError(
1789 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1790 }
1791 // check that rhs has the same distribution as the matrix (now stored
1792 // in Distribution_pt)
1793 if (*this->distribution_pt() != *r.distribution_pt())
1794 {
1795 std::ostringstream error_message;
1796 error_message << "The distribution of the rhs vector and the matrix "
1797 << " must be the same" << std::endl;
1798 throw OomphLibError(
1799 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1800 }
1801
1802 // if the solution is setup make sure it has the same distribution as
1803 // the matrix as well
1804 if (z.built())
1805 {
1806 if (*this->distribution_pt() != *z.distribution_pt())
1807 {
1808 std::ostringstream error_message;
1809 error_message << "The distribution of the solution vector is setup "
1810 << "there it must be the same as the matrix."
1811 << std::endl;
1812 throw OomphLibError(error_message.str(),
1813 OOMPH_CURRENT_FUNCTION,
1814 OOMPH_EXCEPTION_LOCATION);
1815 }
1816 }
1817#endif
1818
1819 // Switch off any timings for the solve
1820 Output_info = false;
1821
1822 // perform hypre_solve
1823 hypre_solve(r, z);
1824
1825 // Add to cumulative solve time
1826 double t_end = TimingHelpers::timer();
1827 Cumulative_preconditioner_solve_time += (t_end - t_start);
1829 My_cumulative_preconditioner_solve_time += (t_end - t_start);
1830 if (Context_string != "")
1831 {
1834 }
1835 }
1836
1837
1838 //===================================================================
1839 /// clean_up_memory() deletes any existing Hypre solver and
1840 /// Hypre matrix
1841 //====================================================================
1843 {
1845 }
1846
1847} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1060
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:1008
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1072
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1096
bool built() const
access function to the Built flag - indicates whether the matrix has been build - i....
Definition: matrices.h:1210
double * value()
Access to C-style value array.
Definition: matrices.h:1084
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1002
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.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
Abstract base class for matrices of doubles – adds abstract interfaces for solving,...
Definition: matrices.h:261
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
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.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
double norm() const
compute the 2 norm of this vector
bool built() const
double * values_pt()
access function to the underlying values
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
bool Output_info
Flag is true to output info and results of timings.
Definition: hypre_solver.h:305
bool Delete_input_data
Internal flag which is true when hypre_setup or hypre_solve can delete input matrix.
Definition: hypre_solver.h:445
unsigned AMG_coarsening
AMG coarsening strategy. Coarsening types include: 0 = CLJP (parallel coarsening using independent se...
Definition: hypre_solver.h:386
int Euclid_level
Euclid level parameter for ILU(k) factorization.
Definition: hypre_solver.h:418
double AMGEuclidSmoother_drop_tol
Definition: hypre_solver.h:283
bool Euclid_using_BJ
Flag to determine if Block Jacobi is used instead of PILU.
Definition: hypre_solver.h:415
bool Euclid_rowScale
Flag to switch on Euclid row scaling.
Definition: hypre_solver.h:409
double AMG_strength
Connection strength threshold parameter for BoomerAMG.
Definition: hypre_solver.h:369
LinearAlgebraDistribution * Hypre_distribution_pt
the distribution for this helpers-
Definition: hypre_solver.h:490
int ParaSails_nlevel
ParaSails nlevel parameter.
Definition: hypre_solver.h:397
double ParaSails_filter
ParaSails filter parameter.
Definition: hypre_solver.h:403
bool AMGEuclidSmoother_use_row_scaling
Definition: hypre_solver.h:280
unsigned AMGEuclidSmoother_level
Definition: hypre_solver.h:282
HYPRE_Solver Solver
The Hypre solver used in solve(...), resolve(...) or preconditioner_solve(...). [This is a C structur...
Definition: hypre_solver.h:477
int ParaSails_symmetry
ParaSails symmetry flag, used to inform ParaSails of Symmetry of definitenss of problem and type of P...
Definition: hypre_solver.h:394
HYPRE_ParCSRMatrix Matrix_par
The Hypre_ParCSRMatrix version of the matrix used in solve(...), resolve(...) or preconditioner_solve...
Definition: hypre_solver.h:473
unsigned Internal_preconditioner
Preconditioner method flag used with Hypre's PCG, GMRES and BiCGStab in solve(...) or resolve(....
Definition: hypre_solver.h:320
void hypre_solve(const DoubleVector &rhs, DoubleVector &solution)
Helper function performs a solve if any solver exists.
unsigned AMGEuclidSmoother_print_level
Definition: hypre_solver.h:284
bool AMGEuclidSmoother_use_block_jacobi
Definition: hypre_solver.h:279
unsigned AMG_complex_smoother
Complex smoothing methods used in BoomerAMG. Relaxation types are: 6 = Schwarz 7 = Pilut 8 = ParaSail...
Definition: hypre_solver.h:360
double ParaSails_thresh
ParaSails thresh parameter.
Definition: hypre_solver.h:400
double AMG_truncation
Interpolation truncation factor for BoomerAMG.
Definition: hypre_solver.h:372
unsigned Krylov_print_level
Used to set the Hypre printing level for the Krylov subspace solvers.
Definition: hypre_solver.h:429
unsigned AMG_smoother_iterations
The number of smoother iterations to apply.
Definition: hypre_solver.h:363
unsigned existing_solver()
Function to return value of which solver (if any) is currently stored.
Definition: hypre_solver.h:267
void hypre_clean_up_memory()
Function deletes all solver data.
bool Hypre_error_messages
Flag to determine if non-zero values of the Hypre error flag plus Hypre error messages are output to ...
Definition: hypre_solver.h:441
double Tolerance
Tolerance used to terminate solver.
Definition: hypre_solver.h:311
bool AMG_using_simple_smoothing
Flag to determine whether simple smoothers (determined by the AMG_simple_smoother flag) or complex sm...
Definition: hypre_solver.h:338
unsigned Existing_solver
Used to keep track of which solver (if any) is currently stored.
Definition: hypre_solver.h:484
unsigned Max_iter
Maximum number of iterations used in solver.
Definition: hypre_solver.h:308
double AMG_damping
Damping factor for BoomerAMG smoothed Jacobi or hybrid SOR.
Definition: hypre_solver.h:366
unsigned AMG_simple_smoother
Simple smoothing methods used in BoomerAMG. Relaxation types include: 0 = Jacobi 1 = Gauss-Seidel,...
Definition: hypre_solver.h:351
unsigned AMG_print_level
Used to set the Hypre printing level for AMG 0: no printout 1: print setup information 2: print solve...
Definition: hypre_solver.h:327
HYPRE_IJMatrix Matrix_ij
The Hypre_IJMatrix version of the matrix used in solve(...), resolve(...) or preconditioner_solve(....
Definition: hypre_solver.h:469
unsigned Euclid_print_level
Flag to set the level of printing from Euclid when the Euclid destructor is called 0: no printing (de...
Definition: hypre_solver.h:425
unsigned AMG_max_levels
Maximum number of levels used in AMG.
Definition: hypre_solver.h:330
unsigned Hypre_method
Hypre method flag. Valid values are specified in enumeration.
Definition: hypre_solver.h:314
void hypre_matrix_setup(CRDoubleMatrix *matrix_pt)
Function which sets values of First_global_row, Last_global_row and other partitioning data and creat...
double AMG_max_row_sum
Parameter to identify diagonally dominant parts of the matrix in AMG.
Definition: hypre_solver.h:333
double Euclid_droptol
Euclid drop tolerance for ILU(k) and ILUT factorization.
Definition: hypre_solver.h:406
unsigned Existing_preconditioner
Used to keep track of which preconditioner (if any) is currently stored.
Definition: hypre_solver.h:487
bool Euclid_using_ILUT
Flag to determine if ILUT (if true) or ILU(k) is used in Euclid.
Definition: hypre_solver.h:412
void hypre_solver_setup()
Sets up the data required for to use as an oomph-lib LinearSolver or Preconditioner....
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: hypre_solver.h:826
unsigned & amg_simple_smoother()
Access function to AMG_simple_smoother flag.
Definition: hypre_solver.h:983
static void reset_cumulative_solve_times()
Reset cumulative solve times.
static std::map< std::string, unsigned > Context_based_cumulative_npreconditioner_solve
Static unsigned that accumulates the number of preconditioner solves of all instantiations of this cl...
Definition: hypre_solver.h:904
std::string Context_string
String can be used to provide context for annotation.
void set_amg_iterations(const unsigned &amg_iterations)
Function to set the number of times to apply BoomerAMG.
Definition: hypre_solver.h:964
void clean_up_memory()
Function deletes all solver data.
static double Cumulative_preconditioner_solve_time
Static double that accumulates the preconditioner solve time of all instantiations of this class....
Definition: hypre_solver.h:882
void setup()
Function to set up a preconditioner for the linear system defined by matrix_pt. This function is requ...
static std::map< std::string, double > Context_based_cumulative_solve_time
Static double that accumulates the preconditioner solve time of all instantiations of this class,...
Definition: hypre_solver.h:889
bool Delete_matrix
Hypre copies matrix data from oomph-lib's CRDoubleMatrix or DistributedCRDoubleMatrix into its own da...
static std::map< std::string, unsigned > Context_based_nrow
Static unsigned that stores nrow for the most recent instantiations of this class,...
Definition: hypre_solver.h:911
void amg_using_simple_smoothing()
Function to select use of 'simple' AMG smoothers as controlled by the flag AMG_simple_smoother.
Definition: hypre_solver.h:977
double My_cumulative_preconditioner_solve_time
Private double that accumulates the preconditioner solve time of thi instantiation of this class....
double & amg_damping()
Access function to AMG_damping parameter.
unsigned & amg_coarsening()
Access function to AMG_coarsening flag.
static void report_cumulative_solve_times()
Report cumulative solve times of all instantiations of this class.
double & amg_strength()
Access function to AMG_strength.
static unsigned Cumulative_npreconditioner_solve
Static unsigned that accumulates the number of preconditioner solves of all instantiations of this cl...
Definition: hypre_solver.h:896
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Function applies solver to vector r for preconditioning. This requires a call to setup(....
unsigned & hypre_method()
Access function to Hypre_method flag – specified via enumeration.
Definition: hypre_solver.h:945
void solve(Problem *const &problem_pt, DoubleVector &solution)
Function which uses problem_pt's get_jacobian(...) function to generate a linear system which is then...
void clean_up_memory()
Function deletes all solver data.
bool Delete_matrix
Hypre copies matrix data from oomph-lib's CRDoubleMatrix or DistributedCRDoubleMatrix into its own da...
Definition: hypre_solver.h:806
void resolve(const DoubleVector &rhs, DoubleVector &solution)
Function to resolve a linear system using the existing solver data, allowing a solve with a new right...
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
void build(const OomphCommunicator *const comm_pt, const unsigned &first_row, const unsigned &nrow_local, const unsigned &nrow=0)
Sets the distribution. Takes first_row, nrow_local and nrow as arguments. If nrow is not provided or ...
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:77
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:73
static bool mpi_has_been_initialised()
return true if MPI has been initialised
static OomphCommunicator * communicator_pt()
access to global communicator. This is the oomph-lib equivalent of MPI_COMM_WORLD
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
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....
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
Definition: problem.cc:3890
int hypre__global_error
void create_HYPRE_Matrix(CRDoubleMatrix *oomph_matrix, HYPRE_IJMatrix &hypre_ij_matrix, HYPRE_ParCSRMatrix &hypre_par_matrix, LinearAlgebraDistribution *dist_pt)
Helper function to create a serial HYPRE_IJMatrix and HYPRE_ParCSRMatrix from a CRDoubleMatrix NOTE: ...
void create_HYPRE_Vector(const DoubleVector &oomph_vec, const LinearAlgebraDistribution *dist_pt, HYPRE_IJVector &hypre_ij_vector, HYPRE_ParVector &hypre_par_vector)
Helper function to create a HYPRE_IJVector and HYPRE_ParVector.
int check_HYPRE_error_flag(std::ostringstream &message)
Helper function to check the Hypre error flag, return the message associated with any error,...
unsigned AMG_coarsening
Default AMG coarsening strategy. Coarsening types include: 0 = CLJP (parallel coarsening using indepe...
double AMG_truncation
AMG interpolation truncation factor.
double AMG_strength
Default for AMG strength (0.25 recommended for 2D problems; larger (0.5-0.75, say) for 3D.
void euclid_settings_helper(const bool &use_block_jacobi, const bool &use_row_scaling, const bool &use_ilut, const int &level, const double &drop_tol, const int &print_level, HYPRE_Solver &euclid_object)
Helper function to set Euclid options using a command line like array.
void set_defaults_for_navier_stokes_momentum_block(HyprePreconditioner *hypre_preconditioner_pt)
Set default parameters for use as preconditioner in for momentum block in Navier-Stokes problem.
Definition: hypre_solver.cc:85
void set_defaults_for_3D_poisson_problem(HyprePreconditioner *hypre_preconditioner_pt)
Set default parameters for use as preconditioner in 3D Poisson-type problem.
Definition: hypre_solver.cc:72
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...