Basic examples

The preferred method for defining boundary conditions (besides using the GUI) is to assign zone-based definitions in the cs_user_boundary_conditions_setup function (cs_user_boundary_conditions.c).

As an alternative, the same definitions may be defined in cs_user_finalize_setup_wrapper (cs_user_parameters), which is called immediately after cs_user_boundary_conditions_setup. Choosing one or the other should be based on the best balance for reability and maintainability,

Simple constant-valued boundary conditions.

The following examples illustrate how different types of boundary conditions can be assigned to given variables at specified zones. In each case, an input array is provided. Values are copied internally, so the array does not need to continue existing after the boundary condition is defined.

Constant value or flux

Imposed values (Dirichlet conditions) can be defined using an array of the same dimension as the variable, as shown here for a velocity vector (1.5, 0, 0) at a zone named "inlet":

    cs_equation_param_t  *eqp = cs_equation_param_by_name("velocity");

    cs_real_t inlet_velocity[] = {1.5, 0, 0};

    cs_equation_add_bc_by_value(eqp,
                                CS_PARAM_BC_DIRICHLET,
                                "inlet",           // zone name
                                inlet_velocity);

Simple Neuman conditions can be defined in a simular manner, replacing CS_PARAM_BC_DIRICHLET with CS_PARAM_BC_NEUMANN.

Constant exchange coefficient

For exchange-condition-type boundary conditions, where a given exchange coefficient and exterior value are provided, a Robin boundary condition can be used. In this case, a boundary condition of the form: can be used by providing the values (alpha, u0, g), where alpha is is the exchange coefficient mulitiplied by -1 (contribution fro mboundary to domain sign convention), and u0 is the external value.

In the following example, and exchange coefficient of value 10.5 and an external value of 0.1 is appplied to a scalar (named "scalar1") on a zone named "exchanger_wall":

    /* Robin BC: K du/dn + alpha*(u - u0) = g
       with alpha = -10.5 and u0 = 0.1 */

    cs_real_t robin_values[3] = {-10.5, 0.1, 0};

    cs_equation_param_t  *eqp = cs_equation_param_by_name("scalar1");

    cs_equation_add_bc_by_value(eqp,
                                CS_PARAM_BC_ROBIN,
                                "exchanger_wall",         // zone name
                                robin_values);

Common boundary conditions types with local values.

The following examples illustrate how slightly more advanced boundary conditions can be assigned with user-defined callback functions. Here, we use "static" (i.e. local) functions that must be defined in the same file (i.e. cs_user_boundary_conditions.c) prior to being referenced in cs_user_boundary_conditions_define. It is of course possible to use functions with global linkage, which can be defined in any file, as log as the matching function prototypes are defined or included at the calling point.

Note the as many callback functions as required may be used, and can be named in any manner chosen by the user (though of course using a consistent naming scheme is recommended for readability and maintainability).

Analytic functions

Analytic functions of type cs_analytic_func_t are one of several callback function types that may be used to define boundary conditions. For each zone, they are called for a set of quadrature points at which values are evaluated. For legacy schemes, the associated quadrature points are the face centers of the selected zone. For CDO-based solvers, the set of points depends on the quadrature type.

In this example, analytic functions are used to define the velocity and a scalar. Here, the following functions are first defined, for the velocity:

static void
_vel_profile(cs_real_t           time,
             cs_lnum_t           n_elts,
             const cs_lnum_t    *elt_ids,
             const cs_real_t    *coords,
             bool                dense_output,
             void               *input,
             cs_real_t          *val)
{
  CS_UNUSED(time);
  CS_UNUSED(input);

  const cs_real_3_t *elt_coords = (const cs_real_3_t *)coords;

  cs_real_3_t  *v = (cs_real_3_t *)val;

  for (cs_lnum_t i = 0; i < n_elts; i++) {

    const cs_lnum_t  elt_id = (elt_ids == NULL) ? i : elt_ids[i];
    const cs_lnum_t  j = dense_output ? i : elt_id;

    cs_real_t  y = elt_coords[elt_id][1] - 0.5;
    cs_real_t  z = elt_coords[elt_id][2] - 0.5;
    cs_real_t r = sqrt(y*y + z*z);

    v[j][0] = 1.5 * (1 - r);
    v[j][1] = 0;
    v[j][2] = 0;

  }
}

