Definition of the variables to post-process

Additional post-processing variables

For the mesh parts defined using the GUI or in cs_user_postprocess.c, the cs_user_postprocess_values function of the cs_user_postprocess.c file may be used to specify the variables to post-process (called for each postprocess output mesh, at every active time step of an associated writer).

The output of a given variable is generated by means of a call to the cs_post_write_var for cell or face values, cs_post_write_vertex_var for vertex values, cs_post_write_particle_values for particle or trajectory values, and cs_post_write_probe_values for probe or profile values.

The examples of post-processing given below use meshes defined in the examples for cs_user_postprocess_meshes above.

Output a given field or component with a given mesh

To output specific field values on a given postprocessing mesh, that field may be attached (associated) to that mesh.

In the following example, the pressure is output on mesh 4, for all associated writers:

  cs_post_mesh_attach_field(4,
                            CS_POST_WRITER_ALL_ASSOCIATED,
                            cs_field_id_by_name("pressure"),
                            -1);

In this second example, the velocity's z component is ouput on mesh 4, only for writer 1:

  cs_post_mesh_attach_field(4,
                            1,
                            CS_F_(vel)->id,
                            2);

If the mesh has the auto_variables option enabled, and a field is already selected for "standard" postprocessing (i.e. selected for output in the GUI, or post_vis key value includes the CS_POST_ON_LOCATION bit), it will not be output an additional time, so this option should be safely usable even with limited precautions.

Output of the turbulent kinetic energy for the Rij-Epsilon model on the volume mesh

One can define, compute and post-process the turbulent kinetic energy for the Rij-Epsilon as shown in the following example:

  if (cat_id == CS_POST_MESH_VOLUME) { /* filter: only for volume
                                          postprocessing mesh */
 
    if (cs_glob_turb_model->itytur == 3) {
 
      cs_real_t *s_cell;
      BFT_MALLOC(s_cell, n_cells, cs_real_t);
 
      const cs_real_6_t *cvar_r = (const cs_real_6_t *)(CS_F_(rij)->val);
      for (cs_lnum_t i = 0; i < n_cells; i++) {
        cs_lnum_t cell_id = cell_list[i];
        s_cell[i] = 0.5* (  cvar_r[cell_id][0]
                          + cvar_r[cell_id][1]
                          + cvar_r[cell_id][2]);
      }
 
      cs_post_write_var(mesh_id,
                        CS_POST_WRITER_ALL_ASSOCIATED,  /* writer id filter */
                        "Turb energy",                  /* var_name */
                        1,                              /* var_dim */
                        true,                           /* interlace, */
                        false,                          /* use_parent */
                        CS_POST_TYPE_cs_real_t,         /* var_type */
                        s_cell,                         /* cel_vals */
                        NULL,                           /* i_face_vals */
                        NULL,                           /* b_face_vals */
                        ts);
 
      BFT_FREE(s_cell);
 
    }
 
  }

Output of a variable on a surface mesh

