7.0
general documentation
Input of calculation parameters (C functions in cs_user_parameters.c)

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 as well as choices of physical models should be done here, if not already done through the GUI.

Choose a turbulent model among the available models

/* Example: Chose a turbulence model
* CS_TURB_NONE: no turbulence model (laminar flow)
* CS_TURB_MIXING_LENGTH: mixing length model
* CS_TURB_K_EPSILON: standard k-epsilon model
* CS_TURB_K_EPSILON_LIN_PROD: k-epsilon model with
* Linear Production (LP) correction
* CS_TURB_K_EPSILON_LS: Launder-Sharma low Re k-epsilon model
* CS_TURB_K_EPSILON_QUAD: Baglietto et al. low Re
* CS_TURB_RIJ_EPSILON_LRR: Rij-epsilon (LRR)
* CS_TURB_RIJ_EPSILON_SSG: Rij-epsilon (SSG)
* CS_TURB_RIJ_EPSILON_EBRSM: Rij-epsilon (EBRSM)
* CS_TURB_LES_SMAGO_CONST: LES (constant Smagorinsky model)
* CS_TURB_LES_SMAGO_DYN: LES ("classical" dynamic Smagorisky model)
* CS_TURB_LES_WALE: LES (WALE)
* CS_TURB_V2F_PHI: v2f phi-model
* CS_TURB_V2F_BL_V2K: v2f BL-v2-k
* CS_TURB_K_OMEGA: k-omega SST
* CS_TURB_SPALART_ALLMARAS: Spalart-Allmaras model */
cs_turb_model_t * cs_get_glob_turb_model(void)
Provide write access to turbulence model structure.
Definition: cs_turbulence_model.c:1432
@ CS_TURB_K_EPSILON_LIN_PROD
Definition: cs_turbulence_model.h:57
Turbulence model general options descriptor.
Definition: cs_turbulence_model.h:114
int iturb
Definition: cs_turbulence_model.h:116

Coupled solver for Rij components (when iturb=30, 31 or 32): 0 to switch off, 1 to switch on

/* Example: Coupled solver for Rij components (when iturb=30, 31 or 32)
* 0: switch off
* 1: switch on (default)
*/
rans_model->irijco = 1;
/* Advanced re-initialization for EBRSM or k-omega models
- 0: switch off
- 1: switch on (default)
*/
rans_model->reinit_turb = 1;
/* Turbulent diffusion model for second moment closure (CS_TURB_RIJ_*)
- 0: scalar diffusivity (Shir model)
- 1: tensorial diffusivity (Daly and Harlow model, default model)
*/
rans_model->idirsm = 0;
cs_turb_rans_model_t * cs_get_glob_turb_rans_model(void)
Provide access to cs_glob_turb_rans_model.
Definition: cs_turbulence_model.c:1540
RANS turbulence model descriptor.
Definition: cs_turbulence_model.h:173
int idirsm
Definition: cs_turbulence_model.h:185
int irijco
Definition: cs_turbulence_model.h:211
int reinit_turb
Definition: cs_turbulence_model.h:208

To set a thermal model (0: none, 1: temperature, 2: entyhalpy, 3: total energy)

/* Example: Choose a thermal model
*
* 0: none
* 1: temperature
* 2: enthalpy
* 3: total energy (only for compressible module)
*
* For temperature, the temperature scale may be set later using itpscl
* (1 for Kelvin, 2 for Celsius).
*
* Warning: When using specific physics, this value is
* set automatically by the physics model.
*
*/
cs_thermal_model_t * cs_get_glob_thermal_model(void)
Definition: cs_thermal_model.c:233
@ CS_THERMAL_MODEL_TEMPERATURE
Definition: cs_thermal_model.h:56
Thermal model descriptor.
Definition: cs_thermal_model.h:74
int itherm
Definition: cs_thermal_model.h:76

Volume of Fluid model with mass transfer Merkle model to take into account vaporization / condensation

cs_vof_parameters_t * cs_get_glob_vof_parameters(void)
Definition: cs_vof.c:410
#define CS_VOF_ENABLED
Definition: cs_vof.h:58
#define CS_VOF_MERKLE_MASS_TRANSFER
Definition: cs_vof.h:64
unsigned vof_model
Definition: cs_vof.h:80
VOF model parameters. Void fraction variable tracks fluid 2.
Definition: cs_vof.h:78

To activate ALE (Arbitrary Lagrangian Eulerian) method (CS_ALE_NONE: switch off, CS_ALE_LEGACY: legacy solver, CS_ALE_CDO: CDO solver)

/* Example: activate ALE (Arbitrary Lagrangian Eulerian) method
* CS_ALE_NONE: switch off
* CS_ALE_LEGACY: legacy solver
* CS_ALE_CDO: CDO solver
*/
int cs_glob_ale
Definition: cs_ale.c:81
@ CS_ALE_LEGACY
Definition: cs_ale.h:56

The user can add a scalar to be solved

/* Example: add 2 scalar variables ("species" in the GUI nomenclature).
*
* Note that at this (very early) stage of the data setup, fields are
* not defined yet. Associated fields will be defined later (after
* model-defined fields) in the same order as that used here, and
* after user-defined variables defined throught the GUI, if used.
*
* Currently, only 1d (scalar) fields are handled.
*
* parameters for cs_parameters_add_variable():
* name <-- name of variable and associated field
* dim <-- variable dimension
*/
cs_parameters_add_variable("species_1", 1);
void cs_parameters_add_variable(const char *name, int dim)
Define a user variable.
Definition: cs_parameters.c:1137

After adding a scalar, the user can add the variance of this scalar

/* Example: add the variance of a user variable.
*
* parameters for cs_parameters_add_variable_variance():
* name <-- name of variance and associated field
* variable_name <-- name of associated variable
*/
"species_1");
void cs_parameters_add_variable_variance(const char *name, const char *variable_name)
Define a user variable which is a variance of another variable.
Definition: cs_parameters.c:1178