And for the scalar:

static void
_scalar_inlet_profile(cs_real_t           time,
                      cs_lnum_t           n_elts,
                      const cs_lnum_t    *elt_ids,
                      const cs_real_t    *coords,
                      bool                dense_output,
                      void               *input,
                      cs_real_t          *val)
{
  CS_UNUSED(time);

  const cs_real_3_t *elt_coords = (const cs_real_3_t *)coords;

  const cs_field_t *f = input;  /* field pointer passed as input
                                   upon assignment */

  if (strcmp(f->name, "scalar1") == 0) {

    for (cs_lnum_t i = 0; i < n_elts; i++) {

      const cs_lnum_t  elt_id = (elt_ids == NULL) ? i : elt_ids[i];
      const cs_lnum_t  j = dense_output ? i : elt_id;

      cs_real_t  y = elt_coords[elt_id][1] - 0.5;
      cs_real_t  z = elt_coords[elt_id][2] - 0.5;
      cs_real_t r = sqrt(y*y + z*z);

      val[j] = 4.0 * (1 - r);

    }
  }
}

They are then associated to the matching zone in cs_user_boundary_conditions_define:

    cs_equation_param_t  *eqp = cs_equation_param_by_name("velocity");

    cs_equation_add_bc_by_analytic(eqp,
                                   CS_PARAM_BC_DIRICHLET,
                                   "inlet",           // zone name
                                   _vel_profile,      // callback function
                                   NULL);             // input structure

and

    cs_field_t  *f = cs_field_by_name("scalar1");
    cs_equation_param_t  *eqp = cs_equation_param_by_name("scalar1");

    cs_equation_add_bc_by_analytic(eqp,
                                   CS_PARAM_BC_DIRICHLET,
                                   "inlet",                 // zone name
                                   _scalar_inlet_profile,   // callback function
                                   f);                      // input structure

Boundary-face-based functions.

Degree-of-freedom-based callbacks of type cs_dof_func_t can be used to define boundary conditions at known locations, such as boundary faces for a given zone, independently of the quadrature actually used. For legacy schemes, boundary faces represent the natural boundary value locations, so using analytic or dof-base functions is equivalent (though point coordinates are not provided in dof-function arguments, so must be accessed through the mesh quantities if needed). For CDO schemes, values are automatically interpolated from the provided boundary face centers to the actual quadrature points.

In the following example, a function based on boundary face values of type cs_dof_func_t is used to define a Robin boundary condition based on an exchange coefficient and external value. As Robin conditions are expressed in the form: for each location pont, a tuple of values (alpha, u0, g) is provided with interleaved storage. This can be seen in the definition of the associated function:

static void
_scalar_exchange_profile(cs_lnum_t         n_elts,
                         const cs_lnum_t  *elt_ids,
                         bool              dense_output,
                         void             *input,
                         cs_real_t        *retval)
{
  CS_UNUSED(input);

  const cs_lnum_t dim = 1;
  const cs_lnum_t stride = 1 + dim + dim*dim;

  /* Face center coordinates */

  const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
  const cs_real_3_t *f_cog = (const cs_real_3_t *)mq->b_face_cog;

  /* Exchange coefficient first, Dirichlet values second. */

  for (cs_lnum_t i = 0; i < n_elts; i++) {
    const cs_lnum_t  face_id = (elt_ids == NULL) ? i : elt_ids[i];
    const cs_lnum_t  j = dense_output ? i : face_id;
    retval[j*stride]   = - 10.5;
    retval[j*stride+1] = 25. + 0.1*f_cog[face_id][0];
    retval[j*stride*3+2] = 0;
  }
}

