Scalar balance on full domain

Scalar balance on full domain

This is an example of cs_user_extra_operations which performs a scalar balance on the full computational domain. It is possible to customize the output to extract the contribution of some boundary zones of interest.

Define local variables

 
  cs_lnum_t n_faces;
  cs_lnum_t *face_list;
 
  cs_lnum_t cell_id, cell_id1, cell_id2, face_id;
  int nt_cur = domain->time_step->nt_cur;
 
  const cs_mesh_t *m = domain->mesh;
  const cs_mesh_quantities_t *mq = domain->mesh_quantities;
 
  const cs_lnum_t n_cells = m->n_cells;
  const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
  const cs_lnum_t n_i_faces = m->n_i_faces;
  const cs_lnum_t n_b_faces = m->n_b_faces;
 
  const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)m->i_face_cells;
  const cs_lnum_t *b_face_cells = (const cs_lnum_t *)m->b_face_cells;
 
  const cs_real_t *cell_vol = mq->cell_vol;
  const cs_real_3_t *diipb = (const cs_real_3_t *)mq->diipb;
  const cs_real_t *b_face_surf = (const cs_real_t *)mq->b_face_surf;
 

Get the physical fields

 
  cs_lnum_t n_faces;
  cs_lnum_t *face_list;
 
  cs_lnum_t cell_id, cell_id1, cell_id2, face_id;
  int nt_cur = domain->time_step->nt_cur;
 
  const cs_mesh_t *m = domain->mesh;
  const cs_mesh_quantities_t *mq = domain->mesh_quantities;
 
  const cs_lnum_t n_cells = m->n_cells;
  const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
  const cs_lnum_t n_i_faces = m->n_i_faces;
  const cs_lnum_t n_b_faces = m->n_b_faces;
 
  const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)m->i_face_cells;
  const cs_lnum_t *b_face_cells = (const cs_lnum_t *)m->b_face_cells;
 
  const cs_real_t *cell_vol = mq->cell_vol;
  const cs_real_3_t *diipb = (const cs_real_3_t *)mq->diipb;
  const cs_real_t *b_face_surf = (const cs_real_t *)mq->b_face_surf;
 

Initialization step

 
  double vol_balance = 0.;
  double div_balance = 0.;
  double a_wall_balance = 0.;
  double h_wall_balance = 0.;
  double sym_balance = 0.;
  double in_balance = 0.;
  double out_balance = 0.;
  double mass_i_balance = 0.;
  double mass_o_balance = 0.;
  double tot_balance = 0.;
 
  /* If the scalar enthalpy is not computed, return */
  if (h == NULL)
    return;
 
  /* Boundary condition coefficient for h */
  const cs_real_t *a_H = h->bc_coeffs->a;
  const cs_real_t *b_H = h->bc_coeffs->b;
  const cs_real_t *af_H = h->bc_coeffs->af;
  const cs_real_t *bf_H = h->bc_coeffs->bf;
 
  /* Convective mass fluxes for inner and boundary faces */
  int iflmas = cs_field_get_key_int(h, cs_field_key_id("inner_mass_flux_id"));
  const cs_real_t *i_mass_flux = cs_field_by_id(iflmas)->val;
 
  int iflmab = cs_field_get_key_int(h, cs_field_key_id("boundary_mass_flux_id"));
  const cs_real_t *b_mass_flux = cs_field_by_id(iflmab)->val;
 
  /* Allocate temporary array */
  cs_real_t *h_reconstructed;
  BFT_MALLOC(h_reconstructed, n_b_faces, cs_real_t);
 
  /* Reconstructed value */
  if (false) {
    cs_real_3_t *grad;
    BFT_MALLOC(grad, n_cells_ext, cs_real_3_t);
 
    cs_field_gradient_scalar(h,
                             true, /* use_previous_t */
                             1, /* inc */
                             true, /* _recompute_cocg */
                             grad);
 
    for (face_id = 0; face_id < n_b_faces; face_id++) {
      cell_id = b_face_cells[face_id]; // associated boundary cell
      h_reconstructed[face_id] = h->val[cell_id]
                               + grad[cell_id][0]*diipb[face_id][0]
                               + grad[cell_id][1]*diipb[face_id][1]
                               + grad[cell_id][2]*diipb[face_id][2];
    }
 
    BFT_FREE(grad);
 
  /* Non-reconstructed value */
  } else {
    for (face_id = 0; face_id < n_b_faces; face_id++) {
      cell_id = b_face_cells[face_id]; // associated boundary cell
      h_reconstructed[face_id] = h->val[cell_id];
    }
  }
 

