Setting porosity values: basic example

Local definitions and initialization

Mesh quantities

It may be useful to access some mesh adjacencies and quantities, in which case local pointers allow for more readable code:

  cs_mesh_t *m = domain->mesh;
  cs_mesh_quantities_t *mq = domain->mesh_quantities;

  const cs_lnum_2_t *i_face_cells
    = (const cs_lnum_2_t *)m->i_face_cells;

  const cs_real_3_t *restrict i_face_cog
    = (const cs_real_3_t *restrict)mq->i_face_cog;
  const cs_real_3_t *restrict b_face_cog
    = (const cs_real_3_t *restrict)mq->b_face_cog;
  const cs_real_3_t *restrict cell_cen
    = (const cs_real_3_t *restrict)mq->cell_cen;
  const cs_real_3_t *restrict i_face_normal
    = (const cs_real_3_t *restrict)mq->i_face_normal;
  cs_real_3_t *restrict i_f_face_normal
    = (cs_real_3_t *restrict)mq->i_f_face_normal;
  const cs_real_3_t *restrict b_face_normal
    = (const cs_real_3_t *restrict)mq->b_face_normal;
  cs_real_3_t *restrict b_f_face_normal
    = (cs_real_3_t *restrict)mq->b_f_face_normal;

  const cs_real_t *i_f_face_surf = mq->i_f_face_surf;
  const cs_real_t *i_face_surf = mq->i_face_surf;

Associated properties

Accessing cell porosity property values is required so values may be set:

  cs_real_t *cpro_porosi = cs_field_by_name("porosity")->val;

Example: define porosity by geometric zones

Individual cell porosity values can be assigned to each cell, so they may be based on groups, geometric criteria, or any other time-independent functions:

  for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {

    cs_real_t x = cell_cen[cell_id][0];

    if (x < 20.)
      cpro_porosi[cell_id] = 1.;
    else
      cpro_porosi[cell_id] = 0.5;
  }

  /* synchronize for use in fluid face factor calculation below */
  cs_halo_type_t halo_type = CS_HALO_STANDARD;
  cs_field_synchronize(cs_field_by_name("porosity"),
                       halo_type);

Matching face equivalent surfaces should also be assigned in a corresponding manner, for interior faces:

  for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {

    cs_real_t x = i_face_cog[face_id][0];

    cs_real_t face_porosity = 1.;
    if (x > 19.9999)
      face_porosity = 0.5;

    for (int i = 0; i < 3; i++)
      i_f_face_normal[face_id][i] = face_porosity * i_face_normal[face_id][i];

    mq->i_f_face_surf[face_id] = cs_math_3_norm(i_f_face_normal[face_id]);

  }

and for boundary faces:

  for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {

    cs_real_t x = b_face_cog[face_id][0];

    cs_real_t face_porosity = 1.;
    if (x > 19.9999)
      face_porosity = 0.5;

    for (int i = 0; i < 3; i++)
      b_f_face_normal[face_id][i] = face_porosity * b_face_normal[face_id][i];

    mq->b_f_face_surf[face_id] = cs_math_3_norm(b_f_face_normal[face_id]);

  }