Add a user property defined on a mesh location (cells, interior faces, boundary faces or vertices).

/* Example: add a user property defined on boundary faces.
*
* parameters for cs_parameters_add_property():
* name <-- name of property and associated field
* dim <-- property dimension
* location_id <-- id of associated mesh location, which must be one of:
* CS_MESH_LOCATION_CELLS
* CS_MESH_LOCATION_INTERIOR_FACES
* CS_MESH_LOCATION_BOUNDARY_FACES
* CS_MESH_LOCATION_VERTICES
*/
cs_parameters_add_property("user_b_property_1",
1,
@ CS_MESH_LOCATION_BOUNDARY_FACES
Definition: cs_mesh_location.h:65
void cs_parameters_add_property(const char *name, int dim, int location_id)
Define a user property.
Definition: cs_parameters.c:1215

General options

Choose a time step option

/* Time step type */
cs_time_step_options_t * cs_get_glob_time_step_options(void)
Provide read/write access to cs_glob_time_step_options.
Definition: cs_time_step.c:447
@ CS_TIME_STEP_CONSTANT
Definition: cs_time_step.h:55
time step options descriptor
Definition: cs_time_step.h:92
cs_time_step_type_t idtvar
Definition: cs_time_step.h:98

Choose a reference time step

/* Reference time step dt_ref
The example given below is probably not adapted to your case. */
cs_real_t dt_ref = 0.005;
domain->time_step->dt_ref = dt_ref;
double cs_real_t
Floating-point value.
Definition: cs_defs.h:304

To set a duration

/* Duration
nt_max is absolute number of the last time step required
if we have already run 10 time steps and want to
run 10 more, nt_max must be set to 10 + 10 = 20 */
domain->time_step->nt_max = (int) (40./dt_ref);

For example, to change the log (run_solver.log) verbosity of all the variables:

{
int n_fields = cs_field_n_fields();
for (int f_id = 0; f_id < n_fields; f_id++) {
if (f->type & CS_FIELD_VARIABLE) {
eqp->verbosity = 2;
}
}
}
cs_field_t * cs_field_by_id(int id)
Return a pointer to a field based on its id.
Definition: cs_field.c:2314
int cs_field_n_fields(void)
Return the number of defined fields.
Definition: cs_field.c:1528
cs_equation_param_t * cs_field_get_equation_param(cs_field_t *f)
Access a field's equation parameters.
Definition: cs_field_default.c:233
#define CS_FIELD_VARIABLE
Definition: cs_field.h:63
Set of parameters to handle an unsteady convection-diffusion-reaction equation with term sources.
Definition: cs_equation_param.h:202
int verbosity
Definition: cs_equation_param.h:221
Field descriptor.
Definition: cs_field.h:125
int type
Definition: cs_field.h:130

For example, to change limiters for the convective scheme for a given scalar:

{
/* retrieve scalar field by its name */
cs_field_t *sca1 = cs_field_by_name("scalar1");
/* ischcv is the type of convective scheme:
0: second order linear upwind
1: centered
2: pure upwind gradient in SOLU
3: blending SOLU and centered
4: NVD/TVD Scheme */
/* isstpc:
0: slope test enabled
1: slope test disabled (default)
2: continuous limiter ensuring boundedness (beta limiter) enabled */
eqp->ischcv = 1;
eqp->isstpc = 0;
/* Min/Max limiter or NVD/TVD limiters
* then "limiter_choice" keyword must be set:
* CS_NVD_Gamma
* CS_NVD_SMART
* CS_NVD_CUBISTA
* CS_NVD_SUPERBEE
* CS_NVD_MUSCL
* CS_NVD_MINMOD
* CS_NVD_CLAM
* CS_NVD_STOIC
* CS_NVD_OSHER
* CS_NVD_WASEB
* --- VOF scheme ---
* CS_NVD_VOF_HRIC
* CS_NVD_VOF_CICSAM
* CS_NVD_VOF_STACS */
int key_lim_id = cs_field_key_id("limiter_choice");
/* Get the Key for the Sup and Inf for the convective scheme */
int kccmin = cs_field_key_id("min_scalar");
int kccmax = cs_field_key_id("max_scalar");
/* Set the Value for the Sup and Inf of the studied scalar
* for the Gamma beta limiter for the temperature */
}
@ CS_NVD_VOF_CICSAM
Definition: cs_convection_diffusion.h:72
int cs_field_set_key_int(cs_field_t *f, int key_id, int value)
Assign a integer value for a given key to a field.
Definition: cs_field.c:2943
int cs_field_key_id(const char *name)
Return an id associated with a given key name.
Definition: cs_field.c:2497
int cs_field_set_key_double(cs_field_t *f, int key_id, double value)
Assign a floating point value for a given key to a field.
Definition: cs_field.c:3123
cs_field_t * cs_field_by_name(const char *name)
Return a pointer to a field based on its name.
Definition: cs_field.c:2338
@ t
Definition: cs_field_pointer.h:98
#define CS_F_(e)
Macro used to return a field pointer by its enumerated value.
Definition: cs_field_pointer.h:51
int isstpc
Definition: cs_equation_param.h:463
int ischcv
Definition: cs_equation_param.h:461

One can also choose the percentage of upwind blending when using the slope test

{
/* retrieve scalar field by its name */
cs_field_t *sca1 = cs_field_by_name("scalar1");
/* blend_st (can be between 0 and 1):
0: full upwind (default)
1: scheme without upwind */
eqp->blend_st = 0.1;
}
double blend_st
Definition: cs_equation_param.h:473

If one wants to declare a scalar as buoyant (i.e. the density depends on this scalar through a given equation of state) and add it in the velocity-pressure optional inner-iterations, one can set the dedicated keyword:

{
/* retrieve scalar field by its name */
cs_field_t *sca1 = cs_field_by_name("scalar1");
int key_is_buoyant = cs_field_key_id("is_buoyant");
cs_field_set_key_int(sca1, key_is_buoyant, 1);
}

If one wants to activate drift terms on a transported scalar:

{
/* retrieve scalar field by its name */
cs_field_t *sca1 = cs_field_by_name("scalar1");
/* Key id for drift scalar */
int key_drift = cs_field_key_id("drift_scalar_model");
int drift = CS_DRIFT_SCALAR_ON;
if (false)
if (false)
if (false)
if (false)
/* Set the key word "drift_scalar_model" into the field structure */
cs_field_set_key_int(sca1, key_drift, drift);
}
@ CS_DRIFT_SCALAR_ON
Definition: cs_parameters.h:127
@ CS_DRIFT_SCALAR_THERMOPHORESIS
Definition: cs_parameters.h:129
@ CS_DRIFT_SCALAR_ADD_DRIFT_FLUX
Definition: cs_parameters.h:128
@ CS_DRIFT_SCALAR_TURBOPHORESIS
Definition: cs_parameters.h:130
@ CS_DRIFT_SCALAR_CENTRIFUGALFORCE
Definition: cs_parameters.h:132
@ CS_DRIFT_SCALAR_ZERO_BNDY_FLUX
Definition: cs_parameters.h:134

To set model and parameter options for the velocity-pressure coupling and resolution, see the cs_velocity_pressure_model_t and cs_velocity_pressure_param_t structures.

For example, to set the Rhie and Chow multiplication factor:

{
vp_param->arak = 0.;
}
cs_velocity_pressure_param_t * cs_get_glob_velocity_pressure_param(void)
Provide access to cs_glob_velocity_pressure_param.
Definition: cs_velocity_pressure.c:427
Inner velocity/pressure iteration options descriptor.
Definition: cs_velocity_pressure.h:80
double arak
Definition: cs_velocity_pressure.h:128

To set the optional number of velocity-pressure inner iterations:

{
vp_param->nterup = 3;
}
int nterup
Definition: cs_velocity_pressure.h:136

Special fields and activate some automatic post-processings

For example, to force the presence of a boundary temperature field, which may be useful for postprocessing:

{
}
cs_field_t * cs_parameters_add_boundary_temperature(void)
Define a boundary values field for temperature, if applicable.
Definition: cs_parameters.c:1506

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++) {
}
}
cs_field_t * cs_parameters_add_boundary_values(cs_field_t *f)
Define a boundary values field for a variable field.
Definition: cs_parameters.c:1401