Values can also be output on a surface mesh, possibly containing a mix of boundary and internal faces. In the following example, we simply average or project adjacent cell values on faces, but more precise techniques could be used:

  if (strcmp(mesh_name, "pressure_surface") == 0) { /* Restrict to this mesh */
 
    cs_real_t *cvar_p = CS_F_(p)->val; /* pressure */
 
    /* Ensure variable is synchronized in parallel or periodic cases;
       should already have been done before, repeated for safety */
    cs_mesh_sync_var_scal(cvar_p);
 
    const cs_mesh_t *m = cs_glob_mesh;
 
    cs_real_t *s_i_faces = NULL, *s_b_faces = NULL;
 
    /* Interior faces  */
 
    if (n_i_faces > 0) {
      BFT_MALLOC(s_i_faces, n_i_faces, cs_real_t);
 
      for (cs_lnum_t i = 0; i < n_i_faces; i++) {
        cs_lnum_t face_id = i_face_list[i];
        /* Use unweighted mean of adjacent cell values here */
        cs_lnum_t c1 = m->i_face_cells[face_id][0];
        cs_lnum_t c2 = m->i_face_cells[face_id][1];
        s_i_faces[i] = 0.5 * (cvar_p[c1] + cvar_p[c2]);
      }
    }
 
    /* Boundary faces  */
 
    if (n_b_faces > 0) {
      BFT_MALLOC(s_b_faces, n_b_faces, cs_real_t);
 
      for (cs_lnum_t i = 0; i < n_b_faces; i++) {
        cs_lnum_t face_id = b_face_list[i];
        /* Use adjacent cell value here */
        cs_lnum_t cell_id = m->b_face_cells[face_id];
        s_b_faces[i] = cvar_p[cell_id];
      }
    }
 
    cs_post_write_var(mesh_id,
                      CS_POST_WRITER_ALL_ASSOCIATED,  /* writer id filter */
                      "Pressure",                     /* var_name */
                      1,                              /* var_dim */
                      true,                           /* interlace, */
                      false,                          /* use_parent */
                      CS_POST_TYPE_cs_real_t,         /* var_type */
                      NULL,                           /* cel_vals */
                      s_i_faces,                      /* i_face_vals */
                      s_b_faces,                      /* b_face_vals */
                      ts);
 
    BFT_FREE(s_i_faces);
    BFT_FREE(s_b_faces);
  }

Simple output of an existing field or array

For fields or arrays already defined on the full mesh, the "use_parent" option of cs_post_write_var may be used to simply reference the values on the parent (i.e. full) mesh when requesting an output. Note that the example below can also be used with probes or profiles:

  if (   cat_id == CS_POST_MESH_VOLUME
      || cat_id == CS_POST_MESH_PROBES) {
 
    const cs_field_t *f = cs_field_by_name_try("my_field");
 
    if (f != NULL)
      cs_post_write_var(mesh_id,
                        CS_POST_WRITER_ALL_ASSOCIATED,  /* writer id filter */
                        f->name,                        /* var_name */
                        1,                              /* var_dim */
                        true,                           /* interlace, */
                        true,                           /* use_parent */
                        CS_POST_TYPE_cs_real_t,         /* var_type */
                        f->val,                         /* cel_vals */
                        NULL,                           /* i_face_vals */
                        NULL,                           /* b_face_vals */
                        ts);
  }

Single output of time-independent values

Finally, a minor modification f the above example shows how it is possible to output time-independent values to a writer also used for time-dependent fields without requiring multiple outputs of those values:

  if (cat_id == CS_POST_MESH_VOLUME) {
 
    const cs_field_t *f = cs_field_by_name_try("my_const_field");
 
    if (f != NULL) {
      if (ts->nt_cur == ts->nt_prev + 1) { /* before time loop */
 
        cs_time_step_t ts0 = *ts;
        ts0.nt_cur = 1; /* Negative time step value implies time-independent */
 
        cs_post_write_var(mesh_id,
                          CS_POST_WRITER_ALL_ASSOCIATED,  /* writer id filter */
                          f->name,                        /* var_name */
                          1,                              /* var_dim */
                          true,                           /* interlace, */
                          true,                           /* use_parent */
                          CS_POST_TYPE_cs_real_t,         /* var_type */
                          f->val,                         /* cel_vals */
                          NULL,                           /* i_face_vals */
                          NULL,                           /* b_face_vals */
                          &ts0);
 
      }
    }
 
  }

Additional profile variables

The following examples match the advanced profile definitions given in Advanced profile definitions.

