Input of calculation parameters (C functions in cs_user_parameters.cpp)

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.

To modify the default value of parameters which do not appear in the examples provided, code should be placed as follows:

• cs_user_model for numerical and physical options
• cs_user_parameters for numerical and physical options
• cs_user_finalize_setup for input-output related options

Base model related options (cs_user_model)

Definition of user variables or properties as well as choices of physical models should be done in the cs_user_model function, if not already done through the GUI. Gravity and reference frame rotation (for Coriolis) can also be (re)-defined here.

Activate user models (such as atmospheric flow module)

 
  /* Activate Atmospheric flow model
   *  CS_ATMO_OFF:                module not activated
   *  CS_ATMO_CONSTANT_DENSITY:   standard modelling
   *  CS_ATMO_DRY:                dry atmosphere
   *  CS_ATMO_HUMID:              humid atmosphere (experimental)
   * */
  cs_glob_physical_model_flag[CS_ATMOSPHERIC] = CS_ATMO_DRY;
 
  /* Activate compressible model
   *  -1: not active, 0: activated, 1: barotropic version,
   *   2: homogeneous two phase model, 3: by pressure increment */
  cs_glob_physical_model_flag[CS_COMPRESSIBLE] = 3;
 
  /* Activate gas combustion model */
  cs_combustion_gas_model_t *gcm
    = cs_combustion_gas_set_model(CS_COMBUSTION_3PT_PERMEATIC);
 
  /* Soot model for gas combustion
   *  if = -1   module not activated
   *  if =  0   constant soot yield
   *  if =  1   2 equations model of Moss et al. */
  gcm->isoot = 0;
  gcm->xsoot = 0.1;    /* (only if isoot = 0 and soot yield is not
                          defined in the thermochemistry data file) */
  gcm->rosoot = 2000.; /* kg/m3 */
 
  /* Enthalpy-Temperature conversion law */
  gcm->use_janaf = true;
 
  if (gcm->use_janaf)
    cs_combustion_gas_set_thermochemical_data_file("dp_C3P");
  else
    cs_combustion_gas_set_thermochemical_data_file("dp_C3PSJ");
 
  /* Activate pulverized coal combustion model */
  cs_coal_model_set_model(CS_COMBUSTION_COAL_STANDARD);
 
  /* Activate the drift (for coal combustion)
   * 0 (no activation),
   * 1 (transported particle velocity)
   * 2 (limit drop particle velocity) */
  cs_glob_coal_model->idrift = 1;
 
  /* Kinetic model for CO <=> CO2 (for coal combustion)
   * 0:  unused (maximal conversion in turbulent model)
   * 1:  transport of CO2 mass fraction
   * 2:  transport of CO mass fraction */
  cs_glob_coal_model->ieqco2 = 0;
 
  /* Cooling towers model
   * -1: not active
   *  1: Poppe's model
   *  2: Merkel's model */
  cs_glob_physical_model_flag[CS_COOLING_TOWERS] = 1;
 
  /* Electric arcs model
   * -1: not active
   *  2: electric potential and vector potential */
  cs_glob_physical_model_flag[CS_ELECTRIC_ARCS] = 2;
 
  /* Joule effcets model
   * -1: not active
   *  1: real potential
   *  2: complex potential
   *  3: real potential + Transfo
   *  4: complex potential + Transfo */
  cs_glob_physical_model_flag[CS_JOULE_EFFECT] = 2;
 
  /* Radiative transfer module
   * CS_RAD_TRANSFER_NONE: not active
   * CS_RAD_TRANSFER_DOM: discrete ordinates model
   * CS_RAD_TRANSFER_P1: P1 model */
  cs_glob_rad_transfer_params->type = CS_RAD_TRANSFER_DOM;
 
  /*  Activate gas mix model
   *  CS_GAS_MIX_OFF                      gas mix model off
   *  CS_GAS_MIX_AIR_HELIUM               air/helium, helium deduced
   *  CS_GAS_MIX_AIR_HYDROGEN             air/hydrogen, hydrogen deduced
   *  CS_GAS_MIX_AIR_STEAM                air/steam, steam deduced
   *  CS_GAS_MIX_AIR_HELIUM_STEAM         helium/steam, steam deduced
   *  CS_GAS_MIX_AIR_HYDROGEN_STEAM       hydrogen/steam, steam deduced
   *  CS_GAS_MIX_HELIUM_AIR               helium/air, O2 from air deduced
   *  CS_GAS_MIX_USER                     user defined
   * */
  cs_glob_physical_model_flag[CS_GAS_MIX] = CS_GAS_MIX_AIR_HELIUM;
 

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 *turb_model = cs_get_glob_turb_model();
  turb_model->model = CS_TURB_K_EPSILON_LIN_PROD;
 

Coupled solver for Rij components (when Rij is used), and other RANS model specific settings:

  /* Coupled solver for Rij components (CS_TURB_RIJ_*)
   *   0: switch off
   *   1: switch on (default)
   */
 
  cs_turb_rans_model_t *rans_model = cs_get_glob_turb_rans_model();
  rans_model->irijco = 0;
 
  /* 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, default)
     - 1: tensorial diffusivity (Daly and Harlow model)
     */
  rans_model->idirsm = 1;
 
  /* Rotation/curvature correction for eddy-viscosity turbulence models
     (k-epsilon, k-omega) */
 
  rans_model->irccor = 1;
 

To set a thermal model (CS_THERMAL_MODEL_NONE: none, CS_THERMAL_MODEL_TEMPERATURE: temperature, CS_THERMAL_MODEL_ENTHALPY: enthalpy, CS_THERMAL_MODEL_TOTAL_ENERGY: total energy)

 
  /* Example: Choose a thermal model
   *
   * CS_THERMAL_MODEL_NONE        : none
   * CS_THERMAL_MODEL_TEMPERATURE : temperature
   * CS_THERMAL_MODEL_ENTHALPY    : enthalpy
   * CS_THERMAL_MODEL_TOTAL_ENERGY: total energy (only for compressible module)
   *
   *  For temperature, the temperature scale may be set later using the
   *  temperature_scale member.
   *
   *  Warning: When using specific physics, this value is
   *           set automatically by the physics model.
   *
   */
 
  cs_thermal_model_t *thermal_model = cs_get_glob_thermal_model();
  thermal_model->thermal_variable = CS_THERMAL_MODEL_TEMPERATURE;
 

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

  cs_vof_parameters_t *vof_param = cs_get_glob_vof_parameters();
  vof_param->vof_model = CS_VOF_ENABLED | CS_VOF_MERKLE_MASS_TRANSFER;

