1 /*============================================================================
2 * This function is called at the end of each time step, and has a very
3 * general purpose
4 * (i.e. anything that does not have another dedicated user function)
5 *============================================================================*/
6
7 /* VERS */
8
9 /*
10 This file is part of Code_Saturne, a general-purpose CFD tool.
11
12 Copyright (C) 1998-2021 EDF S.A.
13
14 This program is free software; you can redistribute it and/or modify it under
15 the terms of the GNU General Public License as published by the Free Software
16 Foundation; either version 2 of the License, or (at your option) any later
17 version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
22 details.
23
24 You should have received a copy of the GNU General Public License along with
25 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
26 Street, Fifth Floor, Boston, MA 02110-1301, USA.
27 */
28
29 /*----------------------------------------------------------------------------*/
30
31 #include "cs_defs.h"
32
33 /*----------------------------------------------------------------------------
34 * Standard C library headers
35 *----------------------------------------------------------------------------*/
36
37 #include <assert.h>
38 #include <math.h>
39 #include <stdio.h>
40
41 #if defined(HAVE_MPI)
42 #include <mpi.h>
43 #endif
44
45 /*----------------------------------------------------------------------------
46 * PLE library headers
47 *----------------------------------------------------------------------------*/
48
49 #include <ple_coupling.h>
50
51 /*----------------------------------------------------------------------------
52 * Local headers
53 *----------------------------------------------------------------------------*/
54
55 #include "cs_headers.h"
56
57 /*----------------------------------------------------------------------------*/
58
59 BEGIN_C_DECLS
60
61 /*----------------------------------------------------------------------------*/
62 /*!
63 * \file cs_user_extra_operations-scalar_balance.c
64 *
65 * \brief This function is called at the end of each time step, and has a very
66 * general purpose (i.e. anything that does not have another dedicated
67 * user function).
68 *
69 * \param[in, out] domain pointer to a cs_domain_t structure
70 */
71 /*----------------------------------------------------------------------------*/
72
73 /*============================================================================
74 * User function definitions
75 *============================================================================*/
76
77 /*----------------------------------------------------------------------------
78 * Example for scalar balance.
79 *----------------------------------------------------------------------------*/
80
81 void
cs_user_extra_operations(cs_domain_t * domain)82 cs_user_extra_operations(cs_domain_t *domain)
83 {
84
85 /* Local variables */
86
87 /*! [local_variables] */
88
89 cs_lnum_t n_faces;
90 cs_lnum_t *face_list;
91
92 cs_lnum_t cell_id, cell_id1, cell_id2, face_id;
93 int nt_cur = domain->time_step->nt_cur;
94
95 const cs_mesh_t *m = domain->mesh;
96 const cs_mesh_quantities_t *mq = domain->mesh_quantities;
97
98 const cs_lnum_t n_cells = m->n_cells;
99 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
100 const cs_lnum_t n_i_faces = m->n_i_faces;
101 const cs_lnum_t n_b_faces = m->n_b_faces;
102
103 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)m->i_face_cells;
104 const cs_lnum_t *b_face_cells = (const cs_lnum_t *)m->b_face_cells;
105
106 const cs_real_t *cell_vol = mq->cell_vol;
107 const cs_real_3_t *diipb = (const cs_real_3_t *)mq->diipb;
108 const cs_real_t *b_face_surf = (const cs_real_t *)mq->b_face_surf;
109
110 /*! [local_variables] */
111
112 /*! [physical_fields] */
113
114 /* Get physical fields */
115 const cs_real_t *dt = CS_F_(dt)->val;
116 const cs_real_t *rho = CS_F_(rho)->val;
117 const cs_field_t *h = cs_field_by_name_try("enthalpy");
118
119 /*! [physical_fields] */
120
121 /*-------------------------------------------------------------------------
122 * This example computes energy balance relative to enthalpy
123 * We assume that we want to compute balances (convective and diffusive)
124 * at the boundaries of the calculation domain represented below
125 * (with boundaries marked by colors).
126 *
127 * The scalar considered if the enthalpy. We will also use the
128 * specific heat (to obtain balances in Joules)
129 *
130 *
131 * Domain and associated boundary colors:
132 * - 2, 3, 4, 7 : adiabatic walls
133 * - 6 : wall with fixed enthalpy
134 * - 1 : inlet
135 * - 5 : outlet
136 * - 8, 9 : symmetry
137 *-------------------------------------------------------------------------*/
138
139 /* 1. Initialization
140 =================
141
142 --> Local variables
143 ---------------
144
145 vol_balance : volume contribution of unsteady terms
146 div_balance : volume contribution due to to term in div(rho u)
147 a_wall_balance: contribution from adiabatic walls
148 h_wall_balance: contribution from walls with fixed temperature
149 sym_balance : contribution from symmetry boundaries
150 in_balance : contribution from inlets
151 out_balance : contribution from outlets
152 mass_i_balance: contribution from mass injections
153 mass_o_balance: constribution from mass suctions
154 tot_balance : total balance */
155
156 /*! [init] */
157
158 double vol_balance = 0.;
159 double div_balance = 0.;
160 double a_wall_balance = 0.;
161 double h_wall_balance = 0.;
162 double sym_balance = 0.;
163 double in_balance = 0.;
164 double out_balance = 0.;
165 double mass_i_balance = 0.;
166 double mass_o_balance = 0.;
167 double tot_balance = 0.;
168
169 /* If the scalar enthalpy is not computed, return */
170 if (h == NULL)
171 return;
172
173 /* Boundary condition coefficient for h */
174 const cs_real_t *a_H = h->bc_coeffs->a;
175 const cs_real_t *b_H = h->bc_coeffs->b;
176 const cs_real_t *af_H = h->bc_coeffs->af;
177 const cs_real_t *bf_H = h->bc_coeffs->bf;
178
179 /* Convective mass fluxes for inner and boundary faces */
180 int iflmas = cs_field_get_key_int(h, cs_field_key_id("inner_mass_flux_id"));
181 const cs_real_t *i_mass_flux = cs_field_by_id(iflmas)->val;
182
183 int iflmab = cs_field_get_key_int(h, cs_field_key_id("boundary_mass_flux_id"));
184 const cs_real_t *b_mass_flux = cs_field_by_id(iflmab)->val;
185
186 /* Allocate temporary array */
187 cs_real_t *h_reconstructed;
188 BFT_MALLOC(h_reconstructed, n_b_faces, cs_real_t);
189
190 /* Reconstructed value */
191 if (false) {
192 cs_real_3_t *grad;
193 BFT_MALLOC(grad, n_cells_ext, cs_real_3_t);
194
195 cs_field_gradient_scalar(h,
196 true, /* use_previous_t */
197 1, /* inc */
198 true, /* _recompute_cocg */
199 grad);
200
201 for (face_id = 0; face_id < n_b_faces; face_id++) {
202 cell_id = b_face_cells[face_id]; // associated boundary cell
203 h_reconstructed[face_id] = h->val[cell_id]
204 + grad[cell_id][0]*diipb[face_id][0]
205 + grad[cell_id][1]*diipb[face_id][1]
206 + grad[cell_id][2]*diipb[face_id][2];
207 }
208
209 BFT_FREE(grad);
210
211 /* Non-reconstructed value */
212 } else {
213 for (face_id = 0; face_id < n_b_faces; face_id++) {
214 cell_id = b_face_cells[face_id]; // associated boundary cell
215 h_reconstructed[face_id] = h->val[cell_id];
216 }
217 }
218
219 /*! [init] */
220
221 /* 2. Compute the balance at time step n
222 ======================================
223 */
224
225 /*! [computation] */
226
227 /*
228 --> Balance on interior volumes
229 ---------------------------
230 */
231
232 for (cell_id = 0; cell_id < n_cells; cell_id++) {
233 vol_balance += cell_vol[cell_id] * rho[cell_id]
234 * (h->val_pre[cell_id] - h->val[cell_id]);
235 }
236
237 /*
238 --> Balance on all faces (interior and boundary), for div(rho u)
239 ------------------------------------------------------------
240 */
241
242 for (face_id = 0; face_id < n_i_faces; face_id++) {
243
244 cell_id1 = i_face_cells[face_id][0]; // associated boundary cell
245 cell_id2 = i_face_cells[face_id][1]; // associated boundary cell
246
247 /* Contribution to flux from the two cells of the current face
248 (The cell is count only once in parallel by checking that
249 the cell_id is not in the halo) */
250
251 if (cell_id1 < n_cells)
252 div_balance += i_mass_flux[face_id] * dt[cell_id1] * h->val[cell_id1];
253
254 if (cell_id2 < n_cells)
255 div_balance -= i_mass_flux[face_id] * dt[cell_id2] * h->val[cell_id2];
256 }
257
258 for (face_id = 0; face_id < n_b_faces; face_id++) {
259
260 cell_id = b_face_cells[face_id]; // associated boundary cell
261
262 /* Contribution to flux from the current face */
263 div_balance += b_mass_flux[face_id] * dt[cell_id] * h->val[cell_id];
264
265 }
266
267 // TODO mass source terms and mass accumulation term
268 // In case of a mass source term, add contribution from Gamma*Tn+1
269
270 /*
271 --> Balance on boundary faces
272 -------------------------
273
274 We handle different types of boundary faces separately to better
275 analyze the information, but this is not mandatory. */
276
277 /*
278 Compute the contribution from walls with colors 2, 3, 4 and 7
279 (adiabatic here, so flux should be 0)
280 */
281 BFT_MALLOC(face_list, n_b_faces, cs_lnum_t);
282
283 cs_selector_get_b_face_list("2 or 3 or 4 or 7", &n_faces, face_list);
284
285 for (cs_lnum_t i = 0; i < n_faces; i++) {
286
287 face_id = face_list[i];
288 cell_id = b_face_cells[face_id]; // associated boundary cell
289
290 /* Contribution to flux from the current face
291 (diffusion and convection flux, negative if incoming) */
292
293 a_wall_balance += - b_face_surf[face_id] * dt[cell_id]
294 * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
295 - b_mass_flux[face_id] * dt[cell_id]
296 * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
297
298 }
299
300 /*
301 Contribution from walls with color 6
302 (here at fixed enthalpy; the convective flux should be 0)
303 */
304 cs_selector_get_b_face_list("6", &n_faces, face_list);
305
306 for (cs_lnum_t i = 0; i < n_faces; i++) {
307
308 face_id = face_list[i];
309 cell_id = b_face_cells[face_id]; // associated boundary cell
310
311 /* Contribution to flux from the current face
312 (diffusion and convection flux, negative if incoming) */
313
314 h_wall_balance += - b_face_surf[face_id] * dt[cell_id]
315 * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
316 - b_mass_flux[face_id] * dt[cell_id]
317 * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
318
319 }
320
321 /*
322 Contribution from symmetries (should be 0).
323 */
324 cs_selector_get_b_face_list("8 or 9", &n_faces, face_list);
325
326 for (cs_lnum_t i = 0; i < n_faces; i++) {
327
328 face_id = face_list[i];
329 cell_id = b_face_cells[face_id]; // associated boundary cell
330
331 /* Contribution to flux from the current face
332 (diffusion and convection flux, negative if incoming) */
333
334 sym_balance += - b_face_surf[face_id] * dt[cell_id]
335 * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
336 - b_mass_flux[face_id] * dt[cell_id]
337 * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
338
339 }
340
341 /*
342 Contribution from inlet (color 1, diffusion and convection flux)
343 */
344 cs_selector_get_b_face_list("1", &n_faces, face_list);
345
346 for (cs_lnum_t i = 0; i < n_faces; i++) {
347
348 face_id = face_list[i];
349 cell_id = b_face_cells[face_id]; // associated boundary cell
350
351 /* Contribution to flux from the current face
352 (diffusion and convection flux, negative if incoming) */
353
354 in_balance += - b_face_surf[face_id] * dt[cell_id]
355 * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
356 - b_mass_flux[face_id] * dt[cell_id]
357 * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
358
359 }
360
361 /*
362 Contribution from outlet (color 5, diffusion and convection flux)
363 */
364 cs_selector_get_b_face_list("5", &n_faces, face_list);
365
366 for (cs_lnum_t i = 0; i < n_faces; i++) {
367
368 face_id = face_list[i];
369 cell_id = b_face_cells[face_id]; // associated boundary cell
370
371 /* Contribution to flux from the current face
372 (diffusion and convection flux, negative if incoming) */
373
374 out_balance += - b_face_surf[face_id] * dt[cell_id]
375 * (af_H[face_id] + bf_H[face_id] * h_reconstructed[face_id])
376 - b_mass_flux[face_id] * dt[cell_id]
377 * (a_H[face_id] + b_H[face_id] * h_reconstructed[face_id]);
378
379 }
380
381 /* Free memory */
382 BFT_FREE(face_list);
383 BFT_FREE(h_reconstructed);
384
385 /* Sum of values on all ranks (parallel calculations) */
386
387 cs_parall_sum(1, CS_DOUBLE, &vol_balance);
388 cs_parall_sum(1, CS_DOUBLE, &div_balance);
389 cs_parall_sum(1, CS_DOUBLE, &a_wall_balance);
390 cs_parall_sum(1, CS_DOUBLE, &h_wall_balance);
391 cs_parall_sum(1, CS_DOUBLE, &sym_balance);
392 cs_parall_sum(1, CS_DOUBLE, &in_balance);
393 cs_parall_sum(1, CS_DOUBLE, &out_balance);
394 cs_parall_sum(1, CS_DOUBLE, &mass_i_balance);
395 cs_parall_sum(1, CS_DOUBLE, &mass_o_balance);
396
397 /* --> Total balance
398 ------------- */
399
400 /* We add the different contributions calculated above */
401
402 tot_balance = vol_balance + div_balance + a_wall_balance + h_wall_balance
403 + sym_balance + in_balance + out_balance + mass_i_balance
404 + mass_o_balance;
405
406 /*! [computation] */
407
408
409 /* 3. Write the balance at time step n
410 ==================================== */
411
412 /*! [writing] */
413
414 bft_printf("\n ** Enthalpy balance **\n"
415 " ----------------\n"
416 "-----------"
417 "----------------------------------------------------------\n"
418 "bt Iter"
419 " Volume Divergence Adia Wall Fixed_H Wall Symmetry"
420 " Inlet Outlet Inj. Mass. Suc. Mass. Total\n"
421 "bt %6i %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e "
422 "%12.4e %12.4e %12.4e\n"
423 "-----------"
424 "----------------------------------------------------------\n",
425 nt_cur, vol_balance, div_balance, a_wall_balance, h_wall_balance,
426 sym_balance, in_balance, out_balance,
427 mass_i_balance, mass_o_balance, tot_balance);
428
429 /*! [writing] */
430
431 }
432
433 END_C_DECLS
434