1 /*============================================================================
2 * Modelling the thermal wall with a 1D approach.
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <math.h>
34 #include <stdlib.h>
35 #include <stdio.h>
36 #include <string.h>
37
38 /*----------------------------------------------------------------------------
39 * Local headers
40 *----------------------------------------------------------------------------*/
41
42 #include "bft_mem.h"
43 #include "bft_error.h"
44 #include "bft_printf.h"
45
46 #include "cs_base.h"
47 #include "cs_field.h"
48 #include "cs_field_pointer.h"
49 #include "cs_lagr.h"
50 #include "cs_mesh.h"
51 #include "cs_mesh_location.h"
52 #include "cs_parall.h"
53 #include "cs_physical_constants.h"
54 #include "cs_restart.h"
55 #include "cs_restart_default.h"
56
57 /*----------------------------------------------------------------------------
58 * Header for the current file
59 *----------------------------------------------------------------------------*/
60
61 #include "cs_1d_wall_thermal.h"
62
63 /*----------------------------------------------------------------------------*/
64
65 BEGIN_C_DECLS
66
67 /*=============================================================================
68 * Additional doxygen documentation
69 *============================================================================*/
70
71 /*!
72 \file cs_1d_wall_thermal.c
73 Modelling the thermal wall with a 1D approach.
74
75 \struct cs_1d_wall_thermal_local_model_t
76 1D wall thermal module parameters and variables
77 for each coupled face.
78
79 \var cs_1d_wall_thermal_local_model_t::nppt1d
80 Number of discretisation cells in the 1D wall
81 for the nfpt1d boundary faces which are coupled
82 with a 1D wall thermal module.
83 \var cs_1d_wall_thermal_local_model_t::iclt1d
84 Boundary condition type at the external (pseudo)
85 wall:
86 - 1: Dirichlet,
87 - 2: Flux condition.
88 \var cs_1d_wall_thermal_local_model_t::eppt1d
89 Thickness of the 1D wall for the nfpt1d boundary faces
90 which are coupled with a 1D wall thermal module.
91 \var cs_1d_wall_thermal_local_model_t::rgpt1d
92 Geometry of the pseudo wall mesh (refined as a fluid
93 if rgt1d is smaller than 1.
94 \var cs_1d_wall_thermal_local_model_t::tept1d
95 External temperature of the pseudo wall in the Dirichlet case.
96 \var cs_1d_wall_thermal_local_model_t::hept1d
97 External coefficient of transfer in the pseudo wall
98 under Dirichlet conditions, (in \f$W.m^{-2}.K\f$).
99 \var cs_1d_wall_thermal_local_model_t::fept1d
100 External heat flux in the pseudo wall under the flux
101 conditions (in \f$W.m^{-2}\f$, negative value for energy
102 entering the wall).
103 \var cs_1d_wall_thermal_local_model_t::xlmbt1
104 Thermal diffusivity.
105 \var cs_1d_wall_thermal_local_model_t::rcpt1d
106 Volumetric heat capacity rho C_p of the wall uniform
107 in thickness (\f$J.m^{-3}.K^{-1}\f$).
108 \var cs_1d_wall_thermal_local_model_t::dtpt1d
109 Physical time step associated with the solved 1D equation
110 of the pseudo wall (which can be different from the time
111 step of the calculation).
112 \var cs_1d_wall_thermal_local_model_t::z
113 Discretization points coordinates.
114 \var cs_1d_wall_thermal_local_model_t::t
115 Temperature at each point of discretization.
116
117 \struct cs_1d_wall_thermal_t
118
119 \brief 1D wall thermal module descriptor.
120
121 \var cs_1d_wall_thermal_t::nfpt1d
122 Number of boundary faces which are coupled
123 with a 1D wall thermal module
124 Zones of t1d, dimensioned with nfabor
125 Global number of boundary faces which are coupled with
126 a 1D wall thermal module, i.e. sum of nfpt1d over all
127 ranks
128 \var cs_1d_wall_thermal_t::nmxt1d
129 \var cs_1d_wall_thermal_t::izft1d
130 Zones of t1d, dimensioned with nfabor
131 \var cs_1d_wall_thermal_t::ifpt1d
132 Array allowing to mark out the numbers of
133 the nfpt1d boundary faces which are coupled with
134 a 1D wall. The test on \ref ifpt1d implicitly assumes
135 that the array is completed in ascending order
136 (i.e ifpt1d[ii]\f$>\f$ifpt1d[jj] if ii\f$>\f$jj.
137 This will be the case if the coupled faces are defined
138 starting from the unique loop on the boundary faces (as
139 in the example, see \ref us_pt1d). If this is not
140 the case, contact the development team to short circuit
141 the test.
142 \var cs_1d_wall_thermal_t::tppt1d
143 Initialization temperature of the wall (uniform in thickness).
144 During the calculation, the array stores the temperature
145 of the solid at the fluid/solid interface.
146 Other than for re-reading a file, \ref tppt1d is not used.
147 \ref cs_1d_wall_thermal_local_model_t::nppt1d "nppt1d" ,
148 \ref ifpt1d, \ref cs_1d_wall_thermal_local_model_t::rgpt1d "rgpt1d"
149 and \ref cs_1d_wall_thermal_local_model_t::eppt1d "eppt1d" are
150 compared to data from the follow-up file and they must be identical.
151 \var cs_1d_wall_thermal_t::local_models
152 Array of structures containing the 1D wall thermal
153 module parameters and variables for each coupled face.
154 */
155
156 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
157
158 /*============================================================================
159 * Local structure definitions
160 *============================================================================*/
161
162 /*============================================================================
163 * Static global variable
164 *============================================================================*/
165
166 static cs_1d_wall_thermal_t _1d_wall_thermal =
167 {
168 .nfpt1d = 0,
169 .nfpt1t = 0,
170 .nmxt1d = 0,
171 .izft1d = NULL,
172 .ifpt1d = NULL,
173 .tppt1d = NULL,
174 .local_models = NULL
175 };
176
177 const cs_1d_wall_thermal_t *cs_glob_1d_wall_thermal = &_1d_wall_thermal;
178
179 static cs_restart_t *cs_glob_tpar1d_suite = NULL;
180
181 /*============================================================================
182 * Prototypes for functions intended for use only by Fortran wrappers.
183 * (descriptions follow, with function bodies).
184 *============================================================================*/
185
186 void
187 cs_f_1d_wall_thermal_get_pointers(cs_lnum_t **nfpt1d,
188 cs_gnum_t **nfpt1t);
189
190 void
191 cs_f_1d_wall_thermal_get_faces(cs_lnum_t **ifpt1d);
192
193 void
194 cs_f_1d_wall_thermal_get_temp(cs_real_t **tppt1d);
195
196 /*============================================================================
197 * Private function definitions
198 *============================================================================*/
199
200 /*----------------------------------------------------------------------------*/
201 /*!
202 * \brief Allocate the discretization points coordinates array and
203 the temperature at each point of discretization.
204 */
205 /*----------------------------------------------------------------------------*/
206
207 static void
_1d_wall_thermal_local_models_init(void)208 _1d_wall_thermal_local_models_init(void)
209 {
210 cs_lnum_t ii;
211
212 /* Computation of nmxt1d */
213 for (ii = 0; ii < _1d_wall_thermal.nfpt1d ; ii++) {
214 _1d_wall_thermal.nmxt1d = CS_MAX(_1d_wall_thermal.local_models[ii].nppt1d,
215 _1d_wall_thermal.nmxt1d);
216 }
217
218 /* if necessary, sum over all the processors */
219 cs_parall_max(1, CS_INT_TYPE, &_1d_wall_thermal.nmxt1d);
220
221 /* Initialization of the number of discretization points in each structure
222 Computation of the total number of discretization points */
223 cs_lnum_t nb_pts_tot = 0;
224
225 for (ii = 0; ii < _1d_wall_thermal.nfpt1d ; ii++)
226 nb_pts_tot += _1d_wall_thermal.local_models[ii].nppt1d;
227
228 /* Allocate the "t" arrays: Temperature in each point of discretization
229 and the "z" arrays: Coordonnates of each point of discretization */
230
231 if (_1d_wall_thermal.nfpt1d > 0) {
232 BFT_MALLOC(_1d_wall_thermal.local_models->z, 2 * nb_pts_tot, cs_real_t);
233 _1d_wall_thermal.local_models->t = _1d_wall_thermal.local_models->z
234 + nb_pts_tot;
235 }
236
237 for (ii = 1 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
238 _1d_wall_thermal.local_models[ii].z
239 = _1d_wall_thermal.local_models[ii-1].z
240 + _1d_wall_thermal.local_models[ii-1].nppt1d;
241 _1d_wall_thermal.local_models[ii].t
242 = _1d_wall_thermal.local_models[ii-1].t
243 + _1d_wall_thermal.local_models[ii-1].nppt1d;
244 }
245 }
246
247 /*============================================================================
248 * Fortran wrapper function definitions
249 *============================================================================*/
250
251 /*----------------------------------------------------------------------------*/
252 /*! \brief Get pointers to members of the global 1d wall thermal structure.
253 *
254 * This function is intended for use by Fortran wrappers, and
255 * enables mapping to Fortran global pointers.
256 *
257 * parameters:
258 * \param[out] nfpt1d pointer to cs_glob_1d_wall_thermal->nfpt1d
259 * \param[out] nfpt1t pointer to cs_glob_1d_wall_thermal->nfpt1t
260 */
261 /*----------------------------------------------------------------------------*/
262
263 void
cs_f_1d_wall_thermal_get_pointers(cs_lnum_t ** nfpt1d,cs_gnum_t ** nfpt1t)264 cs_f_1d_wall_thermal_get_pointers(cs_lnum_t **nfpt1d,
265 cs_gnum_t **nfpt1t)
266 {
267 *nfpt1d = &(_1d_wall_thermal.nfpt1d);
268 *nfpt1t = &(_1d_wall_thermal.nfpt1t);
269 }
270
271 /*----------------------------------------------------------------------------*/
272 /*!
273 * \brief Return the ifpt1d array for the 1D wall thermal module.
274 */
275 /*----------------------------------------------------------------------------*/
276
277 void
cs_f_1d_wall_thermal_get_faces(cs_lnum_t ** ifpt1d)278 cs_f_1d_wall_thermal_get_faces(cs_lnum_t **ifpt1d)
279 {
280 *ifpt1d = _1d_wall_thermal.ifpt1d;
281 }
282
283 /*----------------------------------------------------------------------------*/
284 /*!
285 * \brief Return the tppt1d array for the 1D wall thermal module.
286 */
287 /*----------------------------------------------------------------------------*/
288
289 void
cs_f_1d_wall_thermal_get_temp(cs_real_t ** tppt1d)290 cs_f_1d_wall_thermal_get_temp(cs_real_t **tppt1d)
291 {
292 *tppt1d = _1d_wall_thermal.tppt1d;
293 }
294
295 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
296
297 /*============================================================================
298 * Public function definitions
299 *============================================================================*/
300
301 /*----------------------------------------------------------------------------*/
302 /*!
303 * \brief Initialize the cs_glob_1d_wall_thermal structure.
304 */
305 /*----------------------------------------------------------------------------*/
306
307 void
cs_1d_wall_thermal_create(void)308 cs_1d_wall_thermal_create(void)
309 {
310 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
311
312 _1d_wall_thermal.nfpt1d = 0;
313 _1d_wall_thermal.nfpt1t = 0;
314 _1d_wall_thermal.nmxt1d = 0;
315
316 /* Allocate the izft1d array */
317 BFT_MALLOC(_1d_wall_thermal.izft1d, n_b_faces, cs_lnum_t);
318
319 for (cs_lnum_t ifac = 0 ; ifac < n_b_faces ; ifac++) {
320 _1d_wall_thermal.izft1d[ifac] = 0;
321 }
322 }
323
324 /*----------------------------------------------------------------------------*/
325 /*!
326 * \brief Allocate the array of structures local_models.
327 */
328 /*----------------------------------------------------------------------------*/
329
330 void
cs_1d_wall_thermal_local_models_create(void)331 cs_1d_wall_thermal_local_models_create(void)
332 {
333 cs_lnum_t ii;
334
335 /* Allocate the ifpt1d array */
336 BFT_MALLOC(_1d_wall_thermal.ifpt1d, _1d_wall_thermal.nfpt1d, cs_lnum_t);
337
338 /* Allocate the tppt1d array */
339 BFT_MALLOC(_1d_wall_thermal.tppt1d, _1d_wall_thermal.nfpt1d, cs_real_t);
340
341 /* Allocate the local_models member of the cs_glob_1d_wall_thermal
342 * structure */
343 BFT_MALLOC(_1d_wall_thermal.local_models,
344 _1d_wall_thermal.nfpt1d,
345 cs_1d_wall_thermal_local_model_t);
346
347 for (ii = 0; ii < _1d_wall_thermal.nfpt1d ; ii++) {
348 _1d_wall_thermal.local_models[ii].nppt1d = -999;
349 _1d_wall_thermal.local_models[ii].iclt1d = 3;
350 _1d_wall_thermal.ifpt1d[ii] = -999;
351 _1d_wall_thermal.local_models[ii].eppt1d = -999.;
352 _1d_wall_thermal.local_models[ii].rgpt1d = -999.;
353 _1d_wall_thermal.tppt1d[ii] = 0.;
354 _1d_wall_thermal.local_models[ii].tept1d = 0.;
355 _1d_wall_thermal.local_models[ii].hept1d = 1.e30;
356 _1d_wall_thermal.local_models[ii].fept1d = 0.;
357 _1d_wall_thermal.local_models[ii].xlmbt1 = -999.;
358 _1d_wall_thermal.local_models[ii].rcpt1d = -999.;
359 _1d_wall_thermal.local_models[ii].dtpt1d = -999.;
360 }
361 }
362
363 /*----------------------------------------------------------------------------*/
364 /*!
365 * \brief Create the 1D mesh for each face and initialize the temperature.
366 */
367 /*----------------------------------------------------------------------------*/
368
369 void
cs_1d_wall_thermal_mesh_create(void)370 cs_1d_wall_thermal_mesh_create(void)
371 {
372 cs_lnum_t ii, kk;
373 cs_real_t m, rr, e, n;
374 cs_real_t *zz;
375
376 /* Allocate the global structure: cs_glob_par1d and the number of
377 discretization points on each face */
378 if (_1d_wall_thermal.nfpt1t > 0)
379 _1d_wall_thermal_local_models_init();
380
381 for (ii = 0; ii < _1d_wall_thermal.nfpt1d; ii++) {
382
383 n = _1d_wall_thermal.local_models[ii].nppt1d;
384 e = _1d_wall_thermal.local_models[ii].eppt1d;
385
386 /* Initialization of the Temperature */
387 for (kk = 0; kk < n; kk++) {
388 (_1d_wall_thermal.local_models[ii].t)[kk] = _1d_wall_thermal.tppt1d[ii];
389 }
390
391 /* Mesh */
392 zz = _1d_wall_thermal.local_models[ii].z;
393 rr = _1d_wall_thermal.local_models[ii].rgpt1d;
394
395 /* Regular */
396 if (fabs(rr-1.0) <= 1.e-6) {
397 zz[0] = e/n/2.;
398 for (kk = 1 ; kk < n ; kk++) {
399 zz[kk] = zz[kk-1] + e/n;
400 }
401 }
402
403 /* Geometric */
404 else {
405 m = e*(1.-rr)/(1.-pow(rr,n));
406 *zz = m/2.;
407 for (kk = 1 ; kk < n ; kk++) {
408 zz[kk] = zz[kk-1]+m/2.;
409 m = m*rr;
410 zz[kk] = zz[kk]+m/2.;
411 }
412 }
413 }
414 }
415
416 /*----------------------------------------------------------------------------*/
417 /*!
418 * \brief Solve the 1D equation for a given face
419 *
420 * \param[in] ii face number
421 * \param[in] tf fluid temperature at the boundarys
422 * \param[in] hf exchange coefficient for the fluid
423 */
424 /*----------------------------------------------------------------------------*/
425
426 void
cs_1d_wall_thermal_solve(cs_lnum_t ii,cs_real_t tf,cs_real_t hf)427 cs_1d_wall_thermal_solve(cs_lnum_t ii,
428 cs_real_t tf,
429 cs_real_t hf)
430 {
431 cs_lnum_t kk;
432 cs_real_t qinc, eps;
433 cs_lnum_t ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
434
435 if (cs_glob_lagr_extra_module->radiative_model >= 1) {
436 /* coupling with radiative module, qinc and qeps != 0 */
437 /* incident radiative flux at the boundary at the boundary */
438 qinc = CS_F_(qinci)->val[ifac];
439
440 /* emissivity */
441 eps = CS_F_(emissivity)->val[ifac];
442 } else {
443 /* without coupling with radiative module */
444 qinc = 0.;
445 eps = 0.;
446 }
447
448 /* thermal diffusivity */
449 cs_real_t xlmbt1 = _1d_wall_thermal.local_models[ii].xlmbt1;
450
451 /* volumetric heat capacity of the wall */
452 cs_real_t rcp = _1d_wall_thermal.local_models[ii].rcpt1d;
453
454 /* exchange coefficient on the exterior wall */
455 cs_real_t hept1d = _1d_wall_thermal.local_models[ii].hept1d;
456
457 /* flux on the exterior wall */
458 cs_real_t fept1d = _1d_wall_thermal.local_models[ii].fept1d;
459
460 /* thickness of the 1D wall */
461 cs_real_t eppt1d = _1d_wall_thermal.local_models[ii].eppt1d;
462
463 /* temperature on the exterior boundary */
464 cs_real_t tept1d = _1d_wall_thermal.local_models[ii].tept1d;
465
466 /* time-step for the solid resolution */
467 cs_real_t dtpt1d = _1d_wall_thermal.local_models[ii].dtpt1d;
468
469 /* type of exterior boundary condition */
470 int iclt1d = _1d_wall_thermal.local_models[ii].iclt1d;
471
472 cs_real_t a1; /* extrapolation coefficient for temperature1 */
473 cs_real_t h2; /* thermal exchange coefficient on T(1) */
474 cs_real_t f3; /* thermal flux on Tfluide */
475 cs_real_t a4; /* extrapolation coefficient for temperature4 */
476
477 cs_real_t h5 = 0.; /* thermal exchange coefficient on T(n) */
478 cs_real_t f6 = 0.; /* thermal flux on Text */
479
480 cs_real_t m;
481
482 cs_real_t _al[128];
483 cs_real_t *al, *bl, *cl, *dl;
484 cs_real_t *zz;
485 cs_lnum_t n;
486
487 n = _1d_wall_thermal.local_models[ii].nppt1d;
488
489 if (n > 32)
490 BFT_MALLOC(al, 4*n, cs_real_t);
491 else
492 al = _al;
493
494 bl = al+n;
495 cl = bl+n;
496 dl = cl+n;
497
498 zz = _1d_wall_thermal.local_models[ii].z;
499
500 /* Build the tri-diagonal matrix */
501
502 /* Boundary conditions on the fluid side: flux conservation */
503 /* flux in the fluid = flux in the solid = f3 + h2*T1 */
504 a1 = 1./hf + zz[0]/xlmbt1;
505 h2 = -1./a1; // TAKE CARE TO THE MINUS !
506 f3 = -h2*tf + qinc;
507
508 /* Boundary conditions on the exterior */
509 /* flux in the fluid = flux in the solid = f6 + h5*T(n-1) */
510
511 /* Dirichlet condition */
512 if (iclt1d == 1) {
513 a4 = 1./hept1d + (eppt1d - zz[n-1])/xlmbt1;
514 h5 = -1./a4;
515 f6 = -h5*tept1d;
516 }
517 /* Forced flux condition */
518 else if (iclt1d == 3) {
519 h5 = 0.;
520 f6 = fept1d;
521 }
522
523 /* Mesh interior points */
524 for (kk = 1 ; kk <= n-1; kk++) {
525 al[kk] = -xlmbt1/(zz[kk]-zz[kk-1]);
526 }
527
528 m = 2*zz[0];
529 for (kk = 1 ; kk <= n-2 ; kk++) {
530 m = 2*(zz[kk]-zz[kk-1])-m;
531 bl[kk] = rcp/dtpt1d*m + xlmbt1/(zz[kk+1]-zz[kk]) + xlmbt1/(zz[kk]-zz[kk-1]);
532 }
533
534 for (kk = 0; kk <= n-2; kk++) {
535 cl[kk] = -xlmbt1/(zz[kk+1]-zz[kk]);
536 }
537
538 m = 2*zz[0];
539 dl[0] = rcp/dtpt1d*m*(_1d_wall_thermal.local_models[ii].t)[0];
540
541 for (kk = 1; kk <= n-1; kk++) {
542 m = 2*(zz[kk]-zz[kk-1])-m;
543 dl[kk] = rcp/dtpt1d*m*(_1d_wall_thermal.local_models[ii].t)[kk];
544 }
545
546 /* Boundary points */
547 /* bl[0] and bl[n-1] are initialized here and set later,
548 in the case where 0 = n-1 */
549 bl[0] = 0.;
550 bl[n-1] = 0.;
551 al[0] = 0.;
552 bl[0] += rcp/dtpt1d*2*zz[0] + xlmbt1/(zz[1]-zz[0]) - h2
553 + eps*cs_physical_constants_stephan
554 * pow(_1d_wall_thermal.local_models[ii].t[0],3.);
555 dl[0] += f3;
556 bl[n-1] += rcp/dtpt1d*2*(_1d_wall_thermal.local_models[ii].eppt1d-zz[n-1])
557 + xlmbt1/(zz[n-1]-zz[n-2]) - h5;
558 cl[n-1] = 0.;
559 dl[n-1] += f6;
560
561 /* System resolution by a Cholesky method ("dual-scan") */
562 for (kk = 1 ; kk <= n-1 ; kk++) {
563 bl[kk] -= al[kk]*cl[kk-1]/bl[kk-1];
564 dl[kk] -= al[kk]*dl[kk-1]/bl[kk-1];
565 }
566
567 _1d_wall_thermal.local_models[ii].t[n-1] = dl[n-1]/bl[n-1];
568
569 for (kk=n-2; kk>=0; kk--) {
570 _1d_wall_thermal.local_models[ii].t[kk] =
571 (dl[kk] - cl[kk]*_1d_wall_thermal.local_models[ii].t[kk+1])/bl[kk];
572 }
573
574 /* Compute the new value of tp */
575 _1d_wall_thermal.tppt1d[ii] = hf + xlmbt1/zz[0];
576 _1d_wall_thermal.tppt1d[ii]
577 = 1./_1d_wall_thermal.tppt1d[ii]
578 *(xlmbt1*_1d_wall_thermal.local_models[ii].t[0]/zz[0]
579 + hf*tf);
580
581 if (al != _al)
582 BFT_FREE(al);
583 }
584
585 /*----------------------------------------------------------------------------*/
586 /*!
587 * \brief Read the restart file of the 1D-wall thermal module.
588 */
589 /*----------------------------------------------------------------------------*/
590
591 void
cs_1d_wall_thermal_read(void)592 cs_1d_wall_thermal_read(void)
593 {
594 bool corresp_cel, corresp_fac, corresp_fbr, corresp_som;
595 cs_lnum_t nbvent;
596 cs_lnum_t ii, jj, ifac, indfac, ierror;
597 cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
598
599 char nomsui[] = "1dwall_module.csc";
600
601 cs_restart_t *suite;
602 cs_mesh_location_type_t support;
603 cs_restart_val_type_t typ_val;
604
605 ierror = CS_RESTART_SUCCESS;
606
607 /* Computation of nmxt1d */
608 for (ii = 0; ii < _1d_wall_thermal.nfpt1d ; ii++) {
609 _1d_wall_thermal.nmxt1d = CS_MAX(_1d_wall_thermal.local_models[ii].nppt1d,
610 _1d_wall_thermal.nmxt1d);
611 }
612
613 /* if necessary, sum over all the processors */
614 cs_parall_max(1, CS_INT_TYPE, &_1d_wall_thermal.nmxt1d);
615
616 /* Open the restart file */
617 cs_glob_tpar1d_suite = cs_restart_create(nomsui,
618 NULL,
619 CS_RESTART_MODE_READ);
620
621 if (cs_glob_tpar1d_suite == NULL)
622 bft_error(__FILE__, __LINE__, 0,
623 _("Abort while opening the 1D-wall thermal module restart file "
624 "in read mode.\n"
625 "Verify the existence and the name of the restart file: %s\n"),
626 nomsui);
627
628 /* Pointer to the global restart structure */
629 suite = cs_glob_tpar1d_suite;
630
631 /* Verification of the associated "support" to the restart file */
632 cs_restart_check_base_location(suite, &corresp_cel, &corresp_fac,
633 &corresp_fbr, &corresp_som);
634
635 /* Only boundary faces are of interest */
636 indfac = (corresp_fbr == true ? 1 : 0);
637 if (indfac == 0)
638 bft_error(__FILE__, __LINE__, 0,
639 _("Abort while reading the 1D-wall thermal module restart file.\n"
640 "The number of boundary faces has been modified\n"
641 "Verify that the restart file corresponds to "
642 "the present study.\n"));
643
644
645 {
646 /* Read the header */
647 cs_lnum_t *tabvar;
648
649 BFT_MALLOC(tabvar, 1, cs_lnum_t);
650
651 nbvent = 1;
652 support = CS_MESH_LOCATION_NONE;
653 typ_val = CS_TYPE_int;
654
655 ierror = cs_restart_read_section(suite,
656 "version_fichier_suite_module_1d",
657 support,
658 nbvent,
659 typ_val,
660 tabvar);
661
662 if (ierror < CS_RESTART_SUCCESS)
663 bft_error(__FILE__, __LINE__, 0,
664 _("WARNING: ABORT WHILE READING THE RESTART FILE\n"
665 "******** 1D-WALL THERMAL MODULE\n"
666 " INCORRECT FILE TYPE\n"
667 "\n"
668 "The file %s does not seem to be a restart file\n"
669 "for the 1D-wall thermal module.\n"
670 "The calculation will not be run.\n"
671 "\n"
672 "Verify that the restart file corresponds to a\n"
673 "restart file for the 1D-wall thermal module.\n"),
674 nomsui);
675
676 BFT_FREE(tabvar);
677 }
678
679 {
680 /* Read the number of discretiaztion points and coherency checks
681 with the input data of USPT1D */
682 cs_lnum_t *tabvar;
683 cs_lnum_t mfpt1d;
684 cs_gnum_t mfpt1t;
685 int iok;
686
687 BFT_MALLOC(tabvar, n_b_faces, cs_lnum_t);
688
689 const char nomrub[] = "nb_pts_discretis";
690 nbvent = 1;
691 support = CS_MESH_LOCATION_BOUNDARY_FACES;
692 typ_val = CS_TYPE_int;
693
694 ierror = cs_restart_read_section(suite,
695 nomrub,
696 support,
697 nbvent,
698 typ_val,
699 tabvar);
700
701 if (ierror < CS_RESTART_SUCCESS)
702 bft_error(__FILE__, __LINE__, 0,
703 _("Problem while reading section in the restart file\n"
704 "for the 1D-wall thermal module:\n"
705 "<%s>\n"
706 "The calculation will not be run.\n"), nomrub);
707
708 /* Coherency checks between the read NFPT1T and the one from USPT1D */
709 mfpt1d = 0;
710 for (ifac = 0; ifac < n_b_faces; ifac++) {
711 if (tabvar[ifac] > 0) mfpt1d++;
712 }
713 mfpt1t = mfpt1d;
714 /* if necessary, sum over all the processors */
715 cs_parall_counter(&mfpt1t, 1);
716 if (mfpt1t != _1d_wall_thermal.nfpt1t)
717 bft_error(__FILE__, __LINE__, 0,
718 _("WARNING: ABORT WHILE READING THE RESTART FILE\n"
719 "******** 1D-WALL THERMAL MODULE\n"
720 " CURRENT AND PREVIOUS DATA ARE DIFFERENT\n"
721 "\n"
722 "The number of faces with 1D thermal module has\n"
723 "been modified.\n"
724 "PREVIOUS: %lu boundary faces (total)\n"
725 "CURRENT: %lu boundary faces (total)\n"
726 "\n"
727 "The calculation will not be run.\n"
728 "\n"
729 "Verify that the restart file corresponds to a\n"
730 "restart file for the 1D-wall thermal module.\n"
731 "Verify uspt1d.\n"), mfpt1t, _1d_wall_thermal.nfpt1t);
732
733 /* Coherency check between read NFPT1D/IFPT1D and the ones from USPT1D
734 One already knows that the number of faces are equal, it is then
735 enough to check that all selected faces in USPT1D were also
736 selected in the previous calculation */
737 iok = 0;
738 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d; ii++){
739 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
740 if (tabvar[ifac] != _1d_wall_thermal.local_models[ii].nppt1d)
741 iok++;
742 }
743 if (iok > 0)
744 bft_error(__FILE__, __LINE__, 0,
745 _("WARNING: ABORT WHILE READING THE RESTART FILE\n"
746 "******** 1D-WALL THERMAL MODULE\n"
747 " CURRENT AND PREVIOUS DATA ARE DIFFERENT\n"
748 "\n"
749 "IFPT1D or NPPT1D has been modified with respect\n"
750 "to the restart file on at least on face with\n"
751 "1D thermal module\n"
752 "\n"
753 "The calculation will not be run.\n"
754 "\n"
755 "Verify that the restart file correspond to\n"
756 "the present study.\n"
757 "Verify uspt1d\n"
758 "(refer to the user manual for the specificities\n"
759 "of the test on IFPT1D)"));
760
761 /* Allocate the cs_glob_par1d structure */
762 _1d_wall_thermal_local_models_init();
763
764 BFT_FREE(tabvar);
765 }
766
767 {
768 /* Read the wall thickness and check the coherency with USPT1D*/
769 cs_real_t *tabvar;
770 cs_lnum_t iok;
771 cs_real_t eppt1d;
772
773 BFT_MALLOC(tabvar, n_b_faces, cs_real_t);
774
775 const char nomrub[] = "epaisseur_paroi";
776 nbvent = 1;
777 support = CS_MESH_LOCATION_BOUNDARY_FACES;
778 typ_val = CS_TYPE_cs_real_t;
779
780 ierror = cs_restart_read_section(suite,
781 nomrub,
782 support,
783 nbvent,
784 typ_val,
785 tabvar);
786
787 if (ierror < CS_RESTART_SUCCESS)
788 bft_error(__FILE__, __LINE__, 0,
789 _("Problem while reading section in the restart file\n"
790 "for the 1D-wall thermal module:\n"
791 "<%s>\n"
792 "The calculation will not be run.\n"), nomrub);
793
794 /* Coherence check between the read EPPT1D and the one from USPT1D */
795 iok = 0;
796 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
797 eppt1d = _1d_wall_thermal.local_models[ii].eppt1d;
798 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
799 if (fabs(tabvar[ifac] - eppt1d)/eppt1d > 1.e-10) iok++;
800 }
801 if (iok > 0)
802 bft_error(__FILE__, __LINE__, 0,
803 _("WARNING: ABORT WHILE READING THE RESTART FILE\n"
804 "******** 1D-WALL THERMAL MODULE\n"
805 " CURRENT AND PREVIOUS DATA ARE DIFFERENT\n"
806 "\n"
807 "The parameter EPPT1D has been modified with respect\n"
808 "to the restart file on at least on face with\n"
809 "1D thermal module\n"
810 "\n"
811 "The calculation will not be run.\n"
812 "\n"
813 "Verify that the restart file corresponds to\n"
814 "the present study.\n"
815 "Verify uspt1d\n"));
816
817 for (ii = 0; ii < _1d_wall_thermal.nfpt1d; ii++) {
818 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
819 _1d_wall_thermal.local_models[ii].eppt1d = tabvar[ifac];
820 }
821
822 BFT_FREE(tabvar);
823 }
824
825 {
826 /* Read the interior boundary temperature */
827 cs_real_t *tabvar;
828
829 BFT_MALLOC(tabvar, n_b_faces, cs_real_t);
830
831 const char nomrub[] = "temperature_bord_int";
832 nbvent = 1;
833 support = CS_MESH_LOCATION_BOUNDARY_FACES;
834 typ_val = CS_TYPE_cs_real_t;
835
836 ierror = cs_restart_read_section(suite,
837 nomrub,
838 support,
839 nbvent,
840 typ_val,
841 tabvar);
842
843 if (ierror < CS_RESTART_SUCCESS)
844 bft_error(__FILE__, __LINE__, 0,
845 _("Problem while reading section in the restart file\n"
846 "for the 1D-wall thermal module:\n"
847 "<%s>\n"
848 "The calculation will not be run.\n"), nomrub);
849
850 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
851 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
852 _1d_wall_thermal.tppt1d[ii] = tabvar[ifac];
853 }
854
855 BFT_FREE(tabvar);
856 }
857
858 { /* Read the 1D-mesh coordinates */
859 cs_lnum_t nptmx;
860 cs_lnum_t iok;
861 cs_real_t *tabvar;
862 cs_real_t zz1, zz2, rrgpt1, rgpt1d;
863
864 nptmx = n_b_faces * _1d_wall_thermal.nmxt1d;
865 BFT_MALLOC(tabvar, nptmx, cs_real_t);
866
867 const char nomrub[] = "coords_maillages_1d";
868 nbvent = _1d_wall_thermal.nmxt1d;
869 support = CS_MESH_LOCATION_BOUNDARY_FACES;
870 typ_val = CS_TYPE_cs_real_t;
871
872 ierror = cs_restart_read_section(suite,
873 nomrub,
874 support,
875 nbvent,
876 typ_val,
877 tabvar);
878
879 if (ierror < CS_RESTART_SUCCESS)
880 bft_error(__FILE__, __LINE__, 0,
881 _("Problem while reading section in the restart file\n"
882 "for the 1D-wall thermal module:\n"
883 "<%s>\n"
884 "The calculation will not be run.\n"), nomrub);
885
886 /* Now one have the cell centers, RGPT1D can be tested */
887 iok = 0;
888 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
889 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
890 if (_1d_wall_thermal.local_models[ii].nppt1d > 1) {
891 zz1 = tabvar[0 + _1d_wall_thermal.nmxt1d*ifac];
892 zz2 = tabvar[1 + _1d_wall_thermal.nmxt1d*ifac];
893 rrgpt1 = (zz2-2.*zz1)/zz1;
894 rgpt1d = _1d_wall_thermal.local_models[ii].rgpt1d;
895 if (fabs(rrgpt1-rgpt1d)/rgpt1d > 1.e-10) iok++;
896 }
897 }
898
899 if (iok > 0)
900 bft_error(__FILE__, __LINE__, 0,
901 _("WARNING: ABORT WHILE READING THE RESTART FILE\n"
902 "******** 1D-WALL THERMAL MODULE\n"
903 " CURRENT AND OLD DATA ARE DIFFERENT\n"
904 "\n"
905 "The parameter RGPT1D has been modified with respect\n"
906 "to the restart file on at least on face with\n"
907 "1D thermal module\n"
908 "\n"
909 "The calculation will not be run.\n"
910 "\n"
911 "Verify that the restart file correspond to\n"
912 "the present study\n"
913 "Verify uspt1d\n"));
914
915 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
916 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
917 /* The array is filled until the number of discretization points of
918 the given face is reached */
919 for (jj = 0; jj < _1d_wall_thermal.local_models[ii].nppt1d ; jj++)
920 _1d_wall_thermal.local_models[ii].z[jj]
921 = tabvar[jj + _1d_wall_thermal.nmxt1d*ifac];
922 }
923
924 BFT_FREE(tabvar);
925 }
926
927 { /* Read the wall temperature */
928 cs_lnum_t nptmx;
929 cs_real_t *tabvar;
930
931 nptmx = n_b_faces * _1d_wall_thermal.nmxt1d;
932 BFT_MALLOC(tabvar, nptmx, cs_real_t);
933
934 const char nomrub[] = "temperature_interne";
935 nbvent = _1d_wall_thermal.nmxt1d;
936 support = CS_MESH_LOCATION_BOUNDARY_FACES;
937 typ_val = CS_TYPE_cs_real_t;
938
939 ierror = cs_restart_read_section(suite,
940 nomrub,
941 support,
942 nbvent,
943 typ_val,
944 tabvar);
945
946 if (ierror < CS_RESTART_SUCCESS) {
947 cs_base_warn(__FILE__,__LINE__);
948 bft_printf(_("Problem while reading the section in the restart file\n"
949 "for the 1D-wall thermal module:\n"
950 "<%s>\n"), nomrub);
951 }
952
953 for (ii = 0; ii < _1d_wall_thermal.nfpt1d ; ii++) {
954 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
955
956 /* The array is filled until the number of discretization points of
957 the given face is reached */
958 for (jj = 0 ; jj < _1d_wall_thermal.local_models[ii].nppt1d ; jj++)
959 _1d_wall_thermal.local_models[ii].t[jj]
960 = tabvar[jj + _1d_wall_thermal.nmxt1d*ifac];
961
962 }
963
964 BFT_FREE(tabvar);
965 }
966
967 cs_restart_read_fields(suite, CS_RESTART_1D_WALL_THERMAL);
968
969 /* Close the restart file and free structures */
970 cs_restart_destroy(&cs_glob_tpar1d_suite);
971 }
972
973 /*----------------------------------------------------------------------------*/
974 /*!
975 * \brief Write the restart file of the 1D-wall thermal module.
976 */
977 /*----------------------------------------------------------------------------*/
978
979 void
cs_1d_wall_thermal_write(void)980 cs_1d_wall_thermal_write(void)
981 {
982 cs_lnum_t nbvent;
983 cs_lnum_t ii, jj, ifac;
984
985 char nomsui[] = "1dwall_module.csc";
986
987 cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
988
989 cs_restart_t *suite;
990 cs_mesh_location_type_t support;
991 cs_restart_val_type_t typ_val;
992
993 /* Open the restart file */
994 cs_glob_tpar1d_suite = cs_restart_create(nomsui,
995 NULL,
996 CS_RESTART_MODE_WRITE);
997
998 if (cs_glob_tpar1d_suite == NULL)
999 bft_error(__FILE__, __LINE__, 0,
1000 _("Abort while opening the 1D-wall thermal module restart "
1001 "file in write mode.\n"
1002 "Verify the existence and the name of the restart file: %s\n"),
1003 nomsui);
1004
1005 /* Pointer to the global restart structure */
1006 suite = cs_glob_tpar1d_suite;
1007
1008 {
1009 /* Write the header */
1010 cs_lnum_t tabvar[1] = {120};
1011
1012 nbvent = 1;
1013 support = CS_MESH_LOCATION_NONE;
1014 typ_val = CS_TYPE_int;
1015
1016 cs_restart_write_section(suite,
1017 "version_fichier_suite_module_1d",
1018 support,
1019 nbvent,
1020 typ_val,
1021 tabvar);
1022 }
1023
1024 {
1025 /* Write the number of discretization points */
1026 cs_lnum_t *tabvar;
1027
1028 BFT_MALLOC(tabvar, n_b_faces, cs_lnum_t);
1029
1030 for (ii = 0 ; ii < n_b_faces ; ii++)
1031 tabvar[ii] = 0;
1032
1033 nbvent = 1;
1034 support = CS_MESH_LOCATION_BOUNDARY_FACES;
1035 typ_val = CS_TYPE_int;
1036
1037 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
1038 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
1039 tabvar[ifac] =_1d_wall_thermal.local_models[ii].nppt1d;
1040 }
1041
1042 cs_restart_write_section(suite,
1043 "nb_pts_discretis",
1044 support,
1045 nbvent,
1046 typ_val,
1047 tabvar);
1048
1049 BFT_FREE(tabvar);
1050 }
1051
1052 {
1053 /* Write the wall thickness */
1054 cs_real_t *tabvar;
1055
1056 BFT_MALLOC(tabvar, n_b_faces, cs_real_t);
1057
1058 for (ii = 0 ; ii < n_b_faces ; ii++)
1059 tabvar[ii] = 0.;
1060
1061 nbvent = 1;
1062 support = CS_MESH_LOCATION_BOUNDARY_FACES;
1063 typ_val = CS_TYPE_cs_real_t;
1064
1065 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
1066 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
1067 tabvar[ifac] = _1d_wall_thermal.local_models[ii].eppt1d;
1068 }
1069
1070 cs_restart_write_section(suite,
1071 "epaisseur_paroi",
1072 support,
1073 nbvent,
1074 typ_val,
1075 tabvar);
1076
1077 BFT_FREE(tabvar);
1078 }
1079
1080 {
1081 /* Write the internal wall-boundary temperature */
1082 cs_real_t *tabvar;
1083
1084 BFT_MALLOC(tabvar, n_b_faces, cs_real_t);
1085
1086 for (ii = 0 ; ii < n_b_faces ; ii++)
1087 tabvar[ii] = 0.0;
1088
1089 nbvent = 1;
1090 support = CS_MESH_LOCATION_BOUNDARY_FACES;
1091 typ_val = CS_TYPE_cs_real_t;
1092
1093 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d; ii++) {
1094 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
1095 tabvar[ifac] = _1d_wall_thermal.tppt1d[ii];
1096 }
1097
1098 cs_restart_write_section(suite,
1099 "temperature_bord_int",
1100 support,
1101 nbvent,
1102 typ_val,
1103 tabvar);
1104
1105 BFT_FREE(tabvar);
1106 }
1107
1108 {
1109 /* Write the 1D-mesh coordinates */
1110 cs_lnum_t nptmx;
1111 cs_real_t *tabvar;
1112
1113 nptmx = n_b_faces * _1d_wall_thermal.nmxt1d;
1114 BFT_MALLOC(tabvar, nptmx, cs_real_t);
1115
1116 for (ii = 0 ; ii < nptmx ; ii++)
1117 tabvar[ii] = 0.;
1118
1119 nbvent = _1d_wall_thermal.nmxt1d;
1120 support = CS_MESH_LOCATION_BOUNDARY_FACES;
1121 typ_val = CS_TYPE_cs_real_t;
1122
1123 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
1124 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
1125
1126 /* The array is filled until the number of discretization points of
1127 the given face is reached (the following cases, up to nmxt1d, are
1128 already set to 0 during the initalization of tabvar */
1129 for (jj = 0 ; jj < _1d_wall_thermal.local_models[ii].nppt1d ; jj++)
1130 tabvar[jj + _1d_wall_thermal.nmxt1d*ifac]
1131 = _1d_wall_thermal.local_models[ii].z[jj];
1132 }
1133
1134 cs_restart_write_section(suite,
1135 "coords_maillages_1d",
1136 support,
1137 nbvent,
1138 typ_val,
1139 tabvar);
1140
1141 BFT_FREE(tabvar);
1142 }
1143
1144 {
1145 /* Write the wall-interior temperature */
1146 cs_lnum_t nptmx;
1147 cs_real_t *tabvar;
1148
1149 nptmx = n_b_faces * _1d_wall_thermal.nmxt1d;
1150 BFT_MALLOC(tabvar, nptmx, cs_real_t);
1151
1152 for (ii = 0 ; ii < nptmx ; ii++)
1153 tabvar[ii] = 0.;
1154
1155 for (ii = 0 ; ii < _1d_wall_thermal.nfpt1d ; ii++) {
1156 ifac = _1d_wall_thermal.ifpt1d[ii] - 1;
1157
1158 /* The array is filled until the number of discretization points of
1159 the given face is reached (the following cases, up to nmxt1d, are
1160 already set to 0 during the initalization of tabvar */
1161 for (jj = 0 ; jj < _1d_wall_thermal.local_models[ii].nppt1d ; jj++)
1162 tabvar[jj + _1d_wall_thermal.nmxt1d*ifac]
1163 = _1d_wall_thermal.local_models[ii].t[jj];
1164
1165 }
1166
1167 cs_restart_write_section(suite,
1168 "temperature_interne",
1169 CS_MESH_LOCATION_BOUNDARY_FACES,
1170 _1d_wall_thermal.nmxt1d,
1171 CS_TYPE_cs_real_t,
1172 tabvar);
1173
1174 BFT_FREE(tabvar);
1175 }
1176
1177 cs_restart_write_fields(suite, CS_RESTART_1D_WALL_THERMAL);
1178
1179 /* Close the restart file and free structures */
1180 cs_restart_destroy(&cs_glob_tpar1d_suite);
1181 }
1182
1183 /*----------------------------------------------------------------------------*/
1184 /*!
1185 * \brief Free the array of structures local_models.
1186 */
1187 /*----------------------------------------------------------------------------*/
1188
1189 void
cs_1d_wall_thermal_free(void)1190 cs_1d_wall_thermal_free(void)
1191 {
1192 if (_1d_wall_thermal.local_models != NULL)
1193 BFT_FREE(_1d_wall_thermal.local_models->z);
1194 BFT_FREE(_1d_wall_thermal.local_models);
1195 BFT_FREE(_1d_wall_thermal.ifpt1d);
1196 BFT_FREE(_1d_wall_thermal.tppt1d);
1197 }
1198
1199 /*----------------------------------------------------------------------------*/
1200 /*!
1201 * \brief Destroy the global 1d wall thermal structure.
1202 */
1203 /*----------------------------------------------------------------------------*/
1204
1205 void
cs_1d_wall_thermal_finalize(void)1206 cs_1d_wall_thermal_finalize(void)
1207 {
1208 BFT_FREE(_1d_wall_thermal.izft1d);
1209 cs_glob_1d_wall_thermal = NULL;
1210 }
1211
1212 /*----------------------------------------------------------------------------*/
1213 /*!
1214 * \brief Provide access to cs_glob_1d_wall_thermal.
1215 */
1216 /*----------------------------------------------------------------------------*/
1217
1218 cs_1d_wall_thermal_t *
cs_get_glob_1d_wall_thermal(void)1219 cs_get_glob_1d_wall_thermal(void)
1220 {
1221 return &_1d_wall_thermal;
1222 }
1223
1224 /*----------------------------------------------------------------------------*/
1225
1226 END_C_DECLS
1227