ALE activation

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
   */
 
  cs_glob_ale = CS_ALE_LEGACY;
 
  /* Number of iterations for fluid initialization.
     cs_glob_ale_n_ini_f is not an absolute iteration number, so for a
     restarted calculation, it corresponds to the number of iterations
     for fuid initialization relative to the first restart iteration.
     In general it should be set to 0 in that case. */
 
  cs_glob_ale_n_ini_f = 75;
 
  /* Maximum number of iterations in case of implicit Fluid Structure Coupling
     with structural calculations (internal and/or external
     (i.e. using code_aster).
     For an explicit FSI scheme, set cs_glob_mobile_structures_n_iter_max = 1;
   */
 
  cs_glob_mobile_structures_n_iter_max = 15;
 
  /* Relative precision of sub-cycling Fluid Structure Coupling algorithm */
 
  cs_glob_mobile_structures_i_eps = 1.e-5;
 

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);
  cs_parameters_add_variable("tracer", 1);
 

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
   */
 
  cs_parameters_add_variable_variance("variance_1",
                                      "species_1");
 

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

Define the gravity direction and acceleration (no gravity forces if zero).

  {
    cs_physical_constants_t *pc = cs_get_glob_physical_constants();
 
    pc->gravity[0] = 0.;
    pc->gravity[1] = 0.;
    pc->gravity[2] = -9.81; /* gravity  (m/s2) in the z direction */
  }

Indicate the computation will use a rotating reference frame (with Coriolis forces) and specify the associated rotation.

  {
    cs_physical_constants_t *pc = cs_get_glob_physical_constants();
 
    cs_real_t omega_r[3] = {0., 0., 1.};    /* rotation vector, in rad/s */
    cs_real_t invariant[3] = {0., 0., 0.};  /* invariant point */
 
    cs_rotation_define(omega_r[0], omega_r[1], omega_r[2],
                       invariant[0], invariant[1], invariant[2]);
 
    pc->icorio = 1;  /* take Coriolis source terms into account */
  }

General options (cs_user_parameters)

Most settings can be specified in cs_user_parameters

Choose a time step option

 
  /* Time step type */
 
  cs_time_step_options_t *time_opt = cs_get_glob_time_step_options();
  time_opt->idtvar = CS_TIME_STEP_CONSTANT;
 

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;
 

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

To run a computation with a frozen velocity field (only for a restart):

  {
    cs_time_scheme_t *t_sc = cs_get_glob_time_scheme();
 
    cs_velocity_pressure_param_t *vp_param
      = cs_get_glob_velocity_pressure_param();
    cs_time_control_t *vp_tc =
      &(vp_param->time_control);
    vp_tc->interval_nt = -1;
  }

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++) {
 
      cs_field_t  *f = cs_field_by_id(f_id);
 
      if (f->type & CS_FIELD_VARIABLE) {
        cs_equation_param_t *eqp = cs_field_get_equation_param(f);
        eqp->verbosity = 2;
      }
 
    }
  }

Change the required linear-solver precision required for each increment of the solution on the solved fields:

  {
    int n_fields = cs_field_n_fields();
 
    for (int f_id = 0; f_id < n_fields; f_id++) {
 
      cs_field_t  *f = cs_field_by_id(f_id);
 
      if (f->type & CS_FIELD_VARIABLE) {
        cs_equation_param_t *eqp = cs_field_get_equation_param(f);
        eqp->epsilo = 1.e-6;
      }
 
    }
  }

To compute a dynamic relaxation for the reconstruction sweeps used to handle non-orthogonalities (for any variable):

• iswdyn = 0: no relaxation
• iswdyn = 1: means that the last increment is relaxed
• iswdyn = 2: (default) means that the last two increments are used to relax.

  {
    cs_equation_param_t *eqp = cs_field_get_equation_param(CS_F_(p));
    eqp->iswdyn = 2;
  }

For difficult cases, a stabilization may be obtained by not reconstructing the convective and diffusive flux for variables of the turbulence model, that is for k-epsilon models:

  {
    cs_equation_param_t *eqp;
 
    eqp = cs_field_get_equation_param(CS_F_(k));
    eqp->ircflu = 0;
 
    eqp = cs_field_get_equation_param(CS_F_(eps));
    eqp->ircflu = 0;
  }

Choice of turbulent flux model u'T' for the scalar T:

• Algebraic Models
  • 0) SGDH
  • 10) GGDH
  • 11) EB-GGDH (Elliptic Blending)
  • 20) AFM
  • 21 EB-AFM (Elliptic Blending)
• Models with transport equations
  • 30) 30 DFM
  • 31) EB-DFM (Elliptic Blending)

  {
    const int kturt = cs_field_key_id("turbulent_flux_model");
 
    const int turb_flux_model = 10;
 
    /* GGDH for thermal scalar */
    cs_field_t *f_t = cs_thermal_model_field();
    if (f_t != nullptr)
      cs_field_set_key_int(f_t, kturt, turb_flux_model);
 
    /* GGDH for all user-defined scalars */
    const int n_fields = cs_field_n_fields();
    const int k_scal = cs_field_key_id("scalar_id");
    const int kscavr = cs_field_key_id("first_moment_id");
 
    for (int f_id = 0; f_id < n_fields; f_id++) {
      cs_field_t *f = cs_field_by_id(f_id);
      if (f->type & (CS_FIELD_VARIABLE | CS_FIELD_USER)) {
        int s_num = cs_field_get_key_int(f, k_scal);
        int iscavr = cs_field_get_key_int(f, kscavr);
        if (s_num > 0 && iscavr <= 0)
          cs_field_set_key_int(f, kturt, turb_flux_model);
      }
    }
  }

The following snippet shows how to blend 1st-order upwind and 2nd-order centered convective schemes: To reduce numerical diffusion, it is recommened to use a second-order scheme in space for convection (the default choice). When this causes stability issues, adding a small percent of upwind may help. For exemple, to add 2 percent of upwind, set field's blencv parameter (representing the portion of upwind) to (1 - 0.02) to all variables, instead of relying on the defaults (usually pure centered or upwind depending on the variable type):

  {
    int n_fields = cs_field_n_fields();
 
    for (int f_id = 0; f_id < n_fields; f_id++) {
 
      cs_field_t *f = cs_field_by_id(f_id);
 
      if (f->type & CS_FIELD_VARIABLE) {
        cs_equation_param_t *eqp = cs_field_get_equation_param(f);
        eqp->blencv = 0.98;
      }
 
    }
  }

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 */
 
    cs_equation_param_t *eqp = cs_field_get_equation_param(sca1);
 
    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");
    cs_field_set_key_int(sca1, key_lim_id, CS_NVD_VOF_CICSAM);
 
    /* 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_field_set_key_double(CS_F_(t), kccmin, 0.);
    cs_field_set_key_double(CS_F_(t), kccmax, 1.);
 
  }

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 */
 
    cs_equation_param_t *eqp = cs_field_get_equation_param(sca1);
    eqp->blend_st = 0.1;
  }

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_coupled_with_vel_p = cs_field_key_id("coupled_with_vel_p");
 
    cs_field_set_key_int(sca1, key_coupled_with_vel_p, 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;
    drift |= CS_DRIFT_SCALAR_ADD_DRIFT_FLUX;
 
    if (false)
      drift |= CS_DRIFT_SCALAR_THERMOPHORESIS;
 
    if (false)
      drift |= CS_DRIFT_SCALAR_TURBOPHORESIS;
 
    if (false)
      drift |= CS_DRIFT_SCALAR_CENTRIFUGALFORCE;
 
    if (false)
      drift |= CS_DRIFT_SCALAR_ZERO_BNDY_FLUX;
 
    /* Set the key word "drift_scalar_model" into the field structure */
    cs_field_set_key_int(sca1, key_drift, drift);
  }

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:

  {
    cs_velocity_pressure_param_t *vp_param
      = cs_get_glob_velocity_pressure_param();
    /* Switch off Rhie & Chow filter */
    vp_param->arak = 0.;
 
    /* Activate RT0 reconstruction for the velocity from
     * mass fluxes */
    vp_param->irevmc = 1;
 
    /* Change the hydrostatic pressure algorithm */
    vp_param->iphydr = 1;
  }

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

  {
    cs_velocity_pressure_param_t *vp_param
      = cs_get_glob_velocity_pressure_param();
    vp_param->nterup = 3;
  }

