8.2
general documentation
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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:332
cs_lnum_2_t
cs_lnum_t cs_lnum_2_t[2]
vector of 2 local mesh-entity ids
Definition: cs_defs.h:340
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:347
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:325
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:97
BFT_FREE
#define BFT_FREE(_ptr)
Free allocated memory.
Definition: bft_mem.h:136
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:2465
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:3213
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:2719
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) */
double b_tot[9] = {
vol_balance, div_balance, a_wall_balance, h_wall_balance,
sym_balance, in_balance, out_balance, mass_i_balance, mass_o_balance
};
cs_parall_sum(9, CS_DOUBLE, b_tot);
vol_balance = b_tot[0], div_balance = b_tot[1];
a_wall_balance = b_tot[2], h_wall_balance = b_tot[3];
sym_balance = b_tot[4], in_balance = b_tot[5], out_balance = b_tot[6];
mass_i_balance = b_tot[7], mass_o_balance = b_tot[8];
/* 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:289
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) */
double b_tot[9] = {
vol_balance, div_balance, a_wall_balance, h_wall_balance,
sym_balance, in_balance, out_balance, mass_i_balance, mass_o_balance
};
cs_parall_sum(9, CS_DOUBLE, b_tot);
vol_balance = b_tot[0], div_balance = b_tot[1];
a_wall_balance = b_tot[2], h_wall_balance = b_tot[3];
sym_balance = b_tot[4], in_balance = b_tot[5], out_balance = b_tot[6];
mass_i_balance = b_tot[7], mass_o_balance = b_tot[8];
/* 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;