programmer's documentation
Input of calculation parameters (C functions)

Introduction

C user functions for definition of model options and calculation parameters. These subroutines are called in all cases.

If the Code_Saturne GUI is used, this file is not required (but may be used to override parameters entered through the GUI, and to set parameters not accessible through the GUI).

Several functions are present in the file, each destined to defined specific parameters.

Base model related options

Definition of user variables or properties should be defined here, if not already done throught the GUI.

Most definitions should be done in the cs_user_parameters function, which is the C equivalent of the usipsu Fortran routine.

For example, to add boundary values for all scalars, the following code can be added:

int n_fields = cs_field_n_fields();
for (int f_id = 0; f_id < n_fields; f_id++) {
}

Linear solver related options

By default, Code_Saturne will use a multigrid algorithm for pressure and iterative solver for other variables. For a given case, checking the setup file resulting from a first calculation will provide more info.

Available solvers include a variety of iterative linear solvers, described in more detail at cs_sles_it_create, and a multigrid solver, whose definition and settings are described at cs_multigrid_create, cs_multigrid_set_coarsening_options, and cs_multigrid_set_solver_options.

Simple options may be set using the GUI, but for more advanced settings are described in this section. It is also recommended to read the documentation of cs_sles.c (which is a solver definition "container"), cs_sles_it.c (iterative solvers, with available types cs_sles_it_type_t), and cs_multigrid.c (which are actual solver implementations). The API provided is extensible, so it is possible for a user to define other solvers or link to external solver libraries using this system, without requiring any modification to non-user source files.

The examples which follow illustrate mostly simple setting changes which may be useful.

Example: distance to wall

By default, the wall distance (active only with turbulence models which require it) is computed with a preconditionned conjugate gradient. The following example shows how to use a multigrid solver for this quantity (useful especially if computed repeatedly, such as for ALE).

cs_multigrid_define(-1, "wall_distance");

Example: user variable

The following example shows how to set the linear solver for a given user variable field so as to use a BiCGStab solver with polynomial preconditioning of degree 1.

cs_field_t *cvar_user_1 = cs_field_by_name_try("user_1");
if (cvar_user_1 != NULL) {
cs_sles_it_define(cvar_user_1->id,
NULL, /* name passed is NULL if field_id > -1 */
1, /* polynomial precond. degree (default 0) */
10000); /* n_max_iter */
}

Changing the verbosity

By default, a linear solver uses the same verbosity as its matching variable, and is not verbose for non-variable quantities. The verbosity may be specifically set for linear system resolution, as shown in the following example:

cs_sles_t *sles_p = cs_sles_find_or_add(CS_F_(p)->id, NULL);

Example: advanced multigrid settings

The following example shows how to set advanced settings for the multigrid solver used for the pressure solution.

3, /* aggregation_limit (default 3) */
0, /* coarsening_type (default 0) */
10, /* n_max_levels (default 25) */
30, /* min_g_cells (default 30) */
0.95, /* P0P1 relaxation (default 0.95) */
20); /* postprocessing (default 0) */
(mg,
CS_SLES_JACOBI, /* descent smoother type (default: CS_SLES_PCG) */
CS_SLES_JACOBI, /* ascent smoother type (default: CS_SLES_PCG) */
CS_SLES_PCG, /* coarse solver type (default: CS_SLES_PCG) */
50, /* n max cycles (default 100) */
5, /* n max iter for descent (default 2) */
5, /* n max iter for asscent (default 10) */
1000, /* n max iter coarse solver (default 10000) */
0, /* polynomial precond. degree descent (default 0) */
0, /* polynomial precond. degree ascent (default 0) */
1, /* polynomial precond. degree coarse (default 0) */
-1.0, /* precision multiplier descent (< 0 forces max iters) */
-1.0, /* precision multiplier ascent (< 0 forces max iters) */
0.1); /* requested precision multiplier coarse (default 1) */

Example: multigrid preconditioner settings

The following example shows how to use multigrid as a preconditioner for a conjugate gradient solver (for the pressure correction), and set advanced settings for that multigrid preconditioner.