The temperature scale may be chosen when solving with temperature (except for predefined physical models), as shown below. When coupling with Syrthes, the same temperature scale should be used. Wen using a radiative model, be careful with conversions, as radiation boundary conditions are always defined in Kelvin.

  {
    cs_thermal_model_t *thm = cs_get_glob_thermal_model();
 
    thm->temperature_scale = CS_TEMPERATURE_SCALE_CELSIUS;
  }

If a user-defined scalar behaves like a temperature (relative to Cp): we neeed to set its "is_temperature" keyword to 1 (this is already done for the actual temperature):

  {
    const int kscacp = cs_field_key_id("is_temperature");
 
    cs_field_set_key_int(cs_field_by_name("t_1"), kscacp, 0);
  }

Fluid properties.

Most fluid properties can be defined as part of the global cs_glob_fluid_properties structure, whose members can be modified as illustrated here:

  {
    cs_fluid_properties_t *fp = cs_get_glob_fluid_properties();
 
    /*
      ro0        : density in kg/m3
      viscl0     : dynamic viscosity in kg/(m s)
      cp0        : specific heat in J/(Kelvin kg)
      t0         : reference temperature in Kelvin
      p0         : total reference pressure in Pascal
                   the calculation is based on a
                   reduced pressure P*=Ptot-ro0*g.(x-xref)
                   (except in compressible case)
      xyzp0(3)   : coordinates of the reference point for
                   the total pressure (where it is equal to p0)
 
      In general, it is not necessary to furnish a reference point xyz0.
      If there are outlets, the code will take the center of the
      reference outlet face.
      On the other hand, if we plan to explicitly fix Dirichlet conditions
      for pressure, it is better to indicate to which reference the
      values relate (for a better resolution of reduced pressure).
 
      Other properties are given by default in all cases.
 
      Nonetheless, we may note that:
 
      In the standard case (no combustion, electric arcs, compressibility):
      ---------------------
      ro0, viscl0 and cp0
          are useful and represent either the fluid properties if they
          are constant, either simple mean values for the initialization
          if properties are variable and defined in cs_user_physical_properties.
      t0  is not useful
      p0  is useful but is not used in an equation of state. p0
          is a reference value for the incompressible solver
          which will serve to set the (possible) domain outlet pressure.
          We may also take it as 0 or as a physical value in Pascals.
 
      With the electric module:
      ------------------------
      ro0, viscl0 and cp0
          are useful but simply represent mean initial values;
          the density, molecular dynamic viscosity, and specific
          heat are necessarily defined as fields (whether they are
          physically variable or not): see cs_user_physical_properties
          for the Joule effect
          module and the electric arcs dp_ELE data file.
      t0  is useful an must be in Kelvin (> 0) but represents a simple
          initialization value.
      p0  is useful bu is not used in the equation of state. p0
          is a reference value for the incompressible solver which
          will be used to calibrate the (possible) outlet pressure
          of the domain. We may take it as zero or as a physical
          value in Pascals.
 
      With gas combustion:
      --------------------
      ro0 is not useful (it is automatically recalculated by the
          law of ideal gases from t0 and p0).
      viscl0 is indispensable: it is the molecular dynamic viscosity,
          assumed constant for the fluid.
      cp0 is indispensable: it is the heat capacity, assumed constant,
          (modelization of source terms involving a local Nusselt in
          the Lagrangian module, reference value allowing the
          calculation of a radiative
          (temperature, exchange coefficient) couple).
      t0  is indispensible and must be in Kelvin (> 0).
      p0  is indispensable and must be in Pascal (> 0).
 
      With pulverized coal:
      ---------------------
      ro0 is not useful (it is automatically recalculated by the
          law of ideal gases from t0 and p0).
      viscl0 is indispensable: it is the molecular dynamic viscosity,
          assumed constant for the fluid (its effect is expected to
          be small compared to turbulent effects).
      cp0 is indispensable: it is the heat capacity, assumed constant,
          (modelization of source terms involving a local Nusselt in
          the coal or Lagrangian module, reference value allowing the
          calculation of a radiative
          (temperature, exchange coefficient) couple).
      t0  is indispensable and must be in Kelvin (> 0).
      p0  is indispensable and must be in Pascal (> 0).
 
      With compressibility:
      ---------------------
      ro0 is not useful, stricto sensu; nonetheless, as experience
          shows that users often use this variable, it is required
          to assign to it a strictly positive value (for example,
          an initial value).
      viscl0 is useful and represents the molecular dynamic viscosity,
          when it is constant, or a value which will be used during
          initializations (or in inlet turbulence conditions,
          depending on the user choice.
      cp0 is indispensable: it is the heat capacity, assumed constant
          in the thermodynamics available by default
      t0  is indispensable and must be in Kelvin (> 0).
      p0  is indispensable and must be in Pascal (> 0).
          With the thermodynamic law available by default,
      t0 and p0 are used for the initialization of the density.
      xyzp0 is not useful because the pressure variable directly
          represents the total pressure.
    */
 
    fp->ro0    = 1.17862;
    fp->viscl0 = 1.83337e-5;
    fp->cp0    = 1017.24;
 
    fp->t0 = 20. + 273.15;
    fp->p0 = 1.01325e5;
 
    /* We only specify XYZ0 if we explicitely fix Dirichlet conditions
       for the pressure. */
 
    fp->xyzp0[0] = 0.;
    fp->xyzp0[1] = 0.;
    fp->xyzp0[2] = 0.;
 
    /* irovar, ivivar, icp: constant or variable density,
                              viscosity/diffusivity, and specific heat
 
     When a specific physics module is active
       (coal, combustion, electric arcs, compressible: see cs_user_model)
       we MUST NOT set variables 'irovar', 'ivivar', and 'icp' here, as
       they are defined automatically.
     Nonetheless, for the compressible case, ivivar may be modified
       in the uscfx2 user subroutine.
 
     When no specific physics module is active, we may specify if the
       density, specific heat, and the molecular viscosity
       are constant (irovar=0, ivivar=0, icp=-1), which is the default
       or variable (irovar=1, ivivar=1, icp=0)
 
     For those properties we choose as variable, the corresponding law
       must be defined in cs_user_physical_properties.
       if they are constant, they take values ro0, viscl0, and cp0.
    */
 
    fp->irovar = 1;
    fp->ivivar = 1;
    fp->icp = -1;
 
    /* Account the thermodynamical pressure variation in time
       (only if cs_glob_velocity_pressure_param->idilat = 3)
 
      By default:
      ----------
      - The thermodynamic pressure (pther) is initialized with p0 = p_atmos.
      - The maximum thermodynamic pressure (pthermax) is initialized with -1
        (no maximum by default, this term is used to model a venting effect when
        a positive value is given by the user).
      - A global leak can be set through a leakage surface sleak with a head
       loss kleak of 2.9 (Idelcick) */
 
    fp->ipthrm = 1;
 
    fp->pthermax= -1.;
 
    fp->sleak = 0.;
    fp->kleak = 2.9;
  }

