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 cell_id, cell_id1, cell_id2, face_id;
Get the physical fields
cs_lnum_t cell_id, cell_id1, cell_id2, face_id;
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.;
return;
if (false) {
true,
1,
grad);
for (face_id = 0; face_id < n_b_faces; face_id++) {
cell_id = b_face_cells[face_id];
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];
}
} else {
for (face_id = 0; face_id < n_b_faces; face_id++) {
cell_id = b_face_cells[face_id];
h_reconstructed[face_id] =
h->val[cell_id];
}
}
Computation step
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]);
}
for (face_id = 0; face_id < n_i_faces; face_id++) {
cell_id1 = i_face_cells[face_id][0];
cell_id2 = i_face_cells[face_id][1];
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];
div_balance += b_mass_flux[face_id] *
dt[cell_id] *
h->val[cell_id];
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
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
};
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];
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
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]);
}
for (face_id = 0; face_id < n_i_faces; face_id++) {
cell_id1 = i_face_cells[face_id][0];
cell_id2 = i_face_cells[face_id][1];
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];
div_balance += b_mass_flux[face_id] *
dt[cell_id] *
h->val[cell_id];
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
face_id = face_list[i];
cell_id = b_face_cells[face_id];
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]);
}
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
};
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];
tot_balance = vol_balance + div_balance + a_wall_balance + h_wall_balance
+ sym_balance + in_balance + out_balance + mass_i_balance
+ mass_o_balance;