NULL,
-1,
10000);
assert(strcmp(cs_sles_pc_get_type(cs_sles_it_get_pc(c)), "multigrid") == 0);
(mg,
CS_SLES_P_GAUSS_SEIDEL, /* descent smoother (CS_SLES_P_GAUSS_SEIDEL) */
CS_SLES_P_GAUSS_SEIDEL, /* ascent smoother (CS_SLES_P_GAUSS_SEIDEL) */
CS_SLES_PCG, /* coarse solver (CS_SLES_P_GAUSS_SEIDEL) */
1, /* n max cycles (default 1) */
1, /* n max iter for descent (default 1) */
1, /* n max iter for asscent (default 1) */
500, /* n max iter coarse solver (default 1) */
0, /* polynomial precond. degree descent (default) */
0, /* polynomial precond. degree ascent (default) */
0, /* polynomial precond. degree coarse (default 0) */
-1.0, /* precision multiplier descent (< 0 forces max iters) */
-1.0, /* precision multiplier ascent (< 0 forces max iters) */
1.0); /* requested precision multiplier coarse (default 1) */

Multigrid parallel settings

In parallel, grids may optionally be merged across neigboring ranks when their local size becomes too small. This tends to deepen the grid hierarchy, as some parallel rank boundaries are removed. Depending on the architecture and processor/network performance ratios, this may increase or decrease performance.

cs_grid_set_merge_options(4, /* # of ranks merged at a time */
300, /* mean # of cells under which we merge */
500, /* global # of cells under which we merge */
1); /* # of ranks under which we do not merge */

Example: DOM radiation settings

For DOM radiation models, 1 solver is assigned for each direction this allows using a specific ordering for each direction for the default Block Gauss-Seidel solver.

The example below shows how to set a non-default linear solver for DOM radiation. Here, we assume a quadrature with 32 directions is used (if more solvers than directions are specified, the extra definitions will be unused, but this causes no further issues).

for (int i = 0; i < 32; i++) {
char name[16];
sprintf(name, "radiation_%03d", i+1);
name,
0, /* poly_degree */
1000); /* n_max_iter */
}

Plotting solver convergence

The following example shows how to activate convergence plotting for built-in iterative or multigrid solvers.

const cs_field_t *f = CS_F_(p);
cs_sles_t *sles_p = cs_sles_find_or_add(f->id, NULL);
bool use_iteration = true; /* use iteration or wall clock time for axis */
if (strcmp(cs_sles_get_type(sles_p), "cs_sles_it_t") == 0) {
cs_sles_it_set_plot_options(c, f->name, use_iteration);
}
else if (strcmp(cs_sles_get_type(sles_p), "cs_multigrid_t") == 0) {
cs_multigrid_set_plot_options(c, f->name, use_iteration);
}

Plots will appear as CSV (comma-separated value) files in the monitoring subdirectory.

Using PETSc

The following example shows how to setup a solver to use the PETSc library, if the code was built with PETSc support.

General options (those passed to PETSc through command line options) may be defined directly in cs_user_linear_solvers, for example:

/* Initialization must be called before setting options;
it does not need to be called before calling
cs_sles_petsc_define(), as this is handled automatically. */
PETSC_COMM_WORLD = cs_glob_mpi_comm;
PetscInitializeNoArguments();
/* See the PETSc documentation for the options database */
PetscOptionsSetValue("-ksp_type", "cg");
PetscOptionsSetValue("-pc_type", "jacobi");

A specific system may be set up to use PETsc, as is shown here for the pressure variable:

NULL,
MATSHELL,
_petsc_p_setup_hook,
NULL);

The basic matrix format to be used by PETSc is defined at this stage, using a PETSc MatType string (see PETSc documentation). Further options may be defined in a setup hook function, as follows:

static void
_petsc_p_setup_hook(const void *context,
KSP ksp)
{
PC pc;
KSPSetType(ksp, KSPCG); /* Preconditioned Conjugate Gradient */
KSPSetNormType(ksp, KSP_NORM_UNPRECONDITIONED); /* Try to have "true" norm */
KSPGetPC(ksp, &pc);
PCSetType(pc, PCJACOBI); /* Jacobi (diagonal) preconditioning */
}

If no additional settings are required, the matching parameter in cs_sles_petsc_define may be set to NULL.

Using PETSc with GAMG

To use PETSc's GAMG (geometric-algebraic multigrid) preconditioner, the following general options should be set:

/* Initialization must be called before setting options;
it does not need to be called before calling
cs_sles_petsc_define(), as this is handled automatically. */
PETSC_COMM_WORLD = cs_glob_mpi_comm;
PetscInitializeNoArguments();
/* See the PETSc documentation for the options database */
PetscOptionsSetValue("-ksp_type", "cg");
PetscOptionsSetValue("-pc_type", "gamg");
PetscOptionsSetValue("-pc_gamg_agg_nsmooths", "1");
PetscOptionsSetValue("-mg_levels_ksp_type", "richardson");
PetscOptionsSetValue("-mg_levels_pc_type", "sor");
PetscOptionsSetValue("-mg_levels_ksp_max_it", "1");
PetscOptionsSetValue("-pc_gamg_threshold", "0.02");
PetscOptionsSetValue("-pc_gamg_reuse_interpolation", "TRUE");
PetscOptionsSetValue("-pc_gamg_square_graph", "4");

Setting GAMG-preconditioned PCG for the pressure may be done as in the previous option, ensuring here that a matrix structure native to PETSc is used (SEQAIJ, MPIAIJ, or some other AIJ variant), so all required matrix operations are available:

NULL,
MATMPIAIJ,
_petsc_p_setup_hook_gamg,
NULL);

With the associated setup hook function:

static void
_petsc_p_setup_hook_gamg(const void *context,
KSP ksp)
{
PC pc;
KSPSetType(ksp, KSPCG); /* Preconditioned Conjugate Gradient */
KSPGetPC(ksp, &pc);
PCSetType(pc, PCGAMG); /* GAMG (geometric-algebraic multigrid)
preconditioning */
}

Using PETSc with HYPRE

To use HYPRE's Boomer AMG as a PETSc preconditioner, the following general options should be set:

/* Initialization must be called before setting options;
it does not need to be called before calling
cs_sles_petsc_define(), as this is handled automatically. */
PETSC_COMM_WORLD = cs_glob_mpi_comm;
PetscInitializeNoArguments();
/* See the PETSc documentation for the options database */
PetscOptionsSetValue("-ksp_type", "cg");
PetscOptionsSetValue("-pc_type", "hypre");
PetscOptionsSetValue("-pc_hypre_type","boomeramg");
PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type","HMIS");
PetscOptionsSetValue("-pc_hypre_boomeramg_interp_type","ext+i-cc");
PetscOptionsSetValue("-pc_hypre_boomeramg_agg_nl","2");
PetscOptionsSetValue("-pc_hypre_boomeramg_P_max","4");
PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold","0.5");
PetscOptionsSetValue("-pc_hypre_boomeramg_no_CF","");

Setting BoomerAMG-preconditioned PCG for the pressure may be done as in the previous option, ensuring here that a matrix structure native to PETSc is used (SEQAIJ, MPIAIJ, or some other AIJ variant), so all required matrix operations are available:

NULL,
MATMPIAIJ,
_petsc_p_setup_hook_bamg,
NULL);

With the associated setup hook function:

static void
_petsc_p_setup_hook_bamg(const void *context,
KSP ksp)
{
PC pc;
KSPSetType(ksp, KSPCG); /* Preconditioned Conjugate Gradient */
KSPGetPC(ksp, &pc);
PCSetType(pc, PCHYPRE); /* HYPRE BoomerAMG preconditioning */
}

Additional PETSc features

Many additional features are possible with PETSc; for example, the following setup hook also outputs a view of the matrix, depending on an environment variable, CS_USER_PETSC_MAT_VIEW, which may take values DEFAULT, DRAW_WORLD, or DRAW:

static void
_petsc_p_setup_hook_view(const void *context,
KSP ksp)
{
PC pc;
const char *p = getenv("CS_USER_PETSC_MAT_VIEW");
if (p != NULL) {
/* Get system and preconditioner matrixes */
Mat a, pa;
KSPGetOperators(ksp, &a, &pa);
/* Output matrix in several ways depending on
CS_USER_PETSC_MAT_VIEW environment variable */
if (strcmp(p, "DEFAULT") == 0)
MatView(a, PETSC_VIEWER_DEFAULT);
else if (strcmp(p, "DRAW_WORLD") == 0)
MatView(a, PETSC_VIEWER_DRAW_WORLD);
else if (strcmp(p, "DRAW") == 0) {
PetscViewer viewer;
PetscDraw draw;
PetscViewerDrawOpen(PETSC_COMM_WORLD, NULL, "PETSc View",
0, 0, 600, 600, &viewer);
PetscViewerDrawGetDraw(viewer, 0, &draw);
PetscViewerDrawSetPause(viewer, -1);
MatView(a, viewer);
PetscDrawPause(draw);
PetscViewerDestroy(&viewer);
}
}
}

