1 /*============================================================================
2 * Functions associated to ALE formulation
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 /*----------------------------------------------------------------------------
34 * Local headers
35 *----------------------------------------------------------------------------*/
36
37 #include "bft_mem.h"
38 #include "bft_printf.h"
39
40 #include "cs_base.h"
41 #include "cs_boundary_conditions.h"
42 #include "cs_boundary_zone.h"
43 #include "cs_cell_to_vertex.h"
44 #include "cs_cdo_quantities.h"
45 #include "cs_cdo_connect.h"
46 #include "cs_cdo_main.h"
47 #include "cs_convection_diffusion.h"
48 #include "cs_domain.h"
49 #include "cs_domain_setup.h"
50 #include "cs_equation.h"
51 #include "cs_equation_iterative_solve.h"
52 #include "cs_face_viscosity.h"
53 #include "cs_field.h"
54 #include "cs_field_pointer.h"
55 #include "cs_field_operator.h"
56 #include "cs_gui_mobile_mesh.h"
57 #include "cs_interface.h"
58 #include "cs_log.h"
59 #include "cs_physical_constants.h"
60 #include "cs_math.h"
61 #include "cs_mesh.h"
62 #include "cs_mesh_quantities.h"
63 #include "cs_mesh_bad_cells.h"
64 #include "cs_parall.h"
65 #include "cs_time_step.h"
66
67 /*----------------------------------------------------------------------------
68 * Header for the current file
69 *----------------------------------------------------------------------------*/
70
71 #include "cs_ale.h"
72
73 /*----------------------------------------------------------------------------*/
74
75 BEGIN_C_DECLS
76
77 /*============================================================================
78 * Global variables
79 *============================================================================*/
80
81 int cs_glob_ale = 0;
82
83 /*============================================================================
84 * Fortran wrapper function definitions
85 *============================================================================*/
86
87 /*----------------------------------------------------------------------------
88 * Get pointer to cs_glob_ale
89 *----------------------------------------------------------------------------*/
90
91 void
cs_f_ale_get_pointers(int ** iale)92 cs_f_ale_get_pointers(int **iale)
93 {
94 *iale = &cs_glob_ale;
95 }
96
97 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
98
99 /*============================================================================
100 * Type definitions
101 *============================================================================*/
102
103 typedef struct {
104
105 /* Array of values to set the boundary conditions (array allocated to all
106 mesh vertices since there is no dedicated numbering to identify easily
107 vertices lying on the boundary */
108 cs_real_t *vtx_values;
109
110 /* Pointer on a list of arrays of selected vertices for each zone associated
111 to a definition by array \ref CS_XDEF_BY_ARRAY for
112 CS_BOUNDARY_ALE_IMPOSED_VEL and CS_BOUNDARY_ALE_IMPOSED_DISP (definitions
113 by value or by a sliding condition are excluded. The position of the array
114 in the list is given implicitly and is related to the order in which the
115 boundary conditions are defined. (cf. \ref cs_ale_setup_boundaries
116 definition for more details). The aim of this list of arrays is to speed-up
117 the settings of the boundary conditions by avoiding doing several times the
118 same enforcement. */
119
120 int n_selections; /* Number of selections */
121 cs_lnum_t *n_vertices; /* Number of vertices in each selections */
122 cs_lnum_t **vtx_select; /* List of vertices for each selection */
123
124 } cs_ale_cdo_bc_t;
125
126 /*============================================================================
127 * Static global variables
128 *============================================================================*/
129
130 static cs_real_3_t *_vtx_coord0 = NULL;
131 static cs_ale_cdo_bc_t *_cdo_bc = NULL;
132
133 static bool cs_ale_active = false;
134
135 /*----------------------------------------------------------------------------
136 * Deprecated ALE boundary condition types
137 *----------------------------------------------------------------------------*/
138 enum {
139 CS_ALE_FIXED = 1,
140 CS_ALE_SLIDING = 2,
141 CS_ALE_IMPOSED_VEL = 3
142 };
143
144 /*============================================================================
145 * Private function definitions
146 *============================================================================*/
147
148 /*----------------------------------------------------------------------------*/
149 /*!
150 * \brief Update the values of mesh vertices lying on a free surface boundary
151 * condition.
152 *
153 * \param[in] domain pointer to a \ref cs_domain_t structure
154 * \param[in] z pointer to a \ref cs_zone_t structure
155 * \param[in] select_id id in the list of selected vertices
156 */
157 /*----------------------------------------------------------------------------*/
158
159 static void
_free_surface(const cs_domain_t * domain,const cs_zone_t * z,const int select_id)160 _free_surface(const cs_domain_t *domain,
161 const cs_zone_t *z,
162 const int select_id)
163 {
164 const cs_real_t *grav = cs_glob_physical_constants->gravity;
165 const cs_mesh_t *m = domain->mesh;
166 const cs_real_3_t *restrict vtx_coord
167 = (const cs_real_3_t *restrict)m->vtx_coord;
168 const cs_mesh_quantities_t *mq = domain->mesh_quantities;
169 const cs_real_3_t *restrict b_face_normal
170 = (const cs_real_3_t *restrict)mq->b_face_normal;
171 const cs_real_3_t *restrict b_face_cog
172 = (const cs_real_3_t *restrict)mq->b_face_cog;
173
174 /* Boundary mass flux */
175 int iflmab = cs_field_get_key_int(CS_F_(vel),
176 cs_field_key_id("boundary_mass_flux_id"));
177 const cs_real_t *b_mass_flux = cs_field_by_id(iflmab)->val;
178
179 /* Transform face flux to vertex displacement */
180 cs_real_3_t *_mesh_vel = NULL;
181
182 /* Dual surface associated to vertices */
183 cs_real_t *_v_surf = NULL;
184
185 /* Squared sum of partial surface associated to a vertex */
186 cs_real_t *_is_loc_min = NULL;
187 cs_real_t *_is_loc_max = NULL;
188
189 cs_real_3_t normal;
190 /* Normal direction is given by the gravity */
191 cs_math_3_normalise(grav, normal);
192
193 const cs_real_t invdt = 1./domain->time_step->dt_ref; /* JB: dt[0] ? */
194
195 BFT_MALLOC(_mesh_vel, m->n_vertices, cs_real_3_t);
196 BFT_MALLOC(_v_surf, m->n_vertices, cs_real_t);
197 BFT_MALLOC(_is_loc_min, m->n_vertices, cs_real_t);
198 BFT_MALLOC(_is_loc_max, m->n_vertices, cs_real_t);
199
200 for (cs_lnum_t v_id = 0; v_id < m->n_vertices; v_id++) {
201 _mesh_vel[v_id][0] = 0;
202 _mesh_vel[v_id][1] = 0;
203 _mesh_vel[v_id][2] = 0;
204 _v_surf[v_id] = 0;
205 _is_loc_min[v_id] = 1;
206 _is_loc_max[v_id] = 1;
207 }
208
209 /* First Loop over boundary faces
210 * to compute if there is a local min and max elevation */
211 for (cs_lnum_t elt_id = 0; elt_id < z->n_elts; elt_id++) {
212
213 const cs_lnum_t face_id = z->elt_ids[elt_id];
214
215 const cs_lnum_t s = m->b_face_vtx_idx[face_id];
216 const cs_lnum_t e = m->b_face_vtx_idx[face_id+1];
217
218 /* Compute the portion of surface associated to v_id_1 */
219 for (cs_lnum_t k = s; k < e; k++) {
220
221 const cs_lnum_t v_id = m->b_face_vtx_lst[k];
222
223 cs_real_3_t v_cog = {
224 b_face_cog[face_id][0] - vtx_coord[v_id][0],
225 b_face_cog[face_id][1] - vtx_coord[v_id][1],
226 b_face_cog[face_id][2] - vtx_coord[v_id][2],
227 };
228
229 /* g . (x_f - x_N) S_fN */
230 cs_real_t dz_fn = cs_math_3_dot_product(normal, v_cog);
231
232 /* v1 is higher than x_f */
233 if (dz_fn > 0.)
234 _is_loc_min[v_id] = 0.; /* not a min */
235
236 /* x_f is higher than v1 */
237 if (dz_fn < 0.)
238 _is_loc_max[v_id] = 0.; /* not a max */
239
240 }
241
242 } /* Loop on selected border faces */
243
244 /* Handle parallelism */
245 if (m->vtx_interfaces != NULL) {
246
247 cs_interface_set_min(m->vtx_interfaces,
248 m->n_vertices,
249 1,
250 true,
251 CS_REAL_TYPE,
252 _is_loc_min);
253
254 cs_interface_set_min(m->vtx_interfaces,
255 m->n_vertices,
256 1,
257 true,
258 CS_REAL_TYPE,
259 _is_loc_max);
260 }
261
262 cs_gnum_t _f_count_filter = 0;
263 cs_gnum_t _f_n_elts = z->n_elts;
264
265 /* Loop over boundary faces */
266 for (cs_lnum_t elt_id = 0; elt_id < z->n_elts; elt_id++) {
267
268 const cs_lnum_t face_id = z->elt_ids[elt_id];
269
270 const cs_real_t mf_inv_rho_n_dot_s = b_mass_flux[face_id] /
271 (cs_math_3_dot_product(normal, b_face_normal[face_id])
272 * CS_F_(rho_b)->val[face_id]);
273
274 cs_real_3_t f_vel = {
275 normal[0] * mf_inv_rho_n_dot_s,
276 normal[1] * mf_inv_rho_n_dot_s,
277 normal[2] * mf_inv_rho_n_dot_s
278 };
279
280 const cs_lnum_t s = m->b_face_vtx_idx[face_id];
281 const cs_lnum_t e = m->b_face_vtx_idx[face_id+1];
282
283 cs_real_t f_has_min = 0;
284 cs_real_t f_has_max = 0;
285
286 /* First loop to determine if filtering is needed:
287 * - if face has a local min and a local max
288 * - if the slope is greater than pi/4
289 * */
290 int f_need_filter = 0;
291
292 for (cs_lnum_t k = s; k < e; k++) {
293 const cs_lnum_t k_1 = (k+1 < e) ? k+1 : k+1-(e-s);
294 const cs_lnum_t v_id0 = m->b_face_vtx_lst[k];
295 const cs_lnum_t v_id1 = m->b_face_vtx_lst[k_1];
296 /* Edge to CoG vector */
297 cs_real_3_t e_cog = {
298 0.5 * (b_face_cog[face_id][0] - vtx_coord[v_id0][0])
299 + 0.5 * (b_face_cog[face_id][0] - vtx_coord[v_id1][0]),
300 0.5 * (b_face_cog[face_id][1] - vtx_coord[v_id0][1])
301 + 0.5 * (b_face_cog[face_id][1] - vtx_coord[v_id1][1]),
302 0.5 * (b_face_cog[face_id][2] - vtx_coord[v_id0][2])
303 + 0.5 * (b_face_cog[face_id][2] - vtx_coord[v_id1][2])
304 };
305 cs_real_t dz = CS_ABS(cs_math_3_dot_product(normal, e_cog));
306 cs_real_3_t e_cog_hor;
307 cs_math_3_orthogonal_projection(normal, e_cog, e_cog_hor);
308 cs_real_t dx = cs_math_3_norm(e_cog_hor);
309 /* Too high slope */
310 if (dz > dx)
311 f_need_filter = 1;
312
313 f_has_max = CS_MAX(f_has_max, _is_loc_max[v_id0]);
314 f_has_min = CS_MAX(f_has_min, _is_loc_min[v_id0]);
315 }
316
317 if ((f_has_max > 0.5 && f_has_min > 0.5) || f_need_filter == 1) {
318 f_need_filter = 1;
319 _f_count_filter++;
320 }
321
322 /* Compute the portion of surface associated to v_id_1 */
323 for (cs_lnum_t k = s; k < e; k++) {
324
325 const cs_lnum_t k_1 = (k+1 < e) ? k+1 : k+1-(e-s);
326 const cs_lnum_t k_2 = (k+2 < e) ? k+2 : k+2-(e-s);
327 const cs_lnum_t v_id0 = m->b_face_vtx_lst[k];
328 const cs_lnum_t v_id1 = m->b_face_vtx_lst[k_1];
329 const cs_lnum_t v_id2 = m->b_face_vtx_lst[k_2];
330
331 cs_real_3_t v0v1 = {
332 vtx_coord[v_id1][0] - vtx_coord[v_id0][0],
333 vtx_coord[v_id1][1] - vtx_coord[v_id0][1],
334 vtx_coord[v_id1][2] - vtx_coord[v_id0][2]
335 };
336 cs_real_3_t v1v2 = {
337 vtx_coord[v_id2][0] - vtx_coord[v_id1][0],
338 vtx_coord[v_id2][1] - vtx_coord[v_id1][1],
339 vtx_coord[v_id2][2] - vtx_coord[v_id1][2]
340 };
341 cs_real_3_t v1_cog = {
342 b_face_cog[face_id][0] - vtx_coord[v_id1][0],
343 b_face_cog[face_id][1] - vtx_coord[v_id1][1],
344 b_face_cog[face_id][2] - vtx_coord[v_id1][2]
345 };
346
347 /* Portion of the surface associated to the vertex
348 * projected in the normal direction */
349 cs_real_t portion_surf = -0.25 * (
350 cs_math_3_triple_product(v0v1, v1_cog, normal)
351 + cs_math_3_triple_product(v1v2, v1_cog, normal));
352
353 _v_surf[v_id1] += portion_surf;
354
355 /* g . (x_f - x_N) S_fN */
356 cs_real_t dz_fn = cs_math_3_dot_product(normal, v1_cog);
357
358 for (int i = 0; i < 3; i++)
359 _mesh_vel[v_id1][i] += ( f_vel[i]
360 + invdt * f_need_filter * dz_fn * normal[i])
361 * portion_surf;
362 }
363
364 } /* Loop on selected border faces */
365
366 cs_var_cal_opt_t var_cal_opt;
367 cs_field_get_key_struct(CS_F_(mesh_u),
368 cs_field_key_id("var_cal_opt"),
369 &var_cal_opt);
370 if (var_cal_opt.verbosity >= 1) {
371 cs_parall_sum(1, CS_GNUM_TYPE, &_f_count_filter);
372 cs_parall_sum(1, CS_GNUM_TYPE, &_f_n_elts);
373 bft_printf("Free surface condition %d: %f percents of limited face\n",
374 select_id, (cs_real_t) _f_count_filter / (cs_real_t) _f_n_elts);
375 }
376
377
378 /* Handle parallelism */
379 if (m->vtx_interfaces != NULL) {
380 cs_interface_set_sum(m->vtx_interfaces,
381 m->n_vertices,
382 3,
383 true,
384 CS_REAL_TYPE,
385 _mesh_vel);
386
387 cs_interface_set_sum(m->vtx_interfaces,
388 m->n_vertices,
389 1,
390 true,
391 CS_REAL_TYPE,
392 _v_surf);
393 }
394
395 /* Loop on selected border vertices */
396 for (cs_lnum_t i = 0; i < _cdo_bc->n_vertices[select_id]; i++) {
397
398 const cs_lnum_t v_id = _cdo_bc->vtx_select[select_id][i];
399 const double invsurf = 1./_v_surf[v_id];
400
401 cs_real_t *_val = _cdo_bc->vtx_values + 3*v_id;
402
403 for (int k = 0; k < 3; k++) {
404 _mesh_vel[v_id][k] *= invsurf;
405 _val[k] = _mesh_vel[v_id][k];
406 }
407
408 } /* Loop on selected vertices */
409
410 /* Free memory */
411 BFT_FREE(_mesh_vel);
412 BFT_FREE(_v_surf);
413 BFT_FREE(_is_loc_min);
414 BFT_FREE(_is_loc_max);
415 }
416
417 /*----------------------------------------------------------------------------*/
418 /*!
419 * \brief Add a new selection of vertices related to a definition by array
420 * Update the cs_ale_cdo_bc_t structure.
421 *
422 * \param[in] mesh pointer to a \ref cs_mesh_t structure
423 * \param[in] z pointer to a cs_zone_t structure
424 * \param[in, out] vtag tag array on vertices
425 */
426 /*----------------------------------------------------------------------------*/
427
428 static void
_update_bc_list(const cs_mesh_t * mesh,const cs_zone_t * z,bool vtag[])429 _update_bc_list(const cs_mesh_t *mesh,
430 const cs_zone_t *z,
431 bool vtag[])
432 {
433 const cs_lnum_t *bf2v_idx = mesh->b_face_vtx_idx;
434 const cs_lnum_t *bf2v_lst = mesh->b_face_vtx_lst;
435 const cs_lnum_t n_vertices = mesh->n_vertices;
436 const int id = _cdo_bc->n_selections;
437
438 cs_lnum_t counter = 0;
439
440 _cdo_bc->n_selections++;
441 BFT_REALLOC(_cdo_bc->n_vertices, _cdo_bc->n_selections, cs_lnum_t);
442 BFT_REALLOC(_cdo_bc->vtx_select, _cdo_bc->n_selections, cs_lnum_t *);
443
444 /* Reset vtag */
445 memset(vtag, 0, n_vertices*sizeof(bool));
446
447 /* Count the number of vertices to select */
448 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
449
450 const cs_lnum_t bf_id = z->elt_ids[i];
451 const cs_lnum_t *idx = bf2v_idx + bf_id;
452 const cs_lnum_t *lst = bf2v_lst + idx[0];
453
454 /* Loop on face vertices */
455 for (cs_lnum_t j = 0; j < idx[1]-idx[0]; j++) {
456 cs_lnum_t v_id = lst[j];
457 if (!vtag[v_id]) { /* Not already selected */
458 vtag[v_id] = true;
459 counter++;
460 }
461 }
462
463 } /* Loop on selected border faces */
464
465 _cdo_bc->n_vertices[id] = counter;
466 BFT_MALLOC(_cdo_bc->vtx_select[id], counter, cs_lnum_t);
467
468 /* Fill the list of selected vertices */
469 memset(vtag, 0, n_vertices*sizeof(bool));
470 counter = 0;
471 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
472
473 const cs_lnum_t bf_id = z->elt_ids[i];
474 const cs_lnum_t *idx = bf2v_idx + bf_id;
475 const cs_lnum_t *lst = bf2v_lst + idx[0];
476
477 /* Loop on face vertices */
478 for (cs_lnum_t j = 0; j < idx[1]-idx[0]; j++) {
479 cs_lnum_t v_id = lst[j];
480 if (!vtag[v_id]) { /* Not already selected */
481 vtag[v_id] = true;
482 _cdo_bc->vtx_select[id][counter] = v_id;
483 counter++;
484 }
485 }
486
487 } /* Loop on selected border faces */
488
489 assert(counter == _cdo_bc->n_vertices[id]);
490 }
491
492 /*----------------------------------------------------------------------------*/
493 /*!
494 * \brief This function updates all ALE BCs for CDO except free surface.
495 * These BCs are required for the fluid BCs and therefore needs to be
496 * updated before.
497 *
498 * \param[in] domain domain quantities
499 * \param[out] b_fluid_vel boundary fluid velocity
500 */
501 /*----------------------------------------------------------------------------*/
502
503 static void
_update_bcs(const cs_domain_t * domain,int ale_bc_type[],cs_real_3_t b_fluid_vel[])504 _update_bcs(const cs_domain_t *domain,
505 int ale_bc_type[],
506 cs_real_3_t b_fluid_vel[])
507 {
508 const cs_mesh_t *m = domain->mesh;
509 const cs_mesh_quantities_t *mq = domain->mesh_quantities;
510
511 const cs_real_3_t *b_face_normal = (const cs_real_3_t *)mq->b_face_normal;
512 const cs_real_3_t *restrict vtx_coord
513 = (const cs_real_3_t *restrict)m->vtx_coord;
514 const cs_real_3_t *restrict b_face_cog
515 = (const cs_real_3_t *restrict)mq->b_face_cog;
516
517 cs_field_t *f_displ = cs_field_by_name("mesh_displacement");
518
519 /* Only a part of the boundaries has to be updated */
520 int select_id = 0;
521 for (int b_id = 0; b_id < domain->ale_boundaries->n_boundaries; b_id++) {
522
523 const int z_id = domain->ale_boundaries->zone_ids[b_id];
524 const cs_zone_t *z = cs_boundary_zone_by_id(z_id);
525
526 switch(domain->ale_boundaries->types[b_id]) {
527
528 case CS_BOUNDARY_ALE_IMPOSED_VEL:
529 {
530 /* Retrieve the velocities to enforce on faces */
531 cs_real_t *bc_vals = cs_gui_mobile_mesh_get_fixed_velocity(z->name);
532
533 assert(select_id < _cdo_bc->n_selections);
534
535 #if 0 //TODO when we will have meg on vertices
536 /* Loop on selected border vertices */
537 for (cs_lnum_t i = 0; i < _cdo_bc->n_vertices[select_id]; i++) {
538
539 cs_real_t *_val
540 = _cdo_bc->vtx_values + 3*_cdo_bc->vtx_select[select_id][i];
541 _val[0] = vel[0];
542 _val[1] = vel[1];
543 _val[2] = vel[2];
544 }
545 #endif
546 /* Dual surface associated to vertices */
547 cs_real_t *_v_surf = NULL;
548
549 /* Transform face flux to vertex velocities */
550 cs_real_3_t *_mesh_vel = NULL;
551
552 BFT_MALLOC(_v_surf, m->n_vertices, cs_real_t);
553 BFT_MALLOC(_mesh_vel, m->n_vertices, cs_real_3_t);
554
555 /* Initialize */
556 for (cs_lnum_t v_id = 0; v_id < m->n_vertices; v_id++) {
557 _v_surf[v_id] = 0;
558 _mesh_vel[v_id][0] = 0;
559 _mesh_vel[v_id][1] = 0;
560 _mesh_vel[v_id][2] = 0;
561 }
562
563 /* First Loop over boundary faces
564 * to compute face portion associated to the vertex
565 * and add the portion of face velocity to the vertex velocity */
566 for (cs_lnum_t elt_id = 0; elt_id < z->n_elts; elt_id++) {
567
568 const cs_lnum_t face_id = z->elt_ids[elt_id];
569
570 cs_real_3_t normal;
571 cs_math_3_normalise(b_face_normal[face_id], normal);
572
573 const cs_lnum_t s = m->b_face_vtx_idx[face_id];
574 const cs_lnum_t e = m->b_face_vtx_idx[face_id+1];
575
576
577 /* Compute the portion of surface associated to v_id_1 */
578 for (cs_lnum_t k = s; k < e; k++) {
579
580 const cs_lnum_t k_1 = (k+1 < e) ? k+1 : k+1-(e-s);
581 const cs_lnum_t k_2 = (k+2 < e) ? k+2 : k+2-(e-s);
582 const cs_lnum_t v_id0 = m->b_face_vtx_lst[k];
583 const cs_lnum_t v_id1 = m->b_face_vtx_lst[k_1];
584 const cs_lnum_t v_id2 = m->b_face_vtx_lst[k_2];
585
586 cs_real_3_t v0v1 = {
587 vtx_coord[v_id1][0] - vtx_coord[v_id0][0],
588 vtx_coord[v_id1][1] - vtx_coord[v_id0][1],
589 vtx_coord[v_id1][2] - vtx_coord[v_id0][2]
590 };
591 cs_real_3_t v1v2 = {
592 vtx_coord[v_id2][0] - vtx_coord[v_id1][0],
593 vtx_coord[v_id2][1] - vtx_coord[v_id1][1],
594 vtx_coord[v_id2][2] - vtx_coord[v_id1][2]
595 };
596 cs_real_3_t v1_cog = {
597 b_face_cog[face_id][0] - vtx_coord[v_id1][0],
598 b_face_cog[face_id][1] - vtx_coord[v_id1][1],
599 b_face_cog[face_id][2] - vtx_coord[v_id1][2]
600 };
601
602 /* Portion of the surface associated to the vertex
603 * projected in the normal direction */
604 cs_real_t portion_surf = 0.25 * (
605 cs_math_3_triple_product(v0v1, v1_cog, normal)
606 + cs_math_3_triple_product(v1v2, v1_cog, normal));
607
608 _v_surf[v_id1] += portion_surf;
609
610 /* Portion of the face velocity is added to the node velocity
611 * Warning: bc_vals is not interleaved */
612 for (int i = 0; i < 3; i++)
613 _mesh_vel[v_id1][i] += bc_vals[elt_id + i * z->n_elts] * portion_surf;
614
615 }
616
617 } /* Loop on selected border faces */
618
619 /* Handle parallelism */
620 if (m->vtx_interfaces != NULL) {
621 cs_interface_set_sum(m->vtx_interfaces,
622 m->n_vertices,
623 3,
624 true,
625 CS_REAL_TYPE,
626 _mesh_vel);
627
628 cs_interface_set_sum(m->vtx_interfaces,
629 m->n_vertices,
630 1,
631 true,
632 CS_REAL_TYPE,
633 _v_surf);
634 }
635
636 /* Loop on selected border vertices */
637 for (cs_lnum_t i = 0; i < _cdo_bc->n_vertices[select_id]; i++) {
638
639 const cs_lnum_t v_id = _cdo_bc->vtx_select[select_id][i];
640 const double invsurf = 1./_v_surf[v_id];
641
642 cs_real_t *_val = _cdo_bc->vtx_values + 3*v_id;
643
644 for (int k = 0; k < 3; k++) {
645 _mesh_vel[v_id][k] *= invsurf;
646 _val[k] = _mesh_vel[v_id][k];
647 }
648 }
649
650 /* Loop over boundary faces and impose node mesh velocity */
651 for (cs_lnum_t elt_id = 0; elt_id < z->n_elts; elt_id++) {
652 const cs_lnum_t face_id = z->elt_ids[elt_id];
653
654 /* fluid velocity BC */
655 ale_bc_type[face_id] = CS_ALE_IMPOSED_VEL;
656 for (int d = 0; d < 3; d++)
657 b_fluid_vel[face_id][d] = bc_vals[elt_id + d * z->n_elts];
658 }
659
660 /* Free memory */
661 BFT_FREE(bc_vals);
662 BFT_FREE(_mesh_vel);
663 BFT_FREE(_v_surf);
664
665 select_id++;
666 }
667 break;
668
669 case CS_BOUNDARY_ALE_IMPOSED_DISP:
670 {
671 const cs_real_3_t *restrict disale
672 = (const cs_real_3_t *restrict)(f_displ->val);
673 const cs_real_t invdt = 1./domain->time_step->dt_ref; /* JB: dt[0] ? */
674
675 assert(select_id < _cdo_bc->n_selections);
676
677 /* Loop on selected border vertices */
678 for (cs_lnum_t i = 0; i < _cdo_bc->n_vertices[select_id]; i++) {
679
680 const cs_lnum_t v_id = _cdo_bc->vtx_select[select_id][i];
681 const cs_real_t *_dpl = disale[v_id];
682 const cs_real_t *restrict _xyz = vtx_coord[v_id];
683 const cs_real_t *restrict _xyz0 = _vtx_coord0[v_id];
684
685 cs_real_t *_val = _cdo_bc->vtx_values + 3*v_id;
686
687 for (int k = 0; k < 3; k++)
688 _val[k] = (_dpl[k] + _xyz0[k] - _xyz[k])* invdt;
689
690 } /* Loop on selected vertices */
691
692 /* Loop over boundary faces for the fluid velocity */
693 for (cs_lnum_t elt_id = 0; elt_id < z->n_elts; elt_id++) {
694 const cs_lnum_t face_id = z->elt_ids[elt_id];
695
696 ale_bc_type[face_id] = CS_ALE_IMPOSED_VEL;
697
698 cs_real_3_t normal;
699 /* Normal direction is given by the gravity */
700 cs_math_3_normalise(b_face_normal[face_id], normal);
701 const cs_real_t dsurf = 1. / mq->b_face_surf[face_id];
702
703 b_fluid_vel[face_id][0] = 0.;
704 b_fluid_vel[face_id][1] = 0.;
705 b_fluid_vel[face_id][2] = 0.;
706
707 const cs_lnum_t s = m->b_face_vtx_idx[face_id];
708 const cs_lnum_t e = m->b_face_vtx_idx[face_id+1];
709
710 /* Compute the portion of surface associated to v_id_1 */
711 for (cs_lnum_t k = s; k < e; k++) {
712
713 const cs_lnum_t k_1 = (k+1 < e) ? k+1 : k+1-(e-s);
714 const cs_lnum_t k_2 = (k+2 < e) ? k+2 : k+2-(e-s);
715 const cs_lnum_t v_id0 = m->b_face_vtx_lst[k];
716 const cs_lnum_t v_id1 = m->b_face_vtx_lst[k_1];
717 const cs_lnum_t v_id2 = m->b_face_vtx_lst[k_2];
718
719 cs_real_3_t v0v1 = {
720 vtx_coord[v_id1][0] - vtx_coord[v_id0][0],
721 vtx_coord[v_id1][1] - vtx_coord[v_id0][1],
722 vtx_coord[v_id1][2] - vtx_coord[v_id0][2],
723 };
724 cs_real_3_t v1v2 = {
725 vtx_coord[v_id2][0] - vtx_coord[v_id1][0],
726 vtx_coord[v_id2][1] - vtx_coord[v_id1][1],
727 vtx_coord[v_id2][2] - vtx_coord[v_id1][2],
728 };
729 cs_real_3_t v1_cog = {
730 b_face_cog[face_id][0] - vtx_coord[v_id1][0],
731 b_face_cog[face_id][1] - vtx_coord[v_id1][1],
732 b_face_cog[face_id][2] - vtx_coord[v_id1][2],
733 };
734
735 /* Portion of the surface associated to the vertex
736 * projected in the normal direction */
737 cs_real_t portion_surf = 0.25 * dsurf * (
738 cs_math_3_triple_product(v0v1, v1_cog, normal)
739 + cs_math_3_triple_product(v1v2, v1_cog, normal));
740
741 cs_real_t *_val = _cdo_bc->vtx_values + 3*v_id1;
742
743 b_fluid_vel[face_id][0] += portion_surf * _val[0];
744 b_fluid_vel[face_id][1] += portion_surf * _val[1];
745 b_fluid_vel[face_id][2] += portion_surf * _val[1];
746 }
747
748 }
749
750 select_id++;
751 }
752 break;
753
754 /* Treated elsewhere, only increased selected_id */
755 case CS_BOUNDARY_ALE_FREE_SURFACE:
756 assert(select_id < _cdo_bc->n_selections);
757 select_id++;
758 break;
759
760 default:
761 break; /* Nothing to do */
762 }
763
764 } /* Loop on ALE boundaries */
765 }
766
767 /*----------------------------------------------------------------------------*/
768 /*!
769 * \brief This function updates free surface ALE BCs for CDO.
770 * These BCs are required after the fluid is solved.
771 *
772 * \param[in] domain domain quantities
773 */
774 /*----------------------------------------------------------------------------*/
775
776 static void
_update_bcs_free_surface(const cs_domain_t * domain)777 _update_bcs_free_surface(const cs_domain_t *domain)
778 {
779 /* Only a part of the boundaries has to be updated */
780 int select_id = 0;
781 for (int b_id = 0; b_id < domain->ale_boundaries->n_boundaries; b_id++) {
782
783 const int z_id = domain->ale_boundaries->zone_ids[b_id];
784 const cs_zone_t *z = cs_boundary_zone_by_id(z_id);
785
786 switch(domain->ale_boundaries->types[b_id]) {
787
788 /* Treated elsewhere, only increased selected_id */
789 case CS_BOUNDARY_ALE_IMPOSED_VEL:
790 select_id++;
791 break;
792
793 /* Treated elsewhere, only increased selected_id */
794 case CS_BOUNDARY_ALE_IMPOSED_DISP:
795 select_id++;
796 break;
797
798 case CS_BOUNDARY_ALE_FREE_SURFACE:
799 assert(select_id < _cdo_bc->n_selections);
800 _free_surface(domain, z, select_id);
801 select_id++;
802 break;
803
804 default:
805 break; /* Nothing to do */
806 }
807
808 } /* Loop on ALE boundaries */
809 }
810
811 /*----------------------------------------------------------------------------*/
812 /*!
813 * \brief This subroutine performs the solving of a Poisson equation
814 * on the mesh velocity for ALE module. It also updates the mesh
815 * displacement so that it can be used to update mass fluxes (due to
816 * mesh displacement).
817 *
818 * \param[in] domain domain quantities
819 * \param[in] impale Indicator for fixed node displacement
820 */
821 /*----------------------------------------------------------------------------*/
822
823 static void
_ale_solve_poisson_cdo(const cs_domain_t * domain,const int impale[])824 _ale_solve_poisson_cdo(const cs_domain_t *domain,
825 const int impale[])
826 {
827 const cs_mesh_t *m = domain->mesh;
828
829 /* Build and solve equation on the mesh velocity */
830 cs_equation_t *eq = cs_equation_by_name("mesh_velocity");
831
832 /* Update the values of boundary mesh vertices */
833 _update_bcs_free_surface(domain);
834
835 /* Compute the Poisson equation on the original mesh */
836 cs_real_3_t *vtx_coord = (cs_real_3_t *)m->vtx_coord;
837 cs_real_3_t *vtx_coord0 = (cs_real_3_t *)(cs_field_by_name("vtx_coord0")->val);
838 const cs_lnum_t n_vertices = m->n_vertices;
839 cs_mesh_quantities_t *mq = domain->mesh_quantities;
840
841 /* Back to original mesh */
842 for (cs_lnum_t v_id = 0; v_id < n_vertices; v_id++) {
843 for (cs_lnum_t idim = 0; idim < 3; idim++) {
844 vtx_coord[v_id][idim] = vtx_coord0[v_id][idim];
845 }
846 }
847
848 cs_ale_update_mesh_quantities(&(mq->min_vol), &(mq->max_vol), &(mq->tot_vol));
849
850 /* Solve the algebraic system */
851 cs_equation_solve_steady_state(m, eq);
852
853 /* Retrieving fields */
854 cs_field_t *f_displ = cs_field_by_name("mesh_displacement");
855 cs_real_3_t *disale = (cs_real_3_t *)(f_displ->val);
856 cs_real_3_t *disala = (cs_real_3_t *)(f_displ->val_pre);
857 cs_real_3_t *m_vel = (cs_real_3_t *)(cs_field_by_name("mesh_velocity")->val);
858
859 /* Back to mesh at time n */
860 for (cs_lnum_t v_id = 0; v_id < n_vertices; v_id++) {
861 for (cs_lnum_t idim = 0; idim < 3; idim++) {
862 vtx_coord[v_id][idim] = vtx_coord0[v_id][idim] + disale[v_id][idim];
863 }
864 }
865
866 cs_ale_update_mesh_quantities(&(mq->min_vol), &(mq->max_vol), &(mq->tot_vol));
867
868
869 for (cs_lnum_t v = 0; v < m->n_vertices; v++) {
870 if (impale[v] == 0) {
871 for (int k = 0; k < 3; k++)
872 disale[v][k] = disala[v][k] + m_vel[v][k] * domain->time_step->dt_ref;
873 }
874 }
875 }
876
877 /*----------------------------------------------------------------------------*/
878 /*!
879 * \brief Solve a Poisson equation on the mesh velocity in ALE framework.
880 *
881 * It also updates the mesh displacement
882 * so that it can be used to update mass fluxes (due to mesh displacement).
883 *
884 * \param[in] domain domain quantities
885 * \param[in] iterns Navier-Stokes iteration number
886 * \param[in] impale Indicator for fixed node displacement
887 * \param[in] ale_bc_type Type of boundary for ALE
888 */
889 /*----------------------------------------------------------------------------*/
890
891 static void
_ale_solve_poisson_legacy(const cs_domain_t * domain,const int iterns,const int * impale,const int * ale_bc_type)892 _ale_solve_poisson_legacy(const cs_domain_t *domain,
893 const int iterns,
894 const int *impale,
895 const int *ale_bc_type)
896 {
897 const cs_mesh_t *m = domain->mesh;
898 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
899 const cs_lnum_t n_vertices = m->n_vertices;
900 const cs_lnum_t n_i_faces = m->n_i_faces;
901 const cs_lnum_t n_b_faces = m->n_b_faces;
902 const cs_lnum_t *b_face_cells = (const cs_lnum_t *)m->b_face_cells;
903 const cs_mesh_quantities_t *mq = domain->mesh_quantities;
904 const cs_real_t *b_dist = (const cs_real_t *)mq->b_dist;
905 const cs_real_3_t *b_face_normal = (const cs_real_3_t *)mq->b_face_normal;
906 const cs_real_t *grav = cs_glob_physical_constants->gravity;
907 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
908
909 /* The mass flux is necessary to call cs_equation_iterative_solve_vector
910 but not used (iconv = 0), except for the free surface, where it is used
911 as a boundary condition */
912
913 const int kimasf = cs_field_key_id("inner_mass_flux_id");
914 const int kbmasf = cs_field_key_id("boundary_mass_flux_id");
915 const cs_real_t *i_massflux
916 = cs_field_by_id(cs_field_get_key_int(CS_F_(vel), kimasf))->val;
917 const cs_real_t *b_massflux
918 = cs_field_by_id(cs_field_get_key_int(CS_F_(vel), kbmasf))->val;
919
920 /* 1. Initialization */
921 cs_real_3_t rinfiv = { cs_math_infinite_r,
922 cs_math_infinite_r,
923 cs_math_infinite_r };
924
925 cs_real_3_t *smbr = NULL;
926 cs_real_33_t *fimp = NULL;
927 BFT_MALLOC(smbr, n_cells_ext, cs_real_3_t);
928 BFT_MALLOC(fimp, n_cells_ext, cs_real_33_t);
929
930 cs_real_3_t *mshvel = (cs_real_3_t *)CS_F_(mesh_u)->val;
931 cs_real_3_t *mshvela = (cs_real_3_t *)CS_F_(mesh_u)->val_pre;
932
933 cs_field_t *f_displ = cs_field_by_name("mesh_displacement");
934 cs_real_3_t *disale = (cs_real_3_t *)(f_displ->val);
935 cs_real_3_t *disala = (cs_real_3_t *)(f_displ->val_pre);
936
937 cs_var_cal_opt_t var_cal_opt;
938 cs_field_get_key_struct(CS_F_(mesh_u), key_cal_opt_id, &var_cal_opt);
939
940 if (var_cal_opt.verbosity >= 1)
941 bft_printf("\n ** SOLVING MESH VELOCITY\n"
942 " ---------------------\n");
943
944 /* We compute the boundary condition on the mesh velocity at the free surface
945 * from the new mass flux. */
946
947 /* Density at the boundary */
948 cs_real_t *brom = CS_F_(rho_b)->val;
949
950 cs_field_bc_coeffs_t *bc_coeffs = CS_F_(mesh_u)->bc_coeffs;
951
952 cs_real_3_t *bc_a = (cs_real_3_t *)bc_coeffs->a;
953 cs_real_3_t *bc_af = (cs_real_3_t *)bc_coeffs->af;
954 cs_real_33_t *bc_b = (cs_real_33_t *)bc_coeffs->b;
955 cs_real_33_t *bc_bf = (cs_real_33_t *)bc_coeffs->bf;
956
957 int idftnp = var_cal_opt.idften;
958
959 /* The mesh moves in the direction of the gravity in case of free-surface */
960 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
961 if (ale_bc_type[face_id] == CS_FREE_SURFACE) {
962 cs_lnum_t cell_id = b_face_cells[face_id];
963 cs_real_t distbf = b_dist[face_id];
964
965 cs_real_6_t hintt = {0., 0., 0., 0., 0., 0.};
966 if (idftnp & CS_ISOTROPIC_DIFFUSION) {
967 for (int isou = 0; isou < 3; isou++)
968 hintt[isou] = CS_F_(vism)->val[cell_id] / distbf;
969 } else if (idftnp & CS_ANISOTROPIC_LEFT_DIFFUSION) {
970 for (int isou = 0; isou < 6; isou++)
971 hintt[isou] = CS_F_(vism)->val[6*cell_id+isou] / distbf;
972 }
973
974 cs_real_t prosrf = cs_math_3_dot_product(grav, b_face_normal[face_id]);
975
976 cs_real_3_t pimpv;
977 for (int i = 0; i < 3; i++)
978 pimpv[i] = grav[i]*b_massflux[face_id]/(brom[face_id]*prosrf);
979
980 cs_boundary_conditions_set_dirichlet_vector_aniso(bc_a[face_id],
981 bc_af[face_id],
982 bc_b[face_id],
983 bc_bf[face_id],
984 pimpv,
985 hintt,
986 rinfiv);
987 }
988 }
989
990 /* 2. Solving of the mesh velocity equation */
991
992 if (var_cal_opt.verbosity >= 1)
993 bft_printf("\n\n SOLVING VARIABLE %s\n\n",
994 CS_F_(mesh_u)->name);
995
996 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++) {
997 for (int isou = 0; isou < 3; isou++) {
998 smbr[cell_id][isou] = 0.;
999 for (int jsou = 0; jsou < 3; jsou++)
1000 fimp[cell_id][jsou][isou] = 0.;
1001 }
1002 }
1003
1004 cs_real_t *i_visc = NULL, *b_visc = NULL;
1005
1006 BFT_MALLOC(b_visc, n_b_faces, cs_real_t);
1007
1008 if (idftnp & CS_ISOTROPIC_DIFFUSION) {
1009
1010 BFT_MALLOC(i_visc, n_i_faces, cs_real_t);
1011
1012 cs_face_viscosity(m,
1013 mq,
1014 var_cal_opt.imvisf,
1015 CS_F_(vism)->val,
1016 i_visc,
1017 b_visc);
1018
1019 }
1020 else if (idftnp & CS_ANISOTROPIC_LEFT_DIFFUSION) {
1021
1022 BFT_MALLOC(i_visc, 9*n_i_faces, cs_real_t);
1023
1024 cs_face_anisotropic_viscosity_vector(m,
1025 mq,
1026 var_cal_opt.imvisf,
1027 (cs_real_6_t *)CS_F_(vism)->val,
1028 (cs_real_33_t *)i_visc,
1029 b_visc);
1030 }
1031
1032 var_cal_opt.relaxv = 1.;
1033 var_cal_opt.thetav = 1.;
1034 var_cal_opt.istat = -1;
1035 var_cal_opt.idifft = -1;
1036
1037 cs_equation_iterative_solve_vector(cs_glob_time_step_options->idtvar,
1038 iterns,
1039 CS_F_(mesh_u)->id,
1040 CS_F_(mesh_u)->name,
1041 0, /* ivisep */
1042 0, /* iescap */
1043 &var_cal_opt,
1044 (const cs_real_3_t *)mshvela,
1045 (const cs_real_3_t *)mshvela,
1046 (const cs_real_3_t *)bc_coeffs->a,
1047 (const cs_real_33_t *)bc_coeffs->b,
1048 (const cs_real_3_t *)bc_coeffs->af,
1049 (const cs_real_33_t *)bc_coeffs->bf,
1050 i_massflux,
1051 b_massflux,
1052 i_visc,
1053 b_visc,
1054 i_visc,
1055 b_visc,
1056 NULL, /* i_secvis */
1057 NULL, /* b_secvis */
1058 NULL, /* viscel */
1059 NULL, /* weighf */
1060 NULL, /* weighb */
1061 0, /* icvflv */
1062 NULL, /* icvfli */
1063 (cs_real_33_t *)fimp,
1064 smbr,
1065 mshvel,
1066 NULL); /* eswork */
1067
1068 /* Free memory */
1069 BFT_FREE(smbr);
1070 BFT_FREE(fimp);
1071 BFT_FREE(i_visc);
1072 BFT_FREE(b_visc);
1073
1074 /* 3. Update nodes displacement */
1075
1076 cs_real_3_t *dproj;
1077 cs_real_33_t *gradm;
1078
1079 /* Allocate a temporary array */
1080 BFT_MALLOC(dproj, n_vertices, cs_real_3_t);
1081 BFT_MALLOC(gradm, n_cells_ext, cs_real_33_t);
1082
1083 bool use_previous_t = false;
1084 int inc = 1;
1085
1086 cs_field_gradient_vector(CS_F_(mesh_u),
1087 use_previous_t,
1088 inc,
1089 gradm);
1090
1091 cs_ale_project_displacement(ale_bc_type,
1092 (const cs_real_3_t *)mshvel,
1093 (const cs_real_33_t *)gradm,
1094 (const cs_real_3_t *)bc_coeffs->a,
1095 (const cs_real_33_t *)bc_coeffs->b,
1096 (const cs_real_t *)CS_F_(dt)->val,
1097 dproj);
1098
1099 /* FIXME : warning if nterup > 1, use itrale ? */
1100 /* Update mesh displacement only where it is not
1101 imposed by the user (ie when impale <> 1) */
1102 for (cs_lnum_t v = 0; v < n_vertices; v++) {
1103 if (impale[v] == 0) {
1104 for (int isou = 0; isou < 3; isou++)
1105 disale[v][isou] = disala[v][isou] + dproj[v][isou];
1106 }
1107 }
1108
1109 /* Free memory */
1110 BFT_FREE(dproj);
1111 BFT_FREE(gradm);
1112 }
1113
1114 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1115
1116 /*============================================================================
1117 * Public function definitions
1118 *============================================================================*/
1119
1120 /*----------------------------------------------------------------------------*/
1121 /*!
1122 * \brief Compute cell and face centers of gravity, cell volumes
1123 * and update bad cells.
1124 *
1125 * \param[out] min_vol Minimum cell volume
1126 * \param[out] max_vol Maximum cell volume
1127 * \param[out] tot_vol Total cell volume
1128 */
1129 /*----------------------------------------------------------------------------*/
1130
1131 void
cs_ale_update_mesh_quantities(cs_real_t * min_vol,cs_real_t * max_vol,cs_real_t * tot_vol)1132 cs_ale_update_mesh_quantities(cs_real_t *min_vol,
1133 cs_real_t *max_vol,
1134 cs_real_t *tot_vol)
1135 {
1136 cs_mesh_t *m = cs_glob_mesh;
1137 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1138
1139 cs_gradient_free_quantities();
1140 cs_cell_to_vertex_free();
1141 cs_mesh_quantities_compute(m, mq);
1142 cs_mesh_bad_cells_detect(m, mq);
1143
1144 *min_vol = mq->min_vol;
1145 *max_vol = mq->max_vol;
1146 *tot_vol = mq->tot_vol;
1147 }
1148
1149 /*----------------------------------------------------------------------------*/
1150 /*!
1151 * \brief Project the displacement on mesh vertices (solved on cell center).
1152 *
1153 * \param[in] ale_bc_type Type of boundary for ALE
1154 * \param[in] meshv Mesh velocity
1155 * \param[in] gradm Mesh velocity gradient
1156 * (du_i/dx_j : gradv[][i][j])
1157 * \param[in] claale Boundary conditions A
1158 * \param[in] clbale Boundary conditions B
1159 * \param[in] dt Time step
1160 * \param[out] disp_proj Displacement projected on vertices
1161 */
1162 /*----------------------------------------------------------------------------*/
1163
1164 void
cs_ale_project_displacement(const int ale_bc_type[],const cs_real_3_t * meshv,const cs_real_33_t gradm[],const cs_real_3_t * claale,const cs_real_33_t * clbale,const cs_real_t * dt,cs_real_3_t * disp_proj)1165 cs_ale_project_displacement(const int ale_bc_type[],
1166 const cs_real_3_t *meshv,
1167 const cs_real_33_t gradm[],
1168 const cs_real_3_t *claale,
1169 const cs_real_33_t *clbale,
1170 const cs_real_t *dt,
1171 cs_real_3_t *disp_proj)
1172 {
1173 bool *vtx_interior_indicator = NULL;
1174 cs_real_t *vtx_counter = NULL;
1175 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1176
1177 const cs_mesh_t *m = cs_glob_mesh;
1178 const cs_lnum_t n_vertices = m->n_vertices;
1179 const cs_lnum_t n_cells = m->n_cells;
1180 const cs_lnum_t n_b_faces = m->n_b_faces;
1181 const cs_lnum_t n_i_faces = m->n_i_faces;
1182 const int dim = m->dim;
1183 const cs_real_3_t *restrict vtx_coord
1184 = (const cs_real_3_t *restrict)m->vtx_coord;
1185 const cs_real_3_t *restrict cell_cen
1186 = (const cs_real_3_t *restrict)mq->cell_cen;
1187 const cs_real_3_t *restrict b_face_cog
1188 = (const cs_real_3_t *restrict)mq->b_face_cog;
1189
1190 BFT_MALLOC(vtx_counter, n_vertices, cs_real_t);
1191 BFT_MALLOC(vtx_interior_indicator, n_vertices, bool);
1192
1193 for (cs_lnum_t v_id = 0; v_id < n_vertices; v_id++) {
1194
1195 vtx_counter[v_id] = 0.;
1196 vtx_interior_indicator[v_id] = true;
1197
1198 for (int i = 0; i < dim; i++)
1199 disp_proj[v_id][i] = 0.;
1200
1201 }
1202
1203 /* All nodes wich belong to a boundary face where the
1204 displacement is imposed (that is all faces except sliding BCs)
1205 are boundary nodes, the others are interior nodes. */
1206
1207 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1208
1209 if (ale_bc_type[face_id] != CS_ALE_SLIDING) {
1210
1211 for (cs_lnum_t j = m->b_face_vtx_idx[face_id];
1212 j < m->b_face_vtx_idx[face_id+1];
1213 j++) {
1214
1215 const cs_lnum_t vtx_id = m->b_face_vtx_lst[j];
1216 vtx_interior_indicator[vtx_id] = false;
1217
1218 } /* End of loop on vertices of the face */
1219
1220 } /* Not sliding */
1221
1222 } /* End of loop on border faces */
1223
1224 /* Interior face and nodes treatment */
1225
1226 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1227
1228 const cs_lnum_t cell_id1 = m->i_face_cells[face_id][0];
1229 const cs_lnum_t cell_id2 = m->i_face_cells[face_id][1];
1230 const cs_real_t dvol1 = 1. / mq->cell_vol[cell_id1];
1231 const cs_real_t dvol2 = 1. / mq->cell_vol[cell_id2];
1232 const cs_real_t dt_dvol1 = dt[cell_id1] * dvol1;
1233 const cs_real_t dt_dvol2 = dt[cell_id2] * dvol2;
1234
1235 if (cell_id1 < n_cells) { /* Test to take into account face only once */
1236
1237 for (cs_lnum_t j = m->i_face_vtx_idx[face_id];
1238 j < m->i_face_vtx_idx[face_id+1];
1239 j++) {
1240
1241 /* Get the vertex id */
1242
1243 const cs_lnum_t vtx_id = m->i_face_vtx_lst[j];
1244
1245 if (vtx_interior_indicator[vtx_id]) {
1246
1247 /* Get the vector from the cell center to the node */
1248
1249 cs_real_3_t cen1_node, cen2_node;
1250 for (int i = 0; i < 3; i++) {
1251 cen1_node[i] = vtx_coord[vtx_id][i] - cell_cen[cell_id1][i];
1252 cen2_node[i] = vtx_coord[vtx_id][i] - cell_cen[cell_id2][i];
1253 }
1254
1255 for (int i = 0; i < 3; i++) {
1256 disp_proj[vtx_id][i] +=
1257 dt_dvol1*(meshv[cell_id1][i] + gradm[cell_id1][i][0]*cen1_node[0]
1258 + gradm[cell_id1][i][1]*cen1_node[1]
1259 + gradm[cell_id1][i][2]*cen1_node[2])
1260
1261 + dt_dvol2*(meshv[cell_id2][i] + gradm[cell_id2][i][0]*cen2_node[0]
1262 + gradm[cell_id2][i][1]*cen2_node[1]
1263 + gradm[cell_id2][i][2]*cen2_node[2]);
1264 }
1265
1266 vtx_counter[vtx_id] += dvol1 + dvol2;
1267
1268 } /* End of Interior nodes */
1269
1270 }
1271
1272 }
1273
1274 } /* End of loop on internal faces */
1275
1276 /* Border face treatment.
1277 only border face contribution */
1278
1279 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1280
1281 const cs_lnum_t cell_id = m->b_face_cells[face_id];
1282
1283 for (cs_lnum_t j = m->b_face_vtx_idx[face_id];
1284 j < m->b_face_vtx_idx[face_id+1]
1285 && ale_bc_type[face_id] != CS_ALE_SLIDING; j++) {
1286
1287 const cs_lnum_t vtx_id = m->b_face_vtx_lst[j];
1288
1289 if (!vtx_interior_indicator[vtx_id]) {
1290
1291 /* Get the vector from the face center to the node*/
1292
1293 cs_real_3_t face_node;
1294 for (int i = 0; i < 3; i++)
1295 face_node[i] = vtx_coord[vtx_id][i] - b_face_cog[face_id][i];
1296
1297 /* 1st order extrapolation of the mesh velocity at the face center
1298 * to the node */
1299
1300 cs_real_3_t vel_node;
1301 for (int i = 0; i < 3; i++)
1302 vel_node[i] = claale[face_id][i]
1303 + gradm[cell_id][i][0]*face_node[0]
1304 + gradm[cell_id][i][1]*face_node[1]
1305 + gradm[cell_id][i][2]*face_node[2];
1306
1307 const cs_real_t dsurf = 1./mq->b_face_surf[face_id];
1308
1309 for (int i = 0; i < 3; i++)
1310 disp_proj[vtx_id][i] += dsurf * dt[cell_id] *
1311 (vel_node[i] + clbale[face_id][i][0]*meshv[cell_id][0]
1312 + clbale[face_id][i][1]*meshv[cell_id][1]
1313 + clbale[face_id][i][2]*meshv[cell_id][2]);
1314
1315 vtx_counter[vtx_id] += dsurf;
1316
1317 } /* End of boundary nodes */
1318
1319 } /* End of loop on vertices of the face */
1320
1321 } /* End of loop on border faces */
1322
1323 if (m->vtx_interfaces != NULL) {
1324
1325 cs_interface_set_sum(m->vtx_interfaces,
1326 n_vertices,
1327 3,
1328 true,
1329 CS_REAL_TYPE,
1330 disp_proj);
1331
1332 cs_interface_set_sum(m->vtx_interfaces,
1333 n_vertices,
1334 1,
1335 true,
1336 CS_REAL_TYPE,
1337 vtx_counter);
1338 }
1339
1340 for (cs_lnum_t v_id = 0; v_id < n_vertices; v_id++) {
1341 /* Might be null for a vertex belonging to sliding faces only */
1342 if (vtx_counter[v_id] > 0.) {
1343 for (int i = 0; i < dim; i++)
1344 disp_proj[v_id][i] /= vtx_counter[v_id];
1345 }
1346 }
1347
1348 /* If the boundary face IS a sliding face.
1349 We project the displacment paralelly to the face. */
1350
1351 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1352
1353 if (ale_bc_type[face_id] == CS_ALE_SLIDING) {
1354
1355 for (cs_lnum_t j = m->b_face_vtx_idx[face_id];
1356 j < m->b_face_vtx_idx[face_id+1]; j++) {
1357
1358 const cs_lnum_t vtx_id = m->b_face_vtx_lst[j];
1359
1360 disp_proj[vtx_id][0] = clbale[face_id][0][0]*disp_proj[vtx_id][0]
1361 + clbale[face_id][0][1]*disp_proj[vtx_id][1]
1362 + clbale[face_id][0][2]*disp_proj[vtx_id][2];
1363 disp_proj[vtx_id][1] = clbale[face_id][1][0]*disp_proj[vtx_id][0]
1364 + clbale[face_id][1][1]*disp_proj[vtx_id][1]
1365 + clbale[face_id][1][2]*disp_proj[vtx_id][2];
1366 disp_proj[vtx_id][2] = clbale[face_id][2][0]*disp_proj[vtx_id][0]
1367 + clbale[face_id][2][1]*disp_proj[vtx_id][1]
1368 + clbale[face_id][2][2]*disp_proj[vtx_id][2];
1369
1370 } /* End of loop on vertices of the face */
1371
1372 } /* Sliding condition */
1373
1374 } /* End of loop on border faces */
1375
1376 BFT_FREE(vtx_counter);
1377 BFT_FREE(vtx_interior_indicator);
1378 }
1379
1380 /*----------------------------------------------------------------------------*/
1381 /*!
1382 * \brief Update mesh in the ALE framework.
1383 *
1384 * \param[in] itrale number of the current ALE iteration
1385 */
1386 /*----------------------------------------------------------------------------*/
1387
1388 void
cs_ale_update_mesh(const int itrale)1389 cs_ale_update_mesh(const int itrale)
1390 {
1391 const cs_mesh_t *m = cs_glob_mesh;
1392 const int ndim = m->dim;
1393 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1394 const cs_lnum_t n_vertices = m->n_vertices;
1395 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
1396
1397 cs_real_3_t *vtx_coord = (cs_real_3_t *)m->vtx_coord;
1398 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1399 cs_time_step_t *ts = cs_get_glob_time_step();
1400
1401 /* Initialization */
1402 cs_var_cal_opt_t var_cal_opt;
1403 cs_field_get_key_struct(CS_F_(mesh_u), key_cal_opt_id, &var_cal_opt);
1404
1405 if (var_cal_opt.verbosity >= 1)
1406 bft_printf("\n ---------------------------------------------------"
1407 "---------\n\n"
1408 " Update mesh (ALE)\n"
1409 " =================\n\n");
1410
1411 /* Retrieving fields */
1412 cs_field_t *f_displ = cs_field_by_name("mesh_displacement");
1413 cs_real_3_t *disale = (cs_real_3_t *)(f_displ->val);
1414 cs_real_3_t *disala = (cs_real_3_t *)(f_displ->val_pre);
1415 cs_real_3_t *xyzno0 = (cs_real_3_t *)(cs_field_by_name("vtx_coord0")->val);
1416
1417 /* Update geometry */
1418 for (cs_lnum_t v_id = 0; v_id < n_vertices; v_id++) {
1419 for (cs_lnum_t idim = 0; idim < ndim; idim++) {
1420 vtx_coord[v_id][idim] = xyzno0[v_id][idim] + disale[v_id][idim];
1421 disala[v_id][idim] = vtx_coord[v_id][idim] - xyzno0[v_id][idim];
1422 }
1423 }
1424
1425 cs_ale_update_mesh_quantities(&(mq->min_vol), &(mq->max_vol), &(mq->tot_vol));
1426
1427 /* Abort at the end of the current time-step if there is a negative volume */
1428 if (mq->min_vol <= 0.)
1429 ts->nt_max = ts->nt_cur;
1430
1431 /* The mesh velocity is reverted to its initial value if the current time step
1432 is the initialization time step */
1433 if (itrale == 0) {
1434
1435 cs_field_t *f = cs_field_by_name("mesh_velocity");
1436 if (f->location_id == CS_MESH_LOCATION_VERTICES) {
1437
1438 for (cs_lnum_t v_id = 0; v_id < n_vertices; v_id++)
1439 for (int idim = 0; idim < ndim; idim++)
1440 f->val[3*v_id+idim] = f->val_pre[3*v_id+idim];
1441
1442 }
1443 else if (f->location_id == CS_MESH_LOCATION_CELLS) {
1444
1445 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
1446 for (int idim = 0; idim < ndim; idim++)
1447 f->val[3*cell_id+idim] = f->val_pre[3*cell_id+idim];
1448
1449 } /* Field located at cells */
1450
1451 } /* itrale == 0 */
1452
1453 }
1454
1455 /*----------------------------------------------------------------------------*/
1456 /*!
1457 * \brief Update ALE BCs for required for the fluid
1458 *
1459 * \param[out] ale_bc_type type of ALE bcs
1460 * \param[out] b_fluid_vel Fluid velocity at boundary faces
1461 */
1462 /*----------------------------------------------------------------------------*/
1463
1464 void
cs_ale_update_bcs(int * ale_bc_type,cs_real_3_t * b_fluid_vel)1465 cs_ale_update_bcs(int *ale_bc_type,
1466 cs_real_3_t *b_fluid_vel)
1467 {
1468 assert(cs_glob_ale == CS_ALE_CDO);
1469
1470 /* First update ALE BCs expect free surface that will be updated after */
1471 _update_bcs(cs_glob_domain, ale_bc_type, b_fluid_vel);
1472
1473 }
1474
1475 /*----------------------------------------------------------------------------*/
1476 /*!
1477 * \brief Solve a Poisson equation on the mesh velocity in ALE framework.
1478 *
1479 * It also updates the mesh displacement
1480 * so that it can be used to update mass fluxes (due to mesh displacement).
1481 *
1482 * \param[in] iterns Navier-Stokes iteration number
1483 * \param[in] impale Indicator for fixed node displacement
1484 * \param[in] ale_bc_type Type of boundary for ALE
1485 */
1486 /*----------------------------------------------------------------------------*/
1487
1488 void
cs_ale_solve_mesh_velocity(const int iterns,const int * impale,const int * ale_bc_type)1489 cs_ale_solve_mesh_velocity(const int iterns,
1490 const int *impale,
1491 const int *ale_bc_type)
1492 {
1493 if (cs_glob_ale == CS_ALE_LEGACY)
1494 _ale_solve_poisson_legacy(cs_glob_domain, iterns, impale, ale_bc_type);
1495
1496 else if (cs_glob_ale == CS_ALE_CDO)
1497 _ale_solve_poisson_cdo(cs_glob_domain, impale);
1498 }
1499
1500 /*----------------------------------------------------------------------------*/
1501 /*!
1502 * \brief Activate the mesh velocity solving with CDO
1503 */
1504 /*----------------------------------------------------------------------------*/
1505
1506 void
cs_ale_activate(void)1507 cs_ale_activate(void)
1508 {
1509 if (cs_ale_active)
1510 return;
1511
1512 cs_ale_active = true;
1513
1514 cs_domain_set_cdo_mode(cs_glob_domain, CS_DOMAIN_CDO_MODE_WITH_FV);
1515
1516 cs_equation_t *eq
1517 = cs_equation_add("mesh_velocity", /* equation name */
1518 "mesh_velocity", /* associated variable field name */
1519 CS_EQUATION_TYPE_PREDEFINED,
1520 3, /* dimension of the unknown */
1521 CS_PARAM_BC_HMG_NEUMANN); /* default boundary */
1522
1523 cs_equation_param_t *eqp = cs_equation_get_param(eq);
1524
1525 /* System to solve is SPD by construction */
1526 cs_equation_param_set(eqp, CS_EQKEY_ITSOL, "cg");
1527 cs_equation_param_set(eqp, CS_EQKEY_PRECOND, "jacobi");
1528
1529 cs_equation_param_set(eqp, CS_EQKEY_SPACE_SCHEME, "cdo_vb");
1530
1531 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_RESNORM_TYPE, "filtered");
1532
1533 /* BC settings */
1534 cs_equation_param_set(eqp, CS_EQKEY_BC_ENFORCEMENT, "algebraic");
1535
1536 }
1537
1538 /*----------------------------------------------------------------------------*/
1539 /*!
1540 * \brief Test if mesh velocity solving with CDO is activated
1541 *
1542 * \return true if mesh velocity solving with CDO is requested, false otherwise
1543 */
1544 /*----------------------------------------------------------------------------*/
1545
1546 bool
cs_ale_is_activated(void)1547 cs_ale_is_activated(void)
1548 {
1549 if (cs_ale_active)
1550 return true;
1551 else
1552 return false;
1553 }
1554
1555 /*----------------------------------------------------------------------------*/
1556 /*!
1557 * \brief Setup the equations solving the mesh velocity
1558 *
1559 * \param[in, out] domain pointer to a cs_domain_t structure
1560 */
1561 /*----------------------------------------------------------------------------*/
1562
1563 void
cs_ale_init_setup(cs_domain_t * domain)1564 cs_ale_init_setup(cs_domain_t *domain)
1565 {
1566 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
1567
1568 /* Mesh viscosity (iso or ortho)
1569 * TODO declare it before: add in activate, def here... */
1570 int dim = cs_field_by_name("mesh_viscosity")->dim;
1571
1572 cs_property_t *mesh_visc = cs_property_by_name("mesh_viscosity");
1573 if (mesh_visc == NULL) { /* Not already added */
1574 cs_property_type_t type = 0;
1575 if (dim == 1)
1576 type = CS_PROPERTY_ISO;
1577 else if (dim == 3)
1578 type = CS_PROPERTY_ORTHO;
1579 else if (dim == 6)
1580 type = CS_PROPERTY_ANISO_SYM;
1581 else if (dim == 9)
1582 type = CS_PROPERTY_ANISO;
1583 else
1584 bft_error(__FILE__, __LINE__, 0,
1585 "%s: Invalid dimension (=%d) for the mesh viscosity.\n",
1586 __func__, dim);
1587
1588 /* Add and define this property */
1589 mesh_visc = cs_property_add("mesh_viscosity", type);
1590 cs_property_def_by_field(mesh_visc, cs_field_by_name("mesh_viscosity"));
1591
1592 }
1593
1594 cs_var_cal_opt_t var_cal_opt;
1595 cs_field_get_key_struct(CS_F_(mesh_u), key_cal_opt_id, &var_cal_opt);
1596
1597 //FIXME should be done elsewhere
1598 cs_domain_set_output_param(domain,
1599 -1, /* restart frequency: Only at the end */
1600 cs_glob_log_frequency,
1601 var_cal_opt.verbosity);
1602
1603 cs_equation_param_t *eqp = cs_equation_param_by_name("mesh_velocity");
1604 assert(mesh_visc != NULL);
1605 cs_equation_add_diffusion(eqp, mesh_visc);
1606 }
1607
1608 /*----------------------------------------------------------------------------*/
1609 /*!
1610 * \brief Setup the equations solving the mesh velocity
1611 *
1612 * \param[in] domain pointer to a cs_domain_t structure
1613 */
1614 /*----------------------------------------------------------------------------*/
1615
1616 void
cs_ale_setup_boundaries(const cs_domain_t * domain)1617 cs_ale_setup_boundaries(const cs_domain_t *domain)
1618 {
1619 const cs_mesh_t *mesh = domain->mesh;
1620 const cs_lnum_t n_vertices = mesh->n_vertices;
1621
1622 cs_equation_param_t *eqp = cs_equation_param_by_name("mesh_velocity");
1623
1624 if (_cdo_bc == NULL) {
1625
1626 BFT_MALLOC(_cdo_bc, 1, cs_ale_cdo_bc_t);
1627
1628 BFT_MALLOC(_cdo_bc->vtx_values, 3*n_vertices, cs_real_t);
1629 memset(_cdo_bc->vtx_values, 0, 3*n_vertices*sizeof(cs_real_t));
1630
1631 _cdo_bc->n_selections = 0;
1632 _cdo_bc->n_vertices = NULL;
1633 _cdo_bc->vtx_select = NULL;
1634
1635 }
1636
1637 bool *vtag = NULL;
1638 BFT_MALLOC(vtag, n_vertices, bool);
1639
1640 for (int b_id = 0; b_id < domain->ale_boundaries->n_boundaries; b_id++) {
1641
1642 const int z_id = domain->ale_boundaries->zone_ids[b_id];
1643 const cs_zone_t *z = cs_boundary_zone_by_id(z_id);
1644
1645 switch(domain->ale_boundaries->types[b_id]) {
1646
1647 case CS_BOUNDARY_ALE_FIXED:
1648 {
1649 cs_real_t bc_value[3] = {0., 0., 0.};
1650
1651 cs_equation_add_bc_by_value(eqp,
1652 CS_PARAM_BC_HMG_DIRICHLET,
1653 z->name,
1654 bc_value);
1655 }
1656 break;
1657
1658 case CS_BOUNDARY_ALE_IMPOSED_VEL:
1659 cs_equation_add_bc_by_array(eqp,
1660 CS_PARAM_BC_DIRICHLET,
1661 z->name,
1662 cs_flag_primal_vtx,
1663 _cdo_bc->vtx_values,
1664 false, /* Do not transfer ownership */
1665 NULL); /* No index */
1666
1667 /* Add a new list of selected vertices. */
1668 _update_bc_list(mesh, z, vtag);
1669 break;
1670
1671 case CS_BOUNDARY_ALE_IMPOSED_DISP:
1672 cs_equation_add_bc_by_array(eqp,
1673 CS_PARAM_BC_DIRICHLET,
1674 z->name,
1675 cs_flag_primal_vtx,
1676 _cdo_bc->vtx_values,
1677 false, /* Do not transfer ownership */
1678 NULL); /* No index */
1679
1680 /* Add a new list of selected vertices. */
1681 _update_bc_list(mesh, z, vtag);
1682 break;
1683
1684 case CS_BOUNDARY_ALE_FREE_SURFACE:
1685 cs_equation_add_bc_by_array(eqp,
1686 CS_PARAM_BC_DIRICHLET,
1687 z->name,
1688 cs_flag_primal_vtx,
1689 _cdo_bc->vtx_values,
1690 false, /* Do not transfer ownership */
1691 NULL); /* No index */
1692
1693 /* Add a new list of selected vertices. */
1694 _update_bc_list(mesh, z, vtag);
1695 break;
1696
1697 case CS_BOUNDARY_ALE_SLIDING:
1698 cs_equation_add_sliding_condition(eqp,
1699 z->name);
1700 break;
1701
1702 default:
1703 bft_error(__FILE__, __LINE__, 0,
1704 _(" %s: Boundary for ALE not allowed %s."),
1705 __func__, z->name);
1706 }
1707
1708 }
1709
1710 /* Free memory */
1711 BFT_FREE(vtag);
1712 }
1713
1714 /*----------------------------------------------------------------------------*/
1715 /*!
1716 * \brief Finalize the setup stage for the equation of the mesh velocity
1717 *
1718 * \param[in, out] domain pointer to a cs_domain_t structure
1719 */
1720 /*----------------------------------------------------------------------------*/
1721
1722 void
cs_ale_finalize_setup(cs_domain_t * domain)1723 cs_ale_finalize_setup(cs_domain_t *domain)
1724 {
1725 const cs_mesh_t *m = domain->mesh;
1726
1727 if (_vtx_coord0 == NULL) {
1728
1729 BFT_MALLOC(_vtx_coord0, m->n_vertices, cs_real_3_t);
1730
1731 memcpy(_vtx_coord0, m->vtx_coord, m->n_vertices * sizeof(cs_real_3_t));
1732
1733 }
1734
1735 /* Set the boundary conditions */
1736 cs_ale_setup_boundaries(domain);
1737 }
1738
1739 /*----------------------------------------------------------------------------*/
1740 /*!
1741 * \brief Free the main structure related to the ALE mesh velocity solving
1742 */
1743 /*----------------------------------------------------------------------------*/
1744
1745 void
cs_ale_destroy_all(void)1746 cs_ale_destroy_all(void)
1747 {
1748 BFT_FREE(_vtx_coord0);
1749
1750 if (_cdo_bc != NULL) {
1751 BFT_FREE(_cdo_bc->vtx_values);
1752
1753 for (int i = 0; i < _cdo_bc->n_selections; i++)
1754 BFT_FREE(_cdo_bc->vtx_select[i]);
1755 BFT_FREE(_cdo_bc->vtx_select);
1756 BFT_FREE(_cdo_bc->n_vertices);
1757
1758 BFT_FREE(_cdo_bc);
1759 }
1760 return;
1761 }
1762
1763 /*----------------------------------------------------------------------------*/
1764
1765 END_C_DECLS
1766