As usual, the function is then used as a callback for for Robin condition:

    cs_field_t  *f = cs_field_by_name("scalar1");
    cs_equation_param_t  *eqp = cs_equation_param_by_name("scalar1");

    cs_equation_add_bc_by_dof_func(eqp,
                                   CS_PARAM_BC_ROBIN,
                                   "exchanger_wall",         // zone name
                                   cs_flag_boundary_face,    // location flag
                                   _scalar_exchange_profile, // callback function
                                   f);                       // input structure

Scaling by zone surface.

To define even a constant boundary value using a scaling variable not yet known at setup time, an appropriate callback function must be used. This is the case, for division by a zone's surface, which is unknown when calling cs_user_boundary_conditions_setup (where the mesh has not yet been read or built).

As an example, to define a unit flux such that the total flux over a zone has value 1, we can simply divide it by a zone's measure, such as in the following function for a zone named "wall_top":

static void
_w_flux_top(cs_lnum_t         n_elts,
            const cs_lnum_t  *elt_ids,
            bool              dense_output,
            void             *input,
            cs_real_t        *retval)
{
  CS_NO_WARN_IF_UNUSED(input);

  const cs_zone_t *z = cs_boundary_zone_by_name("wall_top");

  /* Get the fluid measure (i.e. surface) of the zone */
  cs_real_t flux = 1. / z->f_measure;

  /* Exchange coefficient first, Dirichlet values second. */

  for (cs_lnum_t i = 0; i < n_elts; i++) {
    const cs_lnum_t  face_id = (elt_ids == NULL) ? i : elt_ids[i];
    const cs_lnum_t  j = dense_output ? i : face_id;
    retval[j] = flux;
  }
}

which we associate to the given zone and variable:

    cs_equation_param_t  *eqp = cs_equation_param_by_name("scalar_1");

    cs_equation_add_bc_by_dof_func(eqp,
                                   CS_PARAM_BC_NEUMANN,
                                   "wall_top",             // zone name
                                   cs_flag_boundary_face,  // location flag
                                   _w_flux_top,            // callback function
                                   NULL);                  // input structure

A similar definition is done for zone "wall_side", with a slightly improved function definition:

static void
_w_flux_side(cs_lnum_t         n_elts,
             const cs_lnum_t  *elt_ids,
             bool              dense_output,
             void             *input,
             cs_real_t        *retval)
{
  const cs_zone_t *z = (const cs_zone_t *)input;

  /* Get the fluid measure (i.e. surface) of the zone */
  cs_real_t flux = 1. / z->f_measure;

  /* Exchange coefficient first, Dirichlet values second. */

  for (cs_lnum_t i = 0; i < n_elts; i++) {
    const cs_lnum_t  face_id = (elt_ids == NULL) ? i : elt_ids[i];
    const cs_lnum_t  j = dense_output ? i : face_id;
    retval[j] = flux;
  }
}

This function is similar to the one used for "wall_top", but assumes the associated input pointer is a pointer to the associated zone, which can thus be obtained directly using a pointer cast.

In this case, the matching zone pointer must be passed as the associated input when assigning the condition:

    const cs_zone_t *z = cs_boundary_zone_by_name("wall_side");
    cs_equation_param_t  *eqp = cs_equation_param_by_name("scalar_1");

    cs_equation_add_bc_by_dof_func(eqp,
                                   CS_PARAM_BC_NEUMANN,
                                   z->name,                // zone name
                                   cs_flag_boundary_face,  // location flag
                                   _w_flux_side,           // callback function
                                   (void *)z);             // input structure

Note that a structure passed as a callback input must still point to valid data when the callback is executed. This is the case for zone pointers, even for time-varying zones (whose data may change, but base structure pointer is unvarying).