1 /*============================================================================
2 * Routines to handle the SLES settings
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <stdlib.h>
35 #include <string.h>
36
37 /*----------------------------------------------------------------------------
38 * Local headers
39 *----------------------------------------------------------------------------*/
40
41 #include <bft_error.h>
42 #include <bft_mem.h>
43 #include <bft_printf.h>
44
45 #include "cs_fp_exception.h"
46 #include "cs_log.h"
47 #include "cs_multigrid.h"
48 #include "cs_sles.h"
49
50 #if defined(HAVE_MUMPS)
51 #include "cs_sles_mumps.h"
52 #endif
53
54 #if defined(HAVE_PETSC)
55 #include "cs_sles_petsc.h"
56 #endif
57
58 /*----------------------------------------------------------------------------
59 * Header for the current file
60 *----------------------------------------------------------------------------*/
61
62 #include "cs_param_sles.h"
63
64 /*----------------------------------------------------------------------------*/
65
66 BEGIN_C_DECLS
67
68 /*============================================================================
69 * Type definitions
70 *============================================================================*/
71
72 /*============================================================================
73 * Local private variables
74 *============================================================================*/
75
76 /*============================================================================
77 * Private function prototypes
78 *============================================================================*/
79
80 /*----------------------------------------------------------------------------*/
81 /*!
82 * \brief Return true if a solver involving the MUMPS library is requested
83 */
84 /*----------------------------------------------------------------------------*/
85
86 static inline bool
_mumps_is_needed(cs_param_itsol_type_t solver)87 _mumps_is_needed(cs_param_itsol_type_t solver)
88 {
89 if (solver == CS_PARAM_ITSOL_MUMPS ||
90 solver == CS_PARAM_ITSOL_MUMPS_LDLT ||
91 solver == CS_PARAM_ITSOL_MUMPS_FLOAT ||
92 solver == CS_PARAM_ITSOL_MUMPS_FLOAT_LDLT)
93 return true;
94 else
95 return false;
96 }
97
98 /*----------------------------------------------------------------------------*/
99 /*!
100 * \brief Return true if the prescribed solver implies a symmetric linear
101 * system
102 */
103 /*----------------------------------------------------------------------------*/
104
105 static inline bool
_system_should_be_sym(cs_param_itsol_type_t solver)106 _system_should_be_sym(cs_param_itsol_type_t solver)
107 {
108 switch (solver) {
109
110 case CS_PARAM_ITSOL_CG:
111 case CS_PARAM_ITSOL_FCG:
112 case CS_PARAM_ITSOL_GKB_CG:
113 case CS_PARAM_ITSOL_GKB_GMRES:
114 case CS_PARAM_ITSOL_MINRES:
115 case CS_PARAM_ITSOL_MUMPS_LDLT:
116 case CS_PARAM_ITSOL_MUMPS_FLOAT_LDLT:
117 return true;
118
119 default:
120 return false; /* Assume that the system is not symmetric */
121 }
122 }
123
124 #if defined(HAVE_PETSC)
125 /*----------------------------------------------------------------------------*/
126 /*!
127 * \brief Set the command line option for PETSc
128 *
129 * \param[in] use_prefix need a prefix
130 * \param[in] prefix optional prefix
131 * \param[in] keyword command keyword
132 * \param[in] keyval command value
133 */
134 /*----------------------------------------------------------------------------*/
135
136 static inline void
_petsc_cmd(bool use_prefix,const char * prefix,const char * keyword,const char * keyval)137 _petsc_cmd(bool use_prefix,
138 const char *prefix,
139 const char *keyword,
140 const char *keyval)
141 {
142 char cmd_line[128];
143
144 if (use_prefix)
145 sprintf(cmd_line, "-%s_%s", prefix, keyword);
146 else
147 sprintf(cmd_line, "-%s", keyword);
148
149 #if PETSC_VERSION_GE(3,7,0)
150 PetscOptionsSetValue(NULL, cmd_line, keyval);
151 #else
152 PetscOptionsSetValue(cmd_line, keyval);
153 #endif
154 }
155
156 /*----------------------------------------------------------------------------*/
157 /*!
158 * \brief Predefined settings for a block ILU(0) with PETSc
159 *
160 * \param[in] prefix prefix name associated to the current SLES
161 */
162 /*----------------------------------------------------------------------------*/
163
164 static inline void
_petsc_bilu0_hook(const char * prefix)165 _petsc_bilu0_hook(const char *prefix)
166 {
167 /* Sanity checks */
168 assert(prefix != NULL);
169
170 _petsc_cmd(true, prefix, "pc_type", "bjacobi");
171 _petsc_cmd(true, prefix, "pc_jacobi_blocks", "1");
172 _petsc_cmd(true, prefix, "sub_ksp_type", "preonly");
173 _petsc_cmd(true, prefix, "sub_pc_type", "ilu");
174 _petsc_cmd(true, prefix, "sub_pc_factor_level", "0");
175 _petsc_cmd(true, prefix, "sub_pc_factor_reuse_ordering", "");
176 /* If one wants to optimize the memory consumption */
177 /* _petsc_cmd(true, prefix, "sub_pc_factor_in_place", ""); */
178 }
179
180 /*----------------------------------------------------------------------------*/
181 /*!
182 * \brief Predefined settings for a block ICC(0) with PETSc
183 *
184 * \param[in] prefix prefix name associated to the current SLES
185 */
186 /*----------------------------------------------------------------------------*/
187
188 static inline void
_petsc_bicc0_hook(const char * prefix)189 _petsc_bicc0_hook(const char *prefix)
190 {
191 /* Sanity checks */
192 assert(prefix != NULL);
193
194 _petsc_cmd(true, prefix, "pc_type", "bjacobi");
195 _petsc_cmd(true, prefix, "pc_jacobi_blocks", "1");
196 _petsc_cmd(true, prefix, "sub_ksp_type", "preonly");
197 _petsc_cmd(true, prefix, "sub_pc_type", "icc");
198 _petsc_cmd(true, prefix, "sub_pc_factor_level", "0");
199 _petsc_cmd(true, prefix, "sub_pc_factor_reuse_ordering", "");
200 /* If one wants to optimize the memory consumption */
201 /* _petsc_cmd(true, prefix, "sub_pc_factor_in_place", ""); */
202 }
203
204 /*----------------------------------------------------------------------------*/
205 /*!
206 * \brief Predefined settings for a block SSOR with PETSc
207 *
208 * \param[in] prefix prefix name associated to the current SLES
209 */
210 /*----------------------------------------------------------------------------*/
211
212 static inline void
_petsc_bssor_hook(const char * prefix)213 _petsc_bssor_hook(const char *prefix)
214 {
215 /* Sanity checks */
216 assert(prefix != NULL);
217
218 _petsc_cmd(true, prefix, "pc_type", "bjacobi");
219 _petsc_cmd(true, prefix, "pc_jacobi_blocks", "1");
220 _petsc_cmd(true, prefix, "sub_ksp_type", "preonly");
221 _petsc_cmd(true, prefix, "sub_pc_type", "sor");
222 _petsc_cmd(true, prefix, "sub_pc_sor_symmetric", "");
223 _petsc_cmd(true, prefix, "sub_pc_sor_local_symmetric", "");
224 _petsc_cmd(true, prefix, "sub_pc_sor_omega", "1.5");
225 }
226
227 /*----------------------------------------------------------------------------*/
228 /*!
229 * \brief Predefined settings for GAMG as a preconditioner even if another
230 * settings have been defined. One assumes that one really wants to use
231 * GAMG (may be HYPRE is not available)
232 *
233 * \param[in] prefix prefix name associated to the current SLES
234 * \param[in] slesp pointer to a set of SLES parameters
235 * \param[in] is_symm the linear system to solve is symmetric
236 * \param[in, out] pc pointer to a PETSc preconditioner
237 */
238 /*----------------------------------------------------------------------------*/
239
240 static inline void
_petsc_pcgamg_hook(const char * prefix,const cs_param_sles_t * slesp,bool is_symm,PC pc)241 _petsc_pcgamg_hook(const char *prefix,
242 const cs_param_sles_t *slesp,
243 bool is_symm,
244 PC pc)
245 {
246 /* Sanity checks */
247 assert(prefix != NULL);
248 assert(slesp != NULL);
249 assert(slesp->precond == CS_PARAM_PRECOND_AMG);
250
251 /* Remark: -pc_gamg_reuse_interpolation
252 *
253 * Reuse prolongation when rebuilding algebraic multigrid
254 * preconditioner. This may negatively affect the convergence rate of the
255 * method on new matrices if the matrix entries change a great deal, but
256 * allows rebuilding the preconditioner quicker. (default=false)
257 */
258
259 _petsc_cmd(true, prefix, "pc_gamg_reuse_interpolation", "true");
260
261 /* Remark: -pc_gamg_sym_graph
262 * Symmetrize the graph before computing the aggregation. Some algorithms
263 * require the graph be symmetric (default=false)
264 */
265
266 _petsc_cmd(true, prefix, "pc_gamg_sym_graph", "true");
267
268 /* Set smoothers (general settings, i.e. not depending on the symmetry or not
269 of the linear system to solve) */
270
271 _petsc_cmd(true, prefix, "mg_levels_ksp_type", "richardson");
272 _petsc_cmd(true, prefix, "mg_levels_ksp_max_it", "1");
273 _petsc_cmd(true, prefix, "mg_levels_ksp_norm_type", "none");
274 _petsc_cmd(true, prefix, "mg_levels_ksp_richardson_scale", "1.0");
275
276 /* Do not build a coarser level if one reaches the following limit */
277 _petsc_cmd(true, prefix, "pc_gamg_coarse_eq_limit", "100");
278
279 /* In parallel computing, migrate data to another rank if the grid has less
280 than 200 rows */
281 if (cs_glob_n_ranks > 1) {
282
283 _petsc_cmd(true, prefix, "pc_gamg_repartition", "true");
284 _petsc_cmd(true, prefix, "pc_gamg_process_eq_limit", "200");
285
286 }
287 else {
288
289 _petsc_cmd(true, prefix, "mg_coarse_ksp_type", "preonly");
290 _petsc_cmd(true, prefix, "mg_coarse_pc_type", "tfs");
291
292 }
293
294 /* Settings depending on the symmetry or not of the linear system to solve */
295
296 if (is_symm) {
297
298 /* Remark: -pc_gamg_square_graph
299 *
300 * Squaring the graph increases the rate of coarsening (aggressive
301 * coarsening) and thereby reduces the complexity of the coarse grids, and
302 * generally results in slower solver converge rates. Reducing coarse grid
303 * complexity reduced the complexity of Galerkin coarse grid construction
304 * considerably. (default = 1)
305 *
306 * Remark: -pc_gamg_threshold
307 *
308 * Increasing the threshold decreases the rate of coarsening. Conversely
309 * reducing the threshold increases the rate of coarsening (aggressive
310 * coarsening) and thereby reduces the complexity of the coarse grids, and
311 * generally results in slower solver converge rates. Reducing coarse grid
312 * complexity reduced the complexity of Galerkin coarse grid construction
313 * considerably. Before coarsening or aggregating the graph, GAMG removes
314 * small values from the graph with this threshold, and thus reducing the
315 * coupling in the graph and a different (perhaps better) coarser set of
316 * points. (default=0.0) */
317
318 _petsc_cmd(true, prefix, "pc_gamg_agg_nsmooths", "2");
319 _petsc_cmd(true, prefix, "pc_gamg_square_graph", "2");
320 _petsc_cmd(true, prefix, "pc_gamg_threshold", "0.08");
321
322 if (cs_glob_n_ranks > 1) {
323
324 _petsc_cmd(true, prefix, "mg_levels_pc_type", "bjacobi");
325 _petsc_cmd(true, prefix, "mg_levels_pc_jacobi_blocks", "1");
326 _petsc_cmd(true, prefix, "mg_levels_sub_ksp_type", "preonly");
327 _petsc_cmd(true, prefix, "mg_levels_sub_pc_type", "sor");
328 _petsc_cmd(true, prefix, "mg_levels_sub_pc_sor_local_symmetric", "");
329 _petsc_cmd(true, prefix, "mg_levels_sub_pc_sor_omega", "1.5");
330
331 }
332 else { /* serial run */
333
334 _petsc_cmd(true, prefix, "mg_levels_pc_type", "sor");
335 _petsc_cmd(true, prefix, "mg_levels_pc_sor_local_symmetric", "");
336 _petsc_cmd(true, prefix, "mg_levels_pc_sor_omega", "1.5");
337
338 }
339
340 }
341 else { /* Not a symmetric linear system */
342
343 /* Number of smoothing steps to use with smooth aggregation (default=1) */
344 _petsc_cmd(true, prefix, "pc_gamg_agg_nsmooths", "0");
345 _petsc_cmd(true, prefix, "pc_gamg_square_graph", "0");
346 _petsc_cmd(true, prefix, "pc_gamg_threshold", "0.06");
347
348 _petsc_cmd(true, prefix, "mg_levels_pc_type", "bjacobi");
349 _petsc_cmd(true, prefix, "mg_levels_pc_bjacobi_blocks", "1");
350 _petsc_cmd(true, prefix, "mg_levels_sub_ksp_type", "preonly");
351 _petsc_cmd(true, prefix, "mg_levels_sub_pc_type", "ilu");
352 _petsc_cmd(true, prefix, "mg_levels_sub_pc_factor_levels", "0");
353
354 }
355
356 /* After command line options, switch to PETSc setup functions */
357
358 PCSetType(pc, PCGAMG);
359 PCGAMGSetType(pc, PCGAMGAGG);
360 PCGAMGSetNSmooths(pc, 1);
361 PCSetUp(pc);
362
363 switch (slesp->amg_type) {
364
365 case CS_PARAM_AMG_PETSC_GAMG_V:
366 case CS_PARAM_AMG_PETSC_PCMG:
367 case CS_PARAM_AMG_HYPRE_BOOMER_V:
368 PCMGSetCycleType(pc, PC_MG_CYCLE_V);
369 break;
370
371 case CS_PARAM_AMG_PETSC_GAMG_W:
372 case CS_PARAM_AMG_HYPRE_BOOMER_W:
373 PCMGSetCycleType(pc, PC_MG_CYCLE_W);
374 break;
375
376 default:
377 bft_error(__FILE__, __LINE__, 0, "%s: Invalid type of AMG for SLES %s\n",
378 __func__, slesp->name);
379 }
380
381 }
382
383 /*----------------------------------------------------------------------------*/
384 /*!
385 * \brief Predefined settings for BoomerAMG in HYPRE as a preconditioner
386 *
387 * \param[in] prefix prefix name associated to the current SLES
388 * \param[in] slesp pointer to a set of SLES parameters
389 * \param[in] is_symm the linear system to solve is symmetric
390 * \param[in, out] pc pointer to a PETSc preconditioner
391 */
392 /*----------------------------------------------------------------------------*/
393
394 static inline void
_petsc_pchypre_hook(const char * prefix,const cs_param_sles_t * slesp,bool is_symm,PC pc)395 _petsc_pchypre_hook(const char *prefix,
396 const cs_param_sles_t *slesp,
397 bool is_symm,
398 PC pc)
399 {
400 CS_UNUSED(is_symm);
401
402 /* Sanity checks */
403
404 assert(prefix != NULL);
405 assert(slesp != NULL);
406 assert(slesp->precond == CS_PARAM_PRECOND_AMG);
407
408 PCSetType(pc, PCHYPRE);
409 PCHYPRESetType(pc, "boomeramg");
410
411 switch (slesp->amg_type) {
412
413 case CS_PARAM_AMG_HYPRE_BOOMER_V:
414 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_cycle_type", "V");
415 break;
416
417 case CS_PARAM_AMG_HYPRE_BOOMER_W:
418 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_cycle_type", "W");
419 break;
420
421 default:
422 bft_error(__FILE__, __LINE__, 0, "%s: Invalid type of AMG for SLES %s\n",
423 __func__, slesp->name);
424 }
425
426 /* From HYPRE documentation: https://hypre.readthedocs.io/en/lastest
427 *
428 * for three-dimensional diffusion problems, it is recommended to choose a
429 * lower complexity coarsening like HMIS or PMIS (coarsening 10 or 8) and
430 * combine it with a distance-two interpolation (interpolation 6 or 7), that
431 * is also truncated to 4 or 5 elements per row. Additional reduction in
432 * complexity and increased scalability can often be achieved using one or
433 * two levels of aggressive coarsening.
434 */
435
436 /* _petsc_cmd(true, prefix, "pc_hypre_boomeramg_grid_sweeps_down","2"); */
437 /* _petsc_cmd(true, prefix, "pc_hypre_boomeramg_grid_sweeps_up","2"); */
438 /* _petsc_cmd(true, prefix, "pc_hypre_boomeramg_smooth_type","Euclid"); */
439
440 /* Remark: fcf-jacobi or l1scaled-jacobi (or chebyshev) as up/down smoothers
441 can be a good choice
442 */
443
444 /*
445 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_relax_type_down","fcf-jacobi");
446 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_relax_type_up","fcf-jacobi");
447 */
448
449 /* Note that the default coarsening is HMIS in HYPRE */
450
451 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_coarsen_type", "HMIS");
452
453 /* Note that the default interpolation is extended+i interpolation truncated
454 * to 4 elements per row. Using 0 means there is no limitation.
455 * good choices are: ext+i-cc, ext+i, FF1
456 */
457
458 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_interp_type", "ext+i-cc");
459 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_P_max","8");
460
461 /* Number of levels (starting from the finest one) on which one applies an
462 aggressive coarsening */
463
464 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_agg_nl","2");
465
466 /* Number of paths for aggressive coarsening (default = 1) */
467
468 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_agg_num_paths","2");
469
470 /* For best performance, it might be necessary to set certain parameters,
471 * which will affect both coarsening and interpolation. One important
472 * parameter is the strong threshold. The default value is 0.25, which
473 * appears to be a good choice for 2-dimensional problems and the low
474 * complexity coarsening algorithms. For 3-dimensional problems a better
475 * choice appears to be 0.5, when using the default coarsening
476 * algorithm. However, the choice of the strength threshold is problem
477 * dependent.
478 */
479
480 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_strong_threshold","0.5");
481 _petsc_cmd(true, prefix, "pc_hypre_boomeramg_no_CF","");
482 }
483
484 /*----------------------------------------------------------------------------*/
485 /*!
486 * \brief Set command line options for PC according to the kind of
487 * preconditionner
488 *
489 * \param[in, out] slesp set of parameters for the linear algebra
490 * \param[in, out] ksp PETSc solver structure
491 */
492 /*----------------------------------------------------------------------------*/
493
494 static void
_petsc_set_pc_type(cs_param_sles_t * slesp,KSP ksp)495 _petsc_set_pc_type(cs_param_sles_t *slesp,
496 KSP ksp)
497 {
498 if (_mumps_is_needed(slesp->solver))
499 return; /* Direct solver: Nothing to do at this stage */
500
501 PC pc;
502 KSPGetPC(ksp, &pc);
503
504 switch (slesp->precond) {
505
506 case CS_PARAM_PRECOND_NONE:
507 PCSetType(pc, PCNONE);
508 break;
509
510 case CS_PARAM_PRECOND_DIAG:
511 PCSetType(pc, PCJACOBI);
512 break;
513
514 case CS_PARAM_PRECOND_BJACOB_ILU0:
515 if (slesp->solver_class == CS_PARAM_SLES_CLASS_HYPRE) {
516 #if defined(PETSC_HAVE_HYPRE)
517 PCSetType(pc, PCHYPRE);
518 PCHYPRESetType(pc, "euclid");
519 _petsc_cmd(true, slesp->name, "pc_euclid_level", "0");
520 #else
521 _petsc_bilu0_hook(slesp->name);
522 #endif
523 }
524 else
525 _petsc_bilu0_hook(slesp->name);
526 break;
527
528 case CS_PARAM_PRECOND_BJACOB_SGS:
529 _petsc_bssor_hook(slesp->name);
530 break;
531
532 case CS_PARAM_PRECOND_SSOR:
533 if (cs_glob_n_ranks > 1) { /* Switch to a block version */
534
535 slesp->precond = CS_PARAM_PRECOND_BJACOB_SGS;
536 cs_base_warn(__FILE__, __LINE__);
537 cs_log_printf(CS_LOG_DEFAULT,
538 " %s: System %s: Modify the requested preconditioner to"
539 " enable a parallel computation with PETSC.\n"
540 " Switch to a block jacobi preconditioner.\n",
541 __func__, slesp->name);
542
543 _petsc_bssor_hook(slesp->name);
544
545 }
546 else { /* Serial computation */
547 PCSetType(pc, PCSOR);
548 PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP);
549 }
550 break;
551
552 case CS_PARAM_PRECOND_ICC0:
553 if (cs_glob_n_ranks > 1) { /* Switch to a block version */
554
555 cs_base_warn(__FILE__, __LINE__);
556 cs_log_printf(CS_LOG_DEFAULT,
557 " %s: System %s: Modify the requested preconditioner to"
558 " enable a parallel computation with PETSC.\n"
559 " Switch to a block jacobi preconditioner.\n",
560 __func__, slesp->name);
561
562 _petsc_bicc0_hook(slesp->name);
563
564 }
565 else {
566 PCSetType(pc, PCICC);
567 PCFactorSetLevels(pc, 0);
568 }
569 break;
570
571 case CS_PARAM_PRECOND_ILU0:
572 if (slesp->solver_class == CS_PARAM_SLES_CLASS_HYPRE) {
573 #if defined(PETSC_HAVE_HYPRE)
574 /* Euclid is a parallel version of the ILU(0) factorisation */
575 PCSetType(pc, PCHYPRE);
576 PCHYPRESetType(pc, "euclid");
577 _petsc_cmd(true, slesp->name, "pc_euclid_level", "0");
578 #else
579 _petsc_bilu0_hook(slesp->name);
580 if (cs_glob_n_ranks > 1) /* Switch to a block version */
581 slesp->precond = CS_PARAM_PRECOND_BJACOB_ILU0;
582 #endif
583 }
584 else {
585 _petsc_bilu0_hook(slesp->name);
586 if (cs_glob_n_ranks > 1) { /* Switch to a block version */
587
588 slesp->precond = CS_PARAM_PRECOND_BJACOB_ILU0;
589 cs_base_warn(__FILE__, __LINE__);
590 cs_log_printf(CS_LOG_DEFAULT,
591 " %s: System %s: Modify the requested preconditioner to"
592 " enable a parallel computation with PETSC.\n"
593 " Switch to a block jacobi preconditioner.\n",
594 __func__, slesp->name);
595 }
596 }
597 break;
598
599 case CS_PARAM_PRECOND_LU:
600 #if defined(PETSC_HAVE_MUMPS)
601 _petsc_cmd(true, slesp->name, "pc_type", "lu");
602 _petsc_cmd(true, slesp->name, "pc_factor_mat_solver_type", "mumps");
603 #else
604 if (cs_glob_n_ranks == 1)
605 _petsc_cmd(true, slesp->name, "pc_type", "lu");
606 else { /* Switch to a block version (sequential in each block) */
607 _petsc_cmd(true, slesp->name, "pc_type", "bjacobi");
608 _petsc_cmd(true, slesp->name, "pc_jacobi_blocks", "1");
609 _petsc_cmd(true, slesp->name, "sub_ksp_type", "preonly");
610 _petsc_cmd(true, slesp->name, "sub_pc_type", "lu");
611 }
612 #endif
613 break;
614
615 case CS_PARAM_PRECOND_AMG:
616 {
617 bool is_symm = _system_should_be_sym(slesp->solver);
618
619 switch (slesp->amg_type) {
620
621 case CS_PARAM_AMG_PETSC_GAMG_V:
622 case CS_PARAM_AMG_PETSC_GAMG_W:
623 case CS_PARAM_AMG_PETSC_PCMG:
624 _petsc_pcgamg_hook(slesp->name, slesp, is_symm, pc);
625 break;
626
627 case CS_PARAM_AMG_HYPRE_BOOMER_V:
628 case CS_PARAM_AMG_HYPRE_BOOMER_W:
629 #if defined(PETSC_HAVE_HYPRE)
630 _petsc_pchypre_hook(slesp->name, slesp, is_symm, pc);
631 #else
632 cs_base_warn(__FILE__, __LINE__);
633 cs_log_printf(CS_LOG_DEFAULT,
634 "%s: Eq. %s: Switch to GAMG since BoomerAMG is not"
635 " available.\n",
636 __func__, slesp->name);
637 _petsc_pcgamg_hook(slesp->name, slesp, is_symm, pc);
638 #endif
639 break;
640
641 default:
642 bft_error(__FILE__, __LINE__, 0,
643 " %s: Eq. %s: Invalid AMG type for the PETSc library.",
644 __func__, slesp->name);
645 break;
646
647 } /* End of switch on the AMG type */
648
649 } /* AMG as preconditioner */
650 break;
651
652 default:
653 bft_error(__FILE__, __LINE__, 0,
654 " %s: Eq. %s: Preconditioner not interfaced with PETSc.",
655 __func__, slesp->name);
656 }
657
658 /* Apply modifications to the PC structure given with command lines.
659 * This setting stands for a first setting and may be overwritten with
660 * parameters stored in the structure cs_param_sles_t
661 * To get the last word use cs_user_sles_petsc_hook() */
662 PCSetFromOptions(pc);
663 PCSetUp(pc);
664 }
665
666 /*----------------------------------------------------------------------------*/
667 /*!
668 * \brief Set PETSc solver
669 *
670 * \param[in] slesp pointer to SLES parameters
671 * \param[in, out] ksp pointer to PETSc KSP context
672 */
673 /*----------------------------------------------------------------------------*/
674
675 static void
_petsc_set_krylov_solver(cs_param_sles_t * slesp,KSP ksp)676 _petsc_set_krylov_solver(cs_param_sles_t *slesp,
677 KSP ksp)
678 {
679 /* No choice otherwise PETSc yields an error */
680 slesp->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
681 KSPSetNormType(ksp, KSP_NORM_UNPRECONDITIONED);
682
683 /* 2) Set the krylov solver */
684 switch (slesp->solver) {
685
686 case CS_PARAM_ITSOL_NONE:
687 KSPSetType(ksp, KSPPREONLY);
688 break;
689
690 case CS_PARAM_ITSOL_BICG: /* Improved Bi-CG stab */
691 KSPSetType(ksp, KSPIBCGS);
692 break;
693
694 case CS_PARAM_ITSOL_BICGSTAB2: /* Preconditioned BiCGstab2 */
695 KSPSetType(ksp, KSPBCGSL);
696 break;
697
698 case CS_PARAM_ITSOL_CG: /* Preconditioned Conjugate Gradient */
699 if (slesp->precond == CS_PARAM_PRECOND_AMG)
700 KSPSetType(ksp, KSPFCG);
701 else
702 KSPSetType(ksp, KSPCG);
703 break;
704
705 case CS_PARAM_ITSOL_FCG: /* Flexible Conjuguate Gradient */
706 KSPSetType(ksp, KSPFCG);
707 break;
708
709 case CS_PARAM_ITSOL_FGMRES: /* Preconditioned flexible GMRES */
710 KSPSetType(ksp, KSPFGMRES);
711 break;
712
713 case CS_PARAM_ITSOL_GCR: /* Generalized Conjuguate Residual */
714 KSPSetType(ksp, KSPGCR);
715 break;
716
717 case CS_PARAM_ITSOL_GMRES: /* Preconditioned GMRES */
718 KSPSetType(ksp, KSPLGMRES);
719 break;
720
721 case CS_PARAM_ITSOL_MINRES: /* Minimal residual */
722 KSPSetType(ksp, KSPMINRES);
723 break;
724
725 case CS_PARAM_ITSOL_MUMPS: /* Direct solver (factorization) */
726 case CS_PARAM_ITSOL_MUMPS_LDLT:
727 #if defined(PETSC_HAVE_MUMPS)
728 KSPSetType(ksp, KSPPREONLY);
729 #else
730 bft_error(__FILE__, __LINE__, 0,
731 " %s: MUMPS not interfaced with this installation of PETSc.",
732 __func__);
733 #endif
734 break;
735
736 default:
737 bft_error(__FILE__, __LINE__, 0,
738 " %s: Iterative solver not interfaced with PETSc.", __func__);
739 }
740
741 /* 3) Additional settings arising from command lines */
742 switch (slesp->solver) {
743
744 case CS_PARAM_ITSOL_GMRES: /* Preconditioned GMRES */
745 _petsc_cmd(true, slesp->name, "ksp_gmres_modifiedgramschmidt", "1");
746 break;
747
748 default:
749 break; /* Nothing to do. Settings performed with another mechanism */
750
751 }
752
753 /* Apply modifications to the KSP structure given with command lines.
754 * This setting stands for a first setting and may be overwritten with
755 * parameters stored in the structure cs_param_sles_t
756 *
757 * Automatic monitoring
758 * PetscOptionsSetValue(NULL, "-ksp_monitor", "");
759 *
760 */
761
762 KSPSetFromOptions(ksp);
763
764 /* Apply settings from the cs_param_sles_t structure */
765 switch (slesp->solver) {
766
767 case CS_PARAM_ITSOL_GMRES: /* Preconditioned GMRES */
768 case CS_PARAM_ITSOL_FGMRES: /* Flexible GMRES */
769 KSPGMRESSetRestart(ksp, slesp->restart);
770 break;
771
772 case CS_PARAM_ITSOL_GCR: /* Preconditioned GCR */
773 KSPGCRSetRestart(ksp, slesp->restart);
774 break;
775
776 #if defined(PETSC_HAVE_MUMPS)
777 case CS_PARAM_ITSOL_MUMPS:
778 {
779 PC pc;
780 KSPGetPC(ksp, &pc);
781 PCSetType(pc, PCLU);
782 PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
783 }
784 break;
785
786 case CS_PARAM_ITSOL_MUMPS_LDLT:
787 {
788 PC pc;
789 KSPGetPC(ksp, &pc);
790
791 /* Retrieve the matrices related to this KSP */
792 Mat a, pa;
793 KSPGetOperators(ksp, &a, &pa);
794
795 MatSetOption(a, MAT_SPD, PETSC_TRUE); /* set MUMPS id%SYM=1 */
796 PCSetType(pc, PCCHOLESKY);
797
798 PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
799 PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
800 }
801 break;
802 #endif
803
804 default:
805 break; /* Nothing else to do */
806 }
807
808 /* Set KSP tolerances */
809 PetscReal rtol, abstol, dtol;
810 PetscInt maxit;
811 KSPGetTolerances(ksp, &rtol, &abstol, &dtol, &maxit);
812 KSPSetTolerances(ksp,
813 slesp->eps, /* relative convergence tolerance */
814 abstol, /* absolute convergence tolerance */
815 dtol, /* divergence tolerance */
816 slesp->n_max_iter); /* max number of iterations */
817
818 }
819
820 /*----------------------------------------------------------------------------*/
821 /*!
822 * \brief Set PETSc solver and preconditioner
823 *
824 * \param[in, out] context pointer to optional (untyped) value or structure
825 * \param[in, out] ksp pointer to PETSc KSP context
826 */
827 /*----------------------------------------------------------------------------*/
828
829 static void
_petsc_setup_hook(void * context,KSP ksp)830 _petsc_setup_hook(void *context,
831 KSP ksp)
832 {
833 cs_param_sles_t *slesp = (cs_param_sles_t *)context;
834
835 cs_fp_exception_disable_trap(); /* Avoid trouble with a too restrictive
836 SIGFPE detection */
837
838 int len = strlen(slesp->name) + 1;
839 char *prefix = NULL;
840 BFT_MALLOC(prefix, len + 1, char);
841 sprintf(prefix, "%s_", slesp->name);
842 prefix[len] = '\0';
843 KSPSetOptionsPrefix(ksp, prefix);
844 BFT_FREE(prefix);
845
846 /* 1) Set the solver */
847 _petsc_set_krylov_solver(slesp, ksp);
848
849 /* 2) Set the preconditioner */
850 _petsc_set_pc_type(slesp, ksp);
851
852 /* 3) User function for additional settings */
853 cs_user_sles_petsc_hook((void *)slesp, ksp);
854
855 /* Dump the setup related to PETSc in a specific file */
856 if (!slesp->setup_done) {
857 KSPSetUp(ksp);
858 cs_sles_petsc_log_setup(ksp);
859 slesp->setup_done = true;
860 }
861
862 cs_fp_exception_restore_trap(); /* Avoid trouble with a too restrictive
863 SIGFPE detection */
864 }
865
866 /*----------------------------------------------------------------------------*/
867 /*!
868 * \brief Common settings for block preconditioning (when a system is split
869 * according to the Cartesian components: x,y,z)
870 *
871 * \param[in] slesp pointer to the SLES parameter settings
872 * \param[in, out] ksp pointer to PETSc KSP structure (solver)
873 */
874 /*----------------------------------------------------------------------------*/
875
876 static void
_petsc_common_block_hook(const cs_param_sles_t * slesp,KSP ksp)877 _petsc_common_block_hook(const cs_param_sles_t *slesp,
878 KSP ksp)
879 {
880 PC pc;
881 KSPGetPC(ksp, &pc);
882 PCSetType(pc, PCFIELDSPLIT);
883
884 switch (slesp->pcd_block_type) {
885 case CS_PARAM_PRECOND_BLOCK_UPPER_TRIANGULAR:
886 case CS_PARAM_PRECOND_BLOCK_LOWER_TRIANGULAR:
887 case CS_PARAM_PRECOND_BLOCK_FULL_UPPER_TRIANGULAR:
888 case CS_PARAM_PRECOND_BLOCK_FULL_LOWER_TRIANGULAR:
889 PCFieldSplitSetType(pc, PC_COMPOSITE_MULTIPLICATIVE);
890 break;
891
892 case CS_PARAM_PRECOND_BLOCK_SYM_GAUSS_SEIDEL:
893 case CS_PARAM_PRECOND_BLOCK_FULL_SYM_GAUSS_SEIDEL:
894 PCFieldSplitSetType(pc, PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE);
895 break;
896
897 case CS_PARAM_PRECOND_BLOCK_DIAG:
898 case CS_PARAM_PRECOND_BLOCK_FULL_DIAG:
899 default:
900 PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
901 break;
902 }
903
904 /* Apply modifications to the KSP structure */
905 PCFieldSplitSetBlockSize(pc, 3);
906
907 PetscInt id = 0;
908 PCFieldSplitSetFields(pc, "x", 1, &id, &id);
909 id = 1;
910 PCFieldSplitSetFields(pc, "y", 1, &id, &id);
911 id = 2;
912 PCFieldSplitSetFields(pc, "z", 1, &id, &id);
913 }
914
915 /*----------------------------------------------------------------------------*/
916 /*!
917 * \brief Function pointer: setup hook for setting PETSc solver and
918 * preconditioner.
919 * Case of multiplicative AMG block preconditioner for a CG with GAMG
920 * as AMG type
921 *
922 * \param[in, out] context pointer to optional (untyped) value or structure
923 * \param[in, out] ksp pointer to PETSc KSP context
924 */
925 /*----------------------------------------------------------------------------*/
926
927 static void
_petsc_amg_block_gamg_hook(void * context,KSP ksp)928 _petsc_amg_block_gamg_hook(void *context,
929 KSP ksp)
930 {
931 cs_param_sles_t *slesp = (cs_param_sles_t *)context;
932
933 cs_fp_exception_disable_trap(); /* Avoid trouble with a too restrictive
934 SIGFPE detection */
935
936 /* prefix will be extended with the fieldsplit */
937 int len = strlen(slesp->name) + 1;
938 int _len = len + strlen("_fieldsplit_x_") + 1;
939 char *prefix = NULL;
940 BFT_MALLOC(prefix, _len + 1, char);
941 sprintf(prefix, "%s_", slesp->name);
942 prefix[len] = '\0';
943 KSPSetOptionsPrefix(ksp, prefix);
944
945 /* Set the solver */
946 _petsc_set_krylov_solver(slesp, ksp);
947
948 /* Common settings to block preconditionner */
949 _petsc_common_block_hook(slesp, ksp);
950
951 PC pc;
952 KSPGetPC(ksp, &pc);
953 PCSetUp(pc);
954
955 PC _pc;
956 PetscInt n_split;
957 KSP *xyz_subksp;
958 const char xyz[3] = "xyz";
959 bool is_symm = _system_should_be_sym(slesp->solver);
960
961 PCFieldSplitGetSubKSP(pc, &n_split, &xyz_subksp);
962 assert(n_split == 3);
963
964 for (PetscInt id = 0; id < n_split; id++) {
965
966 sprintf(prefix, "%s_fieldsplit_%c", slesp->name, xyz[id]);
967
968 _petsc_cmd(true, prefix, "ksp_type","preonly");
969
970 /* Predefined settings when using AMG as a preconditioner */
971 KSP _ksp = xyz_subksp[id];
972 KSPGetPC(_ksp, &_pc);
973
974 _petsc_pcgamg_hook(prefix, slesp, is_symm, _pc);
975
976 PCSetFromOptions(_pc);
977 KSPSetFromOptions(_ksp);
978
979 } /* Loop on block settings */
980
981 BFT_FREE(prefix);
982 PetscFree(xyz_subksp);
983
984 /* User function for additional settings */
985 cs_user_sles_petsc_hook(context, ksp);
986
987 PCSetFromOptions(pc);
988 KSPSetFromOptions(ksp);
989
990 /* Dump the setup related to PETSc in a specific file */
991 if (!slesp->setup_done) {
992 KSPSetUp(ksp);
993 cs_sles_petsc_log_setup(ksp);
994 slesp->setup_done = true;
995 }
996
997 cs_fp_exception_restore_trap(); /* Avoid trouble with a too restrictive
998 SIGFPE detection */
999 }
1000
1001 #if defined(PETSC_HAVE_HYPRE)
1002 /*----------------------------------------------------------------------------*/
1003 /*!
1004 * \brief Function pointer: setup hook for setting PETSc solver and
1005 * preconditioner.
1006 * Case of multiplicative AMG block preconditioner for a CG with boomer
1007 * as AMG type
1008 *
1009 * \param[in, out] context pointer to optional (untyped) value or structure
1010 * \param[in, out] ksp pointer to PETSc KSP context
1011 */
1012 /*----------------------------------------------------------------------------*/
1013
1014 static void
_petsc_amg_block_boomer_hook(void * context,KSP ksp)1015 _petsc_amg_block_boomer_hook(void *context,
1016 KSP ksp)
1017 {
1018 cs_param_sles_t *slesp = (cs_param_sles_t *)context;
1019
1020 cs_fp_exception_disable_trap(); /* Avoid trouble with a too restrictive
1021 SIGFPE detection */
1022
1023 /* prefix will be extended with the fieldsplit */
1024 int len = strlen(slesp->name) + 1;
1025 int _len = len + strlen("_fieldsplit_x_") + 1;
1026 char *prefix = NULL;
1027 BFT_MALLOC(prefix, _len + 1, char);
1028 sprintf(prefix, "%s_", slesp->name);
1029 prefix[len] = '\0';
1030 KSPSetOptionsPrefix(ksp, prefix);
1031
1032 /* Set the solver */
1033 _petsc_set_krylov_solver(slesp, ksp);
1034
1035 /* Common settings to block preconditionner */
1036 _petsc_common_block_hook(slesp, ksp);
1037
1038 /* Predefined settings when using AMG as a preconditioner */
1039 PC pc;
1040 KSPGetPC(ksp, &pc);
1041 PCSetUp(pc);
1042
1043 PC _pc;
1044 PetscInt n_split;
1045 KSP *xyz_subksp;
1046 const char xyz[3] = "xyz";
1047 bool is_symm = _system_should_be_sym(slesp->solver);
1048
1049 PCFieldSplitGetSubKSP(pc, &n_split, &xyz_subksp);
1050 assert(n_split == 3);
1051
1052 for (PetscInt id = 0; id < n_split; id++) {
1053
1054 sprintf(prefix, "%s_fieldsplit_%c", slesp->name, xyz[id]);
1055
1056 _petsc_cmd(true, prefix, "ksp_type","preonly");
1057
1058 /* Predefined settings when using AMG as a preconditioner */
1059 KSP _ksp = xyz_subksp[id];
1060 KSPGetPC(_ksp, &_pc);
1061
1062 _petsc_pchypre_hook(prefix, slesp, is_symm, _pc);
1063
1064 PCSetFromOptions(_pc);
1065 KSPSetFromOptions(_ksp);
1066
1067 } /* Loop on block settings */
1068
1069 BFT_FREE(prefix);
1070 PetscFree(xyz_subksp);
1071
1072 /* User function for additional settings */
1073 cs_user_sles_petsc_hook(context, ksp);
1074
1075 PCSetFromOptions(pc);
1076 KSPSetFromOptions(ksp);
1077
1078 /* Dump the setup related to PETSc in a specific file */
1079 if (!slesp->setup_done) {
1080 KSPSetUp(ksp);
1081 cs_sles_petsc_log_setup(ksp);
1082 slesp->setup_done = true;
1083 }
1084
1085 cs_fp_exception_restore_trap(); /* Avoid trouble with a too restrictive
1086 SIGFPE detection */
1087 }
1088 #endif
1089
1090 /*----------------------------------------------------------------------------*/
1091 /*!
1092 * \brief Function pointer: setup hook for setting PETSc solver and
1093 * preconditioner.
1094 * Case of block preconditioner
1095 *
1096 * \param[in, out] context pointer to optional (untyped) value or structure
1097 * \param[in, out] ksp pointer to PETSc KSP context
1098 */
1099 /*----------------------------------------------------------------------------*/
1100
1101 static void
_petsc_block_hook(void * context,KSP ksp)1102 _petsc_block_hook(void *context,
1103 KSP ksp)
1104 {
1105 cs_param_sles_t *slesp = (cs_param_sles_t *)context;
1106
1107 cs_fp_exception_disable_trap(); /* Avoid trouble with a too restrictive
1108 SIGFPE detection */
1109
1110 /* prefix will be extended with the fieldsplit */
1111 int len = strlen(slesp->name) + 1;
1112 int _len = len + strlen("_fieldsplit_x_") + 1;
1113 char *prefix = NULL;
1114 BFT_MALLOC(prefix, _len + 1, char);
1115 sprintf(prefix, "%s_", slesp->name);
1116 prefix[len] = '\0';
1117 KSPSetOptionsPrefix(ksp, prefix);
1118
1119 /* Set the solver (tolerance and max_it too) */
1120 _petsc_set_krylov_solver(slesp, ksp);
1121
1122 /* Common settings to block preconditionner */
1123 _petsc_common_block_hook(slesp, ksp);
1124
1125 PC pc;
1126 KSPGetPC(ksp, &pc);
1127 PCSetUp(pc);
1128
1129 PC _pc;
1130 KSP *xyz_subksp;
1131 PetscInt n_split;
1132 const char xyz[3] = "xyz";
1133
1134 PCFieldSplitGetSubKSP(pc, &n_split, &xyz_subksp);
1135 assert(n_split == 3);
1136
1137 for (PetscInt id = 0; id < n_split; id++) {
1138
1139 sprintf(prefix, "%s_fieldsplit_%c", slesp->name, xyz[id]);
1140
1141 switch (slesp->precond) {
1142
1143 case CS_PARAM_PRECOND_NONE:
1144 _petsc_cmd(true, prefix, "ksp_type", "richardson");
1145 break;
1146
1147 case CS_PARAM_PRECOND_DIAG:
1148 _petsc_cmd(true, prefix, "ksp_type", "richardson");
1149 _petsc_cmd(true, prefix, "pc_type", "jacobi");
1150 break;
1151
1152 case CS_PARAM_PRECOND_ILU0:
1153 case CS_PARAM_PRECOND_BJACOB_ILU0:
1154 if (slesp->solver_class == CS_PARAM_SLES_CLASS_HYPRE) {
1155 #if defined(PETSC_HAVE_HYPRE)
1156 _petsc_cmd(true, prefix, "ksp_type", "preonly");
1157 _petsc_cmd(true, prefix, "pc_type", "hypre");
1158 _petsc_cmd(true, prefix, "pc_hypre_type", "euclid");
1159 _petsc_cmd(true, prefix, "pc_hypre_euclid_level", "0");
1160 #else
1161 bft_error(__FILE__, __LINE__, 0,
1162 " %s: Invalid option: HYPRE is not installed.", __func__);
1163 #endif
1164 }
1165 else {
1166
1167 _petsc_cmd(true, prefix, "ksp_type", "richardson");
1168 _petsc_bilu0_hook(prefix);
1169
1170 }
1171 break;
1172
1173 case CS_PARAM_PRECOND_ICC0:
1174 _petsc_cmd(true, prefix, "ksp_type", "richardson");
1175 _petsc_bicc0_hook(prefix);
1176 break;
1177
1178 case CS_PARAM_PRECOND_LU:
1179 _petsc_cmd(true, prefix, "ksp_type", "preonly");
1180 #if defined(PETSC_HAVE_MUMPS)
1181 _petsc_cmd(true, prefix, "pc_type", "lu");
1182 _petsc_cmd(true, prefix, "pc_factor_mat_solver_type", "mumps");
1183 #else
1184 if (cs_glob_n_ranks == 1)
1185 _petsc_cmd(true, prefix, "pc_type", "lu");
1186 else { /* Switch to a block version (sequential in each block) */
1187 _petsc_cmd(true, prefix, "pc_type", "bjacobi");
1188 _petsc_cmd(true, prefix, "pc_jacobi_blocks", "1");
1189 _petsc_cmd(true, prefix, "sub_ksp_type", "preonly");
1190 _petsc_cmd(true, prefix, "sub_pc_type", "lu");
1191 }
1192 #endif
1193 break;
1194
1195 case CS_PARAM_PRECOND_SSOR:
1196 case CS_PARAM_PRECOND_BJACOB_SGS:
1197 _petsc_cmd(true, prefix, "ksp_type", "richardson");
1198 _petsc_bssor_hook(prefix);
1199 break;
1200
1201 default:
1202 bft_error(__FILE__, __LINE__, 0,
1203 " %s: Eq. %s: Invalid preconditioner.",__func__, slesp->name);
1204 break;
1205
1206 } /* Switch on preconditioning type */
1207
1208 KSP _ksp = xyz_subksp[id];
1209 KSPGetPC(_ksp, &_pc);
1210 PCSetFromOptions(_pc);
1211 KSPSetUp(_ksp);
1212
1213 } /* Loop on block settings */
1214
1215 BFT_FREE(prefix);
1216 PetscFree(xyz_subksp);
1217
1218 /* User function for additional settings */
1219 cs_user_sles_petsc_hook(context, ksp);
1220
1221 PCSetFromOptions(pc);
1222 KSPSetFromOptions(ksp);
1223
1224 /* Dump the setup related to PETSc in a specific file */
1225 if (!slesp->setup_done) {
1226 KSPSetUp(ksp);
1227 cs_sles_petsc_log_setup(ksp);
1228 slesp->setup_done = true;
1229 }
1230
1231 cs_fp_exception_restore_trap(); /* Avoid trouble with a too restrictive
1232 SIGFPE detection */
1233 }
1234 #endif /* defined(HAVE_PETSC) */
1235
1236 /*----------------------------------------------------------------------------*/
1237 /*!
1238 * \brief Check if the settings are consitent. Can apply minor modifications.
1239 *
1240 * \param[in, out] slesp pointer to a \ref cs_param_sles_t structure
1241 */
1242 /*----------------------------------------------------------------------------*/
1243
1244 static void
_check_settings(cs_param_sles_t * slesp)1245 _check_settings(cs_param_sles_t *slesp)
1246 {
1247
1248 /* Checks related to MUMPS */
1249
1250 if (_mumps_is_needed(slesp->solver)) {
1251
1252 cs_param_sles_class_t ret_class =
1253 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_MUMPS);
1254 if (ret_class == CS_PARAM_SLES_N_CLASSES)
1255 bft_error(__FILE__, __LINE__, 0,
1256 " %s: Error detected while setting the SLES \"%s\"\n"
1257 " MUMPS is not available with your installation.\n"
1258 " Please check your installation settings.\n",
1259 __func__, slesp->name);
1260 else
1261 slesp->solver_class = ret_class;
1262
1263 }
1264 else {
1265 if (slesp->solver_class == CS_PARAM_SLES_CLASS_MUMPS)
1266 bft_error(__FILE__, __LINE__, 0,
1267 " %s: Error detected while setting the SLES \"%s\"\n"
1268 " MUMPS class is not consistent with your settings.\n"
1269 " Please check your installation settings.\n",
1270 __func__, slesp->name);
1271 }
1272
1273 /* Checks related to GCR/GMRES algorithms */
1274
1275 if (slesp->solver == CS_PARAM_ITSOL_GMRES ||
1276 slesp->solver == CS_PARAM_ITSOL_GCR)
1277 if (slesp->restart < 2)
1278 bft_error(__FILE__, __LINE__, 0,
1279 " %s: Error detected while setting the SLES \"%s\"\n"
1280 " The restart interval (=%d) is not big enough.\n"
1281 " Please check your installation settings.\n",
1282 __func__, slesp->name, slesp->restart);
1283
1284 }
1285
1286 /*----------------------------------------------------------------------------*/
1287 /*!
1288 * \brief Set parameters for initializing SLES structures used for the
1289 * resolution of the linear system.
1290 * Case of saturne's own solvers.
1291 *
1292 * \param[in] use_field_id if false use system name
1293 * \param[in, out] slesp pointer to a \ref cs_param_sles_t structure
1294 */
1295 /*----------------------------------------------------------------------------*/
1296
1297 static void
_set_saturne_sles(bool use_field_id,cs_param_sles_t * slesp)1298 _set_saturne_sles(bool use_field_id,
1299 cs_param_sles_t *slesp)
1300 {
1301 assert(slesp != NULL); /* Sanity checks */
1302
1303 const char *sles_name = use_field_id ? NULL : slesp->name;
1304 assert(slesp->field_id > -1 || sles_name != NULL);
1305
1306 /* 1- Define the preconditioner */
1307 /* ========================= */
1308
1309 int poly_degree;
1310 cs_sles_pc_t *pc = NULL;
1311
1312 switch (slesp->precond) {
1313
1314 case CS_PARAM_PRECOND_DIAG:
1315 poly_degree = 0;
1316 break;
1317
1318 case CS_PARAM_PRECOND_POLY1:
1319 poly_degree = 1;
1320 break;
1321
1322 case CS_PARAM_PRECOND_POLY2:
1323 poly_degree = 2;
1324 break;
1325
1326 case CS_PARAM_PRECOND_AMG:
1327 poly_degree = -1;
1328 switch (slesp->amg_type) {
1329
1330 case CS_PARAM_AMG_HOUSE_V:
1331 pc = cs_multigrid_pc_create(CS_MULTIGRID_V_CYCLE);
1332 break;
1333 case CS_PARAM_AMG_HOUSE_K:
1334 if (slesp->solver == CS_PARAM_ITSOL_CG)
1335 slesp->solver = CS_PARAM_ITSOL_FCG;
1336 pc = cs_multigrid_pc_create(CS_MULTIGRID_K_CYCLE);
1337 break;
1338
1339 default:
1340 bft_error(__FILE__, __LINE__, 0,
1341 " %s: System: %s; Invalid AMG type with Code_Saturne solvers.",
1342 __func__, slesp->name);
1343 break;
1344
1345 } /* Switch on the type of AMG */
1346 break; /* AMG as preconditioner */
1347
1348 case CS_PARAM_PRECOND_GKB_CG:
1349 case CS_PARAM_PRECOND_GKB_GMRES:
1350 poly_degree = -1;
1351 break;
1352
1353 case CS_PARAM_PRECOND_NONE:
1354 default:
1355 poly_degree = -1; /* None or other */
1356
1357 } /* Switch on the preconditioner type */
1358
1359 /* 2- Define the iterative solver */
1360 /* =========================== */
1361
1362 cs_sles_it_t *it = NULL;
1363 cs_multigrid_t *mg = NULL;
1364
1365 switch (slesp->solver) {
1366
1367 case CS_PARAM_ITSOL_AMG:
1368 {
1369 switch (slesp->amg_type) {
1370
1371 case CS_PARAM_AMG_HOUSE_V:
1372 mg = cs_multigrid_define(slesp->field_id,
1373 sles_name,
1374 CS_MULTIGRID_V_CYCLE);
1375
1376 /* Advanced setup (default is specified inside the brackets)
1377 * for AMG as solver */
1378 cs_multigrid_set_solver_options
1379 (mg,
1380 CS_SLES_JACOBI, /* descent smoother type (CS_SLES_PCG) */
1381 CS_SLES_JACOBI, /* ascent smoother type (CS_SLES_PCG) */
1382 CS_SLES_PCG, /* coarse solver type (CS_SLES_PCG) */
1383 slesp->n_max_iter, /* n max cycles (100) */
1384 5, /* n max iter for descent (10) */
1385 5, /* n max iter for ascent (10) */
1386 1000, /* n max iter coarse solver (10000) */
1387 0, /* polynomial precond. degree descent (0) */
1388 0, /* polynomial precond. degree ascent (0) */
1389 -1, /* polynomial precond. degree coarse (0) */
1390 1.0, /* precision multiplier descent (< 0 forces max iters) */
1391 1.0, /* precision multiplier ascent (< 0 forces max iters) */
1392 1); /* requested precision multiplier coarse (default 1) */
1393 break;
1394
1395 case CS_PARAM_AMG_HOUSE_K:
1396 mg = cs_multigrid_define(slesp->field_id,
1397 sles_name,
1398 CS_MULTIGRID_K_CYCLE);
1399
1400 cs_multigrid_set_solver_options
1401 (mg,
1402 CS_SLES_P_SYM_GAUSS_SEIDEL, /* descent smoother */
1403 CS_SLES_P_SYM_GAUSS_SEIDEL, /* ascent smoother */
1404 CS_SLES_PCG, /* coarse smoother */
1405 slesp->n_max_iter, /* n_max_cycles */
1406 1, /* n_max_iter_descent, */
1407 1, /* n_max_iter_ascent */
1408 100, /* n_max_iter_coarse */
1409 0, /* poly_degree_descent */
1410 0, /* poly_degree_ascent */
1411 0, /* poly_degree_coarse */
1412 -1.0, /* precision_mult_descent */
1413 -1.0, /* precision_mult_ascent */
1414 1); /* precision_mult_coarse */
1415 break;
1416
1417 default:
1418 bft_error(__FILE__, __LINE__, 0,
1419 " %s; System: %s -- Invalid AMG type with Code_Saturne"
1420 " solvers.", __func__, slesp->name);
1421 break;
1422
1423 } /* End of switch on the AMG type */
1424
1425 }
1426 break; /* AMG as solver */
1427
1428 case CS_PARAM_ITSOL_BICG:
1429 it = cs_sles_it_define(slesp->field_id,
1430 sles_name,
1431 CS_SLES_BICGSTAB,
1432 poly_degree,
1433 slesp->n_max_iter);
1434 break;
1435
1436 case CS_PARAM_ITSOL_BICGSTAB2:
1437 it = cs_sles_it_define(slesp->field_id,
1438 sles_name,
1439 CS_SLES_BICGSTAB2,
1440 poly_degree,
1441 slesp->n_max_iter);
1442 break;
1443
1444 case CS_PARAM_ITSOL_CG:
1445 it = cs_sles_it_define(slesp->field_id,
1446 sles_name,
1447 CS_SLES_PCG,
1448 poly_degree,
1449 slesp->n_max_iter);
1450 break;
1451
1452 case CS_PARAM_ITSOL_CR3:
1453 it = cs_sles_it_define(slesp->field_id,
1454 sles_name,
1455 CS_SLES_PCR3,
1456 poly_degree,
1457 slesp->n_max_iter);
1458 break;
1459
1460 case CS_PARAM_ITSOL_FCG:
1461 it = cs_sles_it_define(slesp->field_id,
1462 sles_name,
1463 CS_SLES_IPCG,
1464 poly_degree,
1465 slesp->n_max_iter);
1466 break;
1467
1468 case CS_PARAM_ITSOL_GAUSS_SEIDEL:
1469 it = cs_sles_it_define(slesp->field_id,
1470 sles_name,
1471 CS_SLES_P_GAUSS_SEIDEL,
1472 -1, /* Not useful to apply a preconditioner */
1473 slesp->n_max_iter);
1474 break;
1475
1476 case CS_PARAM_ITSOL_GCR:
1477 it = cs_sles_it_define(slesp->field_id,
1478 sles_name,
1479 CS_SLES_GCR,
1480 poly_degree,
1481 slesp->n_max_iter);
1482 break;
1483
1484 case CS_PARAM_ITSOL_GKB_CG:
1485 it = cs_sles_it_define(slesp->field_id,
1486 sles_name,
1487 CS_SLES_IPCG, /* Flexible CG */
1488 poly_degree,
1489 slesp->n_max_iter);
1490 break;
1491
1492 case CS_PARAM_ITSOL_GKB_GMRES:
1493 it = cs_sles_it_define(slesp->field_id,
1494 sles_name,
1495 CS_SLES_GMRES, /* Should be a flexible GMRES */
1496 poly_degree,
1497 slesp->n_max_iter);
1498 break;
1499
1500 case CS_PARAM_ITSOL_GMRES:
1501 it = cs_sles_it_define(slesp->field_id,
1502 sles_name,
1503 CS_SLES_GMRES,
1504 poly_degree,
1505 slesp->n_max_iter);
1506 break;
1507
1508 case CS_PARAM_ITSOL_JACOBI:
1509 it = cs_sles_it_define(slesp->field_id,
1510 sles_name,
1511 CS_SLES_JACOBI,
1512 -1, /* Not useful to apply a preconditioner */
1513 slesp->n_max_iter);
1514 break;
1515
1516 case CS_PARAM_ITSOL_SYM_GAUSS_SEIDEL:
1517 it = cs_sles_it_define(slesp->field_id,
1518 sles_name,
1519 CS_SLES_P_SYM_GAUSS_SEIDEL,
1520 -1, /* Not useful to apply a preconditioner */
1521 slesp->n_max_iter);
1522 break;
1523
1524 case CS_PARAM_ITSOL_USER_DEFINED:
1525 it = cs_sles_it_define(slesp->field_id,
1526 sles_name,
1527 CS_SLES_USER_DEFINED,
1528 poly_degree,
1529 slesp->n_max_iter);
1530 break;
1531
1532 default:
1533 bft_error(__FILE__, __LINE__, 0,
1534 " %s: Invalid iterative solver for solving equation %s.\n"
1535 " Please modify your settings.", __func__, slesp->name);
1536 break;
1537
1538 } /* end of switch */
1539
1540 /* Update the preconditioner settings if needed */
1541 if (slesp->precond == CS_PARAM_PRECOND_AMG) {
1542
1543 assert(pc != NULL && it != NULL);
1544
1545 mg = cs_sles_pc_get_context(pc);
1546 cs_sles_it_transfer_pc(it, &pc);
1547
1548 /* If this is a K-cycle multigrid. Change the default settings when used as
1549 preconditioner */
1550 if (slesp->amg_type == CS_PARAM_AMG_HOUSE_K) {
1551
1552 cs_multigrid_set_solver_options
1553 (mg,
1554 CS_SLES_PCG, /* descent smoother */
1555 CS_SLES_PCG, /* ascent smoother */
1556 CS_SLES_PCG, /* coarse solver */
1557 slesp->n_max_iter, /* n_max_cycles */
1558 2, /* n_max_iter_descent, */
1559 2, /* n_max_iter_ascent */
1560 500, /* n_max_iter_coarse */
1561 0, /* poly_degree_descent */
1562 0, /* poly_degree_ascent */
1563 0, /* poly_degree_coarse */
1564 -1.0, /* precision_mult_descent */
1565 -1.0, /* precision_mult_ascent */
1566 1.0); /* precision_mult_coarse */
1567
1568 cs_multigrid_set_coarsening_options(mg,
1569 8, /* aggregation_limit*/
1570 CS_GRID_COARSENING_SPD_PW,
1571 10, /* n_max_levels */
1572 50, /* min_g_cells */
1573 0., /* P0P1 relaxation */
1574 0); /* postprocess */
1575
1576 } /* AMG K-cycle as preconditioner */
1577
1578 } /* AMG preconditioner */
1579
1580 /* Define the level of verbosity for SLES structure */
1581 if (slesp->verbosity > 3) {
1582
1583 cs_sles_t *sles = cs_sles_find_or_add(slesp->field_id, sles_name);
1584 cs_sles_it_t *sles_it = (cs_sles_it_t *)cs_sles_get_context(sles);
1585
1586 /* true = use_iteration instead of wall clock time */
1587 cs_sles_it_set_plot_options(sles_it, slesp->name, true);
1588
1589 }
1590
1591 }
1592
1593 /*----------------------------------------------------------------------------*/
1594 /*!
1595 * \brief Set parameters for initializing SLES structures used for the
1596 * resolution of the linear system.
1597 * Case of MUMPS's own solvers.
1598 *
1599 * \param[in] use_field_id if false use system name
1600 * \param[in, out] slesp pointer to a \ref cs_param_sles_t structure
1601 */
1602 /*----------------------------------------------------------------------------*/
1603
1604 static void
_set_mumps_sles(bool use_field_id,cs_param_sles_t * slesp)1605 _set_mumps_sles(bool use_field_id,
1606 cs_param_sles_t *slesp)
1607 {
1608 assert(slesp != NULL); /* Sanity checks */
1609
1610 const char *sles_name = use_field_id ? NULL : slesp->name;
1611 assert(slesp->field_id > -1 || sles_name != NULL);
1612
1613 #if defined(HAVE_MUMPS)
1614 cs_sles_mumps_define(slesp->field_id,
1615 sles_name,
1616 slesp,
1617 cs_user_sles_mumps_hook,
1618 NULL);
1619 #else
1620 bft_error(__FILE__, __LINE__, 0,
1621 "%s: System: %s\n"
1622 " MUMPS is not supported directly.\n"
1623 " Please check your settings or your code_saturne installation.",
1624 __func__, slesp->name);
1625 #endif
1626 }
1627
1628 /*----------------------------------------------------------------------------*/
1629 /*!
1630 * \brief Set parameters for initializing SLES structures used for the
1631 * resolution of the linear system.
1632 * Case of PETSc and Hypre families of solvers.
1633 *
1634 * \param[in] use_field_id if false use system name
1635 * \param[in, out] slesp pointer to a \ref cs_param_sles_t structure
1636 */
1637 /*----------------------------------------------------------------------------*/
1638
1639 static void
_set_petsc_hypre_sles(bool use_field_id,cs_param_sles_t * slesp)1640 _set_petsc_hypre_sles(bool use_field_id,
1641 cs_param_sles_t *slesp)
1642 {
1643 assert(slesp != NULL); /* Sanity checks */
1644
1645 const char *sles_name = use_field_id ? NULL : slesp->name;
1646 assert(slesp->field_id > -1 || sles_name != NULL);
1647
1648 #if defined(HAVE_PETSC)
1649 cs_sles_petsc_init();
1650
1651 if (slesp->pcd_block_type != CS_PARAM_PRECOND_BLOCK_NONE) {
1652
1653 if (slesp->precond == CS_PARAM_PRECOND_AMG) {
1654
1655 if (slesp->amg_type == CS_PARAM_AMG_PETSC_GAMG_V ||
1656 slesp->amg_type == CS_PARAM_AMG_PETSC_GAMG_W)
1657 cs_sles_petsc_define(slesp->field_id,
1658 sles_name,
1659 MATMPIAIJ,
1660 _petsc_amg_block_gamg_hook,
1661 (void *)slesp);
1662
1663 else if (slesp->amg_type == CS_PARAM_AMG_HYPRE_BOOMER_V ||
1664 slesp->amg_type == CS_PARAM_AMG_HYPRE_BOOMER_W) {
1665 #if defined(PETSC_HAVE_HYPRE)
1666 cs_sles_petsc_define(slesp->field_id,
1667 sles_name,
1668 MATMPIAIJ,
1669 _petsc_amg_block_boomer_hook,
1670 (void *)slesp);
1671 #else
1672 cs_base_warn(__FILE__, __LINE__);
1673 cs_log_printf(CS_LOG_DEFAULT,
1674 " %s: System: %s.\n"
1675 " Boomer is not available. Switch to GAMG solver.",
1676 __func__, slesp->name);
1677 cs_sles_petsc_define(slesp->field_id,
1678 sles_name,
1679 MATMPIAIJ,
1680 _petsc_amg_block_gamg_hook,
1681 (void *)slesp);
1682 #endif
1683 }
1684 else
1685 bft_error(__FILE__, __LINE__, 0,
1686 " %s: System: %s\n"
1687 " No AMG solver available for a block-AMG.",
1688 __func__, slesp->name);
1689 }
1690 else {
1691
1692 cs_sles_petsc_define(slesp->field_id,
1693 sles_name,
1694 MATMPIAIJ,
1695 _petsc_block_hook,
1696 (void *)slesp);
1697
1698 }
1699
1700 }
1701 else { /* No block preconditioner */
1702
1703 cs_sles_petsc_define(slesp->field_id,
1704 sles_name,
1705 MATMPIAIJ,
1706 _petsc_setup_hook,
1707 (void *)slesp);
1708
1709 }
1710 #else
1711 bft_error(__FILE__, __LINE__, 0,
1712 _(" %s: PETSC algorithms used to solve %s are not linked.\n"
1713 " Please install Code_Saturne with PETSc."),
1714 __func__, slesp->name);
1715 #endif /* HAVE_PETSC */
1716 }
1717
1718 /*============================================================================
1719 * Public function prototypes
1720 *============================================================================*/
1721
1722 /*----------------------------------------------------------------------------*/
1723 /*!
1724 * \brief Create a \ref cs_param_sles_t structure and assign a default
1725 * settings
1726 *
1727 * \param[in] field_id id related to to the variable field or -1
1728 * \param[in] system_name name of the system to solve or NULL
1729 *
1730 * \return a pointer to a cs_param_sles_t stucture
1731 */
1732 /*----------------------------------------------------------------------------*/
1733
1734 cs_param_sles_t *
cs_param_sles_create(int field_id,const char * system_name)1735 cs_param_sles_create(int field_id,
1736 const char *system_name)
1737 {
1738 cs_param_sles_t *slesp = NULL;
1739
1740 BFT_MALLOC(slesp, 1, cs_param_sles_t);
1741
1742 slesp->verbosity = 0; /* SLES verbosity */
1743
1744 slesp->field_id = field_id; /* associated field id */
1745
1746 slesp->solver_class = CS_PARAM_SLES_CLASS_CS; /* solver family */
1747
1748 slesp->precond = CS_PARAM_PRECOND_POLY1; /* preconditioner */
1749 slesp->solver = CS_PARAM_ITSOL_GCR; /* iterative solver */
1750 slesp->amg_type = CS_PARAM_AMG_NONE; /* no predefined AMG type */
1751 slesp->pcd_block_type = CS_PARAM_PRECOND_BLOCK_NONE; /* no block by default */
1752
1753 slesp->restart = 15; /* max. iter. before restarting */
1754 slesp->n_max_iter = 10000; /* max. number of iterations */
1755 slesp->eps = 1e-6; /* relative tolerance to stop
1756 an iterative solver */
1757
1758 slesp->resnorm_type = CS_PARAM_RESNORM_FILTERED_RHS;
1759 slesp->setup_done = false;
1760
1761 slesp->name = NULL;
1762 if (system_name != NULL) {
1763 size_t len = strlen(system_name);
1764 BFT_MALLOC(slesp->name, len + 1, char);
1765 strncpy(slesp->name, system_name, len);
1766 slesp->name[len] = '\0';
1767 }
1768
1769 return slesp;
1770 }
1771
1772 /*----------------------------------------------------------------------------*/
1773 /*!
1774 * \brief Free a \ref cs_param_sles_t structure
1775 *
1776 * \param[in, out] slesp pointer to a \cs_param_sles_t structure to free
1777 */
1778 /*----------------------------------------------------------------------------*/
1779
1780 void
cs_param_sles_free(cs_param_sles_t ** p_slesp)1781 cs_param_sles_free(cs_param_sles_t **p_slesp)
1782 {
1783 if (p_slesp == NULL)
1784 return;
1785
1786 cs_param_sles_t *slesp = *p_slesp;
1787
1788 if (slesp == NULL)
1789 return;
1790
1791 BFT_FREE(slesp->name);
1792 BFT_FREE(slesp);
1793 slesp = NULL;
1794 }
1795
1796 /*----------------------------------------------------------------------------*/
1797 /*!
1798 * \brief Log information related to the linear settings stored in the
1799 * structure
1800 *
1801 * \param[in] slesp pointer to a \ref cs_param_sles_log
1802 */
1803 /*----------------------------------------------------------------------------*/
1804
1805 void
cs_param_sles_log(cs_param_sles_t * slesp)1806 cs_param_sles_log(cs_param_sles_t *slesp)
1807 {
1808 if (slesp == NULL)
1809 return;
1810
1811 cs_log_printf(CS_LOG_SETUP, "\n### %s | Linear algebra settings\n",
1812 slesp->name);
1813 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Family:", slesp->name);
1814 if (slesp->solver_class == CS_PARAM_SLES_CLASS_CS)
1815 cs_log_printf(CS_LOG_SETUP, " Code_Saturne\n");
1816 else if (slesp->solver_class == CS_PARAM_SLES_CLASS_MUMPS)
1817 cs_log_printf(CS_LOG_SETUP, " MUMPS\n");
1818 else if (slesp->solver_class == CS_PARAM_SLES_CLASS_HYPRE)
1819 cs_log_printf(CS_LOG_SETUP, " HYPRE\n");
1820 else if (slesp->solver_class == CS_PARAM_SLES_CLASS_PETSC)
1821 cs_log_printf(CS_LOG_SETUP, " PETSc\n");
1822
1823 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Verbosity: %d\n",
1824 slesp->name, slesp->verbosity);
1825 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Field id: %d\n",
1826 slesp->name, slesp->field_id);
1827
1828 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Solver.Name: %s\n",
1829 slesp->name, cs_param_get_solver_name(slesp->solver));
1830 if (slesp->solver == CS_PARAM_ITSOL_AMG)
1831 cs_log_printf(CS_LOG_SETUP, " * %s | SLES AMG.Type: %s\n",
1832 slesp->name, cs_param_get_amg_type_name(slesp->amg_type));
1833
1834 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Solver.Precond: %s\n",
1835 slesp->name, cs_param_get_precond_name(slesp->precond));
1836 if (slesp->precond == CS_PARAM_PRECOND_AMG)
1837 cs_log_printf(CS_LOG_SETUP, " * %s | SLES AMG.Type: %s\n",
1838 slesp->name, cs_param_get_amg_type_name(slesp->amg_type));
1839 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Block.Precond: %s\n",
1840 slesp->name,
1841 cs_param_get_precond_block_name(slesp->pcd_block_type));
1842
1843 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Solver.MaxIter: %d\n",
1844 slesp->name, slesp->n_max_iter);
1845 if (slesp->solver == CS_PARAM_ITSOL_GMRES ||
1846 slesp->solver == CS_PARAM_ITSOL_FGMRES ||
1847 slesp->solver == CS_PARAM_ITSOL_GCR)
1848 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Solver.Restart: %d\n",
1849 slesp->name, slesp->restart);
1850
1851 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Solver.Eps: % -10.6e\n",
1852 slesp->name, slesp->eps);
1853
1854 cs_log_printf(CS_LOG_SETUP, " * %s | SLES Normalization: ",
1855 slesp->name);
1856 switch (slesp->resnorm_type) {
1857 case CS_PARAM_RESNORM_NORM2_RHS:
1858 cs_log_printf(CS_LOG_SETUP, "Euclidean norm of the RHS\n");
1859 break;
1860 case CS_PARAM_RESNORM_WEIGHTED_RHS:
1861 cs_log_printf(CS_LOG_SETUP, "Weighted Euclidean norm of the RHS\n");
1862 break;
1863 case CS_PARAM_RESNORM_FILTERED_RHS:
1864 cs_log_printf(CS_LOG_SETUP, "Filtered Euclidean norm of the RHS\n");
1865 break;
1866 case CS_PARAM_RESNORM_NONE:
1867 default:
1868 cs_log_printf(CS_LOG_SETUP, "None\n");
1869 break;
1870 }
1871 cs_log_printf(CS_LOG_SETUP, "\n");
1872 }
1873
1874 /*----------------------------------------------------------------------------*/
1875 /*!
1876 * \brief Copy a cs_param_sles_t structure from src to dst
1877 *
1878 * \param[in] src reference cs_param_sles_t structure to copy
1879 * \param[in, out] dst copy of the reference at exit
1880 */
1881 /*----------------------------------------------------------------------------*/
1882
1883 void
cs_param_sles_copy_from(cs_param_sles_t * src,cs_param_sles_t * dst)1884 cs_param_sles_copy_from(cs_param_sles_t *src,
1885 cs_param_sles_t *dst)
1886 {
1887 if (dst == NULL)
1888 return;
1889
1890 /* Remark: name is managed at the creation of the structure */
1891
1892 dst->setup_done = src->setup_done;
1893 dst->verbosity = src->verbosity;
1894 dst->field_id = src->field_id;
1895
1896 dst->solver_class = src->solver_class;
1897 dst->precond = src->precond;
1898 dst->solver = src->solver;
1899 dst->amg_type = src->amg_type;
1900 dst->pcd_block_type = src->pcd_block_type;
1901
1902 dst->resnorm_type = src->resnorm_type;
1903 dst->n_max_iter = src->n_max_iter;
1904 dst->eps = src->eps;
1905 }
1906
1907 /*----------------------------------------------------------------------------*/
1908 /*!
1909 * \brief Define cs_sles_t structure in accordance with the settings of a
1910 * cs_param_sles_t structure (SLES = Sparse Linear Equation Solver)
1911 *
1912 * \param[in] use_field_id if false use system name to define a SLES
1913 * \param[in, out] slesp pointer to a cs_param_sles_t structure
1914 *
1915 * \return an error code (-1 if a problem is encountered, 0 otherwise)
1916 */
1917 /*----------------------------------------------------------------------------*/
1918
1919 int
cs_param_sles_set(bool use_field_id,cs_param_sles_t * slesp)1920 cs_param_sles_set(bool use_field_id,
1921 cs_param_sles_t *slesp)
1922 {
1923 if (slesp == NULL)
1924 return 0;
1925
1926 _check_settings(slesp);
1927
1928 switch (slesp->solver_class) {
1929
1930 case CS_PARAM_SLES_CLASS_CS: /* Code_Saturne's own solvers */
1931 /* true = use field_id instead of slesp->name to set the sles */
1932 _set_saturne_sles(use_field_id, slesp);
1933 break;
1934
1935 case CS_PARAM_SLES_CLASS_MUMPS: /* MUMPS sparse direct solvers */
1936 /* true = use field_id instead of slesp->name to set the sles */
1937 _set_mumps_sles(use_field_id, slesp);
1938 break;
1939
1940 case CS_PARAM_SLES_CLASS_PETSC: /* PETSc solvers */
1941 case CS_PARAM_SLES_CLASS_HYPRE: /* HYPRE solvers through PETSc */
1942 /* true = use field_id instead of slesp->name to set the sles */
1943 _set_petsc_hypre_sles(use_field_id, slesp);
1944 break;
1945
1946 default:
1947 return -1;
1948
1949 } /* End of switch on class of solvers */
1950
1951 /* Define the level of verbosity for the SLES structure */
1952 if (slesp->verbosity > 1) {
1953
1954 /* All the previous SLES are defined thanks to the field_id */
1955 cs_sles_t *sles = NULL;
1956 if (use_field_id)
1957 sles = cs_sles_find_or_add(slesp->field_id, NULL);
1958 else
1959 sles = cs_sles_find_or_add(slesp->field_id, slesp->name);
1960
1961 /* Set verbosity */
1962 cs_sles_set_verbosity(sles, slesp->verbosity);
1963
1964 }
1965
1966 return 0;
1967 }
1968
1969 /*----------------------------------------------------------------------------*/
1970 /*!
1971 * \brief Check the availability of a solver library and return the requested
1972 * one if this is possible or an alternative or CS_PARAM_SLES_N_CLASSES
1973 * if no alternative is available.
1974 *
1975 * \param[in] wanted_class requested class of solvers
1976 *
1977 * \return the available solver class related to the requested class
1978 */
1979 /*----------------------------------------------------------------------------*/
1980
1981 cs_param_sles_class_t
cs_param_sles_check_class(cs_param_sles_class_t wanted_class)1982 cs_param_sles_check_class(cs_param_sles_class_t wanted_class)
1983 {
1984 switch (wanted_class) {
1985
1986 case CS_PARAM_SLES_CLASS_CS: /* No issue */
1987 return CS_PARAM_SLES_CLASS_CS;
1988
1989 case CS_PARAM_SLES_CLASS_HYPRE:
1990 #if defined(HAVE_PETSC)
1991 #if defined(PETSC_HAVE_HYPRE)
1992 return CS_PARAM_SLES_CLASS_HYPRE;
1993 #else
1994 cs_base_warn(__FILE__, __LINE__);
1995 bft_printf(" Switch to PETSc library since Hypre is not available");
1996 return CS_PARAM_SLES_CLASS_PETSC; /* Switch to petsc */
1997 #endif
1998 #else
1999 return CS_PARAM_SLES_N_CLASSES;
2000 #endif
2001
2002 case CS_PARAM_SLES_CLASS_PETSC:
2003 #if defined(HAVE_PETSC)
2004 return CS_PARAM_SLES_CLASS_PETSC;
2005 #else
2006 return CS_PARAM_SLES_N_CLASSES;
2007 #endif
2008
2009 case CS_PARAM_SLES_CLASS_MUMPS:
2010 #if defined(HAVE_MUMPS)
2011 return CS_PARAM_SLES_CLASS_MUMPS;
2012 #else
2013 #if defined(HAVE_PETSC)
2014 #if defined(PETSC_HAVE_MUMPS)
2015 cs_base_warn(__FILE__, __LINE__);
2016 bft_printf(" Switch to PETSc library since MUMPS is not available as"
2017 " a stand-alone library\n");
2018 return CS_PARAM_SLES_CLASS_PETSC;
2019 #else
2020 return CS_PARAM_SLES_N_CLASSES;
2021 #endif /* PETSC_HAVE_MUMPS */
2022 #else
2023 return CS_PARAM_SLES_N_CLASSES; /* PETSc without MUMPS */
2024 #endif /* HAVE_PETSC */
2025 return CS_PARAM_SLES_N_CLASSES; /* Neither MUMPS nor PETSc */
2026 #endif
2027
2028 default:
2029 return CS_PARAM_SLES_N_CLASSES;
2030 }
2031 }
2032
2033 #if defined(HAVE_PETSC)
2034 /*----------------------------------------------------------------------------*/
2035 /*!
2036 * \brief Set the command line option for PETSc
2037 *
2038 * \param[in] use_prefix need a prefix
2039 * \param[in] prefix optional prefix
2040 * \param[in] keyword command keyword
2041 * \param[in] keyval command value
2042 */
2043 /*----------------------------------------------------------------------------*/
2044
2045 void
cs_param_sles_petsc_cmd(bool use_prefix,const char * prefix,const char * keyword,const char * keyval)2046 cs_param_sles_petsc_cmd(bool use_prefix,
2047 const char *prefix,
2048 const char *keyword,
2049 const char *keyval)
2050 {
2051 _petsc_cmd(use_prefix, prefix, keyword, keyval);
2052 }
2053 #endif
2054
2055 /*----------------------------------------------------------------------------*/
2056
2057 END_C_DECLS
2058