The first section is common to both profile series:

    const cs_mesh_t *m = cs_glob_mesh;
    const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
 
    const cs_real_3_t *cell_cen = (const cs_real_3_t *)mq->cell_cen;
 
    const cs_real_3_t *vel = (const cs_real_3_t *)CS_F_(vel)->val;
    const cs_real_t uref = cs_glob_turb_ref_values->uref;
    const cs_real_t uref2 = uref*uref;
 
    const char *name = cs_probe_set_get_name(probes);
 
    bool is_profile = false;
    cs_probe_set_get_post_info(probes,
                               NULL,
                               NULL,
                               &is_profile,
                               NULL,
                               NULL,
                               NULL,
                               NULL,
                               NULL);
 
    const cs_time_step_t *ts_post = ts;
    if (is_profile && ts->nt_cur == ts->nt_max)
      ts = NULL;
 
    /* Common variables */
 
    const cs_real_t href = 1.;

For the profiles along fixed x, the following code is used. Note that this code's complexity is mainly due to extracting Reynolds stresses for different turbulence models and options. Specific values are then computed for each colum, in the switch statement:

    if (strncmp(name, "buicesat", strlen("buicesat")) == 0) {
 
      cs_real_t *val;
      BFT_MALLOC(val, n_cells, cs_real_t);
 
      char var_name[64];
 
      /* mean x */
 
      cs_real_t x_sum[] = {0, 0};
      for (cs_lnum_t i = 0; i < n_cells; i++) {
        cs_lnum_t c_id = cell_list[i];
        x_sum[0] += cell_cen[c_id][0];
      }
      x_sum[1] = n_cells;
      cs_parall_sum(2, CS_REAL_TYPE, x_sum);
 
      cs_real_t xpos = x_sum[0]/x_sum[1];
 
      /* Reynolds stresses */
 
      const cs_turb_model_t *turb_mdl = cs_glob_turb_model;
      const cs_turb_rans_model_t *turb_rans_mdl = cs_glob_turb_rans_model;
 
      cs_real_6_t *rij = NULL;
      BFT_MALLOC(rij, n_cells, cs_real_6_t);
 
      if (   turb_mdl->itytur == 2
          || turb_mdl->itytur == 5
          || turb_mdl->itytur == 6) {
 
        cs_field_interpolate_t interpolation_type = CS_FIELD_INTERPOLATE_MEAN;
        cs_post_evm_reynolds_stresses(interpolation_type,
                                      n_cells,
                                      cell_list,
                                      NULL, /* coords */
                                      rij);
 
      }
      else if (turb_mdl->itytur == 3) {
 
        cs_real_6_t *cvar_rij = (cs_real_6_t *)CS_F_(rij)->val;
        for (cs_lnum_t i = 0; i < n_cells; i++) {
          cs_lnum_t c_id = cell_list[i];
          for (cs_lnum_t j = 0; j < 6; j++)
            rij[i][j] = cvar_rij[c_id][j];
        }
 
      }
 
      /* Reynolds stresses invariants:
       * compute xsi and eta invariant of the Lumley triangle */
 
      cs_real_2_t *inv = NULL;
      BFT_MALLOC(inv, n_cells, cs_real_2_t);
 
 
      cs_post_anisotropy_invariant(n_cells,
                                   cell_list,
                                   NULL, /* coords */
                                   inv);
 
      /* Loop on columns */
 
      for (int col = 0; col < 9; col++) {
 
        switch(col) {
 
        case 0:
          {
            strncpy(var_name, "U*10+x/h", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              cs_lnum_t c_id = cell_list[i];
              val[i] = vel[c_id][0]*10 + xpos;
            }
          }
          break;
 
        case 1:
          {
            strncpy(var_name, "Y/H", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              cs_lnum_t c_id = cell_list[i];
              val[i] = mq->cell_cen[c_id*3 + 1] / href;
            }
          }
          break;
 
        case 2:
          {
            strncpy(var_name, "U/Uc", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              cs_lnum_t c_id = cell_list[i];
              val[i] = vel[c_id][0] / uref;
            }
          }
          break;
 
        case 3:
          {
            strncpy(var_name, "uu/Uc^2", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              val[i] = rij[i][0] / uref2;
            }
          }
          break;
 
        case 4:
          {
            strncpy(var_name, "uv/Uc^2", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              val[i] = rij[i][3] / uref2;
            }
          }
          break;
 
        case 5:
          {
            strncpy(var_name, "vv/Uc^2", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              val[i] = rij[i][1] / uref2;
            }
          }
          break;
 
        case 6:
          {
            strncpy(var_name, "X", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              cs_lnum_t c_id = cell_list[i];
              val[i] = cell_cen[c_id][0];
            }
          }
          break;
 
        case 7:
          {
            strncpy(var_name, "eta", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              cs_lnum_t c_id = cell_list[i];
              val[i] = inv[c_id][0];
            }
          }
          break;
 
        case 8:
          {
            strncpy(var_name, "xsi", 64);
            for (cs_lnum_t i = 0; i < n_cells; i++) {
              cs_lnum_t c_id = cell_list[i];
              val[i] = inv[c_id][1];
            }
          }
          break;
 
        }
 
        cs_post_write_probe_values
          (mesh_id,
           CS_POST_WRITER_ALL_ASSOCIATED,  /* writer id filter */
           var_name,                       /* var_name */
           1,                              /* var_dim */
           CS_POST_TYPE_cs_real_t,
           0,                              /* parent location id */
           NULL,                           /* default interpolation */
           NULL,                           /* interpolation input */
           val,
           ts_post);
 
      }
 
      BFT_FREE(inv);
      BFT_FREE(rij);
      BFT_FREE(val);
 
    }