Variable diffusivity and density

Variable diffusivity for sclar fields is handled differenty, as each scalar may be assigned a specific diffusivity.

So to each thermal and user-defined scalar field, we associate a diffusivity_id keyword, which represents the id of an associated diffusivity field when using a variable (local) diffusivity, and -1 (the default) for a constant diffusivity.

• When initially set to 0, the field will be added automatically, and later calls to cs_field_get_key_int(f, kivisl) will return a field's diffusivity field id (where kivisl is the "diffusivity_id" field keyword id).
• When set to an id > 0, the id of an existing, predefined field is given. This may allow sharing a diffusivity between multiple scalars.

Note that for user scalars which represent the variance of another user scalar, the diffusivity is assumed to be the same as that of the parent scalar, so values set here will be ignored.

For non-user scalars from specific physical models (coal, combustion, electric arcs, ...) implicitly defined in the model, the diffusivity should not be modified.

Caution: do not forget to complete cs_user_physical_properties with the law defining the diffusivity if ifcvsl >= 0 has been set here.

  {
    const int kivisl = cs_field_key_id("diffusivity_id");
 
    /* For thermal scalar */
    if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] < 0)
      cs_field_set_key_int(cs_thermal_model_field(), kivisl, -1);
    else
      cs_field_set_key_int(cs_field_by_name("temperature"), kivisl, -1);
 
    /* For user-defined scalars */
    const int n_fields = cs_field_n_fields();
    const int k_scal = cs_field_key_id("scalar_id");
    const int kscavr = cs_field_key_id("first_moment_id");
 
    for (int f_id = 0; f_id < n_fields; f_id++) {
      cs_field_t *f = cs_field_by_id(f_id);
      if (f->type & (CS_FIELD_VARIABLE | CS_FIELD_USER)) {
        int s_num = cs_field_get_key_int(f, k_scal);
        int iscavr = cs_field_get_key_int(f, kscavr);
        if (s_num > 0 && iscavr <= 0)
          cs_field_set_key_int(f, kivisl, 0);
      }
    }
  }

Variable associated density fields may be defined in a similar manner when a scalar's density is variable and different from the bulk density, which is sometimes the case for for fluid-solid computations with passive scalars with a different density in the solid. This must be consitent with continuity equation.

The logic is similar to that of variable diffusivity fields above:

  {
    const int kromsl = cs_field_key_id("density_id");
 
    /* Example user-defined scalar */
    cs_field_t *f = cs_field_by_name("scalar1");
    cs_field_set_key_int(f, kromsl, 0);
  }

For compressible flows, specific physical property settings may be defined if needed:

  {
    cs_fluid_properties_t *fp = cs_get_glob_fluid_properties();
 
    /* Reference molecular thermal conductivity
       visls0 = lambda0 (molecular thermal conductivity, W/(m K)) */
 
    const int kvisl0 = cs_field_key_id("diffusivity_ref");
 
    cs_field_set_key_double(CS_F_(t), kvisl0, 3.e-2);
 
    /* If the molecular thermal conductivity is variable, its values
       must be provided in 'cs_user_physical_properties'. */
 
    /* Reference volumetric molecular viscosity
       viscv0 = kappa0  (volumetric molecular viscosity, kg/(m s)) */
 
    fp->viscv0 = 0;
 
    /* If the volumetric molecular viscosity is variable, its values
       must be provided in 'cs_user_physical_properties' */
 
    /* Molar mass of the gas (kg/mol)
       For example with dry air, xmasml is around 28.8d-3 kg/mol */
 
    fp->xmasmr = 0.028966;
 
    /* Hydrostatic equilibrium at boundaries
       Specify if the hydrostatic equilibrium must be accounted for
       (yes = 1 , no = 0) */
 
    cs_cf_model_t *cf_model = cs_get_glob_cf_model();
 
    cf_model->icfgrp = 1;
  }

By default, the auxiliary restart file is read if present, but it may be useful to deactivate its use when restarting after a preprocessing stage possibly leading to a different number of faces (such as simply joining meshes on a different architecture or optimization level or with different options).

Writing of auxiliary restart files may also be deactivated setting cs_glob_restart_auxiliary::write_auxiliary to 0,

  cs_glob_restart_auxiliary->read_auxiliary = 0;

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

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  *f = cs_field_by_id(f_id);
 
      if (f->type & CS_FIELD_VARIABLE)
        cs_parameters_add_boundary_values(f);
 
    }
  }

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

  {
    /* add previous values handling for boundary temperature */
    cs_field_set_n_time_vals(CS_F_(t_b), 2);
  }

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.

  {
    cs_field_t *f;
 
    f = cs_field_by_name_try("yplus");
    if (f != nullptr)
      cs_parameters_add_property("yplus",
                                 1,
                                 CS_MESH_LOCATION_BOUNDARY_FACES);
 
    f = cs_field_by_name_try("tplus");
    if (f != nullptr)
      cs_parameters_add_property("tplus",
                                 1,
                                 CS_MESH_LOCATION_BOUNDARY_FACES);
 
    f = cs_field_by_name_try("tstar");
    if (f != nullptr)
      cs_parameters_add_property("tstar",
                                 1,
                                 CS_MESH_LOCATION_BOUNDARY_FACES);
  }

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

 
  cs_function_define_q_criterion();
 

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_t  *f = cs_field_by_id(f_id);
 
      if (f->type & CS_FIELD_VARIABLE)
        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:

 
  cs_field_set_key_int(CS_F_(rij), cs_field_key_id("clipping_id"), 1);
  cs_field_set_key_int(CS_F_(eps), cs_field_key_id("clipping_id"), 1);
 

Error indicators may be activated by creating one of the est_error_pre_2, est_error_der_2, est_error_cor_2, or est_error_tot_2 fields.

 
  /* We recommend running a calculation restart on a few time steps
     with the activation of the most interesting error indicators. */
 
  {
    const int field_type = CS_FIELD_INTENSIVE | CS_FIELD_POSTPROCESS;
 
    const char *name[] = {"est_error_cor_2",
                          "est_error_tot_2"};
 
    for (int i = 0; i < 2; i++) {
      cs_field_t *f = cs_field_create(name[i],
                                      field_type,
                                      CS_MESH_LOCATION_CELLS,
                                      1,
                                      false);
      cs_field_set_key_int(f, cs_field_key_id("log"), 1);
      cs_field_set_key_int(f, cs_field_key_id("post_vis"), 1);
    }
  }
 

