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;
int nt_cur = domain->time_step->nt_cur;
double cs_real_t
Floating-point value.
Definition: cs_defs.h:342
cs_rreal_t cs_rreal_3_t[3]
Definition: cs_defs.h:388
cs_lnum_t cs_lnum_2_t[2]
vector of 2 local mesh-entity ids
Definition: cs_defs.h:352
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:335
double precision, dimension(:,:), pointer diipb
vector II' for interior faces for every boundary face, the three components of the vector ....
Definition: mesh.f90:160
Definition: cs_mesh_quantities.h:92
cs_rreal_3_t * diipb
Definition: cs_mesh_quantities.h:132
cs_real_t * b_face_surf
Definition: cs_mesh_quantities.h:119
cs_real_t * cell_vol
Definition: cs_mesh_quantities.h:97
cs_lnum_t * b_face_cells
Definition: cs_mesh.h:112
cs_lnum_t n_i_faces
Definition: cs_mesh.h:98
cs_lnum_t n_b_faces
Definition: cs_mesh.h:99
cs_lnum_t n_cells_with_ghosts
Definition: cs_mesh.h:151
cs_lnum_t n_cells
Definition: cs_mesh.h:97
cs_lnum_2_t * i_face_cells
Definition: cs_mesh.h:111
Get the physical fields
cs_lnum_t cell_id, cell_id1, cell_id2, face_id;
int nt_cur = domain->time_step->nt_cur;
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];
}
}
#define BFT_FREE(_ptr)
Definition: bft_mem.h:90
#define BFT_MALLOC(_ptr, _ni, _type)
Definition: bft_mem.h:58
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:359
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.cpp:3275
cs_field_t * cs_field_by_id(int id)
Return a pointer to a field based on its id.
Definition: cs_field.cpp:2465
int cs_field_key_id(const char *name)
Return an id associated with a given key name.
Definition: cs_field.cpp:2781
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.cpp:470
@ h
Definition: cs_field_pointer.h:93
cs_real_t * val
Definition: cs_field.h:152
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;
@ CS_DOUBLE
Definition: cs_defs.h:305
@ rho
Definition: cs_field_pointer.h:99
@ dt
Definition: cs_field_pointer.h:65
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:147
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.cpp:213
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;