For the profile defined all around a foil, the following code is used to compute the pressure coefficient and output its values:

      const cs_mesh_t *m = cs_glob_mesh;
      const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
      const cs_real_3_t *b_face_cog = (const cs_real_3_t *)mq->b_face_cog;

 
      cs_real_t *val;
      BFT_MALLOC(val, n_b_faces, cs_real_t);
 
      /* x coordinate */
      for (cs_lnum_t i = 0; i < n_b_faces; i++) {
        cs_lnum_t face_id = b_face_list[i];
        val[i] = b_face_cog[face_id][0];
      }
 
      /* post-process x coordinate */
      cs_post_write_probe_values
        (mesh_id,
         CS_POST_WRITER_ALL_ASSOCIATED,  /* writer id filter */
         "X",                           /* var_name */
         1,                              /* var_dim */
         CS_POST_TYPE_cs_real_t,
         0,                              /* parent location id */
         NULL,                           /* default interpolation */
         NULL,                           /* interpolation input */
         val,
         ts);
 
      /* compute pressure on selected boundary faces */
      cs_post_b_pressure(n_b_faces, b_face_list, val);
 
      const cs_fluid_properties_t *phys_pro = cs_get_glob_fluid_properties();
      cs_real_t p0 = phys_pro->p0; /* reference pressure */
      cs_real_t ro0 = phys_pro->ro0; /* reference density */
 
      const cs_real_t uref = cs_glob_turb_ref_values->uref; /*ref. velocity */
      const cs_real_t uref2 = uref*uref;
 
      /* reference values can be set in GUI */
 
      /* 1/(1/2 rho U^2) */
      cs_real_t div_half_ro0_uref2 = 1. / (0.5 * ro0 * uref2);
 
      /* compute CP at each selected boundary face */
      for (cs_lnum_t i = 0; i < n_b_faces; i++)
        val[i] = (val[i] - p0) * div_half_ro0_uref2;
 
      /* post-process CP */
      cs_post_write_probe_values
        (mesh_id,
         CS_POST_WRITER_ALL_ASSOCIATED,  /* writer id filter */
         "CP",                           /* var_name */
         1,                              /* var_dim */
         CS_POST_TYPE_cs_real_t,
         0,                              /* parent location id */
         NULL,                           /* default interpolation */
         NULL,                           /* interpolation input */
         val,
         ts);
 
      BFT_FREE(val);