Input-output related examples

Input-output settings for fields can be applied as soon as a field is defined. As some fields may be defined automatically based on other field's options, it is recommended to define the settings in cs_user_finalize_setup, where all fields should have been defined.

Basic options

Frequency of log output.

  cs_log_iteration_set_interval(1);

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) != nullptr)
    cs_field_set_key_str(CS_F_(cp), cs_field_key_id("label"), "Cp");

Postprocessing output

Activate or deactivate probes output.

  cs_field_set_key_int_bits(CS_F_(vel),
                            cs_field_key_id("post_vis"),
                            CS_POST_MONITOR);

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

  {
    cs_field_t *f = cs_field_by_name_try("rad_energy");
    if (f != nullptr)
      cs_field_set_key_int_bits(f,
                                cs_field_key_id("post_vis"),
                                CS_POST_ON_LOCATION | CS_POST_MONITOR);
 
    f = cs_field_by_name_try("radiative_flux");
    if (f != nullptr)
      cs_field_set_key_int_bits(f,
                                cs_field_key_id("post_vis"),
                                CS_POST_ON_LOCATION | CS_POST_MONITOR);
  }

Base model for CDO/HHO schemes

Activation of CDO/HHO schemes

Two modes are available:

• CS_PARAM_CDO_MODE_ONLY for a usage of CDO, HHO or MAC in stand-lone (i.e. without the legacy Finite Volume approach)
• CS_PARAM_CDO_MODE_WITH_FV for a usage which can share some equations/models solved with CDO/HHO schemes and some other equations/models solved with the legacy Finite Volume approach

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

  {
    cs_param_cdo_mode_set(CS_PARAM_CDO_MODE_ONLY);
  }

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_domain_t  *domain = cs_glob_domain;
    cs_boundary_t  *bdy = domain->boundaries;
 
    /* Choose a boundary by default */
    cs_boundary_set_default(bdy, CS_BOUNDARY_WALL);
 
    /* 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
     */
 
    cs_boundary_add(bdy, CS_BOUNDARY_SYMMETRY, "sym");
  }

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:

  {
    /* Manage the output logging/checkpoint
     * ====================================
     *   - log frequency in the run_solver.log
     *   - level of verbosity (-1: no, 0, 1, higher is for debugging purpose)
     *   - checkpoint/restart frequency
     */
 
    cs_domain_set_output_param(domain,
                               CS_RESTART_INTERVAL_ONLY_AT_END,
                               10,  /* log frequency */
                               2);  /* verbosity  */
  }

Time step with CDO/HHO schemes

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

  {
    /* Time step management
     * ====================
     * 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.
     */
 
    cs_domain_set_time_param(domain,
                             100,     /* nt_max or -1 (automatic) */
                             -1.);    /* t_max or < 0. (automatic) */
 
    /* Define the value of the time step. Two functions are available to do
       this:
       1. cs_domain_def_time_step_by_value(domain, dt_val);
       2. cs_domain_def_time_step_by_func(domain, dt_func);
 
       The second way enables more complex definitions of the time step.
    */
 
    cs_domain_def_time_step_by_value(domain, 1.0);
  }

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 */
    cs_walldistance_activate();
  }

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:

  {
    cs_equation_add_user("my_equation", /* name of the equation */
                         "my_variable", /* name of the associated variable */
                         1,             /* dimension of the variable */
                         CS_BC_SYMMETRY); /* default BC */
 
    /* Add a new user equation.
     *   The default boundary condition has to be chosen among:
     *    CS_BC_HMG_DIRICHLET
     *    CS_BC_SYMMETRY
     */
 
    cs_equation_add_user("AdvDiff.Upw", // equation name
                         "Pot.Upw",     // associated variable field name
                         1,             // dimension of the unknown
                         CS_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_BC_HMG_DIRICHLET); // default boundary
  }

Add user-defined properties with CDO/HHO schemes

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

  {
    /* For an anistropic property (tensor-valued) */
 
    cs_property_add("conductivity",      /* property name */
                    CS_PROPERTY_ANISO);  /* type of material property */
 
    /* For an isotropic property (scalar-valued) */
 
    cs_property_add("rho.cp",            /* property name */
                    CS_PROPERTY_ISO);    /* type of material property */
 
  }

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_set_option(pty, CS_PTYKEY_POST_FOURIER);
  }

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  *adv = cs_advection_field_add_user("adv_field");
 
    assert(adv != nullptr);
  }

If you need to activate advanced options related to advection fields, you can also specify a user-defined advection field as follows:

  {
    /* Add a user-defined advection field named "adv_field"  */
 
    cs_adv_field_t  *adv = cs_advection_field_add_user("adv_field");
 
    assert(adv != nullptr);
  }

When an advection field is defined, it is possible to retrieve it and then set a post-processing operation as follows:

  {
    /* Retrieve an advection field named "adv_field"  */
    cs_adv_field_t  *adv = cs_advection_field_by_name("adv_field");
 
    /* Compute the Courant number (if unsteady simulation) */
 
    cs_advection_field_set_postprocess(adv, CS_ADVECTION_FIELD_POST_COURANT);
  }

Define or modify general numerical and physical user parameters

CDO/HHO schemes

Modifiy the numerical parameters related to a given equation.

  {
    cs_equation_param_t  *eqp = cs_equation_param_by_name("MyEq");
 
    cs_equation_param_set(eqp, CS_EQKEY_SPACE_SCHEME, "cdo_fb");
  }

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
 
  }

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_adv_field_t  *adv = cs_advection_field_by_name("adv_field");
 
    cs_advection_field_def_by_analytic(adv, _define_adv_field, nullptr);
  }

Set up the boundary conditions for an equation

  {
    cs_equation_param_t  *eqp = cs_equation_param_by_name("AdvDiff.Upw");
 
    cs_equation_add_bc_by_analytic(eqp,
                                   CS_BC_DIRICHLET,
                                   "boundary_faces",  // zone name
                                   _define_bcs,       // pointer to the function
                                   nullptr);             // input structure
  }

