1 #ifndef __CS_SLES_IT_H__ 2 #define __CS_SLES_IT_H__ 3 4 /*============================================================================ 5 * Sparse Linear Equation Solvers 6 *============================================================================*/ 7 8 /* 9 This file is part of Code_Saturne, a general-purpose CFD tool. 10 11 Copyright (C) 1998-2021 EDF S.A. 12 13 This program is free software; you can redistribute it and/or modify it under 14 the terms of the GNU General Public License as published by the Free Software 15 Foundation; either version 2 of the License, or (at your option) any later 16 version. 17 18 This program is distributed in the hope that it will be useful, but WITHOUT 19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 20 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 21 details. 22 23 You should have received a copy of the GNU General Public License along with 24 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 25 Street, Fifth Floor, Boston, MA 02110-1301, USA. 26 */ 27 28 /*----------------------------------------------------------------------------*/ 29 30 /*---------------------------------------------------------------------------- 31 * Local headers 32 *----------------------------------------------------------------------------*/ 33 34 #include "cs_base.h" 35 #include "cs_matrix.h" 36 #include "cs_sles.h" 37 #include "cs_sles_pc.h" 38 39 /*----------------------------------------------------------------------------*/ 40 41 BEGIN_C_DECLS 42 43 /*============================================================================ 44 * Macro definitions 45 *============================================================================*/ 46 47 /*============================================================================ 48 * Type definitions 49 *============================================================================*/ 50 51 /*---------------------------------------------------------------------------- 52 * Solver types 53 *----------------------------------------------------------------------------*/ 54 55 typedef enum { 56 57 CS_SLES_PCG, /*!< Preconditioned conjugate gradient */ 58 CS_SLES_FCG, /*!< Preconditions flexible conjugate gradient, 59 described in \cite Notay:2015 */ 60 CS_SLES_IPCG, /*!< Preconditions inexact conjugate gradient */ 61 CS_SLES_JACOBI, /*!< Jacobi */ 62 CS_SLES_BICGSTAB, /*!< Preconditioned BiCGstab 63 (biconjugate gradient stabilized) */ 64 CS_SLES_BICGSTAB2, /*!< Preconditioned BiCGstab2 */ 65 CS_SLES_GCR, /*!< Generalized conjugate residual */ 66 CS_SLES_GMRES, /*!< Preconditioned GMRES 67 (generalized minimal residual) */ 68 CS_SLES_P_GAUSS_SEIDEL, /*!< Process-local Gauss-Seidel */ 69 CS_SLES_P_SYM_GAUSS_SEIDEL, /*!< Process-local symmetric Gauss-Seidel */ 70 CS_SLES_PCR3, /*!< 3-layer conjugate residual */ 71 CS_SLES_USER_DEFINED, /*!< User-defined iterative solver */ 72 73 CS_SLES_N_IT_TYPES, /*!< Number of resolution algorithms 74 excluding smoother only*/ 75 76 CS_SLES_TS_F_GAUSS_SEIDEL, /*!< Truncated forward Gauss-Seidel smoother */ 77 CS_SLES_TS_B_GAUSS_SEIDEL, /*!< Truncated backward Gauss-Seidel smoother */ 78 79 CS_SLES_N_SMOOTHER_TYPES /*!< Number of resolution algorithms 80 including smoother only */ 81 82 } cs_sles_it_type_t; 83 84 /* Iterative linear solver context (opaque) */ 85 86 typedef struct _cs_sles_it_t cs_sles_it_t; 87 88 /* Forward type declarations */ 89 90 typedef struct _cs_sles_it_convergence_t cs_sles_it_convergence_t; 91 92 /*============================================================================ 93 * Global variables 94 *============================================================================*/ 95 96 /* Names for iterative solver types */ 97 98 extern const char *cs_sles_it_type_name[]; 99 100 /*============================================================================= 101 * User function prototypes 102 *============================================================================*/ 103 104 /*---------------------------------------------------------------------------- 105 * Solution of A.vx = Rhs using a user-defined iterative solver 106 * 107 * On entry, vx is considered initialized. 108 * 109 * parameters: 110 * c <-- pointer to solver context info 111 * a <-- matrix 112 * diag_block_size <-- diagonal block size (unused here) 113 * convergence <-- convergence information structure 114 * rhs <-- right hand side 115 * vx <-> system solution 116 * aux_size <-- number of elements in aux_vectors (in bytes) 117 * aux_vectors --- optional working area (allocation otherwise) 118 * 119 * returns: 120 * convergence state 121 *----------------------------------------------------------------------------*/ 122 123 cs_sles_convergence_state_t 124 cs_user_sles_it_solver(cs_sles_it_t *c, 125 const cs_matrix_t *a, 126 cs_lnum_t diag_block_size, 127 cs_sles_it_convergence_t *convergence, 128 const cs_real_t *rhs, 129 cs_real_t *restrict vx, 130 size_t aux_size, 131 void *aux_vectors); 132 133 /*============================================================================= 134 * Public function prototypes 135 *============================================================================*/ 136 137 /*---------------------------------------------------------------------------- 138 * Define and associate an iterative sparse linear system solver 139 * for a given field or equation name. 140 * 141 * If this system did not previously exist, it is added to the list of 142 * "known" systems. Otherwise, its definition is replaced by the one 143 * defined here. 144 * 145 * This is a utility function: if finer control is needed, see 146 * cs_sles_define() and cs_sles_it_create(). 147 * 148 * Note that this function returns a pointer directly to the iterative solver 149 * management structure. This may be used to set further options, 150 * for example using cs_sles_set_plot_options(). If needed, cs_sles_find() 151 * may be used to obtain a pointer to the matching cs_sles_t container. 152 * 153 * parameters: 154 * f_id <-- associated field id, or < 0 155 * name <-- associated name if f_id < 0, or NULL 156 * solver_type <-- type of solver (PCG, Jacobi, ...) 157 * poly_degree <-- preconditioning polynomial degree 158 * (0: diagonal; -1: non-preconditioned) 159 * n_max_iter <-- maximum number of iterations 160 * 161 * returns: 162 * pointer to newly created iterative solver info object. 163 *----------------------------------------------------------------------------*/ 164 165 cs_sles_it_t * 166 cs_sles_it_define(int f_id, 167 const char *name, 168 cs_sles_it_type_t solver_type, 169 int poly_degree, 170 int n_max_iter); 171 172 /*---------------------------------------------------------------------------- 173 * Create iterative sparse linear system solver info and context. 174 * 175 * parameters: 176 * solver_type <-- type of solver (PCG, Jacobi, ...) 177 * poly_degree <-- preconditioning polynomial degree 178 * (0: diagonal; -1: non-preconditioned) 179 * n_max_iter <-- maximum number of iterations 180 * update_stats <-- automatic solver statistics indicator 181 * 182 * returns: 183 * pointer to newly created solver info object. 184 *----------------------------------------------------------------------------*/ 185 186 cs_sles_it_t * 187 cs_sles_it_create(cs_sles_it_type_t solver_type, 188 int poly_degree, 189 int n_max_iter, 190 bool update_stats); 191 192 /*---------------------------------------------------------------------------- 193 * Destroy iterative sparse linear system solver info and context. 194 * 195 * parameters: 196 * context <-> pointer to iterative sparse linear solver info 197 * (actual type: cs_sles_it_t **) 198 *----------------------------------------------------------------------------*/ 199 200 void 201 cs_sles_it_destroy(void **context); 202 203 /*---------------------------------------------------------------------------- 204 * Create iterative sparse linear system solver info and context 205 * based on existing info and context. 206 * 207 * parameters: 208 * context <-- pointer to reference info and context 209 * (actual type: cs_sles_it_t *) 210 * 211 * returns: 212 * pointer to newly created solver info object 213 * (actual type: cs_sles_it_t *) 214 *----------------------------------------------------------------------------*/ 215 216 void * 217 cs_sles_it_copy(const void *context); 218 219 /*---------------------------------------------------------------------------- 220 * Setup iterative sparse linear equation solver. 221 * 222 * parameters: 223 * context <-> pointer to iterative sparse linear solver info 224 * (actual type: cs_sles_it_t *) 225 * name <-- pointer to system name 226 * a <-- associated matrix 227 * verbosity <-- verbosity level 228 *----------------------------------------------------------------------------*/ 229 230 void 231 cs_sles_it_setup(void *context, 232 const char *name, 233 const cs_matrix_t *a, 234 int verbosity); 235 236 /*---------------------------------------------------------------------------- 237 * Call iterative sparse linear equation solver. 238 * 239 * parameters: 240 * context <-> pointer to iterative sparse linear solver info 241 * (actual type: cs_sles_it_t *) 242 * name <-- pointer to system name 243 * a <-- matrix 244 * verbosity <-- verbosity level 245 * precision <-- solver precision 246 * r_norm <-- residue normalization 247 * n_iter --> number of iterations 248 * residue --> residue 249 * rhs <-- right hand side 250 * vx <-> system solution 251 * aux_size <-- number of elements in aux_vectors (in bytes) 252 * aux_vectors --- optional working area (internal allocation if NULL) 253 * 254 * returns: 255 * convergence state 256 *----------------------------------------------------------------------------*/ 257 258 cs_sles_convergence_state_t 259 cs_sles_it_solve(void *context, 260 const char *name, 261 const cs_matrix_t *a, 262 int verbosity, 263 double precision, 264 double r_norm, 265 int *n_iter, 266 double *residue, 267 const cs_real_t *rhs, 268 cs_real_t *vx, 269 size_t aux_size, 270 void *aux_vectors); 271 272 /*---------------------------------------------------------------------------- 273 * Free iterative sparse linear equation solver setup context. 274 * 275 * This function frees resolution-related data, such as 276 * buffers and preconditioning but does not free the whole context, 277 * as info used for logging (especially performance data) is maintained. 278 279 * parameters: 280 * context <-> pointer to iterative sparse linear solver info 281 * (actual type: cs_sles_it_t *) 282 *----------------------------------------------------------------------------*/ 283 284 void 285 cs_sles_it_free(void *context); 286 287 /*---------------------------------------------------------------------------- 288 * Log sparse linear equation solver info. 289 * 290 * parameters: 291 * context <-> pointer to iterative sparse linear solver info 292 * (actual type: cs_sles_it_t *) 293 * log_type <-- log type 294 *----------------------------------------------------------------------------*/ 295 296 void 297 cs_sles_it_log(const void *context, 298 cs_log_t log_type); 299 300 /*---------------------------------------------------------------------------- 301 * Return iterative solver type. 302 * 303 * parameters: 304 * context <-- pointer to iterative solver info and context 305 * 306 * returns: 307 * selected solver type 308 *----------------------------------------------------------------------------*/ 309 310 cs_sles_it_type_t 311 cs_sles_it_get_type(const cs_sles_it_t *context); 312 313 /*---------------------------------------------------------------------------- 314 * Return the initial residue for the previous solve operation with a solver. 315 * 316 * This is useful for convergence tests when this solver is used as 317 * a preconditioning smoother. 318 * 319 * This operation is only valid between calls to cs_sles_it_setup() 320 * (or cs_sles_it_solve()) and cs_sles_it_free(). 321 * It returns -1 otherwise. 322 * 323 * parameters: 324 * context <-- pointer to iterative solver info and context 325 * 326 * returns: 327 * initial residue from last call to \ref cs_sles_solve with this solver 328 *----------------------------------------------------------------------------*/ 329 330 double 331 cs_sles_it_get_last_initial_residue(const cs_sles_it_t *context); 332 333 /*---------------------------------------------------------------------------- 334 * Return a preconditioner context for an iterative sparse linear 335 * equation solver. 336 * 337 * This allows modifying parameters of a non default (Jacobi or polynomial) 338 * preconditioner. 339 * 340 * parameters: 341 * context <-- pointer to iterative solver info and context 342 * 343 * returns: 344 * pointer to preconditoner context 345 *----------------------------------------------------------------------------*/ 346 347 cs_sles_pc_t * 348 cs_sles_it_get_pc(cs_sles_it_t *context); 349 350 /*---------------------------------------------------------------------------- 351 * Assign a preconditioner to an iterative sparse linear equation 352 * solver, transfering its ownership to to solver context. 353 * 354 * This allows assigning a non default (Jacobi or polynomial) preconditioner. 355 * 356 * The input pointer is set to NULL to make it clear the caller does not 357 * own the preconditioner anymore, though the context can be accessed using 358 * cs_sles_it_get_cp(). 359 * 360 * parameters: 361 * context <-> pointer to iterative solver info and context 362 * pc <-> pointer to preconditoner context 363 *----------------------------------------------------------------------------*/ 364 365 void 366 cs_sles_it_transfer_pc(cs_sles_it_t *context, 367 cs_sles_pc_t **pc); 368 369 /*---------------------------------------------------------------------------- 370 * Copy options from one iterative sparse linear system solver info 371 * and context to another. 372 * 373 * Optional plotting contexts are shared between the source and destination 374 * contexts. 375 * 376 * Preconditioner settings are to be handled separately. 377 * 378 * parameters: 379 * src <-- pointer to source info and context 380 * dest <-> pointer to destination info and context 381 *----------------------------------------------------------------------------*/ 382 383 void 384 cs_sles_it_transfer_parameters(const cs_sles_it_t *src, 385 cs_sles_it_t *dest); 386 387 /*---------------------------------------------------------------------------- 388 * Associate a similar info and context object with which some setup 389 * data may be shared. 390 * 391 * This is especially useful for sharing preconditioning data between 392 * similar solver contexts (for example ascending and descending multigrid 393 * smoothers based on the same matrix). 394 * 395 * For preconditioning data to be effectively shared, cs_sles_it_setup() 396 * (or cs_sles_it_solve()) must be called on "shareable" before being 397 * called on "context" (without cs_sles_it_free() being called in between, 398 * of course). 399 * 400 * It is the caller's responsibility to ensure the context is not used 401 * for a cs_sles_it_setup() or cs_sles_it_solve() operation after the 402 * shareable object has been destroyed (normally by cs_sles_it_destroy()). 403 * 404 * parameters: 405 * context <-> pointer to iterative sparse linear system solver info 406 * shareable <-- pointer to iterative solver info and context 407 * whose context may be shared 408 *----------------------------------------------------------------------------*/ 409 410 void 411 cs_sles_it_set_shareable(cs_sles_it_t *context, 412 const cs_sles_it_t *shareable); 413 414 #if defined(HAVE_MPI) 415 416 /*----------------------------------------------------------------------------*/ 417 /*! 418 * \brief Set MPI communicator for global reductions. 419 * 420 * The system is solved only on ranks with a non-NULL communicator or 421 * if the caller communicator has less than 2 ranks. convergence info 422 * is broadcast across the caller communicator. 423 * 424 * \param[in, out] context pointer to iterative solver info and context 425 * \param[in] comm MPI communicator for solving 426 * \param[in] caller_comm MPI communicator of caller 427 */ 428 /*----------------------------------------------------------------------------*/ 429 430 void 431 cs_sles_it_set_mpi_reduce_comm(cs_sles_it_t *context, 432 MPI_Comm comm, 433 MPI_Comm caller_comm); 434 435 #endif /* defined(HAVE_MPI) */ 436 437 /*---------------------------------------------------------------------------- 438 * Assign ordering to iterative solver. 439 * 440 * The solver context takes ownership of the order array (i.e. it will 441 * handle its later deallocation). 442 * 443 * This is useful only for Block Gauss-Seidel. 444 * 445 * parameters: 446 * context <-> pointer to iterative solver info and context 447 * order <-> pointer to ordering array 448 *----------------------------------------------------------------------------*/ 449 450 void 451 cs_sles_it_assign_order(cs_sles_it_t *context, 452 cs_lnum_t **order); 453 454 /*----------------------------------------------------------------------------*/ 455 /*! 456 * \brief Define convergence level under which the fallback to another 457 * solver may be used if applicable. 458 * 459 * Currently, this mechanism is only by default used for BiCGstab and 460 * 3-layer conjugate residual solvers with scalar matrices, which may 461 * fall back to a preconditioned GMRES solver. For those solvers, the 462 * default threshold is \ref CS_SLES_BREAKDOWN, meaning that divergence 463 * (but not breakdown) will lead to the use of the fallback mechanism. 464 * 465 * \param[in, out] context pointer to iterative solver info and context 466 * \param[in] threshold convergence level under which fallback is used 467 */ 468 /*----------------------------------------------------------------------------*/ 469 470 void 471 cs_sles_it_set_fallback_threshold(cs_sles_it_t *context, 472 cs_sles_convergence_state_t threshold); 473 474 /*----------------------------------------------------------------------------*/ 475 /*! 476 * \brief Define the number of iterations to be done before restarting the 477 * solver. Useful only for GCR or GMRES algorithms. 478 * 479 * \param[in, out] context pointer to iterative solver info and context 480 * \param[in] interval convergence level under which fallback is used 481 */ 482 /*----------------------------------------------------------------------------*/ 483 484 void 485 cs_sles_it_set_restart_interval(cs_sles_it_t *context, 486 int interval); 487 488 /*---------------------------------------------------------------------------- 489 * Query mean number of rows under which Conjugate Gradient algorithm 490 * uses the single-reduction variant. 491 * 492 * The single-reduction variant requires only one parallel sum per 493 * iteration (instead of 2), at the cost of additional vector operations, 494 * so it tends to be more expensive when the number of matrix rows per 495 * MPI rank is high, then becomes cheaper when the MPI latency cost becomes 496 * more significant. 497 * 498 * This option is ignored for non-parallel runs, so 0 is returned. 499 * 500 * return: 501 * mean number of rows per active rank under which the 502 * single-reduction variant will be used 503 *----------------------------------------------------------------------------*/ 504 505 cs_lnum_t 506 cs_sles_it_get_pcg_single_reduction(void); 507 508 /*---------------------------------------------------------------------------- 509 * Set mean number of rows under which Conjugate Gradient algorithm 510 * should use the single-reduction variant. 511 * 512 * The single-reduction variant requires only one parallel sum per 513 * iteration (instead of 2), at the cost of additional vector operations, 514 * so it tends to be more expensive when the number of matrix rows per 515 * MPI rank is high, then becomes cheaper when the MPI latency cost becomes 516 * more significant. 517 * 518 * This option is ignored for non-parallel runs. 519 * 520 * parameters: 521 * threshold <-- mean number of rows per active rank under which the 522 * single-reduction variant will be used 523 *----------------------------------------------------------------------------*/ 524 525 void 526 cs_sles_it_set_pcg_single_reduction(cs_lnum_t threshold); 527 528 /*---------------------------------------------------------------------------- 529 * Log the current global settings relative to parallelism. 530 *----------------------------------------------------------------------------*/ 531 532 void 533 cs_sles_it_log_parallel_options(void); 534 535 /*---------------------------------------------------------------------------- 536 * Error handler for iterative sparse linear equation solver. 537 * 538 * In case of divergence or breakdown, this error handler outputs 539 * postprocessing data to assist debugging, then aborts the run. 540 * It does nothing in case the maximum iteration count is reached. 541 * 542 * parameters: 543 * sles <-> pointer to solver object 544 * state <-- convergence state 545 * a <-- matrix 546 * rhs <-- right hand side 547 * vx <-> system solution 548 * 549 * returns: 550 * false (do not attempt new solve) 551 *----------------------------------------------------------------------------*/ 552 553 bool 554 cs_sles_it_error_post_and_abort(cs_sles_t *sles, 555 cs_sles_convergence_state_t state, 556 const cs_matrix_t *a, 557 const cs_real_t *rhs, 558 cs_real_t *vx); 559 560 /*---------------------------------------------------------------------------- 561 * Set plotting options for an iterative sparse linear equation solver. 562 * 563 * parameters: 564 * context <-> pointer to iterative solver info and context 565 * base_name <-- base plot name to activate, NULL otherwise 566 * use_iteration <-- if true, use iteration as time stamp 567 * otherwise, use wall clock time 568 *----------------------------------------------------------------------------*/ 569 570 void 571 cs_sles_it_set_plot_options(cs_sles_it_t *context, 572 const char *base_name, 573 bool use_iteration); 574 575 /*----------------------------------------------------------------------------*/ 576 577 END_C_DECLS 578 579 #endif /* __CS_SLES_IT_H__ */ 580