7.2
general documentation
Calculation of a local Nusselt number

This is an example of cs_user_extra_operations function which computes a Nusselt number. balances on specified zones.

Faces can be selected either using a selector (as is done here) or directly using boundary zone element ids if an appropriate zone has been defined.

const cs_lnum_t n_cells = domain->mesh->n_cells;
const cs_lnum_t n_b_faces = domain->mesh->n_b_faces;
const cs_lnum_t *b_face_cells = domain->mesh->b_face_cells;
const cs_real_3_t *cell_cen
= (const cs_real_3_t *)domain->mesh_quantities->cell_cen;
= (const cs_real_3_t *)domain->mesh_quantities->b_face_cog;
const cs_real_t *volume = domain->mesh_quantities->cell_vol;
FILE *file = NULL;
cs_real_3_t *vel = (cs_real_3_t *)CS_F_(vel)->val;
/* Parameters, case dependant */
const cs_real_t prandtl = 0.71, height = 1., qwall = 1.;
const cs_field_t *f = CS_F_(t);
const cs_field_t *mu = CS_F_(mu);
cs_lnum_t nlelt;
cs_lnum_t *lstelt;
BFT_MALLOC(lstelt, n_b_faces, cs_lnum_t);
("normal[0, -1, 0, 0.1] and box[-1000, -1000, -1000, 1000, 0.01, 1000]",
&nlelt, lstelt);

Allocate a local array for the selected boundary faces.

cs_real_t *treloc;
BFT_MALLOC(treloc, nlelt, cs_real_t);
int iortho = 0;
if (iortho == 0) {
cs_post_field_cell_to_b_face_values(f, nlelt, lstelt, treloc);
}

Compute boundary value without reconstruction

const cs_real_t *coefap = CS_F_(t)->bc_coeffs->a;
const cs_real_t *coefbp = CS_F_(t)->bc_coeffs->b;
for (cs_lnum_t ielt = 0; ielt < nlelt; ielt++) {
cs_lnum_t face_id = lstelt[ielt];
cs_lnum_t c_id = b_face_cells[face_id];
treloc[ielt] = coefap[face_id] + coefbp[face_id] * f->val[c_id];
}
Note
Here, we assign the non-reconstructed value.

Open file to print values and broadcast values to all parallel ranks. Values are ordered by the xabs array values provided to the cs_parall_allgather_ordered_r, so this function is needed even when not running in parallel.

if (cs_glob_rank_id < 1) {
file = fopen("Nusselt.dat", "w");
}
cs_gnum_t neltg = nlelt;
cs_parall_counter(&neltg, 1);
cs_real_t *xabs = NULL, *xabsg = NULL, *xnusselt = NULL;
cs_real_t *treglo = NULL;
BFT_MALLOC(xabs, nlelt, cs_real_t);
BFT_MALLOC(xnusselt, neltg, cs_real_t);
BFT_MALLOC(xabsg, neltg, cs_real_t);
BFT_MALLOC(treglo, neltg, cs_real_t);
for (cs_lnum_t ilelt = 0; ilelt < nlelt; ilelt ++) {
cs_lnum_t face_id = lstelt[ilelt];
xabs[ilelt] = cdgfbo[face_id][0];
}
cs_parall_allgather_ordered_r(nlelt, neltg, 1, xabs, xabs, xabsg);
cs_parall_allgather_ordered_r(nlelt, neltg, 1, xabs, treloc, treglo);
BFT_FREE(xabs);
BFT_FREE(treloc);
BFT_FREE(lstelt);

Compute the bulk temperature, finalize the Nusselt number, print it and free memory not already freed before:

for (cs_gnum_t ilelt = 0; ilelt < neltg; ilelt ++) {
cs_real_t xtbulk = 0;
cs_real_t xubulk = 0;
cs_real_t lambda = 0;
int npoint = 200;
cs_lnum_t c_id_prev = -1;
int irang1 = -1;
int irangv;
for (int ii = 0; ii < npoint; ii++) {
cs_lnum_t c_id;
cs_real_t xyz[3] = {xabsg[ilelt],
(cs_real_t)ii/(cs_real_t)(npoint-1),
0.};
cell_cen,
xyz,
&c_id, &irangv);
if (c_id != c_id_prev || irangv != irang1) {
c_id_prev = c_id;
irang1 = irangv;
if (cs_glob_rank_id == irangv) {
cs_real_t xtb = volume[c_id]*f->val[c_id]*vel[c_id][0];
cs_real_t xub = volume[c_id]*vel[c_id][0];
xtbulk += xtb;
xubulk += xub;
lambda = phys_pro->cp0*mu->val[c_id] / prandtl;
}
}
}
/* Note: only last location used for lambda.
Would be more appropriate (but still not general at all)
with visls0 in definition of lambda and from the derivation
of Nusselt number; use of lambda at the wall (ii = 0)
would seem better */
cs_parall_bcast(irangv, 1, CS_REAL_TYPE, &lambda);
cs_real_t xbulk[2] = {xtbulk, xubulk};
xtbulk = xbulk[0] / xbulk[1];
xubulk = xbulk[1];
cs_real_t tfac = treglo[ilelt];
xnusselt[ilelt] = qwall * 2. * height / lambda / (tfac - xtbulk);
}
for (cs_gnum_t ii = 0; ii < neltg; ii++)
fprintf(file,
"%17.9e %17.9e\n",
xabsg[ii]*10.,
xnusselt[ii]/(0.023*pow(30000., 0.8)*pow(0.71, 0.4)));
if (file != NULL)
fclose(file);
BFT_FREE(xnusselt);
BFT_FREE(xabsg);
BFT_FREE(treglo);