For the last profiles series, values for each column are also computed, requiring a reference pressure based on the mesh point closest to a given point, and computation of tangential stresses, so as to determine drag coefficients.

    else if (   strcmp(name, "buicstr") == 0
             || strcmp(name, "buicinc") == 0) {
 
      const cs_lnum_t *b_face_cells = m->b_face_cells;
      const cs_real_3_t *face_cog = (const cs_real_3_t *)mq->b_face_cog;
      const cs_real_t *distb = mq->b_dist;
      const cs_fluid_properties_t *phys_pro = cs_get_glob_fluid_properties();
      const cs_real_t *pres = CS_F_(p)->val;
 
      cs_real_t div_half_ro0_uref2 = 1. / (0.5 * phys_pro->ro0 * uref2);
 
      cs_real_t *val;
      BFT_MALLOC(val, n_b_faces, cs_real_t);
 
      char var_name[64];
 
      /* Reference pressure */
 
      cs_real_t xyz_ref[3] = {-1.7, -0.5, 0.005};
 
      cs_real_t pref = 0;
      cs_lnum_t pref_id;
      int pref_rank;
 
      cs_geom_closest_point(m->n_cells,
                            (const cs_real_3_t *)(mq->cell_cen),
                            xyz_ref,
                            &pref_id,
                            &pref_rank);
 
      if (pref_rank == cs_glob_rank_id)
        pref = pres[pref_id];
      cs_parall_bcast(pref_rank, 1, CS_REAL_TYPE, &pref);
 
      /* Stresses */
 
      cs_real_3_t *stresses;
      BFT_MALLOC(stresses, n_b_faces, cs_real_3_t);
      cs_post_stress_tangential(n_b_faces, b_face_list, stresses);
 
      /* Loop on columns */
 
      for (int col = 0; col < 5; col++) {
 
        switch(col) {
 
        case 0:
          {
            strncpy(var_name, "X/H", 64);
            for (cs_lnum_t i = 0; i < n_b_faces; i++) {
              cs_lnum_t f_id = b_face_list[i];
              val[i] = face_cog[f_id][0] / href;
            }
          }
          break;
 
        case 1:
          {
            strncpy(var_name, "CP", 64);
            for (cs_lnum_t i = 0; i < n_b_faces; i++) {
              cs_lnum_t f_id = b_face_list[i];
              cs_lnum_t c_id = b_face_cells[f_id];
              val[i] = (pres[c_id] - pref) * div_half_ro0_uref2;
            }
          }
          break;
 
        case 2:
          {
            strncpy(var_name, "CF", 64);
            for (cs_lnum_t i = 0; i < n_b_faces; i++) {
              val[i] = cs_math_3_norm(stresses[i]) * div_half_ro0_uref2;
            }
          }
          break;
 
        case 3:
          {
            strncpy(var_name, "U/UREF", 64);
            for (cs_lnum_t i = 0; i < n_b_faces; i++) {
              /* previous value in val[i] from case2:
                 norm(stresses[i])/(0.5.ro0.uref^2) */
              val[i] = copysign(val[i], stresses[i][0]);
            }
          }
          break;
 
        case 4:
          {
            strncpy(var_name, "YPLUS", 64);
            for (cs_lnum_t i = 0; i < n_b_faces; i++) {
              cs_lnum_t f_id = b_face_list[i];
              cs_lnum_t c_id = b_face_cells[f_id];
              val[i] = sqrt(fabs(vel[c_id][0])*distb[f_id]*phys_pro->viscl0);
            }
          }
          break;
 
        }
 
        cs_post_write_probe_values
          (mesh_id,
           CS_POST_WRITER_ALL_ASSOCIATED,  /* writer id filter */
           var_name,                       /* var_name */
           1,                              /* var_dim */
           CS_POST_TYPE_cs_real_t,
           0,                              /* parent location id */
           NULL,                           /* default interpolation */
           NULL,                           /* interpolation input */
           val,
           ts_post);
 
      }
 
      BFT_FREE(stresses);
      BFT_FREE(val);
    }