Time moment related options

Code_Saturne allows the calculation of temporal means or variances, either of expressions evaluated through a user function, or of expressions of the type $<f_1*f_2...*f_n>$. The variables may be fields or field components. This is done calling either through the GUI, or in the user function cs_user_time_moments. Each temporal mean is declared using either cs_time_moment_define_by_func, or cs_time_moment_define_by_field_ids.

For each time moment, a starting time step or value is defined. If the starting time step number is negative, the time value is used instead.

The moment values are stored as fields, and updated at each time step, using recursive formulas. Before the matching moment computation starting time step, a moment's values are uniformly 0. For visualization an interpretation reasons, only fields of dimension 1, 3, 6, or 9 (scalars, vectors, or tensors of rank 2) are allowed, so moment definitions not matching this constraint should be split.

To count defined moments, use the cs_time_moment_n_moments function, whether from Fortran or C. To access the matching fields, use time_moment_field_id in Fortran, or cs_time_moment_get_field in C.

Examples

Example 1

In the following example, we define a moment for the mean velocity. All components are used (component -1 means all components), so the moment is a vector.

int moment_f_id[] = {CS_F_(u)->id};
int moment_c_id[] = {-1};
int n_fields = 1;
n_fields,
moment_f_id,
moment_c_id,
1000, /* nt_start */
-1, /* t_start */
NULL);

Example 2

In the next example, we multiply the expression by the density. As the density is of dimension 1, and the velocity of dimension 3, the resulting moment is of dimension 3.

{
/* Moment <U> calculated starting from time step 1000. */
int moment_f_id[] = {CS_F_(rho)->id, CS_F_(u)->id};
int moment_c_id[] = {-1, -1};
int n_fields = 2;
n_fields,
moment_f_id,
moment_c_id,
1000, /* nt_start */
-1, /* t_start */
NULL);
}

Example 3

In the next example, we define a product of several field components, all of dimension 1, as we consider only the x and y components of the velocity; for the density, we cas use either component 0 or -1 (all), since the field is scalar.

This moment's computation is also restarted at each time step.

int moment_f_id[] = {CS_F_(rho)->id, CS_F_(u)->id, CS_F_(u)->id};
int moment_c_id[] = {-1, 0, 1};
int n_fields = 3;
n_fields,
moment_f_id,
moment_c_id,
-1, /* nt_start */
20.0, /* t_start */
NULL);

Example 4

This next example illustrates the use of user-defined functions to evaluate expressions. Here, we compute the moment of the sum ot two variables (which obviously cannot be expressed as a product), so we first need to define an appropriate function, matching the signature of a cs_time_moment_data_t function. We can name that function as we choose, so naming for clarity is recommmended. Note that in this case, the input argument is not used. This argument may be useful to pass data to the function, or distinguish between calls to a same function.

Note also that we compute both means and variances here.

static void
_simple_data_sum(const void *input,
cs_real_t *vals)
{
const int location_id = CS_MESH_LOCATION_CELLS;
const cs_lnum_t n_elts = cs_mesh_location_get_n_elts(location_id)[0];
const cs_real_t *s1 = cs_field_by_name("species_1")->val;
const cs_real_t *s2 = cs_field_by_name("species_2")->val;
for (cs_lnum_t i = 0; i < n_elts; i++) {
vals[i] = s1[i] + s2[i];
}
}

In cs_user_time_moments, we can not assign that function to a moments definition:

const char *sum_comp_name[] = {"species_sum_mean", "species_sum_variance"};
for (int i = 0; i < 2; i++) {
1, /* field dimension */
_simple_data_sum, /* data_func */
NULL, /* data_input */
NULL, /* w_data_func */
NULL, /* w_data_input */
m_type[i],
1000, /* nt_start */
-1, /* t_start */
NULL);
}