To add handling (storage) of previous values for a field, the following following code can be added:

{
/* add previous values handling for boundary temperature */
}
void cs_field_set_n_time_vals(cs_field_t *f, int n_time_vals)
Change the number of time values managed by a field.
Definition: cs_field.c:1652
@ t_b
Definition: cs_field_pointer.h:99

Enforce existence of 'tplus' and 'tstar' fields, so that a Nusselt number may be computed using the post_boundary_nusselt subroutine. When postprocessing this quantity is activated, those fields are present, but if we need to compute them in the cs_user_extra_operations user subroutine without postprocessing them, forcing the definition of these fields to save the values computed for the boundary layer is necessary.

{
f = cs_field_by_name_try("yplus");
if (f != NULL)
1,
f = cs_field_by_name_try("tplus");
if (f != NULL)
1,
f = cs_field_by_name_try("tstar");
if (f != NULL)
1,
}
cs_field_t * cs_field_by_name_try(const char *name)
Return a pointer to a field based on its name if present.
Definition: cs_field.c:2364

You can activate the post-processing of the Q-criterion on the whole domain mesh with:

@ CS_POST_UTIL_Q_CRITERION
Definition: cs_post_util.h:48
int cs_glob_post_util_flag[]

Save contribution of slope test for variables in special fields. These fields are automatically created, with postprocessing output enabled, if the matching variable is convected, does not use a pure upwind scheme, and has a slope test (the slope_test_upwind_id key value for a given variable's field is automatically set to the matching postprocessing field's id, or -1 if not applicable).

{
int n_fields = cs_field_n_fields();
for (int f_id = 0; f_id < n_fields; f_id++) {
cs_field_set_key_int(f, cs_field_key_id("slope_test_upwind_id"), 0);
}
}

You can activate the post-processing of clipping on turbulent quantities on the whole domain mesh with:

@ eps
Definition: cs_field_pointer.h:71
@ rij
Definition: cs_field_pointer.h:79

Input-output related examples (usipes)

Basic options

Frequency of log output.

int cs_glob_log_frequency

Change a property's label (here for density, first checking if it is variable). A field's name cannot be changed, but its label, used for logging and postprocessing output, may be redefined.

if (CS_F_(cp) != NULL)
int cs_field_set_key_str(cs_field_t *f, int key_id, const char *str)
Assign a character string for a given key to a field.
Definition: cs_field.c:3241
@ cp
Definition: cs_field_pointer.h:106

Postprocessing output

Activate or deactivate probes output.

cs_field_key_id("post_vis"),
int cs_field_set_key_int_bits(cs_field_t *f, int key_id, int mask)
Set integer bits matching a mask to 1 for a given key for a field.
Definition: cs_field.c:3062
@ vel
Definition: cs_field_pointer.h:68
#define CS_POST_MONITOR
Definition: cs_post.h:59

Probes for Radiative Transfer (Luminance and radiative density flux vector).

cs_field_t *f = cs_field_by_name_try("luminance");
if (f != NULL)
cs_field_key_id("post_vis"),
f = cs_field_by_name_try("radiative_flux");
if (f != NULL)
cs_field_key_id("post_vis"),

Base model for CDO/HHO schemes

Activation of CDO/HHO schemes

Two modes are available CS_DOMAIN_CDO_MODE_ONLY or CS_DOMAIN_CDO_MODE_WITH_FV

CDO/HHO schemes can be activated within this function as follows:

{
}
void cs_domain_set_cdo_mode(cs_domain_t *domain, int mode)
Set the global variable storing the mode of activation to apply to CDO/HHO schemes.
Definition: cs_domain.c:267
cs_domain_t * cs_glob_domain
Definition: cs_domain.c:85
#define CS_DOMAIN_CDO_MODE_ONLY
Definition: cs_domain.h:59
Structure storing the main features of the computational domain and pointers to the main geometrical ...
Definition: cs_domain.h:88

Domain boundary with CDO/HHO schemes

Several types of domain boundaries can be defined. There are gathered in cs_domain_boundary_type_t The definition of the domain boundaries for CDO/HHO schemes can be specified as follows:

{
cs_boundary_t *bdy = domain->boundaries;
/* Choose a boundary by default */
/* Add a new boundary
* >> cs_domain_add_boundary(boundary, type_of_boundary, zone_name);
*
* zone_name is either a predefined one or user-defined one
* type_of_boundary is one of the following keywords:
* CS_BOUNDARY_WALL,
* CS_BOUNDARY_INLET,
* CS_BOUNDARY_OUTLET,
* CS_BOUNDARY_SYMMETRY
*/
}
void cs_boundary_set_default(cs_boundary_t *boundaries, cs_boundary_type_t type)
Set the default boundary related to the given cs_boundary_t structure.
Definition: cs_boundary.c:429
void cs_boundary_add(cs_boundary_t *bdy, cs_boundary_type_t type, const char *zone_name)
Add a new boundary type for a given boundary zone.
Definition: cs_boundary.c:505
@ CS_BOUNDARY_OUTLET
Definition: cs_boundary.h:86
@ CS_BOUNDARY_INLET
Definition: cs_boundary.h:83
@ CS_BOUNDARY_WALL
Definition: cs_boundary.h:80
Structure storing information related to the "physical" boundaries associated with the computational ...
Definition: cs_boundary.h:155
cs_boundary_t * boundaries
Definition: cs_domain.h:102

Output with CDO/HHO schemes

The management of the level and frequency of details written by the code can be specified for CDO/HHO schemes as follows:

{
-1, // restart frequency
10, // output log frequency
2); // verbosity (-1: no, 0, ...)
}
void cs_domain_set_output_param(cs_domain_t *domain, int nt_interval, int nt_list, int verbosity)
Set auxiliary parameters related to the way output is done.
Definition: cs_domain_setup.c:315

Time step with CDO/HHO schemes

The management of the time step with CDO/HHO schemes can be specified as follows:

{
/*
If there is an inconsistency between the max. number of iteration in
time and the final physical time, the first condition encountered stops
the calculation.
*/
100, // nt_max or -1 (automatic)
-1.); // t_max or < 0. (automatic)
/* Define the value of the time step
>> cs_domain_def_time_step_by_value(domain, dt_val);
>> cs_domain_def_time_step_by_func(domain, dt_func);
The second way to define the time step enable complex definitions.
*/
}
void cs_domain_def_time_step_by_value(cs_domain_t *domain, double dt)
Define the value of the time step.
Definition: cs_domain_setup.c:442
void cs_domain_set_time_param(cs_domain_t *domain, int nt_max, double t_max)
Set parameters for unsteady computations: the max number of time steps or the final physical time of ...
Definition: cs_domain_setup.c:342

Wall distance with CDO/HHO schemes

The computation of the wall distance with CDO schemes is performed as follows:

{
/* Activate predefined module as the computation of the wall distance */
}
void cs_walldistance_activate(void)
Activate the future computation of the wall distance.
Definition: cs_walldistance.c:416

Add a user-defined equation with CDO/HHO schemes

The add of a user-defined equation solved by CDO/HHO schemes is specified as follows:

{
/* Add a new user equation.
* The default boundary condition has to be chosen among:
* CS_PARAM_BC_HMG_DIRICHLET
* CS_PARAM_BC_HMG_NEUMANN
*/
cs_equation_add_user("AdvDiff.Upw", // equation name
"Pot.Upw", // associated variable field name
1, // dimension of the unknown
CS_PARAM_BC_HMG_DIRICHLET); // default boundary
cs_equation_add_user("AdvDiff.SG", // equation name
"Pot.SG", // associated variable field name
1, // dimension of the unknown
CS_PARAM_BC_HMG_DIRICHLET); // default boundary
}
cs_equation_t * cs_equation_add_user(const char *eqname, const char *varname, int dim, cs_param_bc_type_t default_bc)
Add a new user equation structure and set a first set of parameters.
Definition: cs_equation.c:1279
@ CS_PARAM_BC_HMG_DIRICHLET
Definition: cs_param_types.h:433

Add user-defined properties with CDO/HHO schemes

The add of a new user-defined property with CDO/HHO schemes is specified as follows:

{
cs_property_add("conductivity", /* property name */
CS_PROPERTY_ANISO); /* type of material property */
cs_property_add("rho.cp", /* property name */
CS_PROPERTY_ISO); /* type of material property */
}
cs_property_t * cs_property_add(const char *name, cs_property_type_t type)
Create and initialize a new property structure.
Definition: cs_property.c:733
#define CS_PROPERTY_ISO
Definition: cs_property.h:73
#define CS_PROPERTY_ANISO
Definition: cs_property.h:84

If you want to compute the Fourier number related to a given property in an unsteady simulation

{
/* Retrieve a property named "conductivity" */
cs_property_t *pty = cs_property_by_name("conductivity");
/* Activate the computation of the Fourier number for this property */
}
cs_property_t * cs_property_by_name(const char *name)
Find the related property definition from its name.
Definition: cs_property.c:843
void cs_property_set_option(cs_property_t *pty, cs_property_key_t key)
Set optional parameters related to a cs_property_t structure.
Definition: cs_property.c:890
@ CS_PTYKEY_POST_FOURIER
Definition: cs_property.h:113
Structure associated to the definition of a property relying on the cs_xdef_t structure.

Add user-defined advection field with CDO/HHO schemes

