7.1
general documentation
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;