Add terms to an equation

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

  {
    /* Retrieve the set of equation parameters of the equation to modify */
 
    cs_equation_param_t  *eqp = cs_equation_param_by_name("AdvDiff.Upw");
 
    /* Add an unsteady term (draw a link between this equation and a
       property) */
 
    cs_property_t  *rhocp = cs_property_by_name("rho.cp");
 
    cs_equation_add_time(eqp, rhocp);
 
    /* Add a diffusion term (draw a link between this equation and a
       property) */
 
    cs_property_t  *conductivity = cs_property_by_name("conductivity");
 
    cs_equation_add_diffusion(eqp, conductivity);
 
    /* Add an advection term (draw a link between this equation and an
       advection field) */
 
    cs_adv_field_t  *adv = cs_advection_field_by_name("adv_field");
 
    cs_equation_add_advection(eqp, adv);
  }

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 more advanced settings, possible using the cs_user_linear_solvers function 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_define(-1, "wall_distance", CS_MULTIGRID_V_CYCLE);

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 != nullptr) {
    cs_sles_it_define(cvar_user_1->id,
                      nullptr,   /* name passed is null if field_id > -1 */
                      CS_SLES_BICGSTAB2,
                      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, nullptr);
    cs_sles_set_verbosity(sles_p, 4);
  }

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, nullptr);
    cs_sles_set_post_output(sles_p, CS_POST_WRITER_DEFAULT);
 
    cs_sles_t *sles_u = cs_sles_find_or_add(CS_F_(vel)->id, nullptr);
    cs_sles_set_post_output(sles_u, CS_POST_WRITER_DEFAULT);
  }

Example: advanced multigrid settings

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

  {
    cs_multigrid_t *mg = cs_multigrid_define(CS_F_(p)->id,
                                             nullptr,
                                             CS_MULTIGRID_V_CYCLE);
 
    cs_multigrid_set_coarsening_options(mg,
                                        3,                             /* aggregation_limit (default 3) */
                                        CS_GRID_COARSENING_DEFAULT,    /* 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) */
 
    cs_multigrid_set_solver_options
      (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.

  {
    cs_sles_it_t *c = cs_sles_it_define(CS_F_(p)->id,
                                        nullptr,
                                        CS_SLES_FCG,
                                        -1,
                                        10000);
    cs_sles_pc_t *pc = cs_multigrid_pc_create(CS_MULTIGRID_V_CYCLE);
    cs_multigrid_t *mg = (cs_multigrid_t *)cs_sles_pc_get_context(pc);
    cs_sles_it_transfer_pc(c, &pc);
 
    assert(strcmp(cs_sles_pc_get_type(cs_sles_it_get_pc(c)), "multigrid") == 0);
 
    cs_multigrid_set_solver_options
      (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) */
 
  }

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_multigrid_t *mg = cs_multigrid_define(CS_F_(p)->id,
                                             nullptr,
                                             CS_MULTIGRID_V_CYCLE);
 
    cs_multigrid_set_merge_options(mg,
                                   4,    /* # of ranks merged at a time */
                                   300,  /* mean # of cells under which we merge */
                                   500); /* global # of cells under which we 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);
      cs_sles_it_define(-1,
                        name,
                        CS_SLES_JACOBI,
                        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, nullptr);
 
    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_t *c = (cs_sles_it_t *)cs_sles_get_context(sles_p);
      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_t *c = (cs_multigrid_t *)cs_sles_get_context(sles_p);
      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 */
#if PETSC_VERSION_GE(3,7,0)
    PetscOptionsSetValue(nullptr, "-ksp_type", "cg");
    PetscOptionsSetValue(nullptr, "-pc_type", "jacobi");
#else
    PetscOptionsSetValue("-ksp_type", "cg");
    PetscOptionsSetValue("-pc_type", "jacobi");
#endif
  }

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

  {
    cs_sles_petsc_define(CS_F_(p)->id,
                         nullptr,
                         MATSHELL,
                         _petsc_p_setup_hook,
                         nullptr);
 
  }

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,
                    void  *ksp_p)
{
  CS_UNUSED(context);
  KSP ksp = (KSP)ksp_p;
  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 */
#if PETSC_VERSION_GE(3,7,0)
    PetscOptionsSetValue(nullptr, "-ksp_type", "cg");
    PetscOptionsSetValue(nullptr, "-pc_type", "gamg");
    PetscOptionsSetValue(nullptr, "-pc_gamg_agg_nsmooths", "1");
    PetscOptionsSetValue(nullptr, "-mg_levels_ksp_type", "richardson");
    PetscOptionsSetValue(nullptr, "-mg_levels_pc_type", "sor");
    PetscOptionsSetValue(nullptr, "-mg_levels_ksp_max_it", "1");
    PetscOptionsSetValue(nullptr, "-pc_gamg_threshold", "0.02");
    PetscOptionsSetValue(nullptr, "-pc_gamg_reuse_interpolation", "TRUE");
    PetscOptionsSetValue(nullptr, "-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:

  {
    cs_sles_petsc_define(CS_F_(p)->id,
                         nullptr,
                         MATMPIAIJ,
                         _petsc_p_setup_hook_gamg,
                         nullptr);
 
  }

With the associated setup hook function:

static void
_petsc_p_setup_hook_gamg(void  *context,
                         void  *ksp_p)
{
  CS_UNUSED(context);
  KSP ksp = (KSP)ksp_p;
  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(nullptr, "-ksp_type", "cg");
    PetscOptionsSetValue(nullptr, "-pc_type", "hypre");
    PetscOptionsSetValue(nullptr, "-pc_hypre_type","boomeramg");
    PetscOptionsSetValue(nullptr, "-pc_hypre_boomeramg_coarsen_type", "HMIS");
    PetscOptionsSetValue(nullptr, "-pc_hypre_boomeramg_interp_type", "ext+i-cc");
    PetscOptionsSetValue(nullptr, "-pc_hypre_boomeramg_agg_nl","2");
    PetscOptionsSetValue(nullptr, "-pc_hypre_boomeramg_P_max","4");
    PetscOptionsSetValue(nullptr, "-pc_hypre_boomeramg_strong_threshold", "0.5");
    PetscOptionsSetValue(nullptr, "-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:

  {
    cs_sles_petsc_define(CS_F_(p)->id,
                         nullptr,
                         MATMPIAIJ,
                         _petsc_p_setup_hook_bamg,
                         nullptr);
 
  }

With the associated setup hook function:

static void
_petsc_p_setup_hook_bamg(void  *context,
                         void  *ksp_p)
{
  CS_UNUSED(context);
  KSP ksp = (KSP)ksp_p;
  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,
                         void  *ksp_p)
{
  CS_UNUSED(context);
 
  KSP ksp = (KSP)ksp_p;
 
  const char *p = getenv("CS_USER_PETSC_MAT_VIEW");
 
  if (p != nullptr) {
 
    /* 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) {
#if defined(HAVE_MPI)
      if (cs_glob_n_ranks > 1)
        MatView(a, PETSC_VIEWER_STDOUT_(cs_glob_mpi_comm));
#endif
      if (cs_glob_n_ranks == 1)
        MatView(a, PETSC_VIEWER_STDOUT_SELF);
    }
 
    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, nullptr, "PETSc View",
                          0, 0, 600, 600, &viewer);
      PetscViewerDrawGetDraw(viewer, 0, &draw);
      PetscViewerDrawSetPause(viewer, -1);
      MatView(a, viewer);
      PetscDrawPause(draw);
 
      PetscViewerDestroy(&viewer);
 
    }
 
  }
}

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, nullptr);
 
    cs_sles_amgx_set_config_file(amgx_p, "PCG_CLASSICAL_V_JACOBI.json");
  }

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 . 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;
    cs_time_moment_define_by_field_ids("u_mean",
                                       n_fields,
                                       moment_f_id,
                                       moment_c_id,
                                       CS_TIME_MOMENT_MEAN,
                                       1000, /* nt_start */
                                       -1,   /* t_start */
                                       CS_TIME_MOMENT_RESTART_AUTO,
                                       nullptr);
  }

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;
    cs_time_moment_define_by_field_ids("u_variance",
                                       n_fields,
                                       moment_f_id,
                                       moment_c_id,
                                       CS_TIME_MOMENT_VARIANCE,
                                       1000, /* nt_start */
                                       -1,   /* t_start */
                                       CS_TIME_MOMENT_RESTART_AUTO,
                                       nullptr);
  }

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;
    cs_time_moment_define_by_field_ids("rho_u_mean",
                                       n_fields,
                                       moment_f_id,
                                       moment_c_id,
                                       CS_TIME_MOMENT_MEAN,
                                       1000, /* nt_start */
                                       -1,   /* t_start */
                                       CS_TIME_MOMENT_RESTART_AUTO,
                                       nullptr);
  }

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;
    cs_time_moment_define_by_field_ids("rho_u_v_mean",
                                       n_fields,
                                       moment_f_id,
                                       moment_c_id,
                                       CS_TIME_MOMENT_MEAN,
                                       -1,     /* nt_start */
                                       20.0,   /* t_start */
                                       CS_TIME_MOMENT_RESTART_RESET,
                                       nullptr);

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)
{
  CS_NO_WARN_IF_UNUSED(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_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 now assign that function to a moments definition:

  {
    const char *sum_comp_name[] = {"species_sum_mean", "species_sum_variance"};
    cs_time_moment_type_t m_type[] = {CS_TIME_MOMENT_MEAN,
                                      CS_TIME_MOMENT_VARIANCE};
 
    for (int i = 0; i < 2; i++) {
      cs_time_moment_define_by_func(sum_comp_name[i],
                                    CS_MESH_LOCATION_CELLS,
                                    1,                      /* field dimension */
                                    true,                   /* intensive*/
                                    _simple_data_sum,       /* data_func */
                                    nullptr,                   /* data_input */
                                    nullptr,                   /* w_data_func */
                                    nullptr,                   /* w_data_input */
                                    m_type[i],
                                    1000,                   /* nt_start */
                                    -1,                     /* t_start */
                                    CS_TIME_MOMENT_RESTART_AUTO,
                                    nullptr);
    }
  }

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)
{
  CS_NO_WARN_IF_UNUSED(input);
 
  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_field_t *f = cs_thermal_model_field();
 
  cs_post_boundary_flux(f->name, n_elts, nullptr, vals);
}
 

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

  {
    const char *sum_comp_name[] = {"thermal_flux_mean", "thermal_flux_variance"};
    cs_time_moment_type_t m_type[] = {CS_TIME_MOMENT_MEAN,
                                      CS_TIME_MOMENT_VARIANCE};
 
    for (int i = 0; i < 2; i++) {
      cs_time_moment_define_by_func(sum_comp_name[i],
                                    CS_MESH_LOCATION_BOUNDARY_FACES,
                                    1,                      /* field dimension */
                                    true,                   /* intensive*/
                                    _boundary_thermal_flux, /* data_func */
                                    nullptr,                   /* data_input */
                                    nullptr,                   /* w_data_func */
                                    nullptr,                   /* w_data_input */
                                    m_type[i],
                                    1,                      /* nt_start */
                                    -1,                     /* t_start */
                                    CS_TIME_MOMENT_RESTART_AUTO,
                                    nullptr);
    }
  }

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
    = (const cs_real_3_t *)cs_glob_mesh_quantities->cell_cen;
 
  const cs_rotation_t *rot = cs_glob_rotation;
 
  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_real_3_t e_th;
    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_real_3_t e_r;
    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],
                                    CS_MESH_LOCATION_CELLS,
                                    1,
                                    true,                        /* intensive*/
                                    _velocity_moment_data,       /* data_func */
                                    &(vel_comp_input[comp_id]),  /* data_input */
                                    nullptr,                        /* w_data_func */
                                    nullptr,                        /* w_data_input */
                                    CS_TIME_MOMENT_MEAN,
                                    74000,                       /* nt_start */
                                    -1,                          /* t_start */
                                    CS_TIME_MOMENT_RESTART_AUTO,
                                    nullptr);
    }
  }