Computation step

 
  /*
    --> Balance on interior volumes
        ---------------------------
  */
 
  for (cell_id = 0; cell_id < n_cells; cell_id++) {
    vol_balance += cell_vol[cell_id] * rho[cell_id]
                 * (h->val_pre[cell_id] - h->val[cell_id]);
  }
 
  /*
    --> Balance on all faces (interior and boundary), for div(rho u)
        ------------------------------------------------------------
   */
 
  for (face_id = 0; face_id < n_i_faces; face_id++) {
 
    cell_id1 = i_face_cells[face_id][0]; // associated boundary cell
    cell_id2 = i_face_cells[face_id][1]; // associated boundary cell
 
    /* Contribution to flux from the two cells of the current face
      (The cell is count only once in parallel by checking that
       the cell_id is not in the halo) */
 
    if (cell_id1 < n_cells)
      div_balance += i_mass_flux[face_id] * dt[cell_id1] * h->val[cell_id1];
 
    if (cell_id2 < n_cells)
      div_balance -= i_mass_flux[face_id] * dt[cell_id2] * h->val[cell_id2];
  }
 
  for (face_id = 0; face_id < n_b_faces; face_id++) {
 
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face */
    div_balance += b_mass_flux[face_id] * dt[cell_id] * h->val[cell_id];
 
  }
 
  // TODO mass source terms and mass accumulation term
  // In case of a mass source term, add contribution from Gamma*Tn+1
 
  /*
    --> Balance on boundary faces
        -------------------------
 
    We handle different types of boundary faces separately to better
    analyze the information, but this is not mandatory. */
 
  /*
    Compute the contribution from walls with colors 2, 3, 4 and 7
    (adiabatic here, so flux should be 0)
  */
  BFT_MALLOC(face_list, n_b_faces, cs_lnum_t);
 
  cs_selector_get_b_face_list("2 or 3 or 4 or 7", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    a_wall_balance += - b_face_surf[face_id] * dt[cell_id]
                        * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                      - b_mass_flux[face_id] * dt[cell_id]
                        * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /*
    Contribution from walls with color 6
    (here at fixed enthalpy; the convective flux should be 0)
  */
  cs_selector_get_b_face_list("6", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    h_wall_balance += - b_face_surf[face_id] * dt[cell_id]
                        * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                      - b_mass_flux[face_id] * dt[cell_id]
                        * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /*
    Contribution from symmetries (should be 0).
  */
  cs_selector_get_b_face_list("8 or 9", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    sym_balance += - b_face_surf[face_id] * dt[cell_id]
                     * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                   - b_mass_flux[face_id] * dt[cell_id]
                     * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /*
    Contribution from inlet (color 1, diffusion and convection flux)
  */
  cs_selector_get_b_face_list("1", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    in_balance += - b_face_surf[face_id] * dt[cell_id]
                    * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                  - b_mass_flux[face_id] * dt[cell_id]
                    * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /*
    Contribution from outlet (color 5, diffusion and convection flux)
  */
  cs_selector_get_b_face_list("5", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    out_balance += - b_face_surf[face_id] * dt[cell_id]
                     * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                   - b_mass_flux[face_id] * dt[cell_id]
                     * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /* Free memory */
  BFT_FREE(face_list);
  BFT_FREE(h_reconstructed);
 
  /* Sum of values on all ranks (parallel calculations) */
 
  cs_parall_sum(1, CS_DOUBLE, &vol_balance);
  cs_parall_sum(1, CS_DOUBLE, &div_balance);
  cs_parall_sum(1, CS_DOUBLE, &a_wall_balance);
  cs_parall_sum(1, CS_DOUBLE, &h_wall_balance);
  cs_parall_sum(1, CS_DOUBLE, &sym_balance);
  cs_parall_sum(1, CS_DOUBLE, &in_balance);
  cs_parall_sum(1, CS_DOUBLE, &out_balance);
  cs_parall_sum(1, CS_DOUBLE, &mass_i_balance);
  cs_parall_sum(1, CS_DOUBLE, &mass_o_balance);
 
  /* --> Total balance
         ------------- */
 
  /* We add the different contributions calculated above */
 
  tot_balance = vol_balance + div_balance + a_wall_balance + h_wall_balance
              + sym_balance + in_balance + out_balance + mass_i_balance
              + mass_o_balance;
 

Write the balance at time step n

 
  /*
    --> Balance on interior volumes
        ---------------------------
  */
 
  for (cell_id = 0; cell_id < n_cells; cell_id++) {
    vol_balance += cell_vol[cell_id] * rho[cell_id]
                 * (h->val_pre[cell_id] - h->val[cell_id]);
  }
 
  /*
    --> Balance on all faces (interior and boundary), for div(rho u)
        ------------------------------------------------------------
   */
 
  for (face_id = 0; face_id < n_i_faces; face_id++) {
 
    cell_id1 = i_face_cells[face_id][0]; // associated boundary cell
    cell_id2 = i_face_cells[face_id][1]; // associated boundary cell
 
    /* Contribution to flux from the two cells of the current face
      (The cell is count only once in parallel by checking that
       the cell_id is not in the halo) */
 
    if (cell_id1 < n_cells)
      div_balance += i_mass_flux[face_id] * dt[cell_id1] * h->val[cell_id1];
 
    if (cell_id2 < n_cells)
      div_balance -= i_mass_flux[face_id] * dt[cell_id2] * h->val[cell_id2];
  }
 
  for (face_id = 0; face_id < n_b_faces; face_id++) {
 
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face */
    div_balance += b_mass_flux[face_id] * dt[cell_id] * h->val[cell_id];
 
  }
 
  // TODO mass source terms and mass accumulation term
  // In case of a mass source term, add contribution from Gamma*Tn+1
 
  /*
    --> Balance on boundary faces
        -------------------------
 
    We handle different types of boundary faces separately to better
    analyze the information, but this is not mandatory. */
 
  /*
    Compute the contribution from walls with colors 2, 3, 4 and 7
    (adiabatic here, so flux should be 0)
  */
  BFT_MALLOC(face_list, n_b_faces, cs_lnum_t);
 
  cs_selector_get_b_face_list("2 or 3 or 4 or 7", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    a_wall_balance += - b_face_surf[face_id] * dt[cell_id]
                        * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                      - b_mass_flux[face_id] * dt[cell_id]
                        * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /*
    Contribution from walls with color 6
    (here at fixed enthalpy; the convective flux should be 0)
  */
  cs_selector_get_b_face_list("6", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    h_wall_balance += - b_face_surf[face_id] * dt[cell_id]
                        * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                      - b_mass_flux[face_id] * dt[cell_id]
                        * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /*
    Contribution from symmetries (should be 0).
  */
  cs_selector_get_b_face_list("8 or 9", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    sym_balance += - b_face_surf[face_id] * dt[cell_id]
                     * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                   - b_mass_flux[face_id] * dt[cell_id]
                     * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /*
    Contribution from inlet (color 1, diffusion and convection flux)
  */
  cs_selector_get_b_face_list("1", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    in_balance += - b_face_surf[face_id] * dt[cell_id]
                    * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                  - b_mass_flux[face_id] * dt[cell_id]
                    * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /*
    Contribution from outlet (color 5, diffusion and convection flux)
  */
  cs_selector_get_b_face_list("5", &n_faces, face_list);
 
  for (cs_lnum_t i = 0; i < n_faces; i++) {
 
    face_id = face_list[i];
    cell_id = b_face_cells[face_id]; // associated boundary cell
 
    /* Contribution to flux from the current face
      (diffusion and convection flux, negative if incoming) */
 
    out_balance += - b_face_surf[face_id] * dt[cell_id]
                     * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
                   - b_mass_flux[face_id] * dt[cell_id]
                     * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
 
  }
 
  /* Free memory */
  BFT_FREE(face_list);
  BFT_FREE(h_reconstructed);
 
  /* Sum of values on all ranks (parallel calculations) */
 
  cs_parall_sum(1, CS_DOUBLE, &vol_balance);
  cs_parall_sum(1, CS_DOUBLE, &div_balance);
  cs_parall_sum(1, CS_DOUBLE, &a_wall_balance);
  cs_parall_sum(1, CS_DOUBLE, &h_wall_balance);
  cs_parall_sum(1, CS_DOUBLE, &sym_balance);
  cs_parall_sum(1, CS_DOUBLE, &in_balance);
  cs_parall_sum(1, CS_DOUBLE, &out_balance);
  cs_parall_sum(1, CS_DOUBLE, &mass_i_balance);
  cs_parall_sum(1, CS_DOUBLE, &mass_o_balance);
 
  /* --> Total balance
         ------------- */
 
  /* We add the different contributions calculated above */
 
  tot_balance = vol_balance + div_balance + a_wall_balance + h_wall_balance
              + sym_balance + in_balance + out_balance + mass_i_balance
              + mass_o_balance;