8.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;
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:319
cs_lnum_2_t
cs_lnum_t cs_lnum_2_t[2]
vector of 2 local mesh-entity ids
Definition: cs_defs.h:327
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:334
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:313
mesh::diipb
double precision, dimension(:,:), pointer diipb
vector II' for interior faces for every boundary face, the three components of the vector ....
Definition: mesh.f90:168
cs_mesh_quantities_t
Definition: cs_mesh_quantities.h:92
cs_mesh_quantities_t::b_face_surf
cs_real_t * b_face_surf
Definition: cs_mesh_quantities.h:118
cs_mesh_quantities_t::cell_vol
cs_real_t * cell_vol
Definition: cs_mesh_quantities.h:97
cs_mesh_quantities_t::diipb
cs_real_t * diipb
Definition: cs_mesh_quantities.h:131
cs_mesh_t
Definition: cs_mesh.h:85
cs_mesh_t::b_face_cells
cs_lnum_t * b_face_cells
Definition: cs_mesh.h:112
cs_mesh_t::n_i_faces
cs_lnum_t n_i_faces
Definition: cs_mesh.h:98
cs_mesh_t::n_b_faces
cs_lnum_t n_b_faces
Definition: cs_mesh.h:99
cs_mesh_t::n_cells_with_ghosts
cs_lnum_t n_cells_with_ghosts
Definition: cs_mesh.h:151
cs_mesh_t::n_cells
cs_lnum_t n_cells
Definition: cs_mesh.h:97
cs_mesh_t::i_face_cells
cs_lnum_2_t * i_face_cells
Definition: cs_mesh.h:111

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 */
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];
}
}
BFT_MALLOC
#define BFT_MALLOC(_ptr, _ni, _type)
Allocate memory for _ni elements of type _type.
Definition: bft_mem.h:62
BFT_FREE
#define BFT_FREE(_ptr)
Free allocated memory.
Definition: bft_mem.h:101
cs_field_by_id
cs_field_t * cs_field_by_id(int id)
Return a pointer to a field based on its id.
Definition: cs_field.c:2320
cs_field_get_key_int
int cs_field_get_key_int(const cs_field_t *f, int key_id)
Return a integer value for a given key associated with a field.
Definition: cs_field.c:3068
cs_field_key_id
int cs_field_key_id(const char *name)
Return an id associated with a given key name.
Definition: cs_field.c:2574
cs_field_gradient_scalar
void cs_field_gradient_scalar(const cs_field_t *f, bool use_previous_t, int inc, cs_real_3_t *restrict grad)
Compute cell gradient of scalar field or component of vector or tensor field.
Definition: cs_field_operator.c:470
h
@ h
Definition: cs_field_pointer.h:91
cs_field_t::val
cs_real_t * val
Definition: cs_field.h:152

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;
CS_DOUBLE
@ CS_DOUBLE
Definition: cs_defs.h:277
rho
@ rho
Definition: cs_field_pointer.h:97
dt
@ dt
Definition: cs_field_pointer.h:65
cs_parall_sum
static void cs_parall_sum(int n, cs_datatype_t datatype, void *val)
Sum values of a given datatype on all default communicator processes.
Definition: cs_parall.h:154
cs_selector_get_b_face_list
void cs_selector_get_b_face_list(const char *criteria, cs_lnum_t *n_b_faces, cs_lnum_t b_face_list[])
Fill a list of boundary faces verifying a given selection criteria.
Definition: cs_selector.c:213

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;