Example 5

In this last example, we compute components of the mean velocity in the case of a rotating mesh. As the mesh orientation changes at each time step, it is necessary to compensate for this rotation when computing the mean, relative to a given mesh position. When using the matching moment, it will also be necessary to account for the mesh position.

Here, the same function will be called for each component, so an input array is defined, with a different key (here a simple character) used for each call.

static void
_velocity_moment_data(const void *input,
cs_real_t *vals)
{
const char key = *((const char *)input);
const int location_id = CS_MESH_LOCATION_CELLS;
const cs_lnum_t n_elts = cs_mesh_location_get_n_elts(location_id)[0];
const cs_real_3_t *vel = (const cs_real_3_t *)(CS_F_(u)->val);
const cs_real_3_t *restrict cell_cen
double omgnrm = fabs(rot->omega);
/* Axial, tangential and radial unit vectors */
cs_real_3_t e_ax = {rot->axis[0], rot->axis[1], rot->axis[2]};
for (cs_lnum_t i = 0; i < n_elts; i++) {
cs_rotation_velocity(rot, cell_cen[i], e_th);
double xnrm = sqrt(cs_math_3_square_norm(e_th));
e_th[0] /= xnrm;
e_th[1] /= xnrm;
e_th[2] /= xnrm;
cs_rotation_coriolis_v(rot, -1., e_th, e_r);
xnrm = sqrt(cs_math_3_square_norm(e_r));
e_r[0] /= xnrm;
e_r[1] /= xnrm;
e_r[2] /= xnrm;
/* Radius */
cs_real_t xr = cs_math_3_dot_product(cell_cen[i], e_r);
/* Axial, tangential and radial components of velocity */
cs_real_t xva = vel[i][0]*e_ax[0] + vel[i][1]*e_ax[1] + vel[i][2]*e_ax[2];
cs_real_t xvt = vel[i][0]*e_th[0] + vel[i][1]*e_th[1] + vel[i][2]*e_th[2];
cs_real_t xvr = vel[i][0]*e_r[0] + vel[i][1]*e_r[1] + vel[i][2]*e_r[2];
/* Entrainment velocity is removed */
xvt -= omgnrm*xr;
/* Store value */
if (key == 'r')
vals[i] = xvr;
else if (key == 't')
vals[i] = xvt;
else if (key == 'a')
vals[i] = xva;
}
}

Note that the input arrays must be accessible when updating moments at each time step, so the array of inputs is declared static in cs_user_time_moments. Fo more complex inputs, we would have an array of inputs here; in this simple case, we could pass a simple call id as the input, casting from point to integer.

const char *vel_comp_name[] = {"Wr_moy", "Wt,moy", "Wa_moy"};
/* Data input must be "static" so it can be used in later calls */
static char vel_comp_input[3] = {'r', 't', 'a'};
for (int comp_id = 0; comp_id < 3; comp_id++) {
cs_time_moment_define_by_func(vel_comp_name[comp_id],
1,
_velocity_moment_data, /* data_func */
&(vel_comp_input[comp_id]), /* data_input */
NULL, /* w_data_func */
NULL, /* w_data_input */
74000, /* nt_start */
-1, /* t_start */
NULL);
}

Fan modeling options

Code_Saturne allows modeling of some circular fans as volume regions, defined by simple geometric characteristics, and modeled as explicit momentum source terms in those regions.

Fan pressure characteristic curves are defined as a 2nd order polynomial, and a torque may also be specified. For correct results, it is important that the mesh match the fan dimensions and placement (thickness, hub, blades, and total radius).

The following example shows how a fan may be defined:

{
const cs_real_t inlet_axis_coords[3] = {0, 0, 0};
const cs_real_t outlet_axis_coords[3] = {0.1, 0, 0};
const cs_real_t fan_radius = 0.7;
const cs_real_t blades_radius = 0.5;
const cs_real_t hub_radius = 0.1;
const cs_real_t pressure_curve_coeffs[3] = {0.6, -0.1, -0.05};
const cs_real_t axial_torque = 0.01;
cs_fan_define(3, /* fan (mesh) dimension */
inlet_axis_coords,
outlet_axis_coords,
fan_radius,
blades_radius,
hub_radius,
pressure_curve_coeffs,
axial_torque);
}