1 /*============================================================================
2 * Sparse Linear Equation Solver driver
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 <stdarg.h>
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <assert.h>
38 #include <math.h>
39
40 #if defined(HAVE_MPI)
41 #include <mpi.h>
42 #endif
43
44 /*----------------------------------------------------------------------------
45 * Local headers
46 *----------------------------------------------------------------------------*/
47
48 #include "bft_mem.h"
49 #include "bft_error.h"
50 #include "bft_printf.h"
51
52 #include "cs_base.h"
53 #include "cs_blas.h"
54 #include "cs_field.h"
55 #include "cs_log.h"
56 #include "cs_halo.h"
57 #include "cs_map.h"
58 #include "cs_mesh.h"
59 #include "cs_mesh_location.h"
60 #include "cs_matrix.h"
61 #include "cs_matrix_default.h"
62 #include "cs_matrix_util.h"
63 #include "cs_parall.h"
64 #include "cs_post.h"
65 #include "cs_timer.h"
66 #include "cs_timer_stats.h"
67 #include "cs_time_step.h"
68
69 /*----------------------------------------------------------------------------
70 * Header for the current file
71 *----------------------------------------------------------------------------*/
72
73 #include "cs_sles.h"
74
75 /*----------------------------------------------------------------------------*/
76
77 BEGIN_C_DECLS
78
79 /*=============================================================================
80 * Additional doxygen documentation
81 *============================================================================*/
82
83 /*!
84 \file cs_sles.c
85
86 \brief Sparse linear equation solver driver.
87
88 The Sparse Linear Equation Solver subsystem is designed to allow
89 both simple usage of built-in solvers, and plugging of solvers
90 from external libraries.
91
92 As the options associated with different solvers may be very varied,
93 this subsystem is based on the use of a series of callback functions
94 which may be associated with a given linear system (defined either
95 by field id for quick and recurrent access, or by unique system
96 name). It is possible to provide a default function so calls for
97 the resolution of systems not specified before can be handled.
98
99 To summarize, the functions here define a linear system solver
100 driver, with the real work being done by functions bound to this model.
101 The main intent is to help manage passing varied user options to the
102 solvers, and handling consistency of logging.
103
104 \enum cs_sles_convergence_state_t
105
106 \brief Convergence status indicator.
107
108 \var CS_SLES_DIVERGED
109 The solver has diverged
110 \var CS_SLES_BREAKDOWN
111 The solver has broken down, and cannot make any more progress
112 \var CS_SLES_MAX_ITERATION
113 Maximum number of iterations has been reached, without reaching
114 the required convergence
115 \var CS_SLES_ITERATING
116 The solver is iterating
117 \var CS_SLES_CONVERGED
118 The solver has converged
119
120 \typedef cs_sles_setup_t
121
122 \brief Function pointer for pre-resolution setup of a linear system's
123 context.
124
125 This setup may include building a multigrid hierarchy, or a preconditioner.
126
127 Use of this type of function is optional: the context is expected to
128 maintain state, so that if a cs_sles_solve_t function is called before a
129 \ref cs_sles_setup_t function, the latter will be called automatically.
130
131 \param[in, out] context pointer to solver context
132 \param[in] name pointer to name of linear system
133 \param[in] a matrix
134 \param[in] verbosity associated verbosity
135
136 \typedef cs_sles_solve_t
137
138 \brief Function pointer for resolution of a linear system.
139
140 If the associated cs_sles_setup_t function has not been called before
141 this function, it will be called automatically.
142
143 The solution context setup by this call (or that of the matching
144 \ref cs_sles_setup_t function) will be maintained until the matching
145 \ref cs_sles_free_t function is called.
146
147 The matrix is not expected to change between successive calls, although
148 the right hand side may. If the matrix changes, the associated
149 \ref cs_sles_setup_t or \ref cs_sles_free_t function must be called between
150 solves.
151
152 The system is considered to have converged when
153 residue/r_norm <= precision, residue being the L2 norm of a.vx-rhs.
154
155 \param[in, out] context pointer to solver context
156 \param[in] name pointer to name of linear system
157 \param[in] a matrix
158 \param[in] verbosity associated verbosity
159 \param[in] precision solver precision
160 \param[in] r_norm residue normalization
161 \param[out] n_iter number of "equivalent" iterations
162 \param[out] residue residue
163 \param[in] rhs right hand side
164 \param[out] vx system solution
165 \param[in] aux_size number of elements in aux_vectors
166 \param aux_vectors optional working area
167 (internal allocation if NULL)
168
169 \return convergence status
170
171 \typedef cs_sles_free_t
172
173 \brief Function pointer for freeing of a linear system's context data.
174
175 Note that this function should free resolution-related data, such as
176 multigrid hierarchy, preconditioning, and any other temporary arrays or
177 objects required for resolution, but should not free the whole context,
178 as info used for logging (especially performance data) should be
179 maintained.
180
181 \param[in, out] context pointer to solver context.
182
183 \typedef cs_sles_log_t
184
185 \brief Function pointer for logging of linear solver history
186 and performance data.
187
188 This function will be called for each solver when \ref cs_sles_finalize
189 is called.
190
191 \param[in] context pointer to solver context
192 \param[in] log_type log type
193
194 \typedef cs_sles_copy_t
195
196 \brief Function pointer for creation of a solver context based on
197 the copy of another.
198
199 The new context copies the settings of the copied context, but not
200 its setup data and logged info, such as performance data.
201
202 This type of function is optional, but enables associating different
203 solvers to related systems (to differentiate logging) while using
204 the same settings by default.
205
206 \param[in] context context to copy
207
208 \return pointer to newly created context
209
210 \typedef cs_sles_destroy_t
211
212 Function pointer for destruction of a linear system solver context.
213
214 This function should free all context data, and will be called for each
215 system when \ref cs_sles_finalize is called.
216
217 \param[in, out] context pointer to solver context
218
219 \typedef cs_sles_error_handler_t
220
221 \brief Function pointer for handling of non-convergence when solving
222 a linear system.
223
224 Such a function is optional, and may be used for a variety of purposes,
225 such as logging, postprocessing, re-trying with different parameters,
226 aborting the run, or any combination thereof.
227
228 An error handler may be associated with a given solver context using
229 \ref cs_sles_set_error_handler, in which case it will be called whenever
230 convergence fails.
231
232 \param[in, out] sles pointer to solver object
233 \param[in] status convergence status
234 \param[in] a matrix
235 \param[in] rhs Right hand side
236 \param[out] vx System solution
237
238 \return true if solve should be re-executed, false otherwise
239
240 \typedef cs_sles_define_t
241
242 \brief Function pointer for the default definition of a sparse
243 linear equation solver
244
245 The function may be associated using \ref cs_sles_set_default_define,
246 so that it may provide a definition that will be used when
247 \ref cs_sles_setup or \ref cs_sles_solve is used for a system for which
248 no matching call to \ref cs_sles_define has been done.
249
250 The function should call \ref cs_sles_define with arguments \p f_id
251 and \p name, and appropriately chosen function pointers.
252
253 A pointer to the matrix of the system to be solved is also provided,
254 so that the corresponding information may be used to better choose
255 defaults.
256
257 \param[in] f_id associated field id, or < 0
258 \param[in] name associated name if f_id < 0, or NULL
259 \param[in] a matrix
260
261 \typedef cs_sles_verbosity_t
262
263 \brief Function pointer for the default definition of a sparse
264 linear equation solver's verbosity
265
266 The function may be associated using \ref cs_sles_set_default_verbosity, so
267 that it may provide a definition that will be used when
268 \ref cs_sles_find_or_add is called.
269
270 \param[in] f_id associated field id, or < 0
271 \param[in] name associated name if f_id < 0, or NULL
272
273 \return default verbosity value
274 */
275
276 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
277
278 /*=============================================================================
279 * Local Macro Definitions
280 *============================================================================*/
281
282 #define EPZERO 1.E-12
283
284 /*=============================================================================
285 * Local Structure Definitions
286 *============================================================================*/
287
288 /* Postprocessing of linear system residuals */
289 /*-------------------------------------------*/
290
291 typedef struct {
292
293 int writer_id; /* writer id in case of
294 postprocessing (0 for none) */
295
296 cs_lnum_t n_rows; /* number of rows */
297 cs_lnum_t block_size; /* size of block */
298
299 cs_real_t *row_residual; /* residual */
300
301 } cs_sles_post_t;
302
303 /* Basic per linear system options and logging */
304 /*---------------------------------------------*/
305
306 struct _cs_sles_t {
307
308 int n_calls; /* Number of setup or solve
309 calls */
310
311 int n_no_op; /* Number of solves with immediate
312 exit */
313
314 int f_id; /* matching field id, or < 0 */
315
316 const char *name; /* name if f_id < 0, or NULL */
317 char *_name; /* private name if f_id < 0,
318 or NULL */
319
320 int verbosity; /* verbosity level */
321
322 int type_id; /* id of solver type */
323 void *context; /* solver context
324 (options, state, logging) */
325
326 cs_sles_setup_t *setup_func; /* solver setup function */
327 cs_sles_solve_t *solve_func; /* solve function */
328 cs_sles_free_t *free_func; /* free setup function */
329
330 cs_sles_log_t *log_func; /* logging function */
331
332 cs_sles_copy_t *copy_func; /* copy function */
333
334 cs_sles_destroy_t *destroy_func; /* destruction function */
335
336 cs_sles_error_handler_t *error_func; /* error handler */
337
338 cs_sles_post_t *post_info; /* postprocessing info */
339
340 };
341
342 /*============================================================================
343 * Global variables
344 *============================================================================*/
345
346 /* Type name map */
347
348 static cs_map_name_to_id_t *_type_name_map = NULL;
349
350 /* Pointers to default definitions */
351
352 static cs_sles_define_t *_cs_sles_define_default = NULL;
353 static cs_sles_verbosity_t *_cs_sles_default_verbosity = NULL;
354
355 /* Current and maximum number of systems respectively defined by field id,
356 by name, or redefined after use */
357
358 static int _cs_sles_n_systems[3] = {0, 0, 0};
359 static int _cs_sles_n_max_systems[3] = {0, 0, 0};
360
361 /* Arrays of settings and status for linear systems */
362 static cs_sles_t **_cs_sles_systems[3] = {NULL, NULL, NULL};
363
364 /* Timer statistics */
365
366 static cs_timer_counter_t _sles_t_tot; /* Total time in linear solvers */
367 static int _sles_stat_id = -1;
368
369 /*============================================================================
370 * Private function definitions
371 *============================================================================*/
372
373 /*----------------------------------------------------------------------------
374 * Create new sparse linear equation solver structure.
375 *
376 * parameters:
377 * f_id <-- associated field id, or < 0
378 * name <-- associated name if f_id < 0, or NULL
379 *
380 * returns:
381 * pointer to associated linear system object.
382 *----------------------------------------------------------------------------*/
383
384 static cs_sles_t *
_sles_create(int f_id,const char * name)385 _sles_create(int f_id,
386 const char *name)
387 {
388 cs_sles_t *sles;
389
390 BFT_MALLOC(sles, 1, cs_sles_t);
391
392 sles->f_id = f_id;
393
394 if (f_id < 0 && name != NULL) {
395 BFT_MALLOC(sles->_name, strlen(name) + 1, char);
396 strcpy(sles->_name, name);
397 }
398 else
399 sles->_name = NULL;
400
401 if (_cs_sles_default_verbosity != NULL)
402 sles->verbosity = _cs_sles_default_verbosity(f_id, name);
403 else
404 sles->verbosity = 0;
405
406 if (_type_name_map == NULL)
407 _type_name_map = cs_map_name_to_id_create();
408 sles->type_id = cs_map_name_to_id(_type_name_map, "<undefined>");
409
410 sles->name = sles->_name;
411
412 sles->context = NULL;
413 sles->setup_func = NULL;
414 sles->solve_func = NULL;
415 sles->free_func = NULL;
416 sles->log_func = NULL;
417 sles->copy_func = NULL;
418 sles->destroy_func = NULL;
419 sles->error_func = NULL;
420
421 sles->n_calls = 0;
422 sles->n_no_op = 0;
423
424 sles->post_info = NULL;
425
426 return sles;
427 }
428
429 /*----------------------------------------------------------------------------
430 * Return pointer to linear system object, based on matching field id.
431 *
432 * If this system did not previously exist, it is added to the list of
433 * "known" systems.
434 *
435 * parameters:
436 * f_id <-- associated field id
437 *
438 * returns:
439 * pointer to linear system object
440 *----------------------------------------------------------------------------*/
441
442 static cs_sles_t *
_find_or_add_system_by_f_id(int f_id)443 _find_or_add_system_by_f_id(int f_id)
444 {
445 assert(f_id >= 0);
446
447 /* Return system id already defined */
448 if (f_id < _cs_sles_n_max_systems[0]) {
449 if (_cs_sles_systems[0][f_id] != NULL)
450 return _cs_sles_systems[0][f_id];
451 }
452
453 /* Increase size of systems array if required */
454 else {
455 int i = _cs_sles_n_max_systems[0];
456
457 if (_cs_sles_n_max_systems[0] == 0)
458 _cs_sles_n_max_systems[0] = 1;
459 while (_cs_sles_n_max_systems[0] <= f_id)
460 _cs_sles_n_max_systems[0] *= 2;
461 BFT_REALLOC(_cs_sles_systems[0],
462 _cs_sles_n_max_systems[0],
463 cs_sles_t *);
464
465 for (int j = i; j < _cs_sles_n_max_systems[0]; j++)
466 _cs_sles_systems[0][j] = NULL;
467 }
468
469 /* If we have not returned yet, we need to add a new system to the array */
470
471 _cs_sles_systems[0][f_id] = _sles_create(f_id, NULL);
472 _cs_sles_n_systems[0] += 1;
473
474 return _cs_sles_systems[0][f_id];
475 }
476
477 /*----------------------------------------------------------------------------
478 * Return pointer to linear system object, based on system name.
479 *
480 * If this system did not previously exist, it is added to the list of
481 * "known" systems.
482 *
483 * parameters:
484 * name <-- system name
485 *
486 * returns:
487 * pointer to linear system object
488 *----------------------------------------------------------------------------*/
489
490 static cs_sles_t *
_find_or_add_system_by_name(const char * name)491 _find_or_add_system_by_name(const char *name)
492 {
493 int ii, start_id, end_id, mid_id;
494 int cmp_ret = 1;
495
496 /* Use binary search to find system */
497
498 start_id = 0;
499 end_id = _cs_sles_n_systems[1] - 1;
500 mid_id = start_id + ((end_id -start_id) / 2);
501
502 while (start_id <= end_id) {
503 cmp_ret = strcmp((_cs_sles_systems[1][mid_id])->name, name);
504 if (cmp_ret < 0)
505 start_id = mid_id + 1;
506 else if (cmp_ret > 0)
507 end_id = mid_id - 1;
508 else
509 break;
510 mid_id = start_id + ((end_id -start_id) / 2);
511 }
512
513 /* If found, return */
514
515 if (cmp_ret == 0)
516 return _cs_sles_systems[1][mid_id];
517
518 /* Reallocate global array if necessary */
519
520 if (_cs_sles_n_systems[1] >= _cs_sles_n_max_systems[1]) {
521 int i = _cs_sles_n_max_systems[1];
522
523 if (_cs_sles_n_max_systems[1] == 0)
524 _cs_sles_n_max_systems[1] = 1;
525 _cs_sles_n_max_systems[1] *= 2;
526 BFT_REALLOC(_cs_sles_systems[1],
527 _cs_sles_n_max_systems[1],
528 cs_sles_t*);
529
530 for (int j = i; j < _cs_sles_n_max_systems[1]; j++)
531 _cs_sles_systems[1][j] = NULL;
532 }
533
534 /* Insert in sorted list */
535
536 for (ii = _cs_sles_n_systems[1]; ii > mid_id; ii--)
537 _cs_sles_systems[1][ii] = _cs_sles_systems[1][ii - 1];
538
539 _cs_sles_systems[1][mid_id] = _sles_create(-1, name);
540 _cs_sles_n_systems[1] += 1;
541
542 return _cs_sles_systems[1][mid_id];
543 }
544
545 /*----------------------------------------------------------------------------
546 * Copy linear system info whose definition has changed.
547 *
548 * This is intended to maintain logging info, and is useful only
549 * if the solver has been called using the older definition.
550 *
551 * parameters:
552 * s <-> pointer to linear system info
553 *----------------------------------------------------------------------------*/
554
555 static void
_save_system_info(cs_sles_t * s)556 _save_system_info(cs_sles_t *s)
557 {
558 assert(s != NULL);
559
560 /* Resize array if needed */
561
562 int i = _cs_sles_n_systems[2];
563
564 if (i >= _cs_sles_n_max_systems[2]) {
565
566 if (_cs_sles_n_max_systems[2] == 0)
567 _cs_sles_n_max_systems[2] = 1;
568 _cs_sles_n_max_systems[2] *=2;
569 BFT_REALLOC(_cs_sles_systems[2],
570 _cs_sles_n_max_systems[2],
571 cs_sles_t *);
572
573 for (int j = i; j < _cs_sles_n_max_systems[2]; j++)
574 _cs_sles_systems[2][j] = NULL;
575
576 }
577
578 /* Ensure no extra data is maintained for old system */
579
580 if (s->free_func != NULL)
581 s->free_func(s->context);
582
583 /* Save other context and options */
584
585 cs_sles_t *s_old;
586 BFT_MALLOC(s_old, 1, cs_sles_t);
587 memcpy(s_old, s, sizeof(cs_sles_t));
588
589 s_old->_name = NULL; /* still points to new name */
590 s->context = NULL; /* old context now only available through s_old */
591
592 _cs_sles_systems[2][i] = s_old;
593
594 _cs_sles_n_systems[2] += 1;
595 }
596
597 /*----------------------------------------------------------------------------
598 * Compute per-cell residual for Ax = b.
599 *
600 * parameters:
601 * n_vals <-- Number of values
602 * a <-- Linear equation matrix
603 * rhs <-- Right hand side
604 * vx <-> Current system solution
605 * res --> Residual
606 *----------------------------------------------------------------------------*/
607
608 static void
_residual(cs_lnum_t n_vals,const cs_matrix_t * a,const cs_real_t rhs[],cs_real_t vx[],cs_real_t res[])609 _residual(cs_lnum_t n_vals,
610 const cs_matrix_t *a,
611 const cs_real_t rhs[],
612 cs_real_t vx[],
613 cs_real_t res[])
614 {
615 cs_matrix_vector_multiply(a, vx, res);
616
617 # pragma omp parallel for if(n_vals > CS_THR_MIN)
618 for (cs_lnum_t ii = 0; ii < n_vals; ii++)
619 res[ii] = fabs(res[ii] - rhs[ii]);
620 }
621
622 /*----------------------------------------------------------------------------
623 * Test if a general sparse linear system needs solving or if the right-hand
624 * side is already zero within convergence criteria.
625 *
626 * The computed residue is also updated;
627 *
628 * parameters:
629 * name <-- name of the associated system
630 * a <-- pointer to matrix
631 * verbosity <-- verbosity level
632 * precision <-- solver precision
633 * r_norm <-- residue normalization
634 * residue <-> residue
635 * vx <-- initial solution
636 * rhs <-- right hand side
637 *
638 * returns:
639 * 1 if solving is required, 0 if the rhs is already zero within tolerance
640 * criteria (precision of residue normalization)
641 *----------------------------------------------------------------------------*/
642
643 static int
_needs_solving(const char * name,const cs_matrix_t * a,int verbosity,double precision,double r_norm,double * residue,const cs_real_t * vx,const cs_real_t * rhs)644 _needs_solving(const char *name,
645 const cs_matrix_t *a,
646 int verbosity,
647 double precision,
648 double r_norm,
649 double *residue,
650 const cs_real_t *vx,
651 const cs_real_t *rhs)
652 {
653 int retval = 1;
654
655 /* Initialize residue, check for immediate return */
656
657 const cs_lnum_t *diag_block_size = cs_matrix_get_diag_block_size(a);
658 const cs_lnum_t _diag_block_size = diag_block_size[1];
659 assert(diag_block_size[0] == diag_block_size[1]);
660
661 const cs_lnum_t n_rows = cs_matrix_get_n_rows(a) * _diag_block_size;
662
663 double r[2] = {
664 cs_dot_xx(n_rows, rhs),
665 cs_dot_xx(n_rows, vx)
666 };
667 cs_parall_sum(2, CS_DOUBLE, r);
668
669 /* If the initial solution is "true" zero (increment mode), we can determine
670 convergence without resorting to a matrix-vector product */
671
672 if (r[1] < 1e-60) {
673
674 double _precision = CS_MIN(EPZERO, precision); /* prefer to err on the side
675 of caution... */
676 *residue = sqrt(r[0]);
677
678 if (r_norm <= EPZERO)
679 retval = 0;
680 else if (*residue/r_norm <= _precision)
681 retval = 0;
682
683 if (retval == 0 && verbosity > 1)
684 bft_printf(_("[%s]:\n"
685 " immediate exit; r_norm = %11.4e, residual = %11.4e\n"),
686 name, r_norm, *residue);
687
688 }
689 else
690 *residue = HUGE_VAL; /* actually unknown, since we did not multiply
691 by A (we might as well enter the solver,
692 and expect to have vx = 0 most of the time) */
693
694 return retval;
695 }
696
697 /*----------------------------------------------------------------------------
698 * Output post-processing data for failed system convergence.
699 *
700 * parameters:
701 * n_vals <-- Size of val and val_type array
702 * val <-> Values to post-process (set to 0 on output if not
703 * normal floating-point values)
704 * val_type --> 0: normal values, 1: infinite, 2: Nan
705 *
706 * returns:
707 * number of non-normal values
708 *----------------------------------------------------------------------------*/
709
710 static size_t
_value_type(size_t n_vals,cs_real_t val[],cs_real_t val_type[])711 _value_type(size_t n_vals,
712 cs_real_t val[],
713 cs_real_t val_type[])
714 {
715 size_t ii;
716 size_t retval = 0;
717
718 #if (__STDC_VERSION__ >= 199901L)
719
720 for (ii = 0; ii < n_vals; ii++) {
721
722 int v_type = fpclassify(val[ii]);
723
724 if (v_type == FP_INFINITE) {
725 val[ii] = 0.;
726 val_type[ii] = 1;
727 retval += 1;
728 }
729
730 else if (v_type == FP_NAN) {
731 val[ii] = 0.;
732 val_type[ii] = 2;
733 retval += 1;
734 }
735
736 else if (val[ii] > 1.e38 || val[ii] < -1.e38) {
737 val[ii] = 0.;
738 val_type[ii] = 1;
739 retval += 1;
740 }
741
742 else
743 val_type[ii] = 0;
744 }
745
746 #else
747
748 for (ii = 0; ii < n_vals; ii++) {
749
750 if (val[ii] != val[ii]) { /* Test for NaN with IEEE 754 arithmetic */
751 val[ii] = 0.;
752 val_type[ii] = 2;
753 retval += 1;
754 }
755
756 else if (val[ii] > 1.e38 || val[ii] < -1.e38) {
757 val[ii] = 0.;
758 val_type[ii] = 1;
759 retval += 1;
760 }
761
762 else
763 val_type[ii] = 0;
764 }
765
766 #endif
767
768 return retval;
769 }
770
771 /*----------------------------------------------------------------------------
772 * Ensure array for postprocessing output of residuals is allocated
773 *
774 * parameters
775 * sles <-> pointer to solver object
776 * a <-- matrix
777 *----------------------------------------------------------------------------*/
778
779 static void
_ensure_alloc_post(cs_sles_t * sles,const cs_matrix_t * a)780 _ensure_alloc_post(cs_sles_t *sles,
781 const cs_matrix_t *a)
782 {
783 if (sles->post_info != NULL) {
784
785 const cs_lnum_t *diag_block_size = cs_matrix_get_diag_block_size(a);
786 const cs_lnum_t n_vals = cs_matrix_get_n_columns(a) * diag_block_size[1];
787
788 sles->post_info->n_rows = cs_matrix_get_n_rows(a);
789 sles->post_info->block_size = diag_block_size[1];
790
791 BFT_REALLOC(sles->post_info->row_residual, n_vals, cs_real_t);
792
793 assert(diag_block_size[0] == diag_block_size[1]); /* for now */
794
795 }
796 }
797
798 /*----------------------------------------------------------------------------
799 * Post process the residual for a given linear equation solver
800 *
801 * parameters:
802 * sles_p <-- void pointer to sparse linear equation solver context
803 * ts <-- time step status structure
804 *----------------------------------------------------------------------------*/
805
806 static void
_post_function(void * sles_p,const cs_time_step_t * ts)807 _post_function(void *sles_p,
808 const cs_time_step_t *ts)
809 {
810 CS_UNUSED(ts);
811
812 cs_sles_t *sles = sles_p;
813
814 cs_sles_post_t *sp = sles->post_info;
815
816 assert(sp != NULL);
817
818 const cs_mesh_t *mesh = cs_glob_mesh;
819
820 int mesh_id = CS_POST_MESH_VOLUME;
821
822 char base_name[32], val_name[32];
823
824 /* Check for mesh location */
825
826 int location_id = 0, flag = 0;
827 if (sp->n_rows == mesh->n_cells)
828 location_id = CS_MESH_LOCATION_CELLS;
829 else if (sp->n_rows == mesh->n_vertices)
830 location_id = CS_MESH_LOCATION_VERTICES;
831
832 flag = location_id;
833 cs_parall_max(1, CS_INT_TYPE, &flag);
834 if (flag == location_id)
835 flag = 0;
836 else
837 flag =1;
838 cs_parall_max(1, CS_INT_TYPE, &flag);
839 if (flag != 0)
840 return;
841
842 strcpy(base_name, "Residual");
843
844 const char *name = cs_sles_get_name(sles);
845
846 if (strlen(name) + strlen(base_name) < 31) {
847 strcpy(val_name, base_name);
848 strcat(val_name, "_");
849 strcat(val_name, name);
850 }
851 else {
852 strncpy(val_name, base_name, 31);
853 val_name[31] = '\0';
854 }
855
856 cs_sles_post_output_var(val_name,
857 mesh_id,
858 location_id,
859 sp->writer_id,
860 sp->block_size,
861 sp->row_residual);
862
863 BFT_FREE(sp->row_residual);
864 }
865
866 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
867
868 /*============================================================================
869 * Public function definitions for Fortran API
870 *============================================================================*/
871
872 /*============================================================================
873 * Public function definitions
874 *============================================================================*/
875
876 /*----------------------------------------------------------------------------*/
877 /*!
878 * \brief Initialize sparse linear equation solver API.
879 */
880 /*----------------------------------------------------------------------------*/
881
882 void
cs_sles_initialize(void)883 cs_sles_initialize(void)
884 {
885 CS_TIMER_COUNTER_INIT(_sles_t_tot);
886
887 int stats_root = cs_timer_stats_id_by_name("operations");
888
889 if (stats_root > -1) {
890 _sles_stat_id = cs_timer_stats_create("operations",
891 "linear_solvers",
892 "linear solvers");
893 }
894 }
895
896 /*----------------------------------------------------------------------------*/
897 /*!
898 * \brief Finalize sparse linear equation solver API.
899 */
900 /*----------------------------------------------------------------------------*/
901
902 void
cs_sles_finalize(void)903 cs_sles_finalize(void)
904 {
905 for (int i = 0; i < 3; i++) {
906
907 for (int j = 0; j < _cs_sles_n_max_systems[i]; j++) {
908
909 if (_cs_sles_systems[i][j] != NULL) {
910 cs_sles_t *sles = _cs_sles_systems[i][j];
911 if (sles->free_func != NULL)
912 sles->free_func(sles->context);
913 if (sles->destroy_func != NULL)
914 sles->destroy_func(&(sles->context));
915 if (sles->post_info != NULL) {
916 BFT_FREE(sles->post_info->row_residual);
917 BFT_FREE(sles->post_info);
918 }
919 BFT_FREE(sles->_name);
920 BFT_FREE(_cs_sles_systems[i][j]);
921 }
922
923 }
924
925 BFT_FREE(_cs_sles_systems[i]);
926 _cs_sles_n_max_systems[i] = 0;
927 _cs_sles_n_systems[i] = 0;
928
929 }
930
931 cs_map_name_to_id_destroy(&_type_name_map);
932 }
933
934 /*----------------------------------------------------------------------------*/
935 /*!
936 * \brief Log sparse linear equation solver info
937 *
938 * \param[in] log_type log type (CS_LOG_SETUP or CS_LOG_PERFORMANCE)
939 */
940 /*----------------------------------------------------------------------------*/
941
942 void
cs_sles_log(cs_log_t log_type)943 cs_sles_log(cs_log_t log_type)
944 {
945 int n_tot_systems = 0;
946 for (int i = 0; i < 3; i++)
947 n_tot_systems += _cs_sles_n_systems[i];
948
949 if (n_tot_systems < 1)
950 return;
951
952 int log_order[] = {2, 0, 1}; /* log previous setups first, then fields,
953 then others */
954
955 if (log_type == CS_LOG_PERFORMANCE)
956 cs_log_printf
957 (log_type,
958 _("\n"
959 "Total elapsed time for linear equation system solvers: %.3f s\n"),
960 _sles_t_tot.nsec*1e-9);
961
962 const char *option_category[]
963 = {N_("Linear solver options modified during run (previous values)"),
964 N_("Linear solver options for fields"),
965 N_("Linear solver options for other systems")};
966
967 const char *perf_category[]
968 = {N_("Linear solver performance with previous options"),
969 N_("Linear solver performance for fields"),
970 N_("Linear solver performance for other systems")};
971
972 for (int i = 0; i < 3; i++) {
973
974 int j = log_order[i];
975
976 /* Print heading (underlined) line */
977
978 if (_cs_sles_n_systems[j] > 0) {
979
980 size_t l = 0;
981 switch(log_type) {
982 case CS_LOG_SETUP:
983 l = cs_log_strlen(_(option_category[i]));
984 cs_log_printf(log_type, "\n%s\n", _(option_category[i]));
985 break;
986 case CS_LOG_PERFORMANCE:
987 l = cs_log_strlen(_(perf_category[i]));
988 cs_log_printf(log_type, "\n%s\n", _(perf_category[i]));
989 break;
990 default:
991 break;
992 }
993
994 char ul[128];
995 l = CS_MIN(l, 127);
996 for (size_t ll = 0; ll < l; ll++)
997 ul[ll] = '-';
998 ul[l] = '\0';
999 cs_log_printf(log_type, "%s\n", ul);
1000
1001 }
1002
1003 /* Logging for each system */
1004
1005 for (int k = 0; k < _cs_sles_n_max_systems[j]; k++) {
1006 cs_sles_t *sles = _cs_sles_systems[j][k];
1007
1008 if (sles == NULL) continue;
1009
1010 if (sles->log_func != NULL) {
1011
1012 const char *name = cs_sles_base_name(sles->f_id, sles->name);
1013
1014 switch(log_type) {
1015
1016 case CS_LOG_SETUP:
1017 if (sles->f_id > -1)
1018 cs_log_printf
1019 (log_type,
1020 _("\n"
1021 "Linear solver options for \"%s\" (field id %d)\n"),
1022 name, sles->f_id);
1023 else
1024 cs_log_printf
1025 (log_type,
1026 _("\n"
1027 "Linear solver options for \"%s\"\n"),
1028 name);
1029 break;
1030
1031 case CS_LOG_PERFORMANCE:
1032 if (sles->f_id > -1)
1033 cs_log_printf
1034 (log_type,
1035 _("\n"
1036 "Summary of resolutions for \"%s\" (field id %d)\n"),
1037 name, sles->f_id);
1038 else
1039 cs_log_printf
1040 (log_type,
1041 _("\n"
1042 "Summary of resolutions for \"%s\"\n"),
1043 name);
1044 break;
1045
1046 default:
1047 break;
1048 }
1049
1050 sles->log_func(sles->context, log_type);
1051
1052 switch(log_type) {
1053
1054 case CS_LOG_SETUP:
1055 cs_log_printf
1056 (log_type,
1057 _(" Verbosity: %d\n"), sles->verbosity);
1058 if (sles->post_info != NULL)
1059 cs_log_printf
1060 (log_type,
1061 _(" Residual postprocessing writer id: %d\n"),
1062 sles->post_info->writer_id);
1063 break;
1064
1065 case CS_LOG_PERFORMANCE:
1066 if (sles->n_no_op > 0)
1067 cs_log_printf
1068 (log_type,
1069 _("\n"
1070 " Number of immediate solve exits: %d\n"), sles->n_no_op);
1071 break;
1072
1073 default:
1074 break;
1075 }
1076
1077 }
1078 }
1079
1080 }
1081
1082 cs_log_printf(log_type, "\n");
1083 cs_log_separator(log_type);
1084 }
1085
1086 /*----------------------------------------------------------------------------*/
1087 /*!
1088 * \brief Return pointer to linear system object, based on matching field id or
1089 * system name.
1090 *
1091 * If this system did not previously exist, NULL is returned.
1092 *
1093 * \param[in] f_id associated field id, or < 0
1094 * \param[in] name associated name if f_id < 0, or NULL
1095 *
1096 * \return pointer to associated linear system object, or NULL
1097 */
1098 /*----------------------------------------------------------------------------*/
1099
1100 cs_sles_t *
cs_sles_find(int f_id,const char * name)1101 cs_sles_find(int f_id,
1102 const char *name)
1103 {
1104 cs_sles_t *retval = NULL;
1105
1106 if (f_id >= 0) {
1107 if (f_id < _cs_sles_n_max_systems[0]) {
1108 if (_cs_sles_systems[0][f_id] != NULL) {
1109 retval = _cs_sles_systems[0][f_id];
1110 /* Check for masked ("pushed") definition */
1111 if (retval->name != NULL)
1112 retval = cs_sles_find(-1, retval->name);
1113 }
1114 }
1115 }
1116
1117 else if (name != NULL) {
1118
1119 int start_id, end_id, mid_id;
1120 int cmp_ret = 1;
1121
1122 /* Use binary search to find system */
1123
1124 start_id = 0;
1125 end_id = _cs_sles_n_systems[1] - 1;
1126 mid_id = start_id + ((end_id -start_id) / 2);
1127
1128 while (start_id <= end_id) {
1129 cmp_ret = strcmp((_cs_sles_systems[1][mid_id])->name, name);
1130 if (cmp_ret < 0)
1131 start_id = mid_id + 1;
1132 else if (cmp_ret > 0)
1133 end_id = mid_id - 1;
1134 else
1135 break;
1136 mid_id = start_id + ((end_id -start_id) / 2);
1137 }
1138
1139 if (cmp_ret == 0)
1140 retval = _cs_sles_systems[1][mid_id];
1141 }
1142
1143 return retval;
1144 }
1145
1146 /*----------------------------------------------------------------------------*/
1147 /*!
1148 * \brief Return pointer to linear system object, based on matching field id or
1149 * system name.
1150 *
1151 * If this system did not previously exist, it is created and added to
1152 * to the list of "known" systems. In this case, it will be usable
1153 * only if cs_sles_define() is called for the same field id and name
1154 * (in which case calling the present function is redundant), or if
1155 * cs_sles_set_sefault_define() has been previously used to define
1156 * the default solver policy.
1157 *
1158 * \param[in] f_id associated field id, or < 0
1159 * \param[in] name associated name if f_id < 0, or NULL
1160 *
1161 * \return pointer to associated linear system object, or NULL
1162 */
1163 /*----------------------------------------------------------------------------*/
1164
1165 cs_sles_t *
cs_sles_find_or_add(int f_id,const char * name)1166 cs_sles_find_or_add(int f_id,
1167 const char *name)
1168 {
1169 cs_sles_t *retval = NULL;
1170
1171 if (f_id >= 0) {
1172 retval = _find_or_add_system_by_f_id(f_id);
1173 /* Check for masked ("pushed") definition */
1174 if (retval->name != NULL)
1175 retval = _find_or_add_system_by_name(retval->name);
1176 }
1177 else
1178 retval = _find_or_add_system_by_name(name);
1179
1180 return retval;
1181 }
1182
1183 /*----------------------------------------------------------------------------*/
1184 /*!
1185 * \brief Temporarily replace field id with name for matching calls
1186 * to \ref cs_sles_setup, \ref cs_sles_solve, \ref cs_sles_free, and other
1187 * operations involving access through a field id.
1188 *
1189 * \deprecated This function is provided to allow some peculiar
1190 * calling sequences, in which \ref cs_equation_iterative_solve_scalar is
1191 * called with a nonzero \c ivar value, but specific solver options
1192 * must still be set.
1193 * In the future, a more consistent mechanism (using a zero \c ivar
1194 * value or designing a cleaner method to handle those exceptional cases)
1195 * is preferred. As such, only a stack depth of 1 is allowed.
1196 *
1197 * \param[in] f_id associated field id, or < 0
1198 * \param[in] name associated name if f_id < 0, or NULL
1199 */
1200 /*----------------------------------------------------------------------------*/
1201
1202 void
cs_sles_push(int f_id,const char * name)1203 cs_sles_push(int f_id,
1204 const char *name)
1205 {
1206 if (f_id < 0)
1207 bft_error
1208 (__FILE__, __LINE__, 0,
1209 "%s must be called only for an actual field, with id >=0, not %d.",
1210 __func__, f_id);
1211
1212 cs_sles_t *retval = cs_sles_find_or_add(f_id, NULL);
1213
1214 if (retval->name != NULL)
1215 bft_error
1216 (__FILE__, __LINE__, 0,
1217 _("cs_sles_push() only allows a stack of depth 1:\n"
1218 " it may not be called multiple times for a given field (id %d)\n"
1219 " without calling cs_sles_pop between those calls."), f_id);
1220 else {
1221 assert(retval->_name == NULL);
1222 BFT_MALLOC(retval->_name, strlen(name) + 1, char);
1223 strcpy(retval->_name, name);
1224 retval->name = retval->_name;
1225 }
1226 }
1227
1228 /*----------------------------------------------------------------------------*/
1229 /*!
1230 * \brief Restore behavior temporarily modified by \ref cs_sles_push.
1231 *
1232 * \deprecated This function matches \ref cs_sles_push, which is deprecated.
1233 *
1234 * \param[in] f_id associated field id, or < 0
1235 */
1236 /*----------------------------------------------------------------------------*/
1237
1238 void
cs_sles_pop(int f_id)1239 cs_sles_pop(int f_id)
1240 {
1241 if (f_id < 0)
1242 bft_error
1243 (__FILE__, __LINE__, 0,
1244 "%s must be called only for an actual field, with id >=0, not %d.",
1245 __func__, f_id);
1246
1247 cs_sles_t *retval = _find_or_add_system_by_f_id(f_id);
1248
1249 retval->name = NULL;
1250 BFT_FREE(retval->_name);
1251 }
1252
1253 /*----------------------------------------------------------------------------*/
1254 /*!
1255 * \brief Define sparse linear equation solver for a given field or
1256 * equation name.
1257 *
1258 * If this system did not previously exist, it is added to the list of
1259 * "known" systems.
1260 *
1261 * The context pointer is used to point to a structure adapted to the function
1262 * pointers given here, and combined with those functions, allows using
1263 * both built-in, external, or user-defined solvers.
1264 *
1265 * It is recommended the context type name provided here directly relate
1266 * to the associated structure type (for example, "cs_sles_it_t" or
1267 * "cs_multigrid_t").
1268 *
1269 * \param[in] f_id associated field id, or < 0
1270 * \param[in] name associated name if f_id < 0, or NULL
1271 * \param[in, out] context pointer to solver context management
1272 * structure (cs_sles subsystem becomes owner)
1273 * \param[in] type_name context structure or object type name
1274 * \param[in] setup_func pointer to system setup function
1275 * \param[in] solve_func pointer to system solution function (also
1276 * calls setup_func if not done yet or free_func
1277 * called since last solve)
1278 * \param[in] free_func pointer function freeing system setup
1279 * \param[in] log_func pointer to system info logging function
1280 (optional, but recommended)
1281 * \param[in] copy_func pointer to settings copy function (optional)
1282 * \param[in] destroy_func pointer to function destroying solver context
1283 * (called with \ref cs_sles_finalize or with a
1284 * new call to this function for the same system)
1285 *
1286 * \return pointer to associated linear system object
1287 */
1288 /*----------------------------------------------------------------------------*/
1289
1290 cs_sles_t *
cs_sles_define(int f_id,const char * name,void * context,const char * type_name,cs_sles_setup_t * setup_func,cs_sles_solve_t * solve_func,cs_sles_free_t * free_func,cs_sles_log_t * log_func,cs_sles_copy_t * copy_func,cs_sles_destroy_t * destroy_func)1291 cs_sles_define(int f_id,
1292 const char *name,
1293 void *context,
1294 const char *type_name,
1295 cs_sles_setup_t *setup_func,
1296 cs_sles_solve_t *solve_func,
1297 cs_sles_free_t *free_func,
1298 cs_sles_log_t *log_func,
1299 cs_sles_copy_t *copy_func,
1300 cs_sles_destroy_t *destroy_func)
1301 {
1302 cs_sles_t * sles = cs_sles_find_or_add(f_id, name);
1303
1304 /* Check if system was previously defined and used,
1305 and save info for future logging in this case */
1306
1307 if (sles->context != NULL) {
1308 if (sles->n_calls > 0 && sles->log_func != NULL)
1309 _save_system_info(sles);
1310 else if (sles->destroy_func != NULL)
1311 sles->destroy_func(&(sles->context));
1312 }
1313
1314 if (type_name != NULL)
1315 sles->type_id = cs_map_name_to_id(_type_name_map, type_name);
1316
1317 /* Now define options */
1318 sles->context = context;
1319 sles->setup_func = setup_func;
1320 sles->solve_func = solve_func;
1321 sles->free_func = free_func;
1322 sles->log_func = log_func;
1323 sles->copy_func = copy_func;
1324 sles->destroy_func = destroy_func;
1325
1326 return sles;
1327 }
1328
1329 /*----------------------------------------------------------------------------*/
1330 /*!
1331 * \brief Set the verbosity for a given linear equation solver.
1332 *
1333 * This verbosity will be used by cs_sles_setup and cs_sles_solve.
1334 *
1335 * By default, the verbosity is set to 0, or the value returned by the
1336 * function set with cs_sles_set_default_define().
1337 *
1338 * \param[in, out] sles pointer to solver object
1339 * \param[in] verbosity verbosity level
1340 */
1341 /*----------------------------------------------------------------------------*/
1342
1343 void
cs_sles_set_verbosity(cs_sles_t * sles,int verbosity)1344 cs_sles_set_verbosity(cs_sles_t *sles,
1345 int verbosity)
1346 {
1347 sles->verbosity = verbosity;
1348 }
1349
1350 /*----------------------------------------------------------------------------*/
1351 /*!
1352 * \brief Get the verbosity for a given linear equation solver.
1353 *
1354 * This verbosity will be used by cs_sles_setup and cs_sles_solve.
1355 *
1356 * \param[in, out] sles pointer to solver object
1357 *
1358 * \return verbosity level
1359 */
1360 /*----------------------------------------------------------------------------*/
1361
1362 int
cs_sles_get_verbosity(cs_sles_t * sles)1363 cs_sles_get_verbosity(cs_sles_t *sles)
1364 {
1365 return sles->verbosity;
1366 }
1367
1368 /*----------------------------------------------------------------------------*/
1369 /*!
1370 * \brief Activate postprocessing output for a given linear equation solver.
1371 *
1372 * This allows the output of the residual at the end of each solution
1373 * series, using a single postprocessing writer.
1374 * By default, no output is activated.
1375 *
1376 * \param[in, out] sles pointer to solver object
1377 * \param[in] writer_id id of the writer
1378 */
1379 /*----------------------------------------------------------------------------*/
1380
1381 void
cs_sles_set_post_output(cs_sles_t * sles,int writer_id)1382 cs_sles_set_post_output(cs_sles_t *sles,
1383 int writer_id)
1384 {
1385 if (sles->n_calls > 0)
1386 return;
1387
1388 if (sles->post_info == NULL)
1389 cs_post_add_time_dep_output(_post_function, (void *)sles);
1390
1391 BFT_REALLOC(sles->post_info, 1, cs_sles_post_t);
1392 sles->post_info->writer_id = writer_id;
1393 sles->post_info->n_rows = 0;
1394 sles->post_info->block_size = 0;
1395 sles->post_info->row_residual = NULL;
1396 }
1397
1398 /*----------------------------------------------------------------------------*/
1399 /*!
1400 * \brief Return the id of the associated writer if postprocessing output
1401 * is active for a given linear equation solver.
1402 *
1403 * \param[in] sles pointer to solver object
1404 *
1405 * \return id od associated writer, or 0
1406 */
1407 /*----------------------------------------------------------------------------*/
1408
1409 int
cs_sles_get_post_output(cs_sles_t * sles)1410 cs_sles_get_post_output(cs_sles_t *sles)
1411 {
1412 int retval = 0;
1413
1414 if (sles->post_info != NULL)
1415 retval = sles->post_info->writer_id;
1416
1417 return retval;
1418 }
1419
1420 /*----------------------------------------------------------------------------*/
1421 /*!
1422 * \brief Return type name of solver context.
1423 *
1424 * The returned string is intended to help determine which type is associated
1425 * with the void * pointer returned by \ref cs_sles_get_context for a given
1426 * solver definition, so as to be able to call additional specific functions
1427 * beyond the generic functions assigned to a cs_sles_t object.
1428 *
1429 * If no type_name string was associated to the solver upon its definition by
1430 * \ref cs_sles_define, or it has not been defined yet, the string returned
1431 * is "<undefined>". It is recommended the type name provided
1432 * \ref cs_sles_define directly relate to the associated structure type
1433 * (for example, "cs_sles_it_t" or "cs_multigrid_t").
1434 *
1435 * \param[in] sles pointer to solver object
1436 *
1437 * \return pointer to linear system solver specific type name
1438 */
1439 /*----------------------------------------------------------------------------*/
1440
1441 const char *
cs_sles_get_type(cs_sles_t * sles)1442 cs_sles_get_type(cs_sles_t *sles)
1443 {
1444 return cs_map_name_to_id_reverse(_type_name_map, sles->type_id);
1445 }
1446
1447 /*----------------------------------------------------------------------------*/
1448 /*!
1449 * \brief Return pointer to solver context structure pointer.
1450 *
1451 * The context structure depends on the type of solver used, which may in
1452 * turn be determined by the string returned by cs_sles_get_type().
1453 * If may be used by appropriate functions specific to that type.
1454 *
1455 * \param[in] sles pointer to solver object
1456 *
1457 * \return pointer to solver-specific linear system info and context
1458 */
1459 /*----------------------------------------------------------------------------*/
1460
1461 void *
cs_sles_get_context(cs_sles_t * sles)1462 cs_sles_get_context(cs_sles_t *sles)
1463 {
1464 return sles->context;
1465 }
1466
1467 /*----------------------------------------------------------------------------*/
1468 /*!
1469 * \brief Return field id associated with a given sparse linear equation solver.
1470 *
1471 * \param[in] sles pointer to solver object
1472 *
1473 * \return associated field id (or -1 if defined by name)
1474 */
1475 /*----------------------------------------------------------------------------*/
1476
1477 int
cs_sles_get_f_id(const cs_sles_t * sles)1478 cs_sles_get_f_id(const cs_sles_t *sles)
1479 {
1480 return sles->f_id;
1481 }
1482
1483 /*----------------------------------------------------------------------------*/
1484 /*!
1485 * \brief Return name associated with a given sparse linear equation solver.
1486 *
1487 * This is simply a utility function which will return its name argument
1488 * if f_id < 0, and the associated field's name or label otherwise.
1489 *
1490 * \param[in] sles pointer to solver object
1491 *
1492 * \return pointer to associated linear system object name
1493 */
1494 /*----------------------------------------------------------------------------*/
1495
1496 const char *
cs_sles_get_name(const cs_sles_t * sles)1497 cs_sles_get_name(const cs_sles_t *sles)
1498 {
1499 return cs_sles_base_name(sles->f_id, sles->name);
1500 }
1501
1502 /*----------------------------------------------------------------------------*/
1503 /*!
1504 * \brief Setup sparse linear equation solver.
1505 *
1506 * Use of this function is optional: if a \ref cs_sles_solve is called
1507 * for the same system before this function is called, the latter will be
1508 * called automatically.
1509 *
1510 * If no options were previously provided for the matching system,
1511 * default options will be used.
1512 *
1513 * \param[in, out] sles pointer to solver object
1514 * \param[in] a matrix
1515 */
1516 /*----------------------------------------------------------------------------*/
1517
1518 void
cs_sles_setup(cs_sles_t * sles,const cs_matrix_t * a)1519 cs_sles_setup(cs_sles_t *sles,
1520 const cs_matrix_t *a)
1521 {
1522 cs_timer_t t0 = cs_timer_time();
1523
1524 if (sles->context == NULL)
1525 _cs_sles_define_default(sles->f_id, sles->name, a);
1526
1527 int t_top_id = cs_timer_stats_switch(_sles_stat_id);
1528
1529 sles->n_calls += 1;
1530
1531 if (sles->setup_func != NULL) {
1532 const char *sles_name = cs_sles_base_name(sles->f_id, sles->name);
1533 sles->setup_func(sles->context, sles_name, a, sles->verbosity);
1534 }
1535
1536 /* Prepare residual postprocessing if required */
1537
1538 if (sles->post_info != NULL) {
1539 _ensure_alloc_post(sles, a);
1540 const cs_lnum_t n_vals
1541 = cs_matrix_get_n_columns(a) * sles->post_info->block_size;
1542 cs_real_t *r = sles->post_info->row_residual;
1543 # pragma omp parallel for if(n_vals > CS_THR_MIN)
1544 for (cs_lnum_t i = 0; i < n_vals; i++)
1545 r[i] = 0;
1546 }
1547
1548 cs_timer_stats_switch(t_top_id);
1549
1550 cs_timer_t t1 = cs_timer_time();
1551 cs_timer_counter_add_diff(&_sles_t_tot, &t0, &t1);
1552 }
1553
1554 /*----------------------------------------------------------------------------*/
1555 /*!
1556 * \brief General sparse linear system resolution.
1557 *
1558 * If no options were previously provided for the matching system,
1559 * default options will be used.
1560 *
1561 * Note that if \ref cs_sles_setup was previously called for this
1562 * system, and \ref cs_sles_free has not been called since, the matrix
1563 * provided should be the same. The optional separation between the
1564 * two stages is intended to allow amortizing the cost of setup
1565 * over multiple solutions.
1566 *
1567 * The system is considered to have converged when
1568 * residue/r_norm <= precision, residue being the L2 norm of a.vx-rhs.
1569 *
1570 * \param[in, out] sles pointer to solver object
1571 * \param[in] a matrix
1572 * \param[in] precision solver precision
1573 * \param[in] r_norm residue normalization
1574 * \param[out] n_iter number of "equivalent" iterations
1575 * \param[out] residue residue
1576 * \param[in] rhs right hand side
1577 * \param[in, out] vx system solution
1578 * \param[in] aux_size size of aux_vectors (in bytes)
1579 * \param aux_vectors optional working area
1580 * (internal allocation if NULL)
1581 *
1582 * \return convergence state
1583 */
1584 /*----------------------------------------------------------------------------*/
1585
1586 cs_sles_convergence_state_t
cs_sles_solve(cs_sles_t * sles,const cs_matrix_t * a,double precision,double r_norm,int * n_iter,double * residue,const cs_real_t * rhs,cs_real_t * vx,size_t aux_size,void * aux_vectors)1587 cs_sles_solve(cs_sles_t *sles,
1588 const cs_matrix_t *a,
1589 double precision,
1590 double r_norm,
1591 int *n_iter,
1592 double *residue,
1593 const cs_real_t *rhs,
1594 cs_real_t *vx,
1595 size_t aux_size,
1596 void *aux_vectors)
1597 {
1598 cs_timer_t t0 = cs_timer_time();
1599
1600 if (sles->context == NULL)
1601 _cs_sles_define_default(sles->f_id, sles->name, a);
1602
1603 int t_top_id = cs_timer_stats_switch(_sles_stat_id);
1604
1605 sles->n_calls += 1;
1606
1607 assert(sles->solve_func != NULL);
1608
1609 const char *sles_name = cs_sles_base_name(sles->f_id, sles->name);
1610
1611 #if 0
1612 /* Dump linear system to file (for experimenting with external tools) */
1613 cs_matrix_dump_linear_system(a, rhs, sles_name);
1614 #endif
1615
1616 cs_sles_convergence_state_t state;
1617
1618 bool do_solve = _needs_solving(sles_name,
1619 a,
1620 sles->verbosity,
1621 precision,
1622 r_norm,
1623 residue,
1624 vx,
1625 rhs);
1626
1627 if (! do_solve) {
1628 sles->n_no_op += 1;
1629 *n_iter = 0;
1630 state = CS_SLES_CONVERGED;
1631 }
1632
1633 while (do_solve) {
1634
1635 state = sles->solve_func(sles->context,
1636 sles_name,
1637 a,
1638 sles->verbosity,
1639 precision,
1640 r_norm,
1641 n_iter,
1642 residue,
1643 rhs,
1644 vx,
1645 aux_size,
1646 aux_vectors);
1647
1648 if (state < CS_SLES_ITERATING && sles->error_func != NULL)
1649 do_solve = sles->error_func(sles,
1650 state,
1651 a,
1652 rhs,
1653 vx);
1654 else
1655 do_solve = false;
1656
1657 }
1658
1659 /* Prepare postprocessing if needed */
1660
1661 if (sles->post_info != NULL) {
1662 _ensure_alloc_post(sles, a);
1663 const cs_lnum_t n_vals
1664 = sles->post_info->n_rows * sles->post_info->block_size;
1665 _residual(n_vals,
1666 a,
1667 rhs,
1668 vx,
1669 sles->post_info->row_residual);
1670 }
1671
1672 /* Check error */
1673
1674 #if 0
1675 {
1676 const cs_lnum_t block_size = cs_matrix_get_diag_block_size(a)[0];
1677 const cs_lnum_t n_vals_ext = cs_matrix_get_n_columns(a) * block_size;
1678 const cs_lnum_t n_vals = cs_matrix_get_n_rows(a) * block_size;
1679
1680 cs_real_t *resr = NULL;
1681 BFT_MALLOC(resr, n_vals_ext, cs_real_t);
1682
1683 _residual(n_vals, a, rhs, vx, resr);
1684 cs_real_t residue = sqrt(cs_gdot(n_vals, resr, resr));
1685
1686 bft_printf("# residue[%s] = %g (%g * required, precision %g, normalization %g)\n",
1687 sles_name, residue, residue/(precision*r_norm), precision, r_norm);
1688
1689 BFT_FREE(resr);
1690 }
1691 #endif
1692
1693 cs_timer_stats_switch(t_top_id);
1694
1695 cs_timer_t t1 = cs_timer_time();
1696 cs_timer_counter_add_diff(&_sles_t_tot, &t0, &t1);
1697
1698 return state;
1699 }
1700
1701 /*----------------------------------------------------------------------------*/
1702 /*!
1703 * \brief Free sparse linear equation solver setup.
1704 *
1705 * This function frees resolution-related data, such as multigrid hierarchy,
1706 * preconditioning, and any other temporary arrays or objects required for
1707 * resolution, but maintains context information such as that used for
1708 * logging (especially performance data).
1709 *
1710 * \param[in, out] sles pointer to solver object
1711 */
1712 /*----------------------------------------------------------------------------*/
1713
1714 void
cs_sles_free(cs_sles_t * sles)1715 cs_sles_free(cs_sles_t *sles)
1716 {
1717 if (sles != NULL) {
1718 if (sles->free_func != NULL)
1719 sles->free_func(sles->context);
1720 }
1721 }
1722
1723 /*----------------------------------------------------------------------------*/
1724 /*!
1725 * \brief Copy the definition of a sparse linear equation solver to another.
1726 *
1727 * The intended use of this function is to allow associating different
1728 * solvers to related systems, so as to differentiate logging, while using
1729 * the same settings by default.
1730 *
1731 * If the source solver does not provide a \ref cs_sles_copy_t function,
1732 * No modification is done to the solver. If the copy function is available,
1733 * the context is copied, as are the matching function pointers.
1734 *
1735 * If previous settings have been defined and used, they are saved as
1736 * per \ref cs_sles_define.
1737 *
1738 * \param[in, out] dest pointer to destination solver object
1739 * \param[in] src pointer to source solver object
1740 *
1741 * \return 0 in case of success, 1 in case of failure
1742 */
1743 /*----------------------------------------------------------------------------*/
1744
1745 int
cs_sles_copy(cs_sles_t * dest,const cs_sles_t * src)1746 cs_sles_copy(cs_sles_t *dest,
1747 const cs_sles_t *src)
1748 {
1749 int retval = 1;
1750
1751 /* If no copy function is available or source does not have a
1752 context yet, we can do nothing */
1753
1754 if (src->copy_func == NULL)
1755 return retval;
1756
1757 /* Check if system was previously defined and used,
1758 and save info for future logging in this case */
1759
1760 if (dest->context != NULL) {
1761 if (dest->n_calls > 0 && dest->log_func != NULL)
1762 _save_system_info(dest);
1763 else if (dest->destroy_func != NULL)
1764 dest->destroy_func(&(dest->context));
1765 }
1766
1767 dest->type_id = src->type_id;
1768 dest->verbosity = src->verbosity;
1769
1770 /* Now define options */
1771 dest->context = src->copy_func(src->context);
1772 dest->setup_func = src->setup_func;
1773 dest->solve_func = src->solve_func;
1774 dest->free_func = src->free_func;
1775 dest->log_func = src->log_func;
1776 dest->copy_func = src->copy_func;
1777 dest->destroy_func = src->destroy_func;
1778
1779 if (dest->context != NULL)
1780 retval = 0;
1781
1782 return retval;
1783 }
1784
1785 /*----------------------------------------------------------------------------*/
1786 /*!
1787 * \brief Associate a convergence error handler to a given sparse linear
1788 * equation solver.
1789 *
1790 * The error will be called whenever convergence fails. To dissassociate
1791 * the error handler, this function may be called with \p handler = NULL.
1792 *
1793 * The association will only be successful if the matching solver
1794 * has already been defined.
1795 *
1796 * \param[in, out] sles pointer to solver object
1797 * \param[in] error_handler_func pointer to convergence error
1798 * handler function
1799 */
1800 /*----------------------------------------------------------------------------*/
1801
1802 void
cs_sles_set_error_handler(cs_sles_t * sles,cs_sles_error_handler_t * error_handler_func)1803 cs_sles_set_error_handler(cs_sles_t *sles,
1804 cs_sles_error_handler_t *error_handler_func)
1805 {
1806 if (sles != NULL)
1807 sles->error_func = error_handler_func;
1808 }
1809
1810 /*----------------------------------------------------------------------------*/
1811 /*!
1812 * \brief Return pointer to default sparse linear solver definition function.
1813 *
1814 * The associated function will be used to provide a definition when
1815 * \ref cs_sles_setup or \ref cs_sles_solve is used for a system for which no
1816 * matching call to \ref cs_sles_define has been done.
1817 *
1818 * \return define_func pointer to default definition function
1819 */
1820 /*----------------------------------------------------------------------------*/
1821
1822 cs_sles_define_t *
cs_sles_get_default_define(void)1823 cs_sles_get_default_define(void)
1824 {
1825 return _cs_sles_define_default;
1826 }
1827
1828 /*----------------------------------------------------------------------------*/
1829 /*!
1830 * \brief Set default sparse linear solver definition function.
1831 *
1832 * The provided function will be used to provide a definition when
1833 * \ref cs_sles_setup or \ref cs_sles_solve is used for a system for which no
1834 * matching call to \ref cs_sles_define has been done.
1835 *
1836 * \param[in] define_func pointer to default definition function
1837 */
1838 /*----------------------------------------------------------------------------*/
1839
1840 void
cs_sles_set_default_define(cs_sles_define_t * define_func)1841 cs_sles_set_default_define(cs_sles_define_t *define_func)
1842 {
1843 _cs_sles_define_default = define_func;
1844 }
1845
1846 /*----------------------------------------------------------------------------*/
1847 /*!
1848 * \brief Set default verbosity definition function.
1849 *
1850 * The provided function will be used to define the verbosity when
1851 * \ref cs_sles_find_or_add is called.
1852 *
1853 * \param[in] verbosity_func pointer to default verbosity function
1854 */
1855 /*----------------------------------------------------------------------------*/
1856
1857 void
cs_sles_set_default_verbosity(cs_sles_verbosity_t * verbosity_func)1858 cs_sles_set_default_verbosity(cs_sles_verbosity_t *verbosity_func)
1859 {
1860 _cs_sles_default_verbosity = verbosity_func;
1861 }
1862
1863 /*----------------------------------------------------------------------------*/
1864 /*!
1865 * \brief Output default post-processing data for failed system convergence.
1866 *
1867 * \param[in] name variable name
1868 * \param[in] mesh_id id of error output mesh, or 0 if none
1869 * \param[in] a linear equation matrix
1870 * \param[in] rhs right hand side
1871 * \param[in, out] vx current system solution
1872 */
1873 /*----------------------------------------------------------------------------*/
1874
1875 void
cs_sles_post_error_output_def(const char * name,int mesh_id,const cs_matrix_t * a,const cs_real_t * rhs,cs_real_t * vx)1876 cs_sles_post_error_output_def(const char *name,
1877 int mesh_id,
1878 const cs_matrix_t *a,
1879 const cs_real_t *rhs,
1880 cs_real_t *vx)
1881 {
1882 if (mesh_id != 0) {
1883
1884 const cs_mesh_t *mesh = cs_glob_mesh;
1885
1886 char base_name[32], val_name[32];
1887
1888 const cs_lnum_t n_cols = cs_matrix_get_n_columns(a);
1889 const cs_lnum_t n_rows = cs_matrix_get_n_rows(a);
1890 const cs_lnum_t *diag_block_size = cs_matrix_get_diag_block_size(a);
1891
1892 /* Check for mesh location */
1893
1894 int location_id = 0, flag = 0;
1895 if (n_rows == mesh->n_cells)
1896 location_id = CS_MESH_LOCATION_CELLS;
1897 else if (n_rows == mesh->n_vertices)
1898 location_id = CS_MESH_LOCATION_VERTICES;
1899 flag = location_id;
1900 cs_parall_max(1, CS_INT_TYPE, &flag);
1901 if (flag == location_id)
1902 flag = 0;
1903 else
1904 flag =1;
1905 cs_parall_max(1, CS_INT_TYPE, &flag);
1906 if (flag != 0)
1907 return;
1908
1909 /* Now generate output */
1910
1911 cs_real_t *val;
1912 BFT_MALLOC(val, n_cols*diag_block_size[1], cs_real_t);
1913
1914 for (int val_id = 0; val_id < 5; val_id++) {
1915
1916 switch(val_id) {
1917
1918 case 0:
1919 strcpy(base_name, "Diag");
1920 cs_matrix_copy_diagonal(a, val);
1921 break;
1922
1923 case 1:
1924 strcpy(base_name, "RHS");
1925 memcpy(val, rhs, n_rows*diag_block_size[1]*sizeof(cs_real_t));
1926 break;
1927
1928 case 2:
1929 strcpy(base_name, "X");
1930 memcpy(val, vx, n_rows*diag_block_size[1]*sizeof(cs_real_t));
1931 break;
1932
1933 case 3:
1934 strcpy(base_name, "Residual");
1935 _residual(n_rows*diag_block_size[1],
1936 a,
1937 rhs,
1938 vx,
1939 val);
1940 break;
1941
1942 case 4:
1943 strcpy(base_name, "Diag_Dom");
1944 cs_matrix_diag_dominance(a, val);
1945 break;
1946
1947 }
1948
1949 if (strlen(name) + strlen(base_name) < 31) {
1950 strcpy(val_name, base_name);
1951 strcat(val_name, "_");
1952 strcat(val_name, name);
1953 }
1954 else {
1955 strncpy(val_name, base_name, 31);
1956 val_name[31] = '\0';
1957 }
1958
1959 assert(diag_block_size[0] == diag_block_size[1]); /* for now */
1960
1961 cs_sles_post_output_var(val_name,
1962 mesh_id,
1963 location_id,
1964 CS_POST_WRITER_ERRORS,
1965 diag_block_size[1],
1966 val);
1967 }
1968
1969 BFT_FREE(val);
1970 }
1971 }
1972
1973 /*----------------------------------------------------------------------------*/
1974 /*!
1975 * \brief Output post-processing variable related to system convergence.
1976 *
1977 * \param[in] name variable name
1978 * \param[in] mesh_id id of error output mesh, or 0 if none
1979 * \param[in] location_id mesh location id (cells or vertices)
1980 * \param[in] writer_id id of specified associated writer, or
1981 * \ref CS_POST_WRITER_ALL_ASSOCIATED for all
1982 * \param[in] diag_block_size block size for diagonal
1983 * \param[in, out] var variable values
1984 */
1985 /*----------------------------------------------------------------------------*/
1986
1987 void
cs_sles_post_output_var(const char * name,int mesh_id,int location_id,int writer_id,int diag_block_size,cs_real_t var[])1988 cs_sles_post_output_var(const char *name,
1989 int mesh_id,
1990 int location_id,
1991 int writer_id,
1992 int diag_block_size,
1993 cs_real_t var[])
1994 {
1995 if (mesh_id != 0) {
1996
1997 int _diag_block_size[4] = {1, 1, 1, 1};
1998
1999 const cs_mesh_t *mesh = cs_glob_mesh;
2000 const cs_time_step_t *ts = cs_glob_time_step;
2001
2002 size_t n_non_norm;
2003 cs_lnum_t n_rows = 0;
2004 if (location_id == CS_MESH_LOCATION_CELLS)
2005 n_rows = mesh->n_cells;
2006 else if (location_id == CS_MESH_LOCATION_VERTICES)
2007 n_rows = mesh->n_vertices;
2008
2009 cs_real_t *val_type;
2010
2011 assert(_diag_block_size[0] == _diag_block_size[1]); /* no padding */
2012
2013 if (diag_block_size > 1) {
2014 int i;
2015 for (i = 0; i < 3; i++) /* will become false if padding is used */
2016 _diag_block_size[i] = diag_block_size;
2017 _diag_block_size[3] = diag_block_size*diag_block_size;
2018 }
2019
2020 BFT_MALLOC(val_type, _diag_block_size[1]*n_rows, cs_real_t);
2021
2022 n_non_norm = _value_type(_diag_block_size[1]*n_rows, var, val_type);
2023
2024 if (location_id == CS_MESH_LOCATION_CELLS)
2025 cs_post_write_var(mesh_id,
2026 writer_id,
2027 name,
2028 _diag_block_size[0],
2029 true, /* interlace */
2030 true, /* use parents */
2031 CS_POST_TYPE_cs_real_t,
2032 var,
2033 NULL,
2034 NULL,
2035 ts);
2036 else if (location_id == CS_MESH_LOCATION_VERTICES)
2037 cs_post_write_vertex_var(mesh_id,
2038 writer_id,
2039 name,
2040 _diag_block_size[0],
2041 true, /* interlace */
2042 true, /* use parents */
2043 CS_POST_TYPE_cs_real_t,
2044 var,
2045 ts);
2046
2047 int flag = (n_non_norm > 0) ? 1 : 0;
2048 cs_parall_max(1, CS_INT_TYPE, &flag);
2049
2050 if (flag > 0) {
2051
2052 char type_name[32];
2053 size_t l = 31 - strlen("_fp_type");
2054
2055 strncpy(type_name, name, l);
2056 type_name[31] = '\0';
2057
2058 strcat(type_name, "_fp_type");
2059
2060 if (location_id == CS_MESH_LOCATION_CELLS)
2061 cs_post_write_var(mesh_id,
2062 writer_id,
2063 type_name,
2064 _diag_block_size[0],
2065 true, /* interlace */
2066 true, /* use parents */
2067 CS_POST_TYPE_cs_real_t,
2068 val_type,
2069 NULL,
2070 NULL,
2071 ts);
2072 else if (location_id == CS_MESH_LOCATION_VERTICES)
2073 cs_post_write_vertex_var(mesh_id,
2074 writer_id,
2075 name,
2076 _diag_block_size[0],
2077 true, /* interlace */
2078 true, /* use parents */
2079 CS_POST_TYPE_cs_real_t,
2080 var,
2081 ts);
2082
2083 }
2084
2085 BFT_FREE(val_type);
2086 }
2087 }
2088
2089 /*----------------------------------------------------------------------------*/
2090 /*!
2091 * \brief Return base name associated to a field id, name couple.
2092 *
2093 * This is simply a utility function which will return its name argument
2094 * if f_id < 0, and the associated field's name or label otherwise.
2095 *
2096 * \param[in] f_id associated field id, or < 0
2097 * \param[in] name associated name if f_id < 0, or NULL
2098 *
2099 * \return pointer to base name associated to the field id, name couple
2100 */
2101 /*----------------------------------------------------------------------------*/
2102
2103 const char *
cs_sles_base_name(int f_id,const char * name)2104 cs_sles_base_name(int f_id,
2105 const char *name)
2106 {
2107 const char *sles_name = name;
2108
2109 if (f_id > -1) {
2110 const cs_field_t *f = cs_field_by_id(f_id);
2111 sles_name = cs_field_get_label(f);
2112 }
2113
2114 return sles_name;
2115 }
2116
2117 /*----------------------------------------------------------------------------*/
2118 /*!
2119 * \brief Return name associated to a field id, name couple.
2120 *
2121 * \param[in] f_id associated field id, or < 0
2122 * \param[in] name associated name if f_id < 0, or NULL
2123 *
2124 * \return pointer to name associated to the field id, name couple
2125 */
2126 /*----------------------------------------------------------------------------*/
2127
2128 const char *
cs_sles_name(int f_id,const char * name)2129 cs_sles_name(int f_id,
2130 const char *name)
2131 {
2132 const cs_sles_t *sles = cs_sles_find_or_add(f_id, name);
2133
2134 if (sles->name != NULL)
2135 return sles->name;
2136 else
2137 return cs_sles_base_name(f_id, name);
2138 }
2139
2140 /*----------------------------------------------------------------------------*/
2141
2142 END_C_DECLS
2143