1
/*============================================================================
3
* This file is part of the Code_Saturne Kernel, element of the
4
* Code_Saturne CFD tool.
6
* Copyright (C) 1998-2009 EDF S.A., France
8
* contact: saturne-support@edf.fr
10
* The Code_Saturne Kernel is free software; you can redistribute it
11
* and/or modify it under the terms of the GNU General Public License
12
* as published by the Free Software Foundation; either version 2 of
13
* the License, or (at your option) any later version.
15
* The Code_Saturne Kernel is distributed in the hope that it will be
16
* useful, but WITHOUT ANY WARRANTY; without even the implied warranty
17
* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18
* GNU General Public License for more details.
20
* You should have received a copy of the GNU General Public License
21
* along with the Code_Saturne Kernel; if not, write to the
22
* Free Software Foundation, Inc.,
23
* 51 Franklin St, Fifth Floor,
24
* Boston, MA 02110-1301 USA
26
*============================================================================*/
28
/*============================================================================
30
*============================================================================*/
32
#if defined(HAVE_CONFIG_H)
33
#include "cs_config.h"
36
/*----------------------------------------------------------------------------
37
* Standard C library headers
38
*----------------------------------------------------------------------------*/
47
#if defined(__STDC_VERSION__) /* size_t */
48
#if (__STDC_VERSION__ == 199901L)
61
/*----------------------------------------------------------------------------
63
*----------------------------------------------------------------------------*/
66
#include <bft_error.h>
67
#include <bft_printf.h>
68
#include <bft_timer.h>
70
/*----------------------------------------------------------------------------
72
*----------------------------------------------------------------------------*/
76
/*----------------------------------------------------------------------------
78
*----------------------------------------------------------------------------*/
83
#include "cs_matrix.h"
85
#include "cs_mesh_quantities.h"
90
/*----------------------------------------------------------------------------
91
* Header for the current file
92
*----------------------------------------------------------------------------*/
94
#include "cs_multigrid.h"
96
/*----------------------------------------------------------------------------*/
100
/*=============================================================================
101
* Local Macro Definitions
102
*============================================================================*/
104
#define EPZERO 1.E-12
105
#define RINFIN 1.E+30
107
#if !defined(HUGE_VAL)
108
#define HUGE_VAL 1.E+12
112
/*=============================================================================
113
* Local Type Definitions
114
*============================================================================*/
116
/*----------------------------------------------------------------------------
117
* Local Structure Definitions
118
*----------------------------------------------------------------------------*/
120
/* Basic per linear system options and logging */
121
/*---------------------------------------------*/
123
typedef struct _cs_multigrid_info_t {
125
char *name; /* System name */
126
cs_sles_type_t type[3]; /* Descent/ascent smoother
129
unsigned n_builds; /* Number of times grids built */
130
unsigned n_solves; /* Number of times system solved */
132
unsigned long long n_levels_tot; /* Total accumulated number of
134
unsigned n_levels[3]; /* Number of grid levels:
137
unsigned n_iterations[3][4]; /* Number of iterations for
140
[finest, coarsest, total,
141
fine grid equivalent] */
142
unsigned long long n_iterations_tot[4]; /* Total accumulated number of
144
[finest, coarsest, total,
145
fine grid equivalent] */
147
double wt_tot[2]; /* Total wall-clock time used:
149
double cpu_tot[2]; /* Total (local) CPU used:
152
} cs_multigrid_info_t;
157
typedef struct _cs_multigrid_t {
159
cs_multigrid_info_t info; /* Multigrid info */
161
int n_levels; /* Current number of grid levels */
162
int n_levels_max; /* Maximum number of grid levels */
164
int n_levels_post; /* Current number of postprocessed levels */
165
int post_cell_max; /* If > 0, activates postprocessing of
166
coarsening, projecting coarse cell
167
numbers (modulo post_cell_max)
170
cs_grid_t **grid_hierarchy; /* Array of grid pointers */
172
int **post_cell_num; /* If post_cell_max > 0, array of
173
(n_levels - 1) arrays of projected
174
coarse cell numbers on the base grid */
178
/*============================================================================
180
*============================================================================*/
182
static int cs_glob_multigrid_n_systems = 0; /* Current number of systems */
183
static int cs_glob_multigrid_n_max_systems = 0; /* Max. number of sytems for
184
cs_glob_mgrid_systems. */
186
/* System info array */
188
static cs_multigrid_t **cs_glob_multigrid_systems = NULL;
190
/*============================================================================
191
* Private function definitions
192
*============================================================================*/
194
/*----------------------------------------------------------------------------
195
* Initialize multigrid info structure.
198
* name <-- system name
199
* info <-- pointer to multigrid info structure
200
*----------------------------------------------------------------------------*/
203
_multigrid_info_init(cs_multigrid_info_t *info,
208
BFT_MALLOC(info->name, strlen(name) + 1, char);
210
strcpy(info->name, name);
212
for (i = 0; i < 3; i++)
213
info->type[i] = CS_SLES_N_TYPES;
218
info->n_levels_tot = 0;
220
for (i = 0; i < 3; i++) {
221
info->n_levels[i] = 0;
222
for (j = 0; j < 4; j++)
223
info->n_iterations[i][j] = 0;
226
for (i = 0; i < 4; i++)
227
info->n_iterations_tot[i] = 0;
229
for (i = 0; i < 2; i++) {
230
info->wt_tot[i] = 0.0;
231
info->cpu_tot[i] = 0.0;
235
/*----------------------------------------------------------------------------
236
* Destroy multigrid info structure.
239
* this_info <-> pointer to linear system info structure pointer
240
*----------------------------------------------------------------------------*/
243
_multigrid_info_unset(cs_multigrid_info_t *this_info)
245
assert(this_info != NULL);
247
BFT_FREE(this_info->name);
250
/*----------------------------------------------------------------------------
251
* Output information regarding multigrid resolution.
254
* this_info <-> pointer to multigrid info structure
255
*----------------------------------------------------------------------------*/
258
_multigrid_info_dump(const cs_multigrid_info_t *this_info)
260
unsigned long long n_builds_denom = CS_MAX(this_info->n_builds, 1);
261
unsigned long long n_solves_denom = CS_MAX(this_info->n_solves, 1);
262
int n_builds = this_info->n_builds;
263
int n_solves = this_info->n_solves;
264
int n_lv_min = this_info->n_levels[1];
265
int n_lv_max = this_info->n_levels[2];
266
int n_lv_mean = (int)(this_info->n_levels_tot / n_builds_denom);
267
int n_it_f_min = this_info->n_iterations[1][0];
268
int n_it_f_max = this_info->n_iterations[2][0];
269
int n_it_c_min = this_info->n_iterations[1][1];
270
int n_it_c_max = this_info->n_iterations[2][1];
271
int n_it_t_min = this_info->n_iterations[1][2];
272
int n_it_t_max = this_info->n_iterations[2][2];
273
int n_it_e_min = this_info->n_iterations[1][3];
274
int n_it_e_max = this_info->n_iterations[2][3];
275
int n_it_f_mean = (int)(this_info->n_iterations_tot[0] / n_solves_denom);
276
int n_it_c_mean = (int)(this_info->n_iterations_tot[1] / n_solves_denom);
277
int n_it_t_mean = (int)(this_info->n_iterations_tot[2] / n_solves_denom);
278
int n_it_e_mean = (int)(this_info->n_iterations_tot[3] / n_solves_denom);
281
"Summary of multigrid for \"%s\":\n\n"),
284
if (this_info->type[0] != CS_SLES_N_TYPES) {
286
const char *descent_smoother_name = cs_sles_type_name[this_info->type[0]];
287
const char *ascent_smoother_name = cs_sles_type_name[this_info->type[1]];
289
if (this_info->type[0] == this_info->type[1])
290
bft_printf(_(" Smoother: %s\n"), _(descent_smoother_name));
292
bft_printf(_(" Descent smoother: %s\n"
293
" Ascent smoother: %s\n"),
294
_(descent_smoother_name), _(ascent_smoother_name));
296
bft_printf(_(" Coarsest level solver: %s\n"),
297
_(cs_sles_type_name[this_info->type[2]]));
301
bft_printf(_(" Number of constructions: %d\n"
302
" Number of resolutions: %d\n"
303
" Number of levels:\n"
307
" Number of iterations:\n"
312
" on coarsest grid:\n"
320
" equivalent (total weighted by number of cells) :\n"
324
" Associated times (construction, resolution)\n"
325
" total elapsed: %12.3f %12.3f\n"),
326
n_builds, n_solves, n_lv_min, n_lv_max, n_lv_mean,
327
n_it_f_min, n_it_f_max, n_it_f_mean,
328
n_it_c_min, n_it_c_max, n_it_c_mean,
329
n_it_t_min, n_it_t_max, n_it_t_mean,
330
n_it_e_min, n_it_e_max, n_it_e_mean,
331
this_info->wt_tot[0], this_info->wt_tot[1]);
333
#if defined(HAVE_MPI)
335
if (cs_glob_n_ranks > 1) {
337
double cpu_min[2], cpu_max[2], cpu_tot[2], cpu_loc[2];
338
cpu_loc[0] = this_info->cpu_tot[0];
339
cpu_loc[1] = this_info->cpu_tot[1];
341
MPI_Allreduce(cpu_loc, cpu_min, 2, MPI_DOUBLE, MPI_MIN,
343
MPI_Allreduce(cpu_loc, cpu_max, 2, MPI_DOUBLE, MPI_MAX,
345
MPI_Allreduce(cpu_loc, cpu_tot, 2, MPI_DOUBLE, MPI_SUM,
348
bft_printf(_(" Min local total CPU time: %12.3f %12.3f\n"
349
" Max local total CPU time: %12.3f %12.3f\n"
350
" Total CPU time: %12.3f %12.3f\n"),
351
cpu_min[0], cpu_min[1], cpu_max[0], cpu_max[1],
352
cpu_tot[0], cpu_tot[1]);
358
if (cs_glob_n_ranks == 1)
359
bft_printf(_(" Total CPU time: %12.3f %12.3f\n"),
360
this_info->cpu_tot[0], this_info->cpu_tot[1]);
363
/*----------------------------------------------------------------------------
364
* Initialize multigrid info structure.
367
* name <-- system name
370
* pointer to newly created multigrid structure
371
*----------------------------------------------------------------------------*/
373
static cs_multigrid_t *
374
_multigrid_create(const char *name)
379
BFT_MALLOC(mg, 1, cs_multigrid_t);
381
_multigrid_info_init(&(mg->info), name);
384
mg->n_levels_max = 10;
386
mg->n_levels_post = 0;
387
mg->post_cell_max = 0;
389
BFT_MALLOC(mg->grid_hierarchy, mg->n_levels_max, cs_grid_t *);
391
for (ii = 0; ii < mg->n_levels_max; ii++)
392
mg->grid_hierarchy[ii] = NULL;
394
mg->post_cell_num = NULL;
399
/*----------------------------------------------------------------------------
400
* Destroy multigrid structure.
403
* mg <-> pointer multigrid structure pointer
404
*----------------------------------------------------------------------------*/
407
_multigrid_destroy(cs_multigrid_t **mg)
410
cs_multigrid_t *_mg = *mg;
414
_multigrid_info_unset(&(_mg->info));
416
for (ii = 0; ii < _mg->n_levels_max; ii++)
417
cs_grid_destroy(_mg->grid_hierarchy + ii);
419
if (_mg->post_cell_max > 0) {
420
for (ii = 0; ii < _mg->n_levels_max - 1; ii++)
421
if (_mg->post_cell_num[ii] != NULL)
422
BFT_FREE(_mg->post_cell_num[ii]);
423
BFT_FREE(_mg->post_cell_num);
426
BFT_FREE(_mg->grid_hierarchy);
431
/*----------------------------------------------------------------------------
432
* Add grid to multigrid structure hierarchy.
435
* mg <-- multigrid structure
436
* grid <-- grid to add
437
*----------------------------------------------------------------------------*/
440
_multigrid_add_level(cs_multigrid_t *mg,
445
/* Reallocate arrays if necessary */
447
if (mg->n_levels == mg->n_levels_max) {
449
if (mg->n_levels_max == 0)
450
mg->n_levels_max = 10;
451
mg->n_levels_max *= 2;
453
BFT_REALLOC(mg->grid_hierarchy, mg->n_levels_max, cs_grid_t *);
454
for (ii = mg->n_levels; ii < mg->n_levels_max; ii++)
455
mg->grid_hierarchy[ii] = NULL;
457
if (mg->post_cell_num != NULL) {
458
BFT_REALLOC(mg->post_cell_num, mg->n_levels_max, int *);
459
for (ii = mg->n_levels; ii < mg->n_levels_max; ii++)
460
mg->post_cell_num[ii] = NULL;
461
if (mg->n_levels > 0)
462
mg->post_cell_num[mg->n_levels - 1] = NULL;
466
mg->grid_hierarchy[mg->n_levels] = grid;
470
/*----------------------------------------------------------------------------
471
* Get multigrid structure's id in list of know systems
474
* mg <-- multigrid structure
477
* id (0 to n-1) of structure in list of know systems, or -1 otherwise
478
*----------------------------------------------------------------------------*/
481
_multigrid_id(const cs_multigrid_t *mg)
485
for (id = 0; id < cs_glob_multigrid_n_systems; id++) {
486
if (mg == cs_glob_multigrid_systems[id])
490
if (id >= cs_glob_multigrid_n_systems)
496
/*----------------------------------------------------------------------------
497
* Return pointer to linear system info.
499
* If this system did not previously exist, it is added to the list of
503
* name <-- system name
504
*----------------------------------------------------------------------------*/
506
static cs_multigrid_t *
507
_find_or_add_system(const char *name)
509
int ii, start_id, end_id, mid_id;
512
/* Use binary search to find system */
515
end_id = cs_glob_multigrid_n_systems - 1;
516
mid_id = start_id + ((end_id -start_id) / 2);
518
while (start_id <= end_id) {
519
cmp_ret = strcmp((cs_glob_multigrid_systems[mid_id])->info.name, name);
521
start_id = mid_id + 1;
522
else if (cmp_ret > 0)
526
mid_id = start_id + ((end_id -start_id) / 2);
529
/* If found, return */
532
return cs_glob_multigrid_systems[mid_id];
534
/* Reallocate global array if necessary */
536
if (cs_glob_multigrid_n_systems >= cs_glob_multigrid_n_max_systems) {
538
if (cs_glob_multigrid_n_max_systems == 0)
539
cs_glob_multigrid_n_max_systems = 10;
541
cs_glob_multigrid_n_max_systems *= 2;
542
BFT_REALLOC(cs_glob_multigrid_systems,
543
cs_glob_multigrid_n_max_systems,
546
for (ii = cs_glob_multigrid_n_systems;
547
ii < cs_glob_multigrid_n_max_systems;
549
cs_glob_multigrid_systems[ii] = NULL;
553
/* Insert in sorted list */
555
for (ii = cs_glob_multigrid_n_systems; ii > mid_id; ii--)
556
cs_glob_multigrid_systems[ii] = cs_glob_multigrid_systems[ii - 1];
558
cs_glob_multigrid_systems[mid_id] = _multigrid_create(name);
559
cs_glob_multigrid_n_systems += 1;
561
return cs_glob_multigrid_systems[mid_id];
564
/*----------------------------------------------------------------------------
565
* Add postprocessing info to multigrid hierarchy
568
* mg <-> multigrid structure
569
* n_base_cells <-- number of cells in base grid
570
*----------------------------------------------------------------------------*/
573
_multigrid_add_post(cs_multigrid_t *mg,
574
fvm_lnum_t n_base_cells)
580
if (mg->post_cell_max < 1)
583
mg->n_levels_post = mg->n_levels - 1;
585
assert(mg->n_levels_post < mg->n_levels_max);
587
/* Reallocate arrays if necessary */
589
if (mg->post_cell_num == NULL) {
590
BFT_MALLOC(mg->post_cell_num, mg->n_levels_max, int *);
591
for (ii = 0; ii < mg->n_levels_max; ii++)
592
mg->post_cell_num[ii] = NULL;
595
for (ii = 0; ii < mg->n_levels_post; ii++) {
596
BFT_REALLOC(mg->post_cell_num[ii], n_base_cells, int);
597
cs_grid_project_cell_num(mg->grid_hierarchy[ii+1],
600
mg->post_cell_num[ii]);
604
/*----------------------------------------------------------------------------
605
* Post process variables associated with Syrthes couplings
608
* hierarchy_id <-- Id of multigrid hierarchy
609
* nt_cur_abs <-- Current time step
610
* t_cur_abs <-- Current time value
611
*----------------------------------------------------------------------------*/
614
_cs_multigrid_post_function(cs_int_t hierarchy_id,
620
char *var_name = NULL;
621
cs_multigrid_t *mg = NULL;
622
const char name_prefix[] = "mg";
623
const char *base_name = NULL;
625
/* Return if necessary structures inconsistant or have been destroyed */
627
if (hierarchy_id < cs_glob_multigrid_n_systems)
628
mg = cs_glob_multigrid_systems[hierarchy_id];
633
if (mg->post_cell_num == NULL || cs_post_mesh_exists(-1) != true)
636
/* Allocate name buffer */
638
base_name = mg->info.name;
639
name_len = strlen(name_prefix) + 1 + strlen(base_name) + 1 + 3 + 1 + 4 + 1;
640
BFT_MALLOC(var_name, name_len, char);
642
/* Loop on grid levels */
644
for (ii = 0; ii < mg->n_levels_post; ii++) {
646
sprintf(var_name, "%s %s %3d %2d",
647
name_prefix, base_name, (int)(ii+1), (int)nt_cur_abs);
649
cs_post_write_var(-1,
657
mg->post_cell_num[ii],
661
BFT_FREE(mg->post_cell_num[ii]);
664
mg->n_levels_post = 0;
669
/*----------------------------------------------------------------------------
670
* Compute dot product, summing result over all ranks.
673
* n_elts <-- Local number of elements
674
* x <-- first vector in s = x.y
675
* y <-- second vector in s = x.y
679
*----------------------------------------------------------------------------*/
682
_dot_product(cs_int_t n_elts,
686
double s = cblas_ddot(n_elts, x, 1, y, 1);
688
#if defined(HAVE_MPI)
690
if (cs_glob_n_ranks > 1) {
692
MPI_Allreduce(&s, &_sum, 1, MPI_DOUBLE, MPI_SUM, cs_glob_mpi_comm);
696
#endif /* defined(HAVE_MPI) */
701
/*----------------------------------------------------------------------------
702
* Test if convergence is attained.
705
* var_name <-- Variable name
706
* n_f_cells <-- Number of cells on fine mesh
707
* n_max_cycles <-- Maximum number of cycles
708
* cycle_id <-- Number of current cycle
710
* verbosity <-- Verbosity level
711
* n_iters <-- Number of iterations
712
* precision <-- Precision limit
713
* r_norm <-- Residue normalization
714
* residue <-> Residue
715
* rhs --> Right-hand side
718
* 1 if converged, 0 if not converged, -1 if not converged and maximum
719
* cycle number reached, -2 if divergence is detected.
720
*----------------------------------------------------------------------------*/
723
_convergence_test(const char *var_name,
734
const char cycle_h_fmt[]
735
= N_(" ---------------------------------------------------\n"
736
" n. | Cumulative iterations | Norm. residual\n"
737
" cycles | on fine mesh | on fine mesh\n"
738
" ---------------------------------------------------\n");
739
const char cycle_t_fmt[]
740
= N_(" ---------------------------------------------------\n");
741
const char cycle_cv_fmt[]
742
= N_(" %4d | %6d | %12.4e\n");
744
const char cycle_fmt[]
745
= N_(" N. cycles: %4d; Fine mesh cumulative iter: %5d; "
746
"Norm. residual %12.4e\n");
748
static double initial_residue = 0.;
750
/* Compute residue */
752
*residue = sqrt(_dot_product(n_f_cells, rhs, rhs));
755
initial_residue = *residue;
757
if (*residue < precision*r_norm) {
760
bft_printf(_(cycle_fmt), cycle_id, n_iters, *residue/r_norm);
761
else if (verbosity > 2) {
762
bft_printf(_(cycle_h_fmt));
763
bft_printf(_(cycle_cv_fmt),
764
cycle_id, n_iters, *residue/r_norm);
765
bft_printf(_(cycle_t_fmt));
770
else if (cycle_id >= n_max_cycles) {
774
bft_printf(_(cycle_fmt), cycle_id, n_iters, *residue/r_norm);
775
else if (verbosity > 1) {
776
bft_printf(_(cycle_h_fmt));
777
bft_printf(_(cycle_fmt),
778
cycle_id, n_iters, *residue/r_norm);
779
bft_printf(_(cycle_t_fmt));
781
bft_printf(_(" @@ Warning: algebraic multigrid for [%s]\n"
783
" Maximum number of cycles (%d) reached.\n"),
784
var_name, n_max_cycles);
793
bft_printf(_(cycle_fmt), cycle_id, n_iters, *residue/r_norm);
795
if (*residue > initial_residue * 10000.0 && *residue > 100.)
798
#if (_CS_STDC_VERSION >= 199901L)
799
if (isnan(*residue) || isinf(*residue))
807
/*----------------------------------------------------------------------------
808
* Handle error output in case if divergence.
810
* Depending on the multigrid level at which divergence is detected,
811
* output matrix hierarchy information, RHS, and variable info,
812
* and finally abort on error.
815
* mg <-- Multigrid system
816
* level <-- Multigrid level at which divergence is detected
817
* symmetric <-- Symmetric coefficients indicator
818
* rotation_mode <-- Halo update option for rotational periodicity
819
* cycle_id <-- Id of currect cycle
820
* initial_residue <-- Initial residue
821
* residue <-- Residue
822
* rhs <-- Right hand side
823
* vx <-- System solution
824
* c_rhs <-- Right hand side for levels > 0
825
* c_vx <-- System solution for levels > 0
826
*----------------------------------------------------------------------------*/
829
_abort_on_divergence(cs_multigrid_t *mg,
832
cs_perio_rota_t rotation_mode,
834
double initial_residue,
836
const cs_real_t rhs[],
841
int mesh_id = cs_post_init_error_writer_cells();
848
cs_real_t *var = NULL;
850
const cs_real_t *_da = NULL, *_xa = NULL;
851
const cs_grid_t *g = mg->grid_hierarchy[0];
852
const fvm_lnum_t n_base_cells = cs_grid_get_n_cells(g);
854
BFT_MALLOC(var, cs_grid_get_n_cells_ext(g), cs_real_t);
856
/* Output info on main level */
858
cs_grid_get_matrix(g, &_da, &_xa, NULL);
860
cs_sles_post_error_output_def(mg->info.name,
869
/* Output diagonal and diagonal dominance for all coarse levels */
871
for (lv_id = 1; lv_id < mg->n_levels; lv_id++) {
873
g = mg->grid_hierarchy[lv_id];
875
cs_grid_get_matrix(g, &_da, NULL, NULL);
877
cs_grid_project_var(g, n_base_cells, _da, var);
878
sprintf(var_name, "Diag_%04d", lv_id);
879
cs_sles_post_error_output_var(var_name, mesh_id, var);
881
cs_grid_project_diag_dom(g, n_base_cells, var);
882
sprintf(var_name, "Diag_Dom_%04d", lv_id);
883
cs_sles_post_error_output_var(var_name, mesh_id, var);
886
/* Output info on current level if > 0 */
891
fvm_lnum_t n_cells = 0;
892
fvm_lnum_t n_cells_ext = 0;
894
cs_real_t *c_res = NULL;
895
cs_matrix_t *_matrix = NULL;
897
g = mg->grid_hierarchy[level];
898
n_cells = cs_grid_get_n_cells(g);
899
n_cells_ext = cs_grid_get_n_cells_ext(g);
901
cs_grid_project_var(g, n_base_cells, c_rhs[level], var);
902
sprintf(var_name, "RHS_%04d", level);
903
cs_sles_post_error_output_var(var_name, mesh_id, var);
905
cs_grid_project_var(g, n_base_cells, c_vx[level], var);
906
sprintf(var_name, "X_%04d", level);
907
cs_sles_post_error_output_var(var_name, mesh_id, var);
909
/* Compute residual */
911
BFT_MALLOC(c_res, n_cells_ext, cs_real_t);
913
cs_grid_get_matrix(g, &_da, &_xa, &_matrix);
914
cs_matrix_set_coefficients(_matrix, symmetric, _da, _xa);
916
cs_matrix_vector_multiply(rotation_mode, _matrix, c_vx[level], c_res);
918
for (ii = 0; ii < n_cells; ii++)
919
c_res[ii] = fabs(c_res[ii] - c_rhs[level][ii]);
921
cs_grid_project_var(g, n_base_cells, c_res, var);
925
sprintf(var_name, "Residual_%04d", level);
926
cs_sles_post_error_output_var(var_name, mesh_id, var);
935
bft_error(__FILE__, __LINE__, 0,
936
_("algebraic multigrid [%s]: divergence after %d cycles:\n"
937
" initial residual: %11.4e; current residual: %11.4e"),
938
_(mg->info.name), cycle_id, initial_residue, residue);
940
bft_error(__FILE__, __LINE__, 0,
941
_("algebraic multigrid [%s]: divergence after %d cycles\n"
942
" during resolution at level %d:\n"
943
" initial residual: %11.4e; current residual: %11.4e"),
944
_(mg->info.name), cycle_id, level, initial_residue, residue);
947
/*----------------------------------------------------------------------------
948
* Sparse linear system resolution using multigrid.
951
* mg <-- Multigrid system
952
* descent_smoother_type <-- Type of smoother for descent (PCG, Jacobi, ...)
953
* ascent_smoother_type <-- Type of smoother for ascent (PCG, Jacobi, ...)
954
* coarse_solver_type <-- Type of solver (PCG, Jacobi, ...)
955
* symmetric <-- Symmetric coefficients indicator
956
* poly_degree <-- Preconditioning polynomial degree (0: diagonal)
957
* rotation_mode <-- Halo update option for rotational periodicity
958
* verbosity <-- Verbosity level
959
* cycle_id <-- Id of currect cycle
960
* n_max_cycles <-- Maximum number of cycles
961
* n_max_iter <-- Maximum number of iterations per grid level
962
* n_max_iter[level * 2] for descent
963
* n_max_iter[level * 2 + 1] for ascent
964
* precision <-- Precision limit
965
* r_norm <-- Residue normalization
966
* n_level_iter <-> Number of iterations per level
967
* residue <-> Residue
968
* rhs <-- Right hand side
969
* vx --> System solution
970
* aux_size <-- Number of elements in aux_vectors
971
* aux_vectors --- Optional working area (allocation otherwise)
974
* 1 if converged, 0 if not converged, -1 if not converged and maximum
975
* cycle number reached, -2 if divergence is detected.
976
*----------------------------------------------------------------------------*/
979
_multigrid_cycle(cs_multigrid_t *mg,
980
cs_sles_type_t descent_smoother_type,
981
cs_sles_type_t ascent_smoother_type,
982
cs_sles_type_t coarse_solver_type,
985
cs_perio_rota_t rotation_mode,
989
const int n_max_iter[],
994
const cs_real_t *rhs,
999
int level, coarsest_level;
1002
int cvg = 0, c_cvg = 0;
1004
size_t alloc_size = 0;
1005
cs_real_t c_precision = precision;
1006
cs_real_t _residue = -1.;
1008
size_t _aux_size = aux_size;
1009
fvm_lnum_t n_cells = 0, n_cells_ext = 0;
1010
cs_real_t r_norm_l = r_norm;
1012
char _var_lv_name[33];
1013
char *var_lv_name = _var_lv_name;
1015
double initial_residue = *residue;
1016
double _initial_residue = 0.;
1018
cs_real_t *_aux_vectors = aux_vectors;
1019
cs_real_t *wr = NULL;
1020
cs_real_t *_rhs_vx_val = NULL;
1021
cs_real_t **_rhs_vx = NULL, **_rhs = NULL, **_vx = NULL;
1022
cs_matrix_t *_matrix;
1023
cs_multigrid_info_t *mg_info = NULL;
1025
const char *var_name = NULL;
1026
const cs_real_t *_rhs_level = NULL;
1027
const cs_real_t *_da = NULL, *_xa = NULL;
1028
const cs_grid_t *f = NULL, *c= NULL;
1030
cs_bool_t end_cycle = false;
1032
/* Initialization */
1034
mg_info = &(mg->info);
1035
var_name = mg_info->name;
1037
coarsest_level = mg->n_levels - 1;
1039
/* In theory, one should increase precision on coarsest mesh,
1040
but in practice, it more efficient to have a lower precision */
1041
/* c_precision = precision * 0.01; */
1042
c_precision = precision;
1044
f = mg->grid_hierarchy[0];
1054
if (strlen(var_name) + 5 > 32)
1055
BFT_MALLOC(var_lv_name, strlen(var_name) + 5 + 1, char);
1057
/* Allocate wr or use working area */
1059
if (aux_size >= (size_t)n_cells_ext) {
1061
_aux_vectors = wr + n_cells_ext;
1062
_aux_size = aux_size - n_cells_ext;
1065
BFT_MALLOC(wr, n_cells_ext, cs_real_t);
1067
/* reserve memory for rhs and vx;
1068
for the finest level, simply point to input and output arrays */
1070
BFT_MALLOC(_rhs_vx, mg->n_levels*2, cs_real_t *);
1072
_vx = _rhs_vx + mg->n_levels;
1074
_rhs[0] = NULL; /* Use _rhs_level when necessary to avoid const warning */
1077
/* Reserve memory for corrections and residues for coarse levels */
1079
if (mg->n_levels > 1) {
1083
for (level = 1; level < mg->n_levels; level++)
1084
alloc_size += cs_grid_get_n_cells_ext(mg->grid_hierarchy[level]);
1086
BFT_MALLOC(_rhs_vx_val, alloc_size*2, cs_real_t);
1088
_rhs[1] = _rhs_vx_val;
1089
_vx[1] = _rhs_vx_val + alloc_size;
1091
for (level = 2; level < mg->n_levels; level++) {
1092
fvm_lnum_t _n_cells_ext_prev
1093
= cs_grid_get_n_cells_ext(mg->grid_hierarchy[level-1]);
1094
_rhs[level] = _rhs[level - 1] + _n_cells_ext_prev;
1095
_vx[level] = _vx[level - 1] + _n_cells_ext_prev;
1103
bft_printf(_(" Multigrid cycle: descent\n"));
1105
for (level = 0; level < coarsest_level; level++) {
1107
int _poly_degree = (level == 0) ? poly_degree : 0;
1109
_rhs_level = (level == 0) ? rhs : _rhs[level];
1111
sprintf(var_lv_name, "%s:%04d", var_name, level);
1113
c = mg->grid_hierarchy[level+1];
1118
bft_printf(_(" level %3d: smoother\n"), level);
1120
cs_grid_get_matrix(f, &_da, &_xa, &_matrix);
1122
_initial_residue = _residue;
1124
c_cvg = cs_sles_solve(var_lv_name,
1125
descent_smoother_type,
1126
false, /* Stats not updated here */
1135
n_max_iter[level*2],
1146
_abort_on_divergence(mg, level,
1147
symmetric, rotation_mode, cycle_id,
1148
_initial_residue, _residue,
1149
rhs, vx, _rhs, _vx);
1151
n_level_iter[level] += n_iter;
1154
TODO: get residue from cs_sles_solve(). This optimisation would
1155
require adding an argument and exercising caution to ensure the
1156
correct sign and meaning of the residue. */
1158
cs_matrix_set_coefficients(_matrix, symmetric, _da, _xa);
1160
cs_matrix_vector_multiply(rotation_mode,
1165
for (ii = 0; ii < n_cells; ii++)
1166
wr[ii] = _rhs_level[ii] - wr[ii];
1168
n_level_iter[level] += 1;
1170
/* Convergence test in beginning of cycle (fine mesh) */
1174
cvg = _convergence_test(var_name,
1185
/* If converged or cycle limit reached, break from descent loop */
1189
_abort_on_divergence(mg, level,
1190
symmetric, rotation_mode, cycle_id,
1191
initial_residue, *residue,
1192
rhs, vx, _rhs, _vx);
1199
/* Prepare for next level */
1201
cs_grid_restrict_cell_var(f, c, wr, _rhs[level+1]);
1213
for (ii = 0; ii < n_cells; ii++) /* Initialize correction */
1214
_vx[level+1][ii] = 0.0;
1216
} /* End of loop on levels (descent) */
1218
if (end_cycle == false) {
1220
/* Resolve coarsest level to convergence */
1221
/*---------------------------------------*/
1224
bft_printf(_(" Resolution on coarsest level\n"));
1226
assert(level = coarsest_level);
1227
assert(c == mg->grid_hierarchy[coarsest_level]);
1229
/* coarsest level == 0 should never happen, but we play it safe */
1230
_rhs_level = (level == 0) ? rhs : _rhs[coarsest_level];
1232
sprintf(var_lv_name, "%s:%04d", var_name, coarsest_level);
1234
cs_grid_get_matrix(c, &_da, &_xa, &_matrix);
1236
_initial_residue = _residue;
1238
c_cvg = cs_sles_solve(var_lv_name,
1240
false, /* Stats not updated here */
1246
0, /* poly_degree */
1249
n_max_iter[level*2],
1260
_abort_on_divergence(mg, level,
1261
symmetric, rotation_mode, cycle_id,
1262
_initial_residue, _residue,
1263
rhs, vx, _rhs, _vx);
1265
n_level_iter[level] += n_iter;
1271
bft_printf(_(" Multigrid cycle: ascent\n"));
1273
for (level = coarsest_level - 1; level > -1; level--) {
1275
cs_real_t *_f_vx = _vx[level];
1277
c = mg->grid_hierarchy[level+1];
1278
f = mg->grid_hierarchy[level];
1288
/* Prolong correction */
1290
cs_grid_prolong_cell_var(c, f, _vx[level+1], wr);
1292
for (ii = 0; ii < n_cells; ii++)
1293
_f_vx[ii] += wr[ii];
1295
/* Smoother pass if level > 0
1296
(smoother not called for finest mesh, as it will be called in
1297
descent phase of the next cycle, before the convergence test). */
1302
bft_printf(_(" level %3d: smoother\n"), level);
1304
sprintf(var_lv_name, "%s:%04d", var_name, level);
1306
cs_grid_get_matrix(f, &_da, &_xa, &_matrix);
1308
_initial_residue = _residue;
1310
c_cvg = cs_sles_solve(var_lv_name,
1311
ascent_smoother_type,
1312
false, /* Stats not updated here */
1318
0, /* poly_degree */
1321
n_max_iter[level*2 + 1],
1332
_abort_on_divergence(mg, level,
1333
symmetric, rotation_mode, cycle_id,
1334
_initial_residue, _residue,
1335
rhs, vx, _rhs, _vx);
1337
n_level_iter[level] += n_iter;
1340
} /* End loop on levels (ascent) */
1342
} /* End of test on end_cycle */
1346
if (var_lv_name != _var_lv_name)
1347
BFT_FREE(var_lv_name);
1349
if (wr != aux_vectors)
1353
BFT_FREE(_rhs_vx_val);
1358
/*----------------------------------------------------------------------------
1359
* Sparse linear system resolution using multigrid.
1362
* var_name <-- Variable name
1363
* descent_smoother_type <-- Type of smoother for descent (PCG, Jacobi, ...)
1364
* ascent_smoother_type <-- Type of smoother for ascent (PCG, Jacobi, ...)
1365
* coarse_solver_type <-- Type of solver (PCG, Jacobi, ...)
1366
* symmetric <-- Symmetric coefficients indicator
1367
* poly_degree <-- Preconditioning polynomial degree (0: diagonal)
1368
* rotation_mode <-- Halo update option for rotational periodicity
1369
* verbosity <-- Verbosity level
1370
* n_max_cycles <-- Maximum number of cycles
1371
* n_max_iter_descent <-- Maximum nb. of iterations for descent phases
1372
* n_max_iter_ascent <-- Maximum nb. of iterations for ascent phases
1373
* n_max_iter_coarse <-- Maximum nb. of iterations for coarsest solution
1374
* precision <-- Precision limit
1375
* r_norm <-- Residue normalization
1376
* n_cycles --> Number of cycles
1377
* n_iter --> Number of iterations
1378
* residue <-> Residue
1379
* rhs <-- Right hand side
1380
* vx --> System solution
1381
* aux_size <-- Number of elements in aux_vectors
1382
* aux_vectors --- Optional working area (allocation otherwise)
1383
*----------------------------------------------------------------------------*/
1386
_multigrid_solve(const char *var_name,
1387
cs_sles_type_t descent_smoother_type,
1388
cs_sles_type_t ascent_smoother_type,
1389
cs_sles_type_t coarse_solver_type,
1390
cs_bool_t symmetric,
1392
cs_perio_rota_t rotation_mode,
1395
int n_max_iter_descent,
1396
int n_max_iter_ascent,
1397
int n_max_iter_coarse,
1403
const cs_real_t *rhs,
1409
unsigned _n_iter[4] = {0, 0, 0, 0};
1411
fvm_lnum_t n_cells = 0;
1413
cs_multigrid_t *mg = NULL;
1414
cs_multigrid_info_t *mg_info = NULL;
1415
double wt_start = 0.0, wt_stop = 0.0;
1416
double cpu_start = 0.0, cpu_stop = 0.0;
1418
wt_start =bft_timer_wtime();
1419
cpu_start =bft_timer_cpu_time();
1420
mg = _find_or_add_system(var_name);
1421
mg_info = &(mg->info);
1423
cs_grid_get_info(mg->grid_hierarchy[0],
1431
/* Initialize number of iterations and residue,
1432
check for immediate return,
1433
solve sparse linear system using multigrid algorithm. */
1438
if (cs_sles_needs_solving(var_name,
1446
int cycle_id = 1, cvg = 0;
1447
double it_count_num = 0.0;
1449
int *n_max_iter = NULL;
1450
int *n_level_iter = NULL;
1451
size_t _aux_size = n_cells * 6;
1452
cs_real_t *_aux_vectors = aux_vectors;
1454
BFT_MALLOC(n_max_iter, mg->n_levels * 2, int);
1455
BFT_MALLOC(n_level_iter, mg->n_levels, int);
1457
if (_aux_size <= aux_size)
1458
BFT_MALLOC(_aux_vectors, _aux_size, cs_real_t);
1460
_aux_size = aux_size;
1462
for (ii = 0; ii < mg->n_levels; ii++) {
1463
n_max_iter[ii*2] = n_max_iter_descent;
1464
n_max_iter[ii*2 + 1] = n_max_iter_ascent;
1465
n_level_iter[ii] = 0;
1467
n_max_iter[(mg->n_levels-1)*2] = n_max_iter_coarse;
1468
n_max_iter[(mg->n_levels-1)*2 + 1] = n_max_iter_coarse;
1470
if (verbosity == 2) /* More detailed headers later if > 2 */
1471
bft_printf(_("Multigrid [%s]:\n"), var_name);
1473
/* Cycle to solution */
1478
bft_printf(_("Multigrid [%s]: cycle %4d\n"),
1479
var_name, cycle_id);
1481
cvg = _multigrid_cycle(mg,
1482
descent_smoother_type,
1483
ascent_smoother_type,
1505
_n_iter[0] = n_level_iter[0];
1506
_n_iter[1] = n_level_iter[mg->n_levels - 1];
1508
for (ii = 0; ii < mg->n_levels; ii++)
1509
_n_iter[2] += n_level_iter[ii];
1511
/* Estimate "equivalent" iterations */
1513
for (ii = 0; ii < mg->n_levels; ii++) {
1514
fvm_gnum_t n_g_cells = cs_grid_get_n_g_cells(mg->grid_hierarchy[ii]);
1515
it_count_num += n_g_cells * n_level_iter[ii];
1516
_n_iter[2] += n_level_iter[ii];
1519
_n_iter[3] = (int)( it_count_num
1520
/ cs_grid_get_n_g_cells(mg->grid_hierarchy[0]));
1521
*n_iter = _n_iter[3];
1523
if (_aux_vectors != aux_vectors)
1524
BFT_FREE(_aux_vectors);
1525
BFT_FREE(n_level_iter);
1526
BFT_FREE(n_max_iter);
1529
/* Update statistics */
1531
wt_stop =bft_timer_wtime();
1532
cpu_stop =bft_timer_cpu_time();
1534
/* Update stats on number of iterations (last, min, max, total) */
1536
mg_info->type[0] = descent_smoother_type;
1537
mg_info->type[1] = ascent_smoother_type;
1538
mg_info->type[2] = coarse_solver_type;
1540
for (ii = 0; ii < 4; ii++)
1541
mg_info->n_iterations[0][ii] = _n_iter[ii];
1543
if (mg_info->n_solves > 0) {
1544
for (ii = 0; ii < 4; ii++) {
1545
if (mg_info->n_iterations[1][ii] > _n_iter[ii])
1546
mg_info->n_iterations[1][ii] = _n_iter[ii];
1547
if (mg_info->n_iterations[2][ii] < _n_iter[ii])
1548
mg_info->n_iterations[2][ii] = _n_iter[ii];
1552
for (ii = 0; ii < 4; ii++) {
1553
mg_info->n_iterations[1][ii] = _n_iter[ii];
1554
mg_info->n_iterations[2][ii] = _n_iter[ii];
1558
for (ii = 0; ii < 4; ii++)
1559
mg_info->n_iterations_tot[ii] += _n_iter[ii];
1561
/* Update number of resolutions and timing data */
1563
mg_info->n_solves += 1;
1565
mg_info->wt_tot[1] += (wt_stop - wt_start);
1566
mg_info->cpu_tot[1] += (cpu_stop - cpu_start);
1569
/*============================================================================
1570
* Public function definitions for Fortran API
1571
*============================================================================*/
1573
/*----------------------------------------------------------------------------
1574
* Build a hierarchy of meshes starting from a fine mesh, for an
1575
* ACM (Additive Corrective Multigrid) method, grouping cells at
1577
*----------------------------------------------------------------------------*/
1579
void CS_PROCF(clmlga, CLMLGA)
1581
const char *cname, /* <-- variable name */
1582
const cs_int_t *lname, /* <-- variable name length */
1583
const cs_int_t *ncelet, /* <-- Number of cells, halo included */
1584
const cs_int_t *ncel, /* <-- Number of local cells */
1585
const cs_int_t *nfac, /* <-- Number of internal faces */
1586
const cs_int_t *isym, /* <-- Symmetry indicator:
1587
1: symmetric; 2: not symmetric */
1588
cs_int_t *iagmax, /* <-> Maximum agglomeration count */
1589
const cs_int_t *nagmax, /* <-- Agglomeration count limit */
1590
const cs_int_t *ncpost, /* <-- If > 0, postprocess coarsening, using
1591
coarse cell numbers modulo ncpost */
1592
const cs_int_t *iwarnp, /* <-- Verbosity level */
1593
const cs_int_t *ngrmax, /* <-- Maximum number of grid levels */
1594
const cs_int_t *ncegrm, /* <-- Maximum local number of cells on
1596
const cs_real_t *dam, /* <-- Matrix diagonal */
1597
const cs_real_t *xam /* <-- Matrix extra-diagonal terms */
1601
double wt_start, wt_stop, cpu_start, cpu_stop;
1603
fvm_lnum_t n_cells = 0;
1604
fvm_lnum_t n_faces = 0;
1605
fvm_gnum_t n_g_cells = 0;
1606
fvm_gnum_t n_g_cells_prev = 0;
1608
cs_int_t grid_lv = 0;
1609
cs_multigrid_t *mg = NULL;
1611
const cs_mesh_t *mesh = cs_glob_mesh;
1612
const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1614
cs_grid_t *g = NULL;
1616
cs_bool_t symmetric = (*isym == 1) ? true : false;
1618
assert(*ncelet >= *ncel);
1621
/* Initialization */
1623
wt_start = bft_timer_wtime();
1624
cpu_start = bft_timer_cpu_time();
1626
var_name = cs_base_string_f_to_c_create(cname, *lname);
1628
mg = _find_or_add_system(var_name);
1631
bft_printf(_("\n Construction of grid hierarchy for \"%s\"\n"),
1634
/* Destroy previous hierarchy if necessary */
1636
if (mg->n_levels > 0) {
1637
for (grid_lv = mg->n_levels - 1; grid_lv > -1; grid_lv--)
1638
cs_grid_destroy(mg->grid_hierarchy + grid_lv);
1642
/* Build coarse grids hierarchy */
1643
/*------------------------------*/
1645
g = cs_grid_create_from_shared(mesh->n_cells,
1646
mesh->n_cells_with_ghosts,
1651
mesh->i_face_numbering,
1658
_multigrid_add_level(mg, g); /* Assign to hierarchy */
1660
n_cells = mesh->n_cells;
1661
n_faces = mesh->n_i_faces;
1662
n_g_cells = mesh->n_g_cells;
1666
n_g_cells_prev = n_g_cells;
1668
/* Recursion test */
1670
if (n_g_cells <= (fvm_gnum_t)(*ncegrm))
1673
else if (grid_lv >= *ngrmax) {
1674
cs_base_warn(__FILE__, __LINE__);
1675
bft_printf(_(" CLMLGA: maximum number of coarse grids (%d)\n"
1676
" reached for \"%s\".\n"),
1677
(int)(*ngrmax), var_name);
1681
/* Build coarser grid from previous grid */
1686
bft_printf(_("\n building level %2d grid\n"), grid_lv);
1688
g = cs_grid_coarsen(g, *iwarnp, *nagmax, iagmax);
1698
_multigrid_add_level(mg, g); /* Assign to hierarchy */
1700
/* Print coarse mesh stats */
1704
#if defined(HAVE_MPI)
1706
if (cs_glob_n_ranks > 1) {
1708
int lcount[2], gcount[2];
1709
int n_c_min, n_c_max, n_f_min, n_f_max;
1711
lcount[0] = n_cells; lcount[1] = n_faces;
1712
MPI_Allreduce(lcount, gcount, 2, MPI_INT, MPI_MAX,
1714
n_c_max = gcount[0]; n_f_max = gcount[1];
1716
lcount[0] = n_cells; lcount[1] = n_faces;
1717
MPI_Allreduce(lcount, gcount, 2, MPI_INT, MPI_MIN,
1719
n_c_min = gcount[0]; n_f_min = gcount[1];
1722
(_(" total min max\n"
1723
" number of cells: %12llu %10d %10d\n"
1724
" number of faces: %10d %10d\n"),
1725
(unsigned long long)n_g_cells, n_c_min, n_c_max, n_f_min, n_f_max);
1730
if (cs_glob_n_ranks == 1)
1731
bft_printf(_(" number of cells: %10d\n"
1732
" number of faces: %10d\n"),
1733
(int)n_cells, (int)n_faces);
1737
/* If too few cells were grouped, we stop at this level */
1739
if ( n_g_cells > (0.8 * n_g_cells_prev)
1740
|| n_g_cells < (1.5 * cs_glob_n_ranks))
1744
cs_base_string_f_to_c_free(&var_name);
1746
/* Print final info */
1750
(_(" number of coarse grids: %d\n"
1751
" number of cells in coarsest grid: %llu\n\n"),
1752
grid_lv, (unsigned long long)n_g_cells);
1754
/* Prepare preprocessing info if necessary */
1758
if (mg->post_cell_max == 0) {
1759
int mg_id = _multigrid_id(mg);
1761
cs_post_add_time_dep_var(_cs_multigrid_post_function, mg_id);
1762
mg->post_cell_max = *ncpost;
1765
_multigrid_add_post(mg, mesh->n_cells);
1771
mg->info.n_levels_tot += grid_lv;
1773
mg->info.n_levels[0] = grid_lv;
1775
if (mg->info.n_builds > 0) {
1776
if (mg->info.n_levels[0] < mg->info.n_levels[1])
1777
mg->info.n_levels[1] = mg->info.n_levels[0];
1778
if (mg->info.n_levels[0] > mg->info.n_levels[2])
1779
mg->info.n_levels[2] = mg->info.n_levels[0];
1782
mg->info.n_levels[1] = mg->info.n_levels[0];
1783
mg->info.n_levels[2] = mg->info.n_levels[0];
1786
mg->info.n_builds += 1;
1790
wt_stop = bft_timer_wtime();
1791
cpu_stop = bft_timer_cpu_time();
1793
mg->info.wt_tot[0] += (wt_stop - wt_start);
1794
mg->info.cpu_tot[0] += (cpu_stop - cpu_start);
1797
/*----------------------------------------------------------------------------
1798
* Destroy a hierarchy of meshes starting from a fine mesh, keeping
1799
* the corresponding system and postprocessing information for future calls.
1800
*----------------------------------------------------------------------------*/
1802
void CS_PROCF(dsmlga, DSMLGA)
1804
const char *cname, /* <-- variable name */
1805
const cs_int_t *lname /* <-- variable name length */
1810
double wt_start, wt_stop, cpu_start, cpu_stop;
1811
cs_multigrid_t *mg = NULL;
1813
/* Initialization */
1815
wt_start =bft_timer_wtime();
1816
cpu_start =bft_timer_cpu_time();
1818
var_name = cs_base_string_f_to_c_create(cname, *lname);
1820
mg = _find_or_add_system(var_name);
1822
cs_base_string_f_to_c_free(&var_name);
1824
/* Destroy grid hierarchy */
1826
if (mg->n_levels > 0) {
1827
for (ii = mg->n_levels - 1; ii > -1; ii--)
1828
cs_grid_destroy(mg->grid_hierarchy + ii);
1834
wt_stop = bft_timer_wtime();
1835
cpu_stop = bft_timer_cpu_time();
1837
mg->info.wt_tot[0] += (wt_stop - wt_start);
1838
mg->info.cpu_tot[0] += (cpu_stop - cpu_start);
1841
/*----------------------------------------------------------------------------
1842
* General sparse linear system resolution
1843
*----------------------------------------------------------------------------*/
1845
void CS_PROCF(resmgr, RESMGR)
1847
const char *cname, /* <-- variable name */
1848
const cs_int_t *lname, /* <-- variable name length */
1849
const cs_int_t *ncelet, /* <-- Number of cells, halo included */
1850
const cs_int_t *ncel, /* <-- Number of local cells */
1851
const cs_int_t *nfac, /* <-- Number of faces */
1852
const cs_int_t *isym, /* <-- Symmetry indicator:
1853
1: symmetric; 2: not symmetric */
1854
const cs_int_t *iresds, /* <-- Descent smoother type:
1855
0: pcg; 1: Jacobi; 2: cg-stab */
1856
const cs_int_t *iresas, /* <-- Ascent smoother type:
1857
0: pcg; 1: Jacobi; 2: cg-stab */
1858
const cs_int_t *ireslp, /* <-- Coarse Resolution type:
1859
0: pcg; 1: Jacobi; 2: cg-stab */
1860
const cs_int_t *ipol, /* <-- Preconditioning polynomial degree
1862
const cs_int_t *ncymxp, /* <-- Max number of cycles */
1863
const cs_int_t *nitmds, /* <-- Max number of iterations for descent */
1864
const cs_int_t *nitmas, /* <-- Max number of iterations for ascent */
1865
const cs_int_t *nitmap, /* <-- Max number of iterations for
1866
coarsest solution */
1867
const cs_int_t *iinvpe, /* <-- Indicator to cancel increments
1868
in rotational periodicity (2) or
1869
to exchange them as scalars (1) */
1870
const cs_int_t *iwarnp, /* <-- Verbosity level */
1871
cs_int_t *ncyclf, /* --> Number of cycles done */
1872
cs_int_t *niterf, /* --> Number of iterations done */
1873
const cs_real_t *epsilp, /* <-- Precision for iterative resolution */
1874
const cs_real_t *rnorm, /* <-- Residue normalization */
1875
cs_real_t *residu, /* --> Final non normalized residue */
1876
const cs_int_t *ifacel, /* <-- Face -> cell connectivity */
1877
const cs_real_t *rhs, /* <-- System right-hand side */
1878
cs_real_t *vx /* <-> System solution */
1882
cs_sles_type_t type[4] = {CS_SLES_PCG,
1887
int _iresds = *iresds;
1888
int _iresas = *iresas;
1889
int _ireslp = *ireslp;
1891
cs_bool_t symmetric = (*isym == 1) ? true : false;
1892
cs_perio_rota_t rotation_mode = CS_PERIO_ROTA_COPY;
1894
assert(*ncelet >= *ncel);
1896
assert(ifacel != NULL);
1899
rotation_mode = CS_PERIO_ROTA_RESET;
1900
else if (*iinvpe == 3)
1901
rotation_mode = CS_PERIO_ROTA_IGNORE;
1903
var_name = cs_base_string_f_to_c_create(cname, *lname);
1905
assert(*iresds > -1 && *iresds < 3);
1906
assert(*iresas > -1 && *iresas < 3);
1907
assert(*ireslp > -1 && *ireslp < 3);
1909
if (_iresds < 0 || _iresds > 2)
1911
if (_iresas < 0 || _iresas > 2)
1913
if (_ireslp < 0 || _ireslp > 2)
1916
_multigrid_solve(var_name,
1938
cs_base_string_f_to_c_free(&var_name);
1941
/*============================================================================
1942
* Public function definitions
1943
*============================================================================*/
1945
/*----------------------------------------------------------------------------
1946
* Initialize multigrid solver API.
1949
* post_cell_max <-- If > 0, activates postprocessing of coarseninsg,
1950
* projecting coarse cell numbers (modulo post_cell_max)
1952
*----------------------------------------------------------------------------*/
1955
cs_multigrid_initialize(void)
1959
/*----------------------------------------------------------------------------
1960
* Finalize multigrid solver API.
1961
*----------------------------------------------------------------------------*/
1964
cs_multigrid_finalize(void)
1968
/* Print system info */
1970
for (ii = 0; ii < cs_glob_multigrid_n_systems; ii++)
1971
_multigrid_info_dump(&((cs_glob_multigrid_systems[ii])->info));
1973
/* Free multigrid structures */
1975
for (ii = 0; ii < cs_glob_multigrid_n_systems; ii++)
1976
_multigrid_destroy(cs_glob_multigrid_systems + ii);
1978
BFT_FREE(cs_glob_multigrid_systems);
1980
cs_glob_multigrid_n_systems = 0;
1981
cs_glob_multigrid_n_max_systems = 0;
1984
/*----------------------------------------------------------------------------*/
added
removed
src/base/cs_multigrid.c