The definition of an advection field allows one to handle flows with a frozen velocity field or the transport of scalar quantities without solving the Navier-Stokes system. The add of a new user-defined advection field with CDO/HHO schemes is specified as follows:

{
/* Add a user-defined advection field named "adv_field" */
}
cs_adv_field_t * cs_advection_field_add_user(const char *name)
Add and initialize a new user-defined advection field structure.
Definition: cs_advection_field.c:451
Definition: cs_advection_field.h:149

If you need to activate options related to advection fields, you can also specify

Settings related to the groundwater flow module with CDO/HHO schemes

Activation of the groundwater flow module. The second argument is either 0 if no option is needed or a list of flags among CS_GWF_GRAVITATION, CS_GWF_RICHARDS_UNSTEADY, CS_GWF_SOIL_PROPERTY_UNSTEADY, CS_GWF_SOIL_ALL_SATURATED

{
cs_flag_t option_flag = 0;
}
unsigned short int cs_flag_t
Definition: cs_defs.h:306
cs_gwf_t * cs_gwf_activate(cs_property_type_t pty_type, cs_flag_t flag)
Initialize the module dedicated to groundwater flows.
Definition: cs_gwf.c:606

Add soils (settings of the soil is performed in cs_user_gwf_setup). Soils have to be added before adding tracers.

{
}
cs_gwf_soil_t * cs_gwf_soil_add(const char *z_name, cs_gwf_soil_hydraulic_model_t model)
Create and add a new cs_gwf_soil_t structure. A first initialization of all members by default is per...
Definition: cs_gwf_soil.c:381
@ CS_GWF_SOIL_SATURATED
Definition: cs_gwf_soil.h:105

Add tracer equations which correspond to a transport equation using the darcean flux as the advection field. This call implies the creation of a new equation related to a new variable field name given as arguments. Advanced tracer equations can be defined using cs_gwf_add_user_tracer.

{
cs_gwf_tracer_model_t model = 0; /* default model */
cs_gwf_add_tracer(model, "Tracer_01","C1");
}
cs_gwf_tracer_t * cs_gwf_add_tracer(cs_gwf_tracer_model_t model, const char *eq_name, const char *var_name)
Add a new equation related to the groundwater flow module This equation is a particular type of unste...
Definition: cs_gwf.c:814
cs_flag_t cs_gwf_tracer_model_t
Definition: cs_gwf_tracer.h:47

General options

CDO/HHO schemes

Modifiy the numerical parameters related to a given equation.

{
/* The modification of the space discretization should be apply first */
/* Modify other parameters than the space discretization */
/* Linear algebra settings */
#if defined(HAVE_PETSC)
#else
#endif
}
cs_equation_param_t * cs_equation_param_by_name(const char *eqname)
Return the cs_equation_param_t structure associated to a cs_equation_t structure thanks to the equati...
Definition: cs_equation.c:447
void cs_equation_set_param(cs_equation_param_t *eqp, cs_equation_key_t key, const char *keyval)
Set a parameter attached to a keyname in a cs_equation_param_t structure.
Definition: cs_equation_param.c:1222
@ CS_EQKEY_SOLVER_FAMILY
Definition: cs_equation_param.h:1100
@ CS_EQKEY_VERBOSITY
Definition: cs_equation_param.h:1104
@ CS_EQKEY_ITSOL_MAX_ITER
Definition: cs_equation_param.h:1095
@ CS_EQKEY_ITSOL
Definition: cs_equation_param.h:1093
@ CS_EQKEY_HODGE_DIFF_COEF
Definition: cs_equation_param.h:1090
@ CS_EQKEY_HODGE_DIFF_ALGO
Definition: cs_equation_param.h:1089
@ CS_EQKEY_ITSOL_RESNORM_TYPE
Definition: cs_equation_param.h:1096
@ CS_EQKEY_ITSOL_EPS
Definition: cs_equation_param.h:1094
@ CS_EQKEY_PRECOND
Definition: cs_equation_param.h:1098
@ CS_EQKEY_ADV_SCHEME
Definition: cs_equation_param.h:1078
@ CS_EQKEY_SPACE_SCHEME
Definition: cs_equation_param.h:1101

Finalize the set-up for CDO/HHO schemes

Set up properties with CDO/HHO schemes

When a property has been added, the second step is to define this property. According to the type of property (isotropic, orthotropic or anisotropic) definitions differ. Here are two examples:

{
cs_property_t *conductivity = cs_property_by_name("conductivity");
cs_real_33_t tensor = {{1.0, 0.5, 0.0}, {0.5, 1.0, 0.5}, {0.0, 0.5, 1.0}};
cs_property_def_aniso_by_value(conductivity, // property structure
"cells", // name of the volume zone
tensor); // values of the property
cs_property_t *rhocp = cs_property_by_name("rho.cp");
cs_real_t iso_val = 2.0;
cs_property_def_iso_by_value(rhocp, // property structure
"cells", // name of the volume zone
iso_val); // value of the property
}
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:324
cs_xdef_t * cs_property_def_aniso_by_value(cs_property_t *pty, const char *zname, cs_real_t tens[3][3])
Define an anisotropic cs_property_t structure by value for entities related to a volume zone.
Definition: cs_property.c:1317
cs_xdef_t * cs_property_def_iso_by_value(cs_property_t *pty, const char *zname, double val)
Define an isotropic cs_property_t structure by value for entities related to a volume zone.
Definition: cs_property.c:1213

Set up user-defined advection field with CDO/HHO schemes

When an advection field has been added, the second step is to define this advection field. Here are is an example of definition using an anlytic function and the activation of optional features:

{
cs_advection_field_def_by_analytic(adv, _define_adv_field, NULL);
/* Activate the post-processing of the related Courant number */
}
void cs_advection_field_def_by_analytic(cs_adv_field_t *adv, cs_analytic_func_t *func, void *input)
Define a cs_adv_field_t structure thanks to an analytic function.
Definition: cs_advection_field.c:766
void cs_advection_field_set_postprocess(cs_adv_field_t *adv, cs_flag_t post_flag)
Set optional post-processings.
Definition: cs_advection_field.c:717
cs_adv_field_t * cs_advection_field_by_name(const char *name)
Search in the array of advection field structures which one has the name given in argument.
Definition: cs_advection_field.c:400
#define CS_ADVECTION_FIELD_POST_COURANT
Definition: cs_advection_field.h:60

Set up the boundary conditions for an equation

{
"boundary_faces", // zone name
_define_bcs, // pointer to the function
NULL); // input structure
}
cs_xdef_t * cs_equation_add_bc_by_analytic(cs_equation_param_t *eqp, const cs_param_bc_type_t bc_type, const char *z_name, cs_analytic_func_t *analytic, void *input)
Define and initialize a new structure to set a boundary condition related to the given equation param...
Definition: cs_equation_param.c:1909
@ CS_PARAM_BC_DIRICHLET
Definition: cs_param_types.h:434

Add terms to an equation

Add terms like diffusion term, advection term, unsteady term, reaction terms or source terms.

{
cs_property_t *rhocp = cs_property_by_name("rho.cp");
cs_property_t *conductivity = cs_property_by_name("conductivity");
/* Activate the unsteady term */
cs_equation_add_time(eqp, rhocp);
/* Activate the diffusion term */
cs_equation_add_diffusion(eqp, conductivity);
/* Activate advection effect */
/* Simple definition with cs_equation_add_source_term_by_val
where the value of the source term is given by m^3
*/
cs_real_t st_val = -0.1;
"cells",
&st_val);
}
void cs_equation_add_advection(cs_equation_param_t *eqp, cs_adv_field_t *adv_field)
Associate a new term related to the advection operator for the equation associated to the given cs_eq...
Definition: cs_equation_param.c:2296
cs_xdef_t * cs_equation_add_source_term_by_val(cs_equation_param_t *eqp, const char *z_name, cs_real_t *val)
Add a new source term by initializing a cs_xdef_t structure. Case of a definition by a constant value...
Definition: cs_equation_param.c:2384
void cs_equation_add_time(cs_equation_param_t *eqp, cs_property_t *property)
Associate a new term related to the time derivative operator for the equation associated to the given...
Definition: cs_equation_param.c:2271
void cs_equation_add_diffusion(cs_equation_param_t *eqp, cs_property_t *property)
Associate a new term related to the Laplacian operator for the equation associated to the given cs_eq...
Definition: cs_equation_param.c:2188
Structure storing medata for defining a quantity in a very flexible way.
Definition: cs_xdef.h:154

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 preconditioned 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_t * cs_multigrid_define(int f_id, const char *name, cs_multigrid_type_t mg_type)
Define and associate a multigrid sparse linear system solver for a given field or equation name.
Definition: cs_multigrid.c:3748
@ CS_MULTIGRID_V_CYCLE
Definition: cs_multigrid.h:59

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 */
}
cs_sles_it_t * cs_sles_it_define(int f_id, const char *name, cs_sles_it_type_t solver_type, int poly_degree, int n_max_iter)
Define and associate an iterative sparse linear system solver for a given field or equation name.
Definition: cs_sles_it.c:3586
@ CS_SLES_BICGSTAB2
Definition: cs_sles_it.h:66
int id
Definition: cs_field.h:129

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);
}
@ p
Definition: cs_field_pointer.h:67
cs_sles_t * cs_sles_find_or_add(int f_id, const char *name)
Return pointer to linear system object, based on matching field id or system name.
Definition: cs_sles.c:1174
void cs_sles_set_verbosity(cs_sles_t *sles, int verbosity)
Set the verbosity for a given linear equation solver.
Definition: cs_sles.c:1352
struct _cs_sles_t cs_sles_t
Definition: cs_sles.h:68

Example: error visualization

The following example shows how to activate local error visualization output (here for velocity and pressure).

{
cs_sles_t *sles_p = cs_sles_find_or_add(CS_F_(p)->id, NULL);
cs_sles_t *sles_u = cs_sles_find_or_add(CS_F_(vel)->id, NULL);
}
#define CS_POST_WRITER_DEFAULT
Definition: cs_post.h:65
void cs_sles_set_post_output(cs_sles_t *sles, int writer_id)
Activate postprocessing output for a given linear equation solver.
Definition: cs_sles.c:1390

Example: advanced multigrid settings

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

{
NULL,
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) */
}
void cs_multigrid_set_coarsening_options(cs_multigrid_t *mg, int aggregation_limit, int coarsening_type, int n_max_levels, cs_gnum_t min_g_rows, double p0p1_relax, int postprocess)
Set multigrid coarsening parameters.
Definition: cs_multigrid.c:4013
void cs_multigrid_set_solver_options(cs_multigrid_t *mg, cs_sles_it_type_t descent_smoother_type, cs_sles_it_type_t ascent_smoother_type, cs_sles_it_type_t coarse_solver_type, int n_max_cycles, int n_max_iter_descent, int n_max_iter_ascent, int n_max_iter_coarse, int poly_degree_descent, int poly_degree_ascent, int poly_degree_coarse, double precision_mult_descent, double precision_mult_ascent, double precision_mult_coarse)
Set multigrid parameters for associated iterative solvers.
Definition: cs_multigrid.c:4066
struct _cs_multigrid_t cs_multigrid_t
Definition: cs_multigrid.h:68
@ CS_SLES_JACOBI
Definition: cs_sles_it.h:63
@ CS_SLES_PCG
Definition: cs_sles_it.h:59

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_SYM_GAUSS_SEIDEL) */
CS_SLES_P_GAUSS_SEIDEL, /* ascent smoother (CS_SLES_P_SYM_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) */
}
cs_sles_pc_t * cs_multigrid_pc_create(cs_multigrid_type_t mg_type)
Create a multigrid preconditioner.
Definition: cs_multigrid.c:4504
void cs_sles_it_transfer_pc(cs_sles_it_t *context, cs_sles_pc_t **pc)
Assign a preconditioner to an iterative sparse linear equation solver, transfering its ownership to t...
Definition: cs_sles_it.c:4350
cs_sles_pc_t * cs_sles_it_get_pc(cs_sles_it_t *context)
Return a preconditioner context for an iterative sparse linear equation solver.
Definition: cs_sles_it.c:4321
struct _cs_sles_it_t cs_sles_it_t
Definition: cs_sles_it.h:86
@ CS_SLES_P_GAUSS_SEIDEL
Definition: cs_sles_it.h:69
void * cs_sles_pc_get_context(cs_sles_pc_t *pc)
Return pointer to preconditioner context structure pointer.
Definition: cs_sles_pc.c:850
const char * cs_sles_pc_get_type(cs_sles_pc_t *pc)
Return type name of preconditioner context.
Definition: cs_sles_pc.c:804
struct _cs_sles_pc_t cs_sles_pc_t
Definition: cs_sles_pc.h:66

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.

{
NULL,
4, /* # of ranks merged at a time */
300, /* mean # of cells under which we merge */
500); /* global # of cells under which we merge */
}
void cs_multigrid_set_merge_options(cs_multigrid_t *mg, int rank_stride, int rows_mean_threshold, cs_gnum_t rows_glob_threshold)
Set global multigrid parameters for parallel grid merging behavior.
Definition: cs_multigrid.c:4833

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);
}
}
void cs_multigrid_set_plot_options(cs_multigrid_t *mg, const char *base_name, bool use_iteration)
Set plotting options for multigrid.
Definition: cs_multigrid.c:4750
const char * cs_sles_get_type(cs_sles_t *sles)
Return type name of solver context.
Definition: cs_sles.c:1450
void * cs_sles_get_context(cs_sles_t *sles)
Return pointer to solver context structure pointer.
Definition: cs_sles.c:1470
void cs_sles_it_set_plot_options(cs_sles_it_t *context, const char *base_name, bool use_iteration)
Set plotting options for an iterative sparse linear equation solver.
Definition: cs_sles_it.c:4688
const char * name
Definition: cs_field.h:127

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 */
#if PETSC_VERSION_GE(3,7,0)
PetscOptionsSetValue(NULL, "-ksp_type", "cg");
PetscOptionsSetValue(NULL, "-pc_type", "jacobi");
#else
PetscOptionsSetValue("-ksp_type", "cg");
PetscOptionsSetValue("-pc_type", "jacobi");
#endif
}
MPI_Comm cs_glob_mpi_comm
Definition: cs_defs.c:179

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);
}
cs_sles_petsc_t * cs_sles_petsc_define(int f_id, const char *name, MatType matrix_type, cs_sles_petsc_setup_hook_t *setup_hook, void *context)
Define and associate a PETSc linear system solver for a given field or equation name.
Definition: cs_sles_petsc.c:533

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(void *context,
KSP ksp)
{
CS_UNUSED(context);
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 */
}
#define CS_UNUSED(x)
Definition: cs_defs.h:478

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 */
#if PETSC_VERSION_GE(3,7,0)
PetscOptionsSetValue(NULL, "-ksp_type", "cg");
PetscOptionsSetValue(NULL, "-pc_type", "gamg");
PetscOptionsSetValue(NULL, "-pc_gamg_agg_nsmooths", "1");
PetscOptionsSetValue(NULL, "-mg_levels_ksp_type", "richardson");
PetscOptionsSetValue(NULL, "-mg_levels_pc_type", "sor");
PetscOptionsSetValue(NULL, "-mg_levels_ksp_max_it", "1");
PetscOptionsSetValue(NULL, "-pc_gamg_threshold", "0.02");
PetscOptionsSetValue(NULL, "-pc_gamg_reuse_interpolation", "TRUE");
PetscOptionsSetValue(NULL, "-pc_gamg_square_graph", "4");
#else
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");
#endif
}

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(void *context,
KSP ksp)
{
CS_UNUSED(context);
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 */
#if PETSC_VERSION_GE(3,7,0)
PetscOptionsSetValue(NULL, "-ksp_type", "cg");
PetscOptionsSetValue(NULL, "-pc_type", "hypre");
PetscOptionsSetValue(NULL, "-pc_hypre_type","boomeramg");
PetscOptionsSetValue(NULL, "-pc_hypre_boomeramg_coarsen_type","HMIS");
PetscOptionsSetValue(NULL, "-pc_hypre_boomeramg_interp_type","ext+i-cc");
PetscOptionsSetValue(NULL, "-pc_hypre_boomeramg_agg_nl","2");
PetscOptionsSetValue(NULL, "-pc_hypre_boomeramg_P_max","4");
PetscOptionsSetValue(NULL, "-pc_hypre_boomeramg_strong_threshold","0.5");
PetscOptionsSetValue(NULL, "-pc_hypre_boomeramg_no_CF","");
#else
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","");
#endif
}

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(void *context,
KSP ksp)
{
CS_UNUSED(context);
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(void *context,
KSP ksp)
{
CS_UNUSED(context);
const char *p = getenv("CS_USER_PETSC_MAT_VIEW");
if (p != NULL) {
/* Get system and preconditioner matrices */
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);
}
}
}
double precision, save a
Definition: cs_fuel_incl.f90:146

Using AmgX

The AmgX library provides advanced solvers targeting NVIDIA GPUs.

As described in its documentation, AmgX solvers can be configured either using a configuration string (containing key/value pairs), or .json formatted files.

The following example shows how to setup a solver to use the AmgX library, if the code was built with AmgX support. In this example, a configuration file (which must be present in the case's DATA directory) is used.

{
cs_sles_amgx_t *amgx_p = cs_sles_amgx_define(CS_F_(p)->id, NULL);
cs_sles_amgx_set_config_file(amgx_p, "PCG_CLASSICAL_V_JACOBI.json");
}
cs_sles_amgx_t * cs_sles_amgx_define(int f_id, const char *name)
Define and associate an AmgX linear system solver for a given field or equation name.
Definition: cs_sles_amgx.c:740
void cs_sles_amgx_set_config_file(void *context, const char *path)
Set the solver configuration file for an AmgX solver.
Definition: cs_sles_amgx.c:1024
struct _cs_sles_amgx_t cs_sles_amgx_t
Definition: cs_sles_amgx.h:60

To set options using a string, the cs_sles_amgx_set_config function may be used intead of cs_sles_amgx_set_config_file. If neither is called, a defaut configuration is used.

The cs_sles_amgx_set_pin_memory and cs_sles_amgx_set_use_device functions may also be called to modify default behavior relative to using pinned memory and running on the device or host.

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.

{
/* Moment <U> calculated starting from time step 1000. */
/* resulting field is a vector */
int moment_f_id[] = {CS_F_(vel)->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);
}
int cs_time_moment_define_by_field_ids(const char *name, int n_fields, const int field_id[], const int component_id[], cs_time_moment_type_t type, int nt_start, double t_start, cs_time_moment_restart_t restart_mode, const char *restart_name)
Define a moment of a product of existing fields components.
Definition: cs_time_moment.c:1522
@ CS_TIME_MOMENT_RESTART_AUTO
Definition: cs_time_moment.h:68
@ CS_TIME_MOMENT_MEAN
Definition: cs_time_moment.h:58

In the next example, we define the variance of the vector velocity. All components are used again (component -1 means all components), so the moment is a tensor.

{
/* Second order moment <UU>-<U><U> calculated starting from time step 1000. */
/* resulting field is a tensor */
int moment_f_id[] = {CS_F_(vel)->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);
}
@ CS_TIME_MOMENT_VARIANCE
Definition: cs_time_moment.h:59

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 <rho U> (vector) calculated starting from time step 1000. */
/* resulting field is a vector */
int moment_f_id[] = {CS_F_(rho)->id, CS_F_(vel)->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);
}
@ rho
Definition: cs_field_pointer.h:103

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_(vel)->id, CS_F_(vel)->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);
@ CS_TIME_MOMENT_RESTART_RESET
Definition: cs_time_moment.h:67

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];
}
}
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:298
const cs_lnum_t * cs_mesh_location_get_n_elts(int id)
Get a mesh location's number of elements.
Definition: cs_mesh_location.c:823
@ CS_MESH_LOCATION_CELLS
Definition: cs_mesh_location.h:63
cs_real_t * val
Definition: cs_field.h:146

In cs_user_time_moments, we can now 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);
}
}
int cs_time_moment_define_by_func(const char *name, int location_id, int dim, cs_time_moment_data_t *data_func, const void *data_input, cs_time_moment_data_t *w_data_func, void *w_data_input, cs_time_moment_type_t type, int nt_start, double t_start, cs_time_moment_restart_t restart_mode, const char *restart_name)
Define a moment whose data values will be computed using a specified function.
Definition: cs_time_moment.c:1579
cs_time_moment_type_t
Definition: cs_time_moment.h:56

Example 5

This next example illustrates the use of another user-defined function to evaluate expressions. Here, we compute the moment of the thermal flux at the boundary. We also that we compute both means and variances here.

static void
_boundary_thermal_flux(const void *input,
cs_real_t *vals)
{
const int location_id = CS_MESH_LOCATION_BOUNDARY_FACES;
const cs_lnum_t n_elts = cs_mesh_location_get_n_elts(location_id)[0];
cs_post_boundary_flux(f->name, n_elts, NULL, vals);
}
void cs_post_boundary_flux(const char *scalar_name, cs_lnum_t n_loc_b_faces, const cs_lnum_t b_face_ids[], cs_real_t b_face_flux[])
Compute scalar flux on a specific boundary region.
Definition: cs_post_util.c:1269
cs_field_t * cs_thermal_model_field(void)
Definition: cs_thermal_model.c:205

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

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

Example 6

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_(vel)->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;
}
}
#define restrict
Definition: cs_defs.h:124
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:317
cs_mesh_quantities_t * cs_glob_mesh_quantities
cs_rotation_t * cs_glob_rotation
cs_real_t * cell_cen
Definition: cs_mesh_quantities.h:91
Subdomain rotation description.
Definition: cs_rotation.h:46

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);
}
}

To activate means for all variables:

for (int f_id = 0; f_id < cs_field_n_fields(); f_id++) {
if (f->type & CS_FIELD_VARIABLE) {
int moment_f_id[] = {f_id};
int moment_c_id[] = {-1};
int n_fields = 1;
const char extension[] = "_mean";
char *mean_name;
BFT_MALLOC(mean_name, strlen(f->name) + 1 + 5, char);
strcpy(mean_name, f->name); /* copy field name into the new var */
strcat(mean_name, extension); /* add the extension */
n_fields,
moment_f_id,
moment_c_id,
10, /* nt_start */
-1, /* t_start */
NULL);
}
}
#define BFT_MALLOC(_ptr, _ni, _type)
Allocate memory for _ni elements of type _type.
Definition: bft_mem.h:62

Fan modelling options

Code_Saturne allows modelling 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 (2D or 3D) */
0, /* 0: Fan mode or 1: wind turbine mode */
inlet_axis_coords,
outlet_axis_coords,
fan_radius,
blades_radius,
hub_radius,
pressure_curve_coeffs,
axial_torque);
}
void cs_fan_define(int fan_dim, int mode, const cs_real_t inlet_axis_coords[3], const cs_real_t outlet_axis_coords[3], cs_real_t fan_radius, cs_real_t blades_radius, cs_real_t hub_radius, const cs_real_t curve_coeffs[3], cs_real_t axial_torque)
Fan definition (added to the ones previously defined)
Definition: cs_fan.c:231