To activate means for all variables:

  for (int f_id = 0; f_id < cs_field_n_fields(); f_id++) {
 
    cs_field_t *f = cs_field_by_id(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 */
 
      cs_time_moment_define_by_field_ids(mean_name,
                                         n_fields,
                                         moment_f_id,
                                         moment_c_id,
                                         CS_TIME_MOMENT_MEAN,
                                         10, /* nt_start */
                                         -1, /* t_start */
                                         CS_TIME_MOMENT_RESTART_AUTO,
                                         nullptr);
    }
  }

Calculation options for the atmospheric module

In the following example, options for the atmospheric module are set in cs_user_model.

Atmospheric options are available in the following structure:

 
  cs_atmo_option_t *at_opt = cs_glob_atmo_option;
 

One can set the microphysics parameterization opitons as follows:

  /*  Microphysics parameterization options */
 
  /* Option for nucleation for humid atmosphere
   *  0: without nucleation
   *  1: Pruppacher and Klett 1997
   *  2: Cohard et al. 1998,1999
   *  3: Abdul-Razzak et al. 1998,2000
   *  logarithmic standard deviation of the log-normal law of the droplet
   *  spectrum
   */
  at_opt->nucleation_model = 3;
 
  /* Option for liquid water content distribution models
   *  1: all or nothing
   *  2: Gaussian distribution
   */
  at_opt->distribution_model = 1;
 
  /*  Option for subgrid models
   *   0: the simplest parameterization (for numerical verifications)
   *   1: Bechtold et al. 1995 (Luc Musson-Genon)
   *   2: Bouzereau et al. 2004
   *   3: Cuijpers and Duynkerke 1993, Deardorff 1976, Sommeria and
   *                 Deardorff 1977
   */
  at_opt->subgrid_model = 0;
 
  /* Sedimentation flag */
  at_opt->sedimentation_model = 1;
 
  /* Deposition flag */
  at_opt->deposition_model = 1;

