1 /*============================================================================
2 * Turbomachinery modeling features.
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 You should have received a copy of the GNU General Public License along with
20 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
21 Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 */
23
24 /*----------------------------------------------------------------------------*/
25
26 #include "cs_defs.h"
27
28 /*----------------------------------------------------------------------------
29 * Standard C library headers
30 *----------------------------------------------------------------------------*/
31
32 #include <assert.h>
33 #include <math.h>
34 #include <string.h>
35
36 /*----------------------------------------------------------------------------
37 * Local headers
38 *----------------------------------------------------------------------------*/
39
40 #include "bft_mem.h"
41 #include "bft_printf.h"
42
43 #include "fvm_selector.h"
44
45 #include "cs_interface.h"
46
47 #include "cs_base.h"
48 #include "cs_boundary_zone.h"
49 #include "cs_coupling.h"
50 #include "cs_cell_to_vertex.h"
51 #include "cs_ext_neighborhood.h"
52 #include "cs_gradient.h"
53 #include "cs_gui.h"
54 #include "cs_gui_mesh.h"
55 #include "cs_gui_output.h"
56 #include "cs_gradient.h"
57 #include "cs_join.h"
58 #include "cs_halo.h"
59 #include "cs_halo_perio.h"
60 #include "cs_matrix_default.h"
61 #include "cs_mesh.h"
62 #include "cs_mesh_adjacencies.h"
63 #include "cs_mesh_coherency.h"
64 #include "cs_mesh_location.h"
65 #include "cs_mesh_quantities.h"
66 #include "cs_mesh_to_builder.h"
67 #include "cs_multigrid.h"
68 #include "cs_parall.h"
69 #include "cs_post.h"
70 #include "cs_preprocess.h"
71 #include "cs_prototypes.h"
72 #include "cs_renumber.h"
73 #include "cs_rotation.h"
74 #include "cs_time_step.h"
75 #include "cs_timer.h"
76 #include "cs_timer_stats.h"
77 #include "cs_restart.h"
78 #include "cs_sat_coupling.h"
79 #include "cs_preprocessor_data.h"
80 #include "cs_volume_zone.h"
81
82 /*----------------------------------------------------------------------------
83 * Header for the current file
84 *----------------------------------------------------------------------------*/
85
86 #include "cs_turbomachinery.h"
87
88 /*----------------------------------------------------------------------------*/
89
90 BEGIN_C_DECLS
91
92 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
93
94 /*============================================================================
95 * Local structure definitions
96 *============================================================================*/
97
98 /* Turbomachinery structure */
99
100 typedef struct {
101
102 cs_turbomachinery_model_t model; /* turbomachinery model type */
103
104 int n_rotors; /* number of rotors */
105
106 int n_couplings; /* number of couplings */
107
108 cs_rotation_t *rotation; /* rotation structures */
109 char **rotor_cells_c; /* rotor cells selection
110 criteria (for each rotor) */
111
112 int n_max_join_tries; /* maximum number of tries
113 for joining differences */
114 double dt_retry; /* time shift multiplier for
115 retry position */
116 double t_cur; /* current time for update */
117
118 cs_mesh_t *reference_mesh; /* reference mesh (before
119 rotation and joining) */
120
121 cs_lnum_t n_b_faces_ref; /* reference number of
122 boundary faces */
123
124 int *cell_rotor_num; /* cell rotation axis number */
125
126 bool active;
127
128 } cs_turbomachinery_t;
129
130 /*============================================================================
131 * Static global variables
132 *============================================================================*/
133
134 cs_turbomachinery_t *_turbomachinery = NULL;
135
136 /*============================================================================
137 * Prototypes for functions intended for use only by Fortran wrappers.
138 * (descriptions follow, with function bodies).
139 *============================================================================*/
140
141 void cs_f_map_turbomachinery_model(int *iturbo,
142 int *ityint);
143
144 void cs_f_map_turbomachinery_rotor(int **irotce);
145
146 /*============================================================================
147 * Private function definitions
148 *============================================================================*/
149
150 /*----------------------------------------------------------------------------
151 * Function for selection of faces with joining changes.
152 *
153 * parameters:
154 * input <-- pointer to input (count[0]: n_previous, count[1]: n_current)
155 * n_faces --> number of selected faces
156 * face_ids --> array of selected face ids (0 to n-1 numbering)
157 *----------------------------------------------------------------------------*/
158
159 static void
_post_error_faces_select(void * input,cs_lnum_t * n_faces,cs_lnum_t ** face_ids)160 _post_error_faces_select(void *input,
161 cs_lnum_t *n_faces,
162 cs_lnum_t **face_ids)
163 {
164 cs_lnum_t face_id;
165
166 cs_lnum_t _n_faces = 0;
167 cs_lnum_t *_face_ids = NULL;
168
169 const cs_lnum_t *count = input;
170
171 BFT_MALLOC(_face_ids, count[1], cs_lnum_t);
172
173 for (face_id = count[0]; face_id < count[1]; face_id++)
174 _face_ids[_n_faces++] = face_id;
175
176 *n_faces = _n_faces;
177 *face_ids = _face_ids;
178 }
179
180 /*----------------------------------------------------------------------------
181 * Function for tagging of coupling mesh.
182 *
183 * We use the rotor number as a tag
184 *
185 * Note: if the context pointer is non-NULL, it must point to valid data
186 * when the selection function is called, so that value or structure
187 * should not be temporary (i.e. local);
188 *
189 * parameters:
190 * context <-> pointer to optional (untyped) value or structure.
191 * mesh <-> nodal mesh which should be tagged
192 * n_points <-- number of points to tag
193 * point_list_base <-- base numbering for point_list
194 * point_list <-- optional indirection for points
195 * point_tag --> point tag values (size: n_tags)
196 *----------------------------------------------------------------------------*/
197
198 static void
_turbomachinery_coupling_tag(void * context,fvm_nodal_t * mesh,cs_lnum_t n_points,cs_lnum_t point_list_base,const cs_lnum_t point_list[],int * point_tag)199 _turbomachinery_coupling_tag(void *context,
200 fvm_nodal_t *mesh,
201 cs_lnum_t n_points,
202 cs_lnum_t point_list_base,
203 const cs_lnum_t point_list[],
204 int *point_tag)
205 {
206 const cs_turbomachinery_t *tbm = context;
207 const cs_mesh_t *m = cs_glob_mesh;
208
209 /* Tag elements (boundary faces) */
210
211 if (mesh != NULL) {
212
213 cs_lnum_t n_elts = fvm_nodal_get_n_entities(mesh, 3);
214
215 const int ent_dim = 3;
216 int *elt_tag;
217 cs_lnum_t *parent_num;
218
219 BFT_MALLOC(elt_tag, n_elts, int);
220 BFT_MALLOC(parent_num, n_elts, cs_lnum_t);
221
222 fvm_nodal_get_parent_num(mesh, ent_dim, parent_num);
223 for (cs_lnum_t i = 0; i < n_elts; i++) {
224 cs_lnum_t c_id = parent_num[i] - 1;
225 elt_tag[i] = tbm->cell_rotor_num[c_id];
226 }
227
228 BFT_FREE(parent_num);
229
230 fvm_nodal_set_tag(mesh, elt_tag, ent_dim);
231
232 BFT_FREE(elt_tag);
233
234 }
235
236 /* Tag vertices */
237
238 if (point_list != NULL) {
239 for (cs_lnum_t i = 0; i < n_points; i++) {
240 cs_lnum_t f_id = point_list[i] - point_list_base;
241 cs_lnum_t c_id = m->b_face_cells[f_id];
242 point_tag[i] = tbm->cell_rotor_num[c_id];
243 }
244 }
245 else {
246 for (cs_lnum_t i = 0; i < n_points; i++) {
247 cs_lnum_t c_id = m->b_face_cells[i];
248 point_tag[i] = tbm->cell_rotor_num[c_id];
249 }
250 }
251 }
252
253 /*----------------------------------------------------------------------------
254 * Create an empty turbomachinery structure
255 *----------------------------------------------------------------------------*/
256
257 static cs_turbomachinery_t *
_turbomachinery_create(void)258 _turbomachinery_create(void)
259 {
260 cs_turbomachinery_t *tbm = NULL;
261
262 BFT_MALLOC(tbm, 1, cs_turbomachinery_t);
263
264 tbm->n_rotors = 0;
265 tbm->rotor_cells_c = NULL;
266
267 BFT_MALLOC(tbm->rotation, 1, cs_rotation_t); /* Null rotation at id 0 */
268 cs_rotation_t *r = tbm->rotation;
269 r->omega = 0;
270 r->angle = 0;
271 for (int i = 0; i < 3; i++) {
272 r->axis[i] = 0;
273 r->invariant[i] = 0;
274 }
275
276 tbm->n_max_join_tries = 5;
277 tbm->t_cur = 0;
278 tbm->dt_retry = 1e-2;
279
280 tbm->reference_mesh = cs_mesh_create();
281 tbm->n_b_faces_ref = -1;
282 tbm->cell_rotor_num = NULL;
283 tbm->model = CS_TURBOMACHINERY_NONE;
284 tbm->n_couplings = 0;
285
286 return tbm;
287 }
288
289 /*----------------------------------------------------------------------------
290 * Compute a matrix/vector product to apply a transformation to a vector.
291 *
292 * parameters:
293 * m[3][4] <-- matrix of the transformation in homogeneous coord.
294 * last line = [0; 0; 0; 1] (Not used here)
295 * c[3] <-> coordinates
296 *----------------------------------------------------------------------------*/
297
298 static inline void
_apply_vector_transfo(double matrix[3][4],cs_real_t c[3])299 _apply_vector_transfo(double matrix[3][4],
300 cs_real_t c[3])
301 {
302 int i, j;
303
304 double c_a[4] = {c[0], c[1], c[2], 1.}; /* homogeneous coords */
305 double c_b[3] = {0, 0, 0};
306
307 for (i = 0; i < 3; i++)
308 for (j = 0; j < 4; j++)
309 c_b[i] += matrix[i][j]*c_a[j];
310
311 for (i = 0; i < 3; i++)
312 c[i] = c_b[i];
313 }
314
315 /*----------------------------------------------------------------------------
316 * Compute a matrix/vector product to apply a rotation to a vector.
317 *
318 * parameters:
319 * m[3][4] <-- matrix of the transformation in homogeneous coord.
320 * last line = [0; 0; 0; 1] (Not used here)
321 * v[3] <-> vector
322 *----------------------------------------------------------------------------*/
323
324 static inline void
_apply_vector_rotation(double m[3][4],cs_real_t v[3])325 _apply_vector_rotation(double m[3][4],
326 cs_real_t v[3])
327 {
328 double t[3] = {v[0], v[1], v[2]};
329
330 for (int i = 0; i < 3; i++)
331 v[i] = m[i][0]*t[0] + m[i][1]*t[1] + m[i][2]*t[2];
332 }
333
334 /*----------------------------------------------------------------------------
335 * Compute a matrix * tensor * Tmatrix product to apply a rotation to a
336 * given symmetric tensor
337 *
338 * parameters:
339 * matrix[3][4] --> transformation matrix in homogeneous coords.
340 * last line = [0; 0; 0; 1] (Not used here)
341 * tensor[6] <-> incoming (6) symmetric tensor
342 *----------------------------------------------------------------------------*/
343
344 static inline void
_apply_sym_tensor_rotation(cs_real_t matrix[3][4],cs_real_t t[6])345 _apply_sym_tensor_rotation(cs_real_t matrix[3][4],
346 cs_real_t t[6])
347 {
348 int i, j, k, l;
349
350 double _t[3][3], _t0[3][3];
351
352 _t0[0][0] = t[0];
353 _t0[1][1] = t[1];
354 _t0[2][2] = t[2];
355 _t0[0][1] = t[3];
356 _t0[1][0] = t[3];
357 _t0[1][2] = t[4];
358 _t0[2][1] = t[4];
359 _t0[0][2] = t[5];
360 _t0[2][0] = t[5];
361
362 for (k = 0; k < 3; k++) {
363 for (j = 0; j < 3; j++) {
364 _t[k][j] = 0.;
365 for (l = 0; l < 3; l++)
366 _t[k][j] += matrix[j][l] * _t0[k][l];
367 }
368 }
369
370 for (i = 0; i < 3; i++) {
371 for (j = 0; j < 3; j++) {
372 _t0[i][j] = 0.;
373 for (k = 0; k < 3; k++)
374 _t0[i][j] += matrix[i][k] * _t[k][j];
375 }
376 }
377
378 t[0] = _t0[0][0];
379 t[1] = _t0[1][1];
380 t[2] = _t0[2][2];
381 t[3] = _t0[0][1];
382 t[3] = _t0[1][0];
383 t[4] = _t0[2][1];
384 t[5] = _t0[2][0];
385 }
386
387 /*----------------------------------------------------------------------------
388 * Compute velocity relative to fixed coordinates at a given point
389 *
390 * parameters:
391 * omega <-- rotation velocity
392 * axis <-- components of rotation axis direction vector (3)
393 * invariant_point <-- components of invariant point (3)
394 * coords <-- coordinates at point
395 * velocity --> resulting relative velocity
396 *---------------------------------------------------------------------------*/
397
398 static void
_relative_velocity(double omega,const double axis[3],const double invariant_point[3],const double coords[3],cs_real_t velocity[3])399 _relative_velocity(double omega,
400 const double axis[3],
401 const double invariant_point[3],
402 const double coords[3],
403 cs_real_t velocity[3])
404 {
405 velocity[0] = (- axis[2] * (coords[1] - invariant_point[1])
406 + axis[1] * (coords[2] - invariant_point[2])) * omega;
407 velocity[1] = ( axis[2] * (coords[0] - invariant_point[0])
408 - axis[0] * (coords[2] - invariant_point[2])) * omega;
409 velocity[2] = (- axis[1] * (coords[0] - invariant_point[0])
410 + axis[0] * (coords[1] - invariant_point[1])) * omega;
411 }
412
413 /*----------------------------------------------------------------------------
414 * Duplicate a cs_mesh_t structure.
415 *
416 * Note that some fields which are recomputable are not copied.
417 *
418 * In case of coupling only (i.e. no joining), only a partial copy of the
419 * mesh is done.
420 *
421 * parameters:
422 * mesh <-- reference mesh
423 * mesh_copy <-> mesh copy
424 *----------------------------------------------------------------------------*/
425
426 static void
_copy_mesh(const cs_mesh_t * mesh,cs_mesh_t * mesh_copy)427 _copy_mesh(const cs_mesh_t *mesh,
428 cs_mesh_t *mesh_copy)
429 {
430 cs_lnum_t n_elts;
431
432 /* General features */
433
434 mesh_copy->dim = mesh->dim;
435 mesh_copy->domain_num = mesh->domain_num;
436 mesh_copy->n_domains = mesh->n_domains;
437
438 /* Local dimensions */
439
440 mesh_copy->n_cells = mesh->n_cells;
441 mesh_copy->n_i_faces = mesh->n_i_faces;
442 mesh_copy->n_b_faces = mesh->n_b_faces;
443 mesh_copy->n_vertices = mesh->n_vertices;
444
445 mesh_copy->i_face_vtx_connect_size = mesh->i_face_vtx_connect_size;
446 mesh_copy->b_face_vtx_connect_size = mesh->b_face_vtx_connect_size;
447
448 /* Local structures */
449
450 BFT_REALLOC(mesh_copy->vtx_coord, (3*mesh->n_vertices), cs_real_t);
451 memcpy(mesh_copy->vtx_coord,
452 mesh->vtx_coord,
453 3*mesh->n_vertices*sizeof(cs_real_t));
454
455 if (cs_glob_n_joinings < 1)
456 return;
457
458 BFT_MALLOC(mesh_copy->i_face_cells, mesh->n_i_faces, cs_lnum_2_t);
459 memcpy(mesh_copy->i_face_cells,
460 mesh->i_face_cells,
461 mesh->n_i_faces*sizeof(cs_lnum_2_t));
462
463 if (mesh->n_b_faces > 0) {
464 BFT_MALLOC(mesh_copy->b_face_cells, mesh->n_b_faces, cs_lnum_t);
465 memcpy(mesh_copy->b_face_cells,
466 mesh->b_face_cells,
467 mesh->n_b_faces*sizeof(cs_lnum_t));
468 }
469
470 BFT_MALLOC(mesh_copy->i_face_vtx_idx, mesh->n_i_faces + 1, cs_lnum_t);
471 memcpy(mesh_copy->i_face_vtx_idx,
472 mesh->i_face_vtx_idx,
473 (mesh->n_i_faces + 1)*sizeof(cs_lnum_t));
474
475 BFT_MALLOC(mesh_copy->i_face_vtx_lst,
476 mesh->i_face_vtx_connect_size,
477 cs_lnum_t);
478 memcpy(mesh_copy->i_face_vtx_lst, mesh->i_face_vtx_lst,
479 mesh->i_face_vtx_connect_size*sizeof(cs_lnum_t));
480
481 BFT_MALLOC(mesh_copy->b_face_vtx_idx, mesh->n_b_faces + 1, cs_lnum_t);
482 memcpy(mesh_copy->b_face_vtx_idx,
483 mesh->b_face_vtx_idx,
484 (mesh->n_b_faces + 1)*sizeof(cs_lnum_t));
485
486 if (mesh->b_face_vtx_connect_size > 0) {
487 BFT_MALLOC(mesh_copy->b_face_vtx_lst,
488 mesh->b_face_vtx_connect_size,
489 cs_lnum_t);
490 memcpy(mesh_copy->b_face_vtx_lst,
491 mesh->b_face_vtx_lst,
492 mesh->b_face_vtx_connect_size*sizeof(cs_lnum_t));
493 }
494
495 /* Global dimension */
496
497 mesh_copy->n_g_cells = mesh->n_g_cells;
498 mesh_copy->n_g_i_faces = mesh->n_g_i_faces;
499 mesh_copy->n_g_b_faces = mesh->n_g_b_faces;
500 mesh_copy->n_g_vertices = mesh->n_g_vertices;
501
502 /* Global numbering */
503
504 if (mesh->global_cell_num != NULL) {
505 BFT_MALLOC(mesh_copy->global_cell_num, mesh->n_cells, cs_gnum_t);
506 memcpy(mesh_copy->global_cell_num,
507 mesh->global_cell_num,
508 mesh->n_cells*sizeof(cs_gnum_t));
509 }
510
511 if (mesh->global_i_face_num != NULL) {
512 BFT_MALLOC(mesh_copy->global_i_face_num, mesh->n_i_faces, cs_gnum_t);
513 memcpy(mesh_copy->global_i_face_num,
514 mesh->global_i_face_num,
515 mesh->n_i_faces*sizeof(cs_gnum_t));
516 }
517
518 if (mesh->global_b_face_num != NULL) {
519 BFT_MALLOC(mesh_copy->global_b_face_num, mesh->n_b_faces, cs_gnum_t);
520 memcpy(mesh_copy->global_b_face_num,
521 mesh->global_b_face_num,
522 mesh->n_b_faces*sizeof(cs_gnum_t));
523 }
524
525 if (mesh->global_vtx_num != NULL) {
526 BFT_MALLOC(mesh_copy->global_vtx_num, mesh->n_vertices, cs_gnum_t);
527 memcpy(mesh_copy->global_vtx_num,
528 mesh->global_vtx_num,
529 mesh->n_vertices*sizeof(cs_gnum_t));
530 }
531
532 /* Parallelism and/or periodic features */
533 /* and extended neighborhood features */
534
535 mesh_copy->n_init_perio = mesh->n_init_perio;
536 mesh_copy->n_transforms = mesh->n_transforms;
537 mesh_copy->have_rotation_perio = mesh->have_rotation_perio;
538
539 mesh_copy->halo_type = mesh->halo_type;
540
541 /* modified in cs_mesh_init_halo, if necessary */
542 /* ------------------------------------------- */
543 /* cs_lnum_t n_cells_with_ghosts; */
544 /* cs_lnum_t n_ghost_cells; */
545
546 mesh_copy->n_cells_with_ghosts = mesh->n_cells_with_ghosts;
547 mesh_copy->n_ghost_cells = mesh->n_ghost_cells;
548
549 /* Re-computable connectivity features */
550 /*-------------------------------------*/
551
552 mesh_copy->n_b_cells = mesh->n_b_cells;
553
554 BFT_MALLOC(mesh_copy->b_cells, mesh->n_b_cells, cs_lnum_t);
555 memcpy(mesh_copy->b_cells, mesh->b_cells, mesh->n_b_cells*sizeof(cs_lnum_t));
556
557 /* Group and family features */
558 /*---------------------------*/
559
560 mesh_copy->n_groups = mesh->n_groups;
561
562 if (mesh->n_groups > 0) {
563 BFT_MALLOC(mesh_copy->group_idx, mesh->n_groups + 1, int);
564 memcpy(mesh_copy->group_idx, mesh->group_idx,
565 (mesh->n_groups + 1)*sizeof(cs_lnum_t));
566 BFT_MALLOC(mesh_copy->group, mesh->group_idx[mesh->n_groups], char);
567 memcpy(mesh_copy->group, mesh->group,
568 mesh->group_idx[mesh->n_groups]*sizeof(char));
569 }
570
571 mesh_copy->n_families = mesh->n_families;
572 mesh_copy->n_max_family_items = mesh->n_max_family_items;
573
574 n_elts = mesh->n_families*mesh->n_max_family_items;
575 if (n_elts > 0) {
576 BFT_MALLOC(mesh_copy->family_item, n_elts , int);
577 memcpy(mesh_copy->family_item, mesh->family_item, n_elts*sizeof(int));
578 }
579
580 BFT_MALLOC(mesh_copy->cell_family, mesh->n_cells_with_ghosts, int);
581 memcpy(mesh_copy->cell_family, mesh->cell_family,
582 mesh->n_cells_with_ghosts*sizeof(int));
583
584 BFT_MALLOC(mesh_copy->i_face_family, mesh->n_i_faces, int);
585 memcpy(mesh_copy->i_face_family, mesh->i_face_family,
586 mesh->n_i_faces*sizeof(int));
587
588 if (mesh->n_b_faces > 0) {
589 BFT_MALLOC(mesh_copy->b_face_family, mesh->n_b_faces, int);
590 memcpy(mesh_copy->b_face_family, mesh->b_face_family,
591 mesh->n_b_faces*sizeof(int));
592 }
593
594 if (mesh->i_face_r_gen != NULL) {
595 BFT_MALLOC(mesh_copy->i_face_r_gen, mesh->n_i_faces, char);
596 memcpy(mesh_copy->i_face_r_gen, mesh->i_face_r_gen,
597 mesh->n_i_faces);
598 }
599 }
600
601 /*----------------------------------------------------------------------------
602 * Update mesh vertex positions
603 *
604 * parameters:
605 * t <-- associated time
606 *----------------------------------------------------------------------------*/
607
608 static void
_update_angle(cs_real_t t)609 _update_angle(cs_real_t t)
610 {
611 cs_turbomachinery_t *tbm = _turbomachinery;
612
613 /* Now update coordinates */
614
615 double dt = t - tbm->t_cur;
616
617 if (dt > 0) {
618 for (int j = 0; j < tbm->n_rotors+1; j++) {
619 cs_rotation_t *r = tbm->rotation + j;
620 r->angle += r->omega * dt;
621 }
622 tbm->t_cur = t;
623 }
624 }
625
626 /*----------------------------------------------------------------------------
627 * Update mesh vertex positions
628 *
629 * parameters:
630 * mesh <-> mesh to update
631 * dt <-- associated time delta (0 for current, unmodified time)
632 *----------------------------------------------------------------------------*/
633
634 static void
_update_geometry(cs_mesh_t * mesh,cs_real_t dt)635 _update_geometry(cs_mesh_t *mesh,
636 cs_real_t dt)
637 {
638 cs_turbomachinery_t *tbm = _turbomachinery;
639
640 cs_lnum_t f_id, v_id;
641
642 cs_lnum_t *vtx_rotor_num = NULL;
643
644 const int *cell_flag = tbm->cell_rotor_num;
645
646 BFT_MALLOC(vtx_rotor_num, mesh->n_vertices, cs_lnum_t);
647
648 for (v_id = 0; v_id < mesh->n_vertices; v_id++)
649 vtx_rotor_num[v_id] = 0;
650
651 /* Mark from interior faces */
652
653 for (f_id = 0; f_id < mesh->n_i_faces; f_id++) {
654 cs_lnum_t c_id_0 = mesh->i_face_cells[f_id][0];
655 cs_lnum_t c_id_1 = mesh->i_face_cells[f_id][1];
656 assert(c_id_0 > -1);
657 assert(c_id_1 > -1);
658 if (c_id_0 < mesh->n_cells && cell_flag[c_id_0] != 0) {
659 for (cs_lnum_t i = mesh->i_face_vtx_idx[f_id];
660 i < mesh->i_face_vtx_idx[f_id+1];
661 i++)
662 vtx_rotor_num[mesh->i_face_vtx_lst[i]] = cell_flag[c_id_0];
663 }
664 else if (c_id_1 < mesh->n_cells && cell_flag[c_id_1] != 0) {
665 for (cs_lnum_t i = mesh->i_face_vtx_idx[f_id];
666 i < mesh->i_face_vtx_idx[f_id+1];
667 i++)
668 vtx_rotor_num[mesh->i_face_vtx_lst[i]] = cell_flag[c_id_1];
669 }
670 }
671
672 /* Mark from boundary faces */
673
674 for (f_id = 0; f_id < mesh->n_b_faces; f_id++) {
675 cs_lnum_t c_id = mesh->b_face_cells[f_id];
676 if (cell_flag[c_id] != 0) {
677 for (cs_lnum_t i = mesh->b_face_vtx_idx[f_id];
678 i < mesh->b_face_vtx_idx[f_id+1];
679 i++)
680 vtx_rotor_num[mesh->b_face_vtx_lst[i]] = cell_flag[c_id];
681 }
682 }
683
684 /* Now update coordinates */
685
686 cs_real_34_t *m;
687
688 BFT_MALLOC(m, tbm->n_rotors+1, cs_real_34_t);
689
690 for (int j = 0; j < tbm->n_rotors+1; j++) {
691 cs_rotation_t *r = tbm->rotation + j;
692
693 cs_rotation_matrix(r->angle + r->omega*dt,
694 r->axis,
695 r->invariant,
696 m[j]);
697 }
698
699 for (v_id = 0; v_id < mesh->n_vertices; v_id++) {
700 if (vtx_rotor_num[v_id] > 0)
701 _apply_vector_transfo(m[vtx_rotor_num[v_id]],
702 &(mesh->vtx_coord[3*v_id]));
703 }
704
705 BFT_FREE(m);
706 BFT_FREE(vtx_rotor_num);
707 }
708
709 /*----------------------------------------------------------------------------
710 * Check that rotors and stators are originally disjoint
711 *
712 * parameters:
713 * mesh <-- mesh to check
714 *----------------------------------------------------------------------------*/
715
716 static void
_check_geometry(cs_mesh_t * mesh)717 _check_geometry(cs_mesh_t *mesh)
718 {
719 cs_turbomachinery_t *tbm = _turbomachinery;
720
721 cs_gnum_t n_errors = 0;
722
723 const int *cell_flag = tbm->cell_rotor_num;
724
725 /* Mark from interior faces */
726
727 for (cs_lnum_t f_id = 0; f_id < mesh->n_i_faces; f_id++) {
728 cs_lnum_t c_id_0 = mesh->i_face_cells[f_id][0];
729 cs_lnum_t c_id_1 = mesh->i_face_cells[f_id][1];
730 if (cell_flag[c_id_1] - cell_flag[c_id_0] != 0)
731 n_errors ++;
732 }
733
734 cs_parall_counter(&n_errors, 1);
735 if (n_errors > 0)
736 bft_error
737 (__FILE__, __LINE__, 0,
738 _("%s: some faces of the initial mesh belong to different\n"
739 "rotor/stator sections.\n"
740 "These sections must be initially disjoint to rotate freely."),
741 __func__);
742 }
743
744 /*----------------------------------------------------------------------------
745 * Select rotor cells
746 *
747 * parameters:
748 * tbm <-> turbomachinery options structure
749 *----------------------------------------------------------------------------*/
750
751 static void
_select_rotor_cells(cs_turbomachinery_t * tbm)752 _select_rotor_cells(cs_turbomachinery_t *tbm)
753 {
754 cs_lnum_t _n_cells = 0;
755 cs_lnum_t *_cell_list = NULL;
756
757 cs_mesh_t *m = cs_glob_mesh;
758
759 assert(tbm->rotor_cells_c != NULL);
760
761 BFT_REALLOC(tbm->cell_rotor_num, m->n_cells_with_ghosts, int);
762
763 for (cs_lnum_t i = 0; i < m->n_cells_with_ghosts; i++)
764 tbm->cell_rotor_num[i] = 0;
765
766 BFT_MALLOC(_cell_list, m->n_cells_with_ghosts, cs_lnum_t);
767
768 for (int r_id = 0; r_id < tbm->n_rotors; r_id++) {
769 cs_selector_get_cell_list(tbm->rotor_cells_c[r_id],
770 &_n_cells,
771 _cell_list);
772 cs_gnum_t _n_g_cells = _n_cells;
773 cs_parall_counter(&_n_g_cells, 1);
774 if (_n_g_cells == 0)
775 bft_error
776 (__FILE__, __LINE__, 0,
777 _("%sd: The rotor %d with cell selection criteria\n"
778 " \"%s\"\n"
779 "does not contain any cell.\n"
780 "This rotor should be removed or its selection criteria modified."),
781 __func__, r_id + 1, tbm->rotor_cells_c[r_id]);
782 for (cs_lnum_t i = 0; i < _n_cells; i++)
783 tbm->cell_rotor_num[_cell_list[i]] = r_id + 1;
784 }
785
786 BFT_FREE(_cell_list);
787
788 if (m->halo != NULL)
789 cs_halo_sync_untyped(m->halo,
790 CS_HALO_EXTENDED,
791 sizeof(int),
792 tbm->cell_rotor_num);
793
794 if (tbm->model == CS_TURBOMACHINERY_TRANSIENT)
795 _check_geometry(m);
796 }
797
798 /*----------------------------------------------------------------------------
799 * Update mesh for unsteady rotor/stator computation when no joining is used.
800 *
801 * parameters:
802 * restart_mode true for restart, false otherwise
803 * t_cur_mob current rotor time
804 * t_elapsed elapsed computation time
805 */
806 /*----------------------------------------------------------------------------*/
807
808 static void
_update_mesh_coupling(double t_cur_mob,double * t_elapsed)809 _update_mesh_coupling(double t_cur_mob,
810 double *t_elapsed)
811 {
812 double t_start, t_end;
813
814 cs_turbomachinery_t *tbm = _turbomachinery;
815
816 int t_stat_id = cs_timer_stats_id_by_name("mesh_processing");
817 int t_top_id = cs_timer_stats_switch(t_stat_id);
818
819 t_start = cs_timer_wtime();
820
821 /* Indicates we are in the framework of turbomachinery */
822
823 tbm->active = true;
824
825 /* Cell and boundary face numberings can be moved from old mesh
826 to new one, as the corresponding parts of the mesh should not change */
827
828 _copy_mesh(tbm->reference_mesh, cs_glob_mesh);
829
830 /* Update geometry, if necessary */
831
832 _update_angle(t_cur_mob);
833
834 if (tbm->n_rotors > 0)
835 _update_geometry(cs_glob_mesh, 0);
836
837 /* Recompute geometric quantities related to the mesh */
838
839 cs_mesh_quantities_compute(cs_glob_mesh, cs_glob_mesh_quantities);
840
841 /* Update linear algebra APIs relative to mesh */
842
843 t_end = cs_timer_wtime();
844
845 *t_elapsed = t_end - t_start;
846
847 cs_timer_stats_switch(t_top_id);
848 }
849
850 /*----------------------------------------------------------------------------
851 * Update mesh for unsteady rotor/stator computation.
852 *
853 * parameters:
854 * restart_mode true for restart, false otherwise
855 * t_cur_mob current rotor time
856 * t_elapsed elapsed computation time
857 */
858 /*----------------------------------------------------------------------------*/
859
860 static void
_update_mesh(bool restart_mode,double t_cur_mob,double * t_elapsed)861 _update_mesh(bool restart_mode,
862 double t_cur_mob,
863 double *t_elapsed)
864 {
865 double t_start, t_end;
866
867 cs_halo_type_t halo_type = cs_glob_mesh->halo_type;
868 cs_turbomachinery_t *tbm = _turbomachinery;
869
870 int t_stat_id = cs_timer_stats_id_by_name("mesh_processing");
871 int t_top_id = cs_timer_stats_switch(t_stat_id);
872
873 t_start = cs_timer_wtime();
874
875 /* Indicates we are in the framework of turbomachinery */
876
877 tbm->active = true;
878
879 /* In case of simple coupling, use simpler update */
880
881 if (cs_glob_n_joinings < 1) {
882 _update_mesh_coupling(t_cur_mob,
883 t_elapsed);
884 return;
885 }
886
887 /* Cell and boundary face numberings can be moved from old mesh
888 to new one, as the corresponding parts of the mesh should not change */
889
890 cs_numbering_t *cell_numbering = NULL;
891
892 if (restart_mode == false) {
893 cell_numbering = cs_glob_mesh->cell_numbering;
894 cs_glob_mesh->cell_numbering = NULL;
895 }
896
897 /* Destroy previous global mesh and related entities */
898
899 cs_mesh_quantities_destroy(cs_glob_mesh_quantities);
900
901 cs_mesh_destroy(cs_glob_mesh);
902
903 /* Create new global mesh and related entities */
904
905 cs_glob_mesh = cs_mesh_create();
906 cs_glob_mesh->verbosity = 0;
907 cs_glob_mesh_builder = cs_mesh_builder_create();
908 cs_glob_mesh_quantities = cs_mesh_quantities_create();
909
910 _update_angle(t_cur_mob);
911
912 if (restart_mode == false) {
913
914 int n_retry = CS_MAX(tbm->n_max_join_tries, 1);
915 cs_lnum_t boundary_changed = 0;
916 double eps_dt = 0.;
917
918 const cs_time_step_t *ts = cs_glob_time_step;
919
920 do {
921
922 n_retry -= 1;
923
924 _copy_mesh(tbm->reference_mesh, cs_glob_mesh);
925
926 /* Update geometry, if necessary */
927
928 if (tbm->n_rotors > 0)
929 _update_geometry(cs_glob_mesh, eps_dt);
930
931 /* Reset the interior faces -> cells connectivity */
932 /* (in order to properly build the halo of the joined mesh) */
933
934 cs_mesh_to_builder_perio_faces(cs_glob_mesh, cs_glob_mesh_builder);
935
936 {
937 int i;
938 cs_lnum_t f_id;
939 cs_lnum_2_t *i_face_cells = (cs_lnum_2_t *)cs_glob_mesh->i_face_cells;
940 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
941 for (f_id = 0; f_id < cs_glob_mesh->n_i_faces; f_id++) {
942 for (i = 0; i < 2; i++) {
943 if (i_face_cells[f_id][i] >= n_cells)
944 i_face_cells[f_id][i] = -1;
945 }
946 }
947 }
948
949 /* Join meshes and build periodicity links */
950
951 cs_join_all(false);
952
953 boundary_changed = 0;
954 if (tbm->n_b_faces_ref > -1) {
955 if (cs_glob_mesh->n_b_faces != tbm->n_b_faces_ref)
956 boundary_changed = 1;
957 }
958 cs_parall_counter_max(&boundary_changed, 1);
959
960 /* Check that joining has not added or removed boundary faces.
961 Postprocess new faces appearing on boundary or inside of mesh
962 (which assumes that joining appends new faces at the end of the mesh)
963 or try again with slightly different rotation angle
964 */
965
966 if (boundary_changed) {
967
968 const char join_err_fmt[]
969 = N_("Error in turbomachinery mesh update:\n"
970 "Number of boundary faces has changed from %llu to %llu.\n"
971 "There are probably unjoined faces, "
972 "due to an insufficiently regular mesh;\n"
973 "adjusting mesh joining parameters might help.");
974
975 cs_gnum_t n_g_b_faces_ref = tbm->n_b_faces_ref;
976 cs_parall_counter(&n_g_b_faces_ref, 1);
977
978 if (n_retry < 1) {
979 const int writer_id = -2;
980 const int writer_ids[] = {writer_id};
981 const int mesh_id = cs_post_get_free_mesh_id();
982 cs_lnum_t b_face_count[] = {tbm->n_b_faces_ref,
983 cs_glob_mesh->n_b_faces};
984 cs_post_init_error_writer();
985 cs_post_define_surface_mesh_by_func(mesh_id,
986 _("Added boundary faces"),
987 NULL,
988 _post_error_faces_select,
989 NULL,
990 b_face_count,
991 false, /* time varying */
992 true, /* add groups if present */
993 false, /* auto variables */
994 1,
995 writer_ids);
996 cs_post_activate_writer(writer_id, 1);
997 cs_post_write_meshes(NULL);
998 bft_error(__FILE__, __LINE__, 0,
999 _(join_err_fmt),
1000 (unsigned long long)n_g_b_faces_ref,
1001 (unsigned long long)cs_glob_mesh->n_g_b_faces);
1002 }
1003 else {
1004 eps_dt += tbm->dt_retry * ts->dt[0];
1005 bft_printf(_(join_err_fmt),
1006 (unsigned long long)n_g_b_faces_ref,
1007 (unsigned long long)cs_glob_mesh->n_g_b_faces);
1008 bft_printf("\nTrying again with eps_dt = %lg\n", eps_dt);
1009
1010 /* Destroy previous global mesh and related entities */
1011
1012 cs_mesh_quantities_destroy(cs_glob_mesh_quantities);
1013
1014 cs_mesh_destroy(cs_glob_mesh);
1015
1016 /* Create new global mesh and related entities */
1017
1018 cs_glob_mesh = cs_mesh_create();
1019 cs_glob_mesh->verbosity = 0;
1020 cs_glob_mesh_builder = cs_mesh_builder_create();
1021 cs_glob_mesh_quantities = cs_mesh_quantities_create();
1022 }
1023 }
1024 } while (boundary_changed && n_retry >= 0);
1025 }
1026 else {
1027
1028 cs_mesh_to_builder_partition(tbm->reference_mesh,
1029 cs_glob_mesh_builder);
1030
1031 cs_preprocessor_data_add_file("restart/mesh", 0, NULL, NULL);
1032
1033 cs_preprocessor_data_read_headers(cs_glob_mesh,
1034 cs_glob_mesh_builder);
1035
1036 if (tbm->reference_mesh->n_g_cells != cs_glob_mesh->n_g_cells)
1037 bft_error
1038 (__FILE__, __LINE__, 0,
1039 _("Error in turbomachinery mesh restart:\n"
1040 " number of cells expected/read: %llu/%llu\n"
1041 "Check your restart directory."),
1042 (unsigned long long)tbm->reference_mesh->n_g_cells,
1043 (unsigned long long)cs_glob_mesh->n_g_cells);
1044
1045 cs_preprocessor_data_read_mesh(cs_glob_mesh,
1046 cs_glob_mesh_builder);
1047 }
1048
1049 tbm->n_b_faces_ref = cs_glob_mesh->n_b_faces;
1050
1051 /* Initialize extended connectivity, ghost cells and other remaining
1052 parallelism-related structures */
1053
1054 cs_mesh_init_halo(cs_glob_mesh, cs_glob_mesh_builder, halo_type,
1055 cs_glob_mesh->verbosity, true);
1056 cs_mesh_update_auxiliary(cs_glob_mesh);
1057
1058 /* Destroy the temporary structure used to build the main mesh */
1059
1060 cs_mesh_builder_destroy(&cs_glob_mesh_builder);
1061
1062 /* Update numberings (cells saved from previous, faces
1063 faces updated; it is important that boundary faces
1064 renumbering produce the same result at each iteration) */
1065
1066 if (restart_mode)
1067 cs_renumber_cells(cs_glob_mesh);
1068 else
1069 cs_glob_mesh->cell_numbering = cell_numbering;
1070
1071 cs_renumber_i_faces(cs_glob_mesh);
1072 cs_renumber_b_faces(cs_glob_mesh);
1073
1074 /* Build group classes */
1075
1076 cs_mesh_init_group_classes(cs_glob_mesh);
1077
1078 /* Print info on mesh */
1079
1080 if (cs_glob_mesh->verbosity > 0)
1081 cs_mesh_print_info(cs_glob_mesh, _("Mesh"));
1082
1083 /* Compute geometric quantities related to the mesh */
1084
1085 cs_mesh_quantities_compute(cs_glob_mesh, cs_glob_mesh_quantities);
1086 cs_mesh_bad_cells_detect(cs_glob_mesh, cs_glob_mesh_quantities);
1087 cs_user_mesh_bad_cells_tag(cs_glob_mesh, cs_glob_mesh_quantities);
1088
1089 cs_ext_neighborhood_reduce(cs_glob_mesh,
1090 cs_glob_mesh_quantities);
1091
1092 /* Initialize selectors and locations for the mesh */
1093
1094 cs_mesh_init_selectors();
1095 cs_mesh_location_build(cs_glob_mesh, -1);
1096 cs_volume_zone_build_all(true);
1097 cs_boundary_zone_build_all(true);
1098
1099 /* Check coherency if debugging */
1100
1101 #if 0
1102 cs_mesh_coherency_check();
1103 #endif
1104
1105 /* Update Fortran mesh sizes and quantities */
1106
1107 cs_preprocess_mesh_update_fortran();
1108
1109 /* Update mapping for accelerated devices */
1110
1111 #if defined(HAVE_ACCEL)
1112 cs_preprocess_mesh_update_device(cs_alloc_mode);
1113 #endif
1114
1115 /* Update rotor cells flag array in case of parallelism and/or periodicity */
1116
1117 if (cs_glob_mesh->halo != NULL) {
1118
1119 const cs_mesh_t *m = cs_glob_mesh;
1120 BFT_REALLOC(tbm->cell_rotor_num,
1121 m->n_cells_with_ghosts,
1122 int);
1123
1124 cs_halo_sync_untyped(m->halo,
1125 CS_HALO_EXTENDED,
1126 sizeof(int),
1127 tbm->cell_rotor_num);
1128
1129 }
1130
1131 cs_gradient_free_quantities();
1132 cs_cell_to_vertex_free();
1133 cs_mesh_adjacencies_update_mesh();
1134
1135 /* Update linear algebra APIs relative to mesh */
1136
1137 cs_matrix_update_mesh();
1138
1139 t_end = cs_timer_wtime();
1140
1141 *t_elapsed = t_end - t_start;
1142
1143 cs_timer_stats_switch(t_top_id);
1144 }
1145
1146 /*============================================================================
1147 * Fortran wrapper function definitions
1148 *============================================================================*/
1149
1150 /*----------------------------------------------------------------------------
1151 * Assist map turbomachinery to fortran module
1152 *
1153 * parameters:
1154 * iturbo <-- turbomachinery type flag
1155 * ityint <-- 0: joining; 1: internal coupling
1156 *----------------------------------------------------------------------------*/
1157
1158 void
cs_f_map_turbomachinery_model(int * iturbo,int * ityint)1159 cs_f_map_turbomachinery_model(int *iturbo,
1160 int *ityint)
1161 {
1162 if (_turbomachinery != NULL)
1163 *iturbo = _turbomachinery->model;
1164 else
1165 *iturbo = CS_TURBOMACHINERY_NONE;
1166
1167 if (_turbomachinery->n_couplings > 0)
1168 *ityint = 1;
1169 else
1170 *ityint = 0;
1171 }
1172
1173 /*----------------------------------------------------------------------------
1174 * Map turbomachinery rotor to fortran module
1175 *
1176 * parameters:
1177 * irotce <-- pointer to cell flag for rotation
1178 *----------------------------------------------------------------------------*/
1179
1180 void
cs_f_map_turbomachinery_rotor(int ** irotce)1181 cs_f_map_turbomachinery_rotor(int **irotce)
1182 {
1183 if (_turbomachinery != NULL)
1184 *irotce = _turbomachinery->cell_rotor_num;
1185 else
1186 *irotce = NULL;
1187 }
1188
1189 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1190
1191 /*============================================================================
1192 * Public function definitions
1193 *============================================================================*/
1194
1195 /*----------------------------------------------------------------------------*/
1196 /*!
1197 * \brief Define rotor/stator model.
1198 */
1199 /*----------------------------------------------------------------------------*/
1200
1201 void
cs_turbomachinery_set_model(cs_turbomachinery_model_t model)1202 cs_turbomachinery_set_model(cs_turbomachinery_model_t model)
1203 {
1204 if ( model == CS_TURBOMACHINERY_NONE
1205 && _turbomachinery != NULL) {
1206 cs_turbomachinery_finalize();
1207 return;
1208 }
1209
1210 else if (_turbomachinery == NULL)
1211 _turbomachinery = _turbomachinery_create();
1212
1213 _turbomachinery->model = model;
1214 }
1215
1216 /*----------------------------------------------------------------------------*/
1217 /*!
1218 * \brief Return rotor/stator model.
1219 */
1220 /*----------------------------------------------------------------------------*/
1221
1222 cs_turbomachinery_model_t
cs_turbomachinery_get_model(void)1223 cs_turbomachinery_get_model(void)
1224 {
1225 if (_turbomachinery == NULL)
1226 return CS_TURBOMACHINERY_NONE;
1227 else
1228 return _turbomachinery->model;
1229 }
1230
1231 /*----------------------------------------------------------------------------*/
1232 /*!
1233 * \brief Define a rotor by its axis and cell selection criteria.
1234 *
1235 * \param[in] cell_criteria cell selection criteria string
1236 * \param[in] rotation_velocity rotation velocity, in radians/second
1237 * \param[in] rotation_axis rotation axis vector
1238 * \param[in] rotation_invariant rotation invariant point
1239 */
1240 /*----------------------------------------------------------------------------*/
1241
1242 void
cs_turbomachinery_add_rotor(const char * cell_criteria,double rotation_velocity,const double rotation_axis[3],const double rotation_invariant[3])1243 cs_turbomachinery_add_rotor(const char *cell_criteria,
1244 double rotation_velocity,
1245 const double rotation_axis[3],
1246 const double rotation_invariant[3])
1247 {
1248 cs_turbomachinery_t *tbm = _turbomachinery;
1249 if (tbm == NULL)
1250 return;
1251
1252 const double *v = rotation_axis;
1253 double len = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
1254
1255 int r_id = tbm->n_rotors;
1256 tbm->n_rotors += 1;
1257
1258 BFT_REALLOC(tbm->rotation, tbm->n_rotors + 1, cs_rotation_t);
1259 cs_rotation_t *r = tbm->rotation + r_id + 1;
1260 r->omega = rotation_velocity;
1261 r->angle = 0;
1262 for (int i = 0; i < 3; i++) {
1263 r->axis[i] = v[i]/len;
1264 r->invariant[i] = rotation_invariant[i];
1265 }
1266
1267 BFT_REALLOC(tbm->rotor_cells_c, tbm->n_rotors, char *);
1268 BFT_MALLOC(tbm->rotor_cells_c[r_id], strlen(cell_criteria) + 1, char);
1269 strcpy(tbm->rotor_cells_c[r_id], cell_criteria);
1270 }
1271
1272 /*----------------------------------------------------------------------------*/
1273 /*!
1274 * \brief Add a cs_join_t structure to the list of rotor/stator joinings.
1275 *
1276 * \param[in] sel_criteria boundary face selection criteria
1277 * \param[in] fraction value of the fraction parameter
1278 * \param[in] plane value of the plane parameter
1279 * \param[in] verbosity level of verbosity required
1280 * \param[in] visualization level of visualization required
1281 *
1282 * \return number (1 to n) associated with new joining
1283 */
1284 /*----------------------------------------------------------------------------*/
1285
1286 int
cs_turbomachinery_join_add(const char * sel_criteria,float fraction,float plane,int verbosity,int visualization)1287 cs_turbomachinery_join_add(const char *sel_criteria,
1288 float fraction,
1289 float plane,
1290 int verbosity,
1291 int visualization)
1292 {
1293 /* Allocate and initialize a cs_join_t structure */
1294
1295 BFT_REALLOC(cs_glob_join_array, cs_glob_n_joinings + 1, cs_join_t *);
1296
1297 cs_glob_join_array[cs_glob_n_joinings]
1298 = cs_join_create(cs_glob_n_joinings + 1,
1299 sel_criteria,
1300 fraction,
1301 plane,
1302 FVM_PERIODICITY_NULL,
1303 NULL,
1304 verbosity,
1305 visualization,
1306 false);
1307
1308 cs_glob_join_count++; /* Store number of joining (without periodic ones) */
1309 cs_glob_n_joinings++;
1310
1311 /* Set mesh modification type */
1312
1313 if (_turbomachinery != NULL) {
1314 if ( _turbomachinery->model == CS_TURBOMACHINERY_TRANSIENT
1315 && cs_glob_mesh->time_dep < CS_MESH_TRANSIENT_CONNECT)
1316 cs_glob_mesh->time_dep = CS_MESH_TRANSIENT_CONNECT;
1317 }
1318
1319 return cs_glob_n_joinings;
1320 }
1321
1322 /*----------------------------------------------------------------------------*/
1323 /*!
1324 * \brief Add a boundary coupling to the list of rotor/stator couplings.
1325 *
1326 * \param[in] sel_criteria boundary face selection criteria
1327 * \param[in] tolerance value of the search tolerance
1328 * \param[in] verbosity level of verbosity required
1329 *
1330 * \return number (1 to n) associated with new coupling
1331 */
1332 /*----------------------------------------------------------------------------*/
1333
1334 int
cs_turbomachinery_coupling_add(const char * sel_criteria,float tolerance,int verbosity)1335 cs_turbomachinery_coupling_add(const char *sel_criteria,
1336 float tolerance,
1337 int verbosity)
1338 {
1339 cs_sat_coupling_add_internal(_turbomachinery_coupling_tag,
1340 _turbomachinery,
1341 sel_criteria,
1342 NULL,
1343 NULL,
1344 "all[]",
1345 tolerance,
1346 verbosity);
1347
1348 _turbomachinery->n_couplings += 1;
1349
1350 if (_turbomachinery != NULL) {
1351 if ( _turbomachinery->model == CS_TURBOMACHINERY_TRANSIENT
1352 && cs_glob_mesh->time_dep < CS_MESH_TRANSIENT_COORDS)
1353 cs_glob_mesh->time_dep = CS_MESH_TRANSIENT_COORDS;
1354 }
1355
1356 return cs_sat_coupling_n_couplings();
1357 }
1358
1359 /*----------------------------------------------------------------------------*/
1360 /*!
1361 * \brief Update mesh for unsteady rotor/stator computation.
1362 *
1363 * \param[in] t_cur_mob current rotor time
1364 * \param[out] t_elapsed elapsed computation time
1365 */
1366 /*----------------------------------------------------------------------------*/
1367
1368 void
cs_turbomachinery_update_mesh(double t_cur_mob,double * t_elapsed)1369 cs_turbomachinery_update_mesh(double t_cur_mob,
1370 double *t_elapsed)
1371 {
1372 _update_mesh(false, t_cur_mob, t_elapsed);
1373 }
1374
1375 /*----------------------------------------------------------------------------*/
1376 /*!
1377 * \brief Update mesh for unsteady rotor/stator computation in case of restart.
1378 *
1379 * Reads mesh from checkpoint when available.
1380 */
1381 /*----------------------------------------------------------------------------*/
1382
1383 void
cs_turbomachinery_restart_mesh(void)1384 cs_turbomachinery_restart_mesh(void)
1385 {
1386 if (cs_turbomachinery_get_model() != CS_TURBOMACHINERY_TRANSIENT)
1387 return;
1388
1389 if (cs_glob_time_step->nt_prev > 0) {
1390 double t_elapsed;
1391 if (cs_file_isreg("restart/mesh"))
1392 _update_mesh(true, cs_glob_time_step->t_cur, &t_elapsed);
1393 else
1394 _update_mesh(false, cs_glob_time_step->t_cur, &t_elapsed);
1395 }
1396 }
1397
1398 /*----------------------------------------------------------------------------*/
1399 /*!
1400 * \brief Definitions for turbomachinery computation.
1401 */
1402 /*----------------------------------------------------------------------------*/
1403
1404 void
cs_turbomachinery_define(void)1405 cs_turbomachinery_define(void)
1406 {
1407 /* Define model; could be moved anywhere before time loop. */
1408
1409 cs_gui_turbomachinery();
1410 cs_user_turbomachinery();
1411
1412 if (_turbomachinery == NULL)
1413 return;
1414
1415 cs_turbomachinery_t *tbm = _turbomachinery;
1416
1417 if (tbm->model == CS_TURBOMACHINERY_NONE)
1418 return;
1419
1420 /* Define rotor(s) */
1421
1422 cs_gui_turbomachinery_rotor();
1423 cs_user_turbomachinery_rotor();
1424 }
1425
1426 /*----------------------------------------------------------------------------*/
1427 /*!
1428 * \brief Initializations for turbomachinery computation.
1429 *
1430 * \note This function should be called after the mesh is built,
1431 * but before cs_post_init_meshes() so that postprocessing meshes are
1432 * updated correctly in the transient case.
1433 */
1434 /*----------------------------------------------------------------------------*/
1435
1436 void
cs_turbomachinery_initialize(void)1437 cs_turbomachinery_initialize(void)
1438 {
1439 if (_turbomachinery == NULL)
1440 return;
1441
1442 cs_turbomachinery_t *tbm = _turbomachinery;
1443
1444 if (tbm->model == CS_TURBOMACHINERY_NONE)
1445 return;
1446
1447 /* Select rotor cells */
1448
1449 _select_rotor_cells(tbm);
1450
1451 /* Build the reference mesh that duplicates the global mesh before joining;
1452 first remove the boundary face numbering, as it will need to be
1453 rebuilt after the first joining */
1454
1455 if (cs_glob_mesh->b_face_numbering != NULL && cs_glob_n_joinings > 0)
1456 cs_numbering_destroy(&(cs_glob_mesh->b_face_numbering));
1457
1458 _copy_mesh(cs_glob_mesh, tbm->reference_mesh);
1459
1460 /* Reorder reference mesh by global number to avoid some issues with
1461 joining, especially in serial mode where global numbers are not
1462 expected to be present at joining stages */
1463
1464 cs_renumber_i_faces_by_gnum(tbm->reference_mesh);
1465 cs_renumber_b_faces_by_gnum(tbm->reference_mesh);
1466
1467 /* Complete the mesh with rotor-stator joining */
1468
1469 if (cs_glob_n_joinings > 0) {
1470 cs_real_t t_elapsed;
1471 cs_turbomachinery_update_mesh(0.0, &t_elapsed);
1472 }
1473
1474 /* Adapt postprocessing options if required;
1475 must be called before cs_post_init_meshes(). */
1476
1477 if (tbm->model == CS_TURBOMACHINERY_TRANSIENT)
1478 cs_post_set_changing_connectivity();
1479
1480 /* Destroy the reference mesh, if required */
1481
1482 if (tbm->model == CS_TURBOMACHINERY_FROZEN) {
1483 cs_mesh_destroy(tbm->reference_mesh);
1484 tbm->reference_mesh = NULL;
1485 }
1486
1487 /* Set global rotations pointer */
1488
1489 cs_glob_rotation = tbm->rotation;
1490 }
1491
1492 /*----------------------------------------------------------------------------*/
1493 /*!
1494 * \brief Free turbomachinery structure.
1495 */
1496 /*----------------------------------------------------------------------------*/
1497
1498 void
cs_turbomachinery_finalize(void)1499 cs_turbomachinery_finalize(void)
1500 {
1501 if (_turbomachinery != NULL) {
1502
1503 cs_turbomachinery_t *tbm = _turbomachinery;
1504
1505 for (int i = tbm->n_rotors -1; i >= 0; i--)
1506 BFT_FREE(tbm->rotor_cells_c[i]);
1507 BFT_FREE(tbm->rotor_cells_c);
1508
1509 BFT_FREE(tbm->rotation);
1510
1511 BFT_FREE(tbm->cell_rotor_num);
1512
1513 if (tbm->reference_mesh != NULL)
1514 cs_mesh_destroy(tbm->reference_mesh);
1515
1516 /* Unset global rotations pointer for safety */
1517 cs_glob_rotation = NULL;
1518 }
1519
1520 BFT_FREE(_turbomachinery);
1521 }
1522
1523 /*----------------------------------------------------------------------------*/
1524 /*!
1525 * \brief Reinitialize interior face-based fields.
1526 */
1527 /*----------------------------------------------------------------------------*/
1528
1529 void
cs_turbomachinery_reinit_i_face_fields(void)1530 cs_turbomachinery_reinit_i_face_fields(void)
1531 {
1532 const int n_fields = cs_field_n_fields();
1533
1534 for (int f_id = 0; f_id < n_fields; f_id++) {
1535
1536 cs_field_t *f = cs_field_by_id(f_id);
1537
1538 if ( cs_mesh_location_get_type(f->location_id)
1539 == CS_MESH_LOCATION_INTERIOR_FACES)
1540 cs_field_allocate_values(f);
1541 }
1542 }
1543
1544 /*----------------------------------------------------------------------------*/
1545 /*!
1546 * \brief Resize cell-based fields.
1547 *
1548 * This function only handles fields owning their values.
1549 */
1550 /*----------------------------------------------------------------------------*/
1551
1552 void
cs_turbomachinery_resize_cell_fields(void)1553 cs_turbomachinery_resize_cell_fields(void)
1554 {
1555 const int n_fields = cs_field_n_fields();
1556
1557 const cs_halo_t *halo = cs_glob_mesh->halo;
1558 const cs_lnum_t *n_elts = cs_mesh_location_get_n_elts(CS_MESH_LOCATION_CELLS);
1559 cs_lnum_t _n_cells = n_elts[2];
1560
1561 for (int f_id = 0; f_id < n_fields; f_id++) {
1562
1563 cs_field_t *f = cs_field_by_id(f_id);
1564
1565 if ( f->location_id == CS_MESH_LOCATION_CELLS
1566 && f->is_owner) {
1567
1568 /* Ghost cell sizes may change, but the number of main cells
1569 is unchanged, so a simple reallocation will do */
1570
1571 for (int kk = 0; kk < f->n_time_vals; kk++) {
1572
1573 BFT_REALLOC(f->vals[kk], _n_cells*f->dim, cs_real_t);
1574
1575 if (halo != NULL) {
1576
1577 cs_halo_sync_untyped(halo,
1578 CS_HALO_EXTENDED,
1579 f->dim*sizeof(cs_real_t),
1580 f->vals[kk]);
1581 if (f->dim == 3)
1582 cs_halo_perio_sync_var_vect(halo,
1583 CS_HALO_EXTENDED,
1584 f->vals[kk],
1585 f->dim);
1586 }
1587 }
1588
1589 f->val = f->vals[0];
1590 if (f->n_time_vals > 1) f->val_pre = f->vals[1];
1591
1592 }
1593
1594 }
1595
1596 }
1597
1598 /*----------------------------------------------------------------------------*/
1599 /*!
1600 * \brief Compute rotation matrix
1601 *
1602 * \param[in] rotor_num rotor number (1 to n numbering)
1603 * \param[in] theta rotation angle, in radians
1604 * \param[out] matrix resulting rotation matrix
1605 */
1606 /*----------------------------------------------------------------------------*/
1607
1608 void
cs_turbomachinery_rotation_matrix(int rotor_num,double theta,cs_real_t matrix[3][4])1609 cs_turbomachinery_rotation_matrix(int rotor_num,
1610 double theta,
1611 cs_real_t matrix[3][4])
1612 {
1613 cs_turbomachinery_t *tbm = _turbomachinery;
1614 const cs_rotation_t *r = tbm->rotation + rotor_num;
1615
1616 cs_rotation_matrix(theta, r->axis, r->invariant, matrix);
1617 }
1618
1619 /*----------------------------------------------------------------------------*/
1620 /*!
1621 * \brief Return number of rotors.
1622 *
1623 * Note that the number of associated rotations is n_rotors + 1, as the
1624 * first rotation id is reserved for the fixed portion of the domain.
1625 *
1626 * \return number of rotors
1627 */
1628 /*----------------------------------------------------------------------------*/
1629
1630 int
cs_turbomachinery_n_rotors(void)1631 cs_turbomachinery_n_rotors(void)
1632 {
1633 int retval = 1;
1634
1635 if (_turbomachinery != NULL)
1636 retval = _turbomachinery->n_rotors;
1637
1638 return retval;
1639 }
1640
1641 /*----------------------------------------------------------------------------*/
1642 /*!
1643 * \brief Return cell rotor number.
1644 *
1645 * Each cell may be associated with a given rotor, or rotation, with 0
1646 * indicating that that cell does not rotate.
1647 *
1648 * \return array defining rotor number associated with each cell
1649 * (0 for none, 1 to n otherwise)
1650 */
1651 /*----------------------------------------------------------------------------*/
1652
1653 const int *
cs_turbomachinery_get_cell_rotor_num(void)1654 cs_turbomachinery_get_cell_rotor_num(void)
1655 {
1656 return _turbomachinery->cell_rotor_num;
1657 }
1658
1659 /*----------------------------------------------------------------------------*/
1660 /*!
1661 * \brief Return rotation velocity
1662 *
1663 * \param[in] rotor_num rotor number (1 to n numbering)
1664 */
1665 /*----------------------------------------------------------------------------*/
1666
1667 double
cs_turbomachinery_get_rotation_velocity(int rotor_num)1668 cs_turbomachinery_get_rotation_velocity(int rotor_num)
1669 {
1670 cs_turbomachinery_t *tbm = _turbomachinery;
1671
1672 return (tbm->rotation + rotor_num)->omega;
1673 }
1674
1675 /*----------------------------------------------------------------------------*/
1676 /*!
1677 * \brief Set rotation velocity
1678 *
1679 * param[in] rotor_num rotor number (1 to n numbering)
1680 * param[in] omega rotation velocity
1681 */
1682 /*----------------------------------------------------------------------------*/
1683
1684 void
cs_turbomachinery_set_rotation_velocity(int rotor_num,double omega)1685 cs_turbomachinery_set_rotation_velocity(int rotor_num,
1686 double omega)
1687 {
1688 cs_turbomachinery_t *tbm = _turbomachinery;
1689
1690 (tbm->rotation + rotor_num)->omega = omega;
1691 }
1692
1693 /*----------------------------------------------------------------------------*/
1694 /*!
1695 * \brief Set turbomachinery joining retry parameters.
1696 *
1697 * When a joing leads to a different number of boundary faces from the
1698 * previous position, the rotor positions may be perturbed by a small
1699 * quantity to try to obtain a better joining.
1700 *
1701 * param[in] n_max_join_retries maximum number of retries before considering
1702 * the joining has failed
1703 * param[in] dt_retry_multiplier time step multiplier for new position retry
1704 */
1705 /*----------------------------------------------------------------------------*/
1706
1707 void
cs_turbomachinery_set_rotation_retry(int n_max_join_retries,double dt_retry_multiplier)1708 cs_turbomachinery_set_rotation_retry(int n_max_join_retries,
1709 double dt_retry_multiplier)
1710 {
1711 cs_turbomachinery_t *tbm = _turbomachinery;
1712
1713 tbm->n_max_join_tries = n_max_join_retries;
1714 tbm->dt_retry = dt_retry_multiplier;
1715 }
1716
1717 /*----------------------------------------------------------------------------*/
1718 /*!
1719 * \brief Build rotation matrices for a given time interval.
1720 *
1721 * The caller is responsible for freeing the array when not needed.
1722 *
1723 * \param[in] dt associated time delta (0 for current, unmodified time)
1724 *
1725 * \return array of rotation matrices.
1726 */
1727 /*----------------------------------------------------------------------------*/
1728
1729 cs_real_34_t *
cs_turbomachinery_get_rotation_matrices(double dt)1730 cs_turbomachinery_get_rotation_matrices(double dt)
1731 {
1732 const cs_turbomachinery_t *tbm = _turbomachinery;
1733 cs_real_34_t *m;
1734
1735 BFT_MALLOC(m, tbm->n_rotors+1, cs_real_34_t);
1736
1737 for (int j = 0; j < tbm->n_rotors+1; j++) {
1738 cs_rotation_t *r = tbm->rotation + j;
1739
1740 cs_rotation_matrix(r->omega*dt,
1741 r->axis,
1742 r->invariant,
1743 m[j]);
1744 }
1745
1746 return m;
1747 }
1748
1749 /*----------------------------------------------------------------------------*/
1750 /*!
1751 * \brief Rotation of vector and tensor fields.
1752 */
1753 /*----------------------------------------------------------------------------*/
1754
1755 void
cs_turbomachinery_rotate_fields(const cs_real_t dt[])1756 cs_turbomachinery_rotate_fields(const cs_real_t dt[])
1757 {
1758 cs_lnum_t i;
1759 cs_real_34_t *m;
1760
1761 cs_turbomachinery_t *tbm = _turbomachinery;
1762 cs_real_t time_step = dt[0];
1763
1764 BFT_MALLOC(m, tbm->n_rotors+1, cs_real_34_t);
1765
1766 for (int j = 0; j < tbm->n_rotors+1; j++) {
1767 cs_rotation_t *r = tbm->rotation + j;
1768 cs_rotation_matrix(r->omega * time_step,
1769 r->axis,
1770 r->invariant,
1771 m[j]);
1772 }
1773
1774 int n_fields = cs_field_n_fields();
1775
1776 for (int f_id = 0; f_id < n_fields; f_id++) {
1777
1778 cs_field_t *f = cs_field_by_id(f_id);
1779
1780 if (! (f->dim > 1 && (f->type & CS_FIELD_VARIABLE)))
1781 continue;
1782
1783 const cs_lnum_t *n_elts = cs_mesh_location_get_n_elts(f->location_id);
1784 cs_lnum_t _n_elts = n_elts[2];
1785
1786 if (f->dim == 3) {
1787 for (i = 0; i < _n_elts; i++)
1788 _apply_vector_rotation(m[tbm->cell_rotor_num[i]], f->val + i*3);
1789 }
1790
1791 else if (f->dim == 6) {
1792 for (i = 0; i < _n_elts; i++)
1793 _apply_sym_tensor_rotation(m[tbm->cell_rotor_num[i]], f->val + i*6);
1794 }
1795
1796 assert(f->dim == 3 || f->dim == 6);
1797
1798 } /* End of loop on fields */
1799
1800 /* Specific handling of Reynolds stresses */
1801
1802 cs_field_t *frij = cs_field_by_name("rij");
1803 if (frij != NULL) {
1804 const cs_lnum_t *n_elts = cs_mesh_location_get_n_elts(frij->location_id);
1805 cs_lnum_t _n_elts = n_elts[2];
1806 cs_real_6_t *cvar_r = (cs_real_6_t *)(frij->val);
1807 for (i = 0; i < _n_elts; i++) {
1808 _apply_sym_tensor_rotation(m[tbm->cell_rotor_num[i]], cvar_r[i]);
1809 }
1810
1811 }
1812
1813 BFT_FREE(m);
1814 }
1815
1816 /*----------------------------------------------------------------------------*/
1817 /*!
1818 * \brief Compute velocity relative to fixed coordinates at a given point
1819 *
1820 * \deprecated
1821 * Use \ref cs_rotation_velocity for more consistent naming of this reference
1822 * frame velocity.
1823 *
1824 * \param[in] rotor_num rotor number (1 to n numbering)
1825 * \param[in] coords point coordinates
1826 * \param[out] velocity velocity relative to fixed coordinates
1827 */
1828 /*----------------------------------------------------------------------------*/
1829
1830 void
cs_turbomachinery_relative_velocity(int rotor_num,const cs_real_t coords[3],cs_real_t velocity[3])1831 cs_turbomachinery_relative_velocity(int rotor_num,
1832 const cs_real_t coords[3],
1833 cs_real_t velocity[3])
1834 {
1835 cs_turbomachinery_t *tbm = _turbomachinery;
1836 const cs_rotation_t *r = tbm->rotation + rotor_num;
1837
1838 _relative_velocity(r->omega,
1839 r->axis,
1840 r->invariant,
1841 coords,
1842 velocity);
1843 }
1844
1845 /*----------------------------------------------------------------------------*/
1846 /*!
1847 * \brief Read turbomachinery metadata from restart file.
1848 *
1849 * The mesh is handled separately.
1850 *
1851 * \param[in, out] r associated restart file pointer
1852 */
1853 /*----------------------------------------------------------------------------*/
1854
1855 void
cs_turbomachinery_restart_read(cs_restart_t * r)1856 cs_turbomachinery_restart_read(cs_restart_t *r)
1857 {
1858 cs_turbomachinery_t *tbm = _turbomachinery;
1859
1860 if (tbm == NULL)
1861 return;
1862
1863 cs_real_t *t_angle;
1864 BFT_MALLOC(t_angle, tbm->n_rotors+2, cs_real_t);
1865
1866 t_angle[0] = tbm->t_cur;
1867 for (int i = 0; i < tbm->n_rotors+1; i++) {
1868 cs_rotation_t *rot = tbm->rotation + i;
1869 t_angle[i+1] = rot->angle;
1870 }
1871
1872 int retcode = cs_restart_read_section(r,
1873 "turbomachinery:rotor_time_and_angle",
1874 CS_MESH_LOCATION_NONE,
1875 tbm->n_rotors+2,
1876 CS_TYPE_cs_real_t,
1877 t_angle);
1878
1879 if (retcode == CS_RESTART_SUCCESS) {
1880 tbm->t_cur = t_angle[0];
1881 for (int i = 0; i < tbm->n_rotors+1; i++) {
1882 cs_rotation_t *rot = tbm->rotation + i;
1883 rot->angle = t_angle[i+1];
1884 }
1885 }
1886
1887 BFT_FREE(t_angle);
1888 }
1889
1890 /*----------------------------------------------------------------------------*/
1891 /*!
1892 * \brief Write turbomachinery metadata to checkpoint file.
1893 *
1894 * The mesh is handled separately.
1895 *
1896 * \param[in, out] r associated restart file pointer
1897 */
1898 /*----------------------------------------------------------------------------*/
1899
1900 void
cs_turbomachinery_restart_write(cs_restart_t * r)1901 cs_turbomachinery_restart_write(cs_restart_t *r)
1902 {
1903 const cs_turbomachinery_t *tbm = _turbomachinery;
1904
1905 if (tbm == NULL)
1906 return;
1907
1908 cs_real_t *t_angle;
1909 BFT_MALLOC(t_angle, tbm->n_rotors+2, cs_real_t);
1910
1911 t_angle[0] = tbm->t_cur;
1912 for (int i = 0; i < tbm->n_rotors+1; i++) {
1913 cs_rotation_t *rot = tbm->rotation + i;
1914 t_angle[i+1] = rot->angle;
1915 }
1916
1917 cs_restart_write_section(r,
1918 "turbomachinery:rotor_time_and_angle",
1919 CS_MESH_LOCATION_NONE,
1920 tbm->n_rotors+2,
1921 CS_TYPE_cs_real_t,
1922 t_angle);
1923
1924 BFT_FREE(t_angle);
1925 }
1926
1927 /*----------------------------------------------------------------------------*/
1928
1929 END_C_DECLS
1930