One can set information relative to meteo profiles:

  /* Read the meteo file (1) or impose directly the input values to compute it
   * in code_saturne (2) */
  at_opt->meteo_profile = 2;
 
  /* set meteo file if meteo_profile == 1 */
  cs_atmo_set_meteo_file_name("meteo");
 
  /* Advanced choice of universal functions among for stable
   *  - CS_ATMO_UNIV_FN_CHENG (default)
   *  - CS_ATMO_UNIV_FN_HOGSTROM
   *  - CS_ATMO_UNIV_FN_BUSINGER
   *  - CS_ATMO_UNIV_FN_HARTOGENSIS
   *  for unstable
   *  - CS_ATMO_UNIV_FN_HOGSTROM (default)
   *  - CS_ATMO_UNIV_FN_BUSINGER
   *  - CS_ATMO_UNIV_FN_CARL
   *
   * */
 
  /* Hartogensis for stable */
  at_opt->meteo_phim_s = CS_ATMO_UNIV_FN_HARTOGENSIS;
  at_opt->meteo_phih_s = CS_ATMO_UNIV_FN_HARTOGENSIS;
 
  /* Inverse LMO length (m^-1) */
  at_opt->meteo_dlmo = 0.;
  /* Large scale roughness (m) */
  at_opt->meteo_z0 = 0.1;
  /* Elevation for reference velocity (m) */
  at_opt->meteo_zref = 10.;
  /* Friction velocity (m/s) */
  at_opt->meteo_ustar0 = 1.;
  /* Velocity direction (degrees from North) */
  at_opt->meteo_angle = 270.;
  /* Reference temperature at ground level  (K) */
  at_opt->meteo_t0 = 288.;
  /* Pressure at sea level (Pa) */
  at_opt->meteo_psea = 1.01325e5;
 
  /* Option to compute ground elevation in the domain */
  at_opt->compute_z_ground = true;
 
  /* Automatic open boundary conditions
   *   1: meteo mass flow rate is imposed with a constant large scale
   *      pressure gradient
   *   2: same plus velocity profile imposed at ingoing faces
   */
  at_opt->open_bcs_treatment = 1;
 
  /* Time of the simulation (for radiative model or chemistry)
   * syear:  starting year
   * squant: starting quantile
   * shour:  starting hour (UTC)
   * smin:   starting minute
   * ssec:   starting second
   */
  at_opt->syear = 2020;
  at_opt->squant = 1;
  at_opt->shour = 1;
  at_opt->smin = 0;
  at_opt->ssec = 0.;
 
  /* Geographic position
   *  longitude: longitude of the domain origin
   *  latitude: latitude of the domain origin
   */
  at_opt->longitude = 0.;
  at_opt->latitude = 45.0;

One can set information relative to large scale imbrication:

  /* activation flag */
  cs_glob_atmo_imbrication->imbrication_flag = false;
  cs_glob_atmo_imbrication->imbrication_verbose = false;
 
  /* flags for activating the cressman interpolation
   * for the boundary conditions */
  cs_glob_atmo_imbrication->cressman_u     = true;
  cs_glob_atmo_imbrication->cressman_v     = true;
  cs_glob_atmo_imbrication->cressman_tke   = true;
  cs_glob_atmo_imbrication->cressman_eps   = true;
  cs_glob_atmo_imbrication->cressman_theta = true;
  cs_glob_atmo_imbrication->cressman_qw    = true;
  cs_glob_atmo_imbrication->cressman_nc    = true;
 
  /* numerical parameters for the cressman interpolation formulas */
  cs_glob_atmo_imbrication->horizontal_influence_radius = 8500.;
  cs_glob_atmo_imbrication->vertical_influence_radius = 100.;

One can set information relative to meteo chemistry:

  /* Chemistry:
   *   model: choice of chemistry resolution scheme
   *     0: no atmospheric chemistry
   *     1: quasi steady equilibrium NOx scheme with 4 species and 5 reactions
   *     2: scheme with 20 species and 34 reactions
   *     3: scheme CB05 with 52 species and 155 reactions
   *     4: user defined scheme
   *      for model = 4, a SPACK file must be provided using
   *        cs_atmo_chemistry_set_spack_file_name("species.spack.dat")
 
   *      the following sources generated by SPACK should be included in the
   *      SRC folder
   *        kinetic.f90, fexchem.f90, jacdchemdc.f90, rates.f90, dratedc.f90
   *        dimensions.f90, LU_decompose.f90, LU_solve.f90
   */
  cs_glob_atmo_chemistry->model = 0;
 
  /* Default file for the chemistry profile is "chemistry" */
  cs_atmo_set_chem_conc_file_name("chem_01_01_2000");
 
  /* Chemistry with photolysis: inclusion (true) or not (false) of photolysis
   * reactions
   * warning: photolysis is not compatible with space-variable time step
   */
  cs_glob_atmo_chemistry->chemistry_with_photolysis = true;
 
  /* chemistry_sep_mode: split (=1) or semi-coupled (=2, pu-sun)
   * resolution of chemistry.
   * Split (=1) mandatory for aerosols.
   * Semi-coupled (=2) by default. */
  cs_glob_atmo_chemistry->chemistry_sep_mode = 1;
 
  /* Aerosol chemistry
   * -----------------*/
 
  /*   aerosol_model: flag to activate aerosol chemistry
   *   CS_ATMO_AEROSOL_OFF: aerosol chemistry deactivated, default
   *   CS_ATMO_AEROSOL_SSH: model automatically set to 4
   *     External library SSH-aerosol is used
   *     Corresponding SPACK files must be provided
   *     The SSH namelist file can be specified using
   *       call atmo_chemistry_set_aerosol_file_name("namelist_coag.ssh")
   *       if no namelist file is specified, "namelist.ssh" is used
   */
  cs_glob_atmo_chemistry->aerosol_model = CS_ATMO_AEROSOL_SSH;
 
  /* Default file for the aerosol profile is "aerosols" */
  cs_atmo_set_aero_conc_file_name("aero_01_01_2001");
 
  /* Frozen gaseous chemistry
   *   false: gaseous chemistry is activated (default)
   *   true: gaseous chemistry is frozen
   */
  cs_glob_atmo_chemistry->frozen_gas_chem = false;

One can set information relative to soil modelling:

  at_opt->soil_model = 1; /* Switch on soil model */
 
  /* Set the number of predefined categories (+1 which is the default one)
   * among:
   *  - CS_ATMO_SOIL_5_CAT
   *  - CS_ATMO_SOIL_7_CAT
   * */
  at_opt->soil_cat= CS_ATMO_SOIL_5_CAT; /* Switch on soil model */
 
  /* Specify the boundary zone which is modeled */
  at_opt->soil_zone_id = cs_boundary_zone_by_name("Sol")->id;

One can set information relative to 1D radiative transfer:

  /* Activate 1-D radiative transfer model */
  at_opt->radiative_model_1d = 1;
 
  /* Specify the number of verticals and the number of levels.
   * The mesh levels can be specified in cs_user_parameters function.
   * */
  at_opt->rad_1d_nvert = 1;
  at_opt->rad_1d_nlevels = 50;
  at_opt->rad_1d_nlevels_max = at_opt->rad_1d_nlevels;
 
  /* Complete 1-D mesh to ztop in case of radiative transfer */
  if (at_opt->rad_1d_nvert > 0) {
    cs_real_t zvmax = 1975.;/* top of the domain */
    cs_real_t ztop = 11000.;/* top of the troposphere */
    for (cs_real_t zzmax = (((int) zvmax)/1000)*1000.;
        zzmax <= (ztop -1000);
        zzmax += 1000.) {
      (at_opt->rad_1d_nlevels_max)++;
    }
 
  }

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