1 /*============================================================================
2 * SYRTHES coupling
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <stdarg.h>
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <assert.h>
38 #include <math.h>
39
40 #if defined(HAVE_MPI)
41 #include <mpi.h>
42 #endif
43
44 /*----------------------------------------------------------------------------
45 * PLE library headers
46 *----------------------------------------------------------------------------*/
47
48 #include <ple_defs.h>
49 #include <ple_coupling.h>
50 #include <ple_locator.h>
51
52 /*----------------------------------------------------------------------------
53 * Local headers
54 *----------------------------------------------------------------------------*/
55
56 #include "bft_mem.h"
57 #include "bft_error.h"
58 #include "bft_printf.h"
59
60 #include "fvm_nodal.h"
61 #include "fvm_nodal_extract.h"
62 #include "fvm_nodal_project.h"
63
64 #if defined(HAVE_MPI)
65 #include "cs_coupling.h"
66 #endif
67
68 #include "cs_cf_thermo.h"
69 #include "cs_coupling.h"
70 #include "cs_field.h"
71 #include "cs_field_pointer.h"
72 #include "cs_ht_convert.h"
73 #include "cs_log.h"
74 #include "cs_mesh.h"
75 #include "cs_mesh_connect.h"
76 #include "cs_parall.h"
77 #include "cs_post.h"
78 #include "cs_parameters.h"
79 #include "cs_prototypes.h"
80 #include "cs_physical_model.h"
81 #include "cs_thermal_model.h"
82 #include "cs_timer_stats.h"
83
84 /*----------------------------------------------------------------------------
85 * Header for the current file
86 *----------------------------------------------------------------------------*/
87
88 #include "cs_syr_coupling.h"
89
90 /*----------------------------------------------------------------------------*/
91
92 BEGIN_C_DECLS
93
94 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
95
96 /*=============================================================================
97 * Local Macro Definitions
98 *============================================================================*/
99
100 const int cs_syr_coupling_tag = 'C'+'S'+'_'+'C'+'O'+'U'+'P'+'L'+'A'+'G'+'E';
101
102 /*=============================================================================
103 * Local Structure Definitions
104 *============================================================================*/
105
106 /* Structure associated with Syrthes coupling */
107
108 typedef struct {
109
110 /* Mesh-related members */
111
112 ple_locator_t *locator; /* Associated locator */
113
114 int elt_dim; /* Element dimension */
115 cs_lnum_t n_elts; /* Number of coupled elements */
116
117 fvm_nodal_t *elts; /* Coupled elements */
118
119 /* Saved arrays for post processing (float for reduced memory use) */
120
121 int post_mesh_id; /* 0 if post-processing is not active,
122 or post-processing mesh_id (< 0) */
123 cs_real_t *solid_temp; /* Solid temperature received
124 from SYRTHES */
125 float *flux; /* Flux (calculated) */
126 float *tfluid_tmp; /* Fluid temperature (points to flux in
127 transient stage where solid_temp and
128 fluid_temp are known, NULL once
129 flux is calculated) */
130
131 /* Saved array for volume coupling. Will be used to build source term. */
132
133 double *hvol; /* Volumetric exchange coefficient. */
134
135 } cs_syr_coupling_ent_t ;
136
137 /* Structure associated with Syrthes coupling */
138
139 typedef struct {
140
141 /* Mesh-related members */
142
143 int dim; /* Coupled mesh dimension */
144 int ref_axis; /* Axis for edge extraction */
145
146 char *syr_name; /* Application name */
147
148 int n_b_locations; /* Numbero of boundary locations */
149 int n_v_locations; /* Numbero of volume locations */
150 int *b_location_ids; /* Boundary location ids */
151 int *v_location_ids; /* Volume location ids */
152
153 cs_syr_coupling_ent_t *faces; /* Wall coupling structure */
154 cs_syr_coupling_ent_t *cells; /* Volume coupling structure */
155
156 bool allow_nearest; /* Allow nearest-neighbor
157 mapping beyond basic matching
158 tolerance */
159 float tolerance; /* Tolerance */
160 int verbosity; /* Verbosity level */
161 int visualization; /* Visualization output flag */
162
163 /* Communication-related members */
164
165 #if defined(HAVE_MPI)
166
167 MPI_Comm comm; /* Associated MPI communicator */
168
169 int n_syr_ranks; /* Number of associated SYRTHES ranks */
170 int syr_root_rank; /* First associated SYRTHES rank */
171
172 #endif
173
174 } cs_syr_coupling_t;
175
176 /*============================================================================
177 * Global variables
178 *============================================================================*/
179
180 static int _syr_n_couplings = 0;
181 static cs_syr_coupling_t **_syr_couplings = NULL;
182
183 /* Start and end (negative) numbers associated with
184 dedicated post processing meshes */
185
186 static int _syr_post_mesh_ext[2] = {0, 1};
187
188 static int _syr_coupling_conservativity = 0; /* No forcing by default */
189 static int _syr_coupling_implicit = 1;
190
191 /*============================================================================
192 * Private function definitions
193 *============================================================================*/
194
195 /*----------------------------------------------------------------------------
196 * Initialize communicator for Syrthes coupling
197 *
198 * parameters:
199 * syr_coupling <-> Syrthes coupling structure
200 * coupling_id <-- id of this coupling (for log file message)
201 *----------------------------------------------------------------------------*/
202
203 static void
_init_comm(cs_syr_coupling_t * syr_coupling,int coupling_id)204 _init_comm(cs_syr_coupling_t *syr_coupling,
205 int coupling_id)
206
207 {
208 #if defined(HAVE_MPI)
209
210 int mpi_flag = 0;
211 int local_range[2] = {-1, -1};
212 int distant_range[2] = {-1, -1};
213
214 MPI_Initialized(&mpi_flag);
215
216 if (mpi_flag == 0)
217 return;
218
219 bft_printf(_(" SYRTHES coupling %d: initializing MPI communication ... "),
220 coupling_id);
221 bft_printf_flush();
222
223 ple_coupling_mpi_intracomm_create(MPI_COMM_WORLD,
224 cs_glob_mpi_comm,
225 syr_coupling->syr_root_rank,
226 &(syr_coupling->comm),
227 local_range,
228 distant_range);
229
230 bft_printf(_("[ok]\n"));
231 bft_printf(_(" Local ranks = [%d..%d], distant ranks = [%d..%d].\n\n"),
232 local_range[0], local_range[1] - 1,
233 distant_range[0], distant_range[1] - 1);
234 bft_printf_flush();
235
236 syr_coupling->n_syr_ranks = distant_range[1] - distant_range[0];
237 syr_coupling->syr_root_rank = distant_range[0];
238
239 #endif
240 }
241
242 /*----------------------------------------------------------------------------
243 * Free communicator for Syrthes coupling
244 *
245 * parameters:
246 * syr_coupling <-> Syrthes coupling structure
247 *----------------------------------------------------------------------------*/
248
249 static void
_finalize_comm(cs_syr_coupling_t * syr_coupling)250 _finalize_comm(cs_syr_coupling_t *syr_coupling)
251 {
252 #if defined(HAVE_MPI)
253
254 if (syr_coupling == NULL)
255 return;
256
257 if (syr_coupling->comm != MPI_COMM_NULL) {
258 MPI_Comm_free(&(syr_coupling->comm));
259 syr_coupling->comm = MPI_COMM_NULL;
260 }
261
262 #endif
263 }
264
265 /*----------------------------------------------------------------------------
266 * Exchange synchronization messages between Code_Saturne and SYRTHES.
267 *
268 * parameters:
269 * syr_coupling <-- Syrthes coupling structure
270 * op_name_send <-- operation name to send, or NULL. Only the 32
271 * first characters are sent if the nae is longer.
272 * op_name_recv <-- operation name to receive, or NULL (size: 33)
273 *----------------------------------------------------------------------------*/
274
275 static void
_exchange_sync(cs_syr_coupling_t * syr_coupling,const char * op_name_send,char * op_name_recv)276 _exchange_sync(cs_syr_coupling_t *syr_coupling,
277 const char *op_name_send,
278 char *op_name_recv)
279 {
280 #if defined(HAVE_MPI)
281
282 if (cs_glob_rank_id < 1) {
283
284 MPI_Status status;
285
286 if (op_name_send != NULL) {
287
288 char _op_name_send[33];
289 strncpy(_op_name_send, op_name_send, 32);
290 _op_name_send[32] = '\0';
291
292 /* Exchange command messages */
293 if (op_name_recv != NULL) {
294 MPI_Sendrecv(_op_name_send, 32, MPI_CHAR,
295 syr_coupling->syr_root_rank, cs_syr_coupling_tag,
296 op_name_recv, 32, MPI_CHAR,
297 syr_coupling->syr_root_rank, cs_syr_coupling_tag,
298 syr_coupling->comm, &status);
299 }
300
301 else
302 MPI_Send(_op_name_send, 32, MPI_CHAR,
303 syr_coupling->syr_root_rank, cs_syr_coupling_tag,
304 syr_coupling->comm);
305
306 }
307 else if (op_name_recv != NULL) {
308 MPI_Recv(op_name_recv, 32, MPI_CHAR,
309 syr_coupling->syr_root_rank, cs_syr_coupling_tag,
310 syr_coupling->comm, &status);
311 }
312
313 }
314
315 if (op_name_recv != NULL && cs_glob_rank_id > -1) {
316 MPI_Bcast(op_name_recv, 32, MPI_CHAR, 0, cs_glob_mpi_comm);
317 op_name_recv[32] = '\0';
318 }
319
320 #endif
321 }
322
323 /*----------------------------------------------------------------------------
324 * Post process variables associated with Syrthes couplings
325 *
326 * parameters:
327 * coupling <-- Void pointer to SYRTHES coupling structure
328 * ts <-- time step status structure, or NULL
329 *----------------------------------------------------------------------------*/
330
331 static void
_cs_syr_coupling_post_function(void * coupling,const cs_time_step_t * ts)332 _cs_syr_coupling_post_function(void *coupling,
333 const cs_time_step_t *ts)
334 {
335 int type_id;
336
337 const cs_syr_coupling_t *syr_coupling = coupling;
338 cs_syr_coupling_ent_t *coupling_ent = NULL;
339
340 for (type_id = 0; type_id < 2; type_id++) {
341
342 if (type_id == 0)
343 coupling_ent = syr_coupling->faces;
344 else
345 coupling_ent = syr_coupling->cells;
346
347 if (coupling_ent != NULL) {
348
349 if (coupling_ent->post_mesh_id != 0) {
350
351 const cs_real_t *cell_temp = NULL;
352 const cs_real_t *face_temp = NULL;
353
354 if (type_id == 0)
355 face_temp = coupling_ent->solid_temp;
356 else
357 cell_temp = coupling_ent->solid_temp;
358
359 cs_post_write_var(coupling_ent->post_mesh_id,
360 CS_POST_WRITER_ALL_ASSOCIATED,
361 _("Solid T"),
362 1,
363 false,
364 false,
365 CS_POST_TYPE_cs_real_t,
366 cell_temp,
367 NULL,
368 face_temp,
369 ts);
370
371 if (type_id == 1)
372 cs_post_write_var(coupling_ent->post_mesh_id,
373 CS_POST_WRITER_ALL_ASSOCIATED,
374 _("Solid heat flux"),
375 1,
376 false,
377 false,
378 CS_POST_TYPE_float,
379 coupling_ent->flux,
380 NULL,
381 NULL,
382 ts);
383
384 }
385 }
386 }
387 }
388
389 /*----------------------------------------------------------------------------
390 * Initialize post-processing of a SYRTHES coupling
391 *
392 * parameters:
393 * syr_coupling <-- partially initialized SYRTHES coupling structure
394 * coupling_ent <-- associated coupling mesh entity
395 *----------------------------------------------------------------------------*/
396
397 static void
_post_init(cs_syr_coupling_t * syr_coupling,cs_syr_coupling_ent_t * coupling_ent)398 _post_init(cs_syr_coupling_t *syr_coupling,
399 cs_syr_coupling_ent_t *coupling_ent)
400 {
401 int dim_shift = 0;
402 int coupling_id = -1;
403
404 const int writer_id = -1;
405 const int writer_ids[] = {writer_id};
406
407 /* Determine coupling id */
408
409 for (coupling_id = 0;
410 ( coupling_id < _syr_n_couplings
411 && _syr_couplings[coupling_id] != syr_coupling);
412 coupling_id++);
413
414 /* Exit silently if associated writer is not available */
415
416 if (cs_post_writer_exists(writer_id) != true)
417 return;
418
419 int t_top_id
420 = cs_timer_stats_switch(cs_timer_stats_id_by_name("postprocessing_stage"));
421
422 /* Initialize post processing flag */
423
424 coupling_ent->post_mesh_id = cs_post_get_free_mesh_id();
425
426 /* Allocate arrays if not already present */
427
428 if (coupling_ent->n_elts > 0) {
429 if (coupling_ent->solid_temp == NULL) /* surface coupling */
430 BFT_MALLOC(coupling_ent->solid_temp, coupling_ent->n_elts, cs_real_t);
431 if (coupling_ent->elt_dim == syr_coupling->dim) { /* volume coupling */
432 if (coupling_ent->flux == NULL)
433 BFT_MALLOC(coupling_ent->flux, coupling_ent->n_elts, float);
434 }
435 }
436 coupling_ent->tfluid_tmp = NULL;
437
438 /* Associate external mesh description with post processing subsystem */
439
440 if (syr_coupling->dim == 2)
441 dim_shift = 1;
442
443 cs_post_define_existing_mesh(coupling_ent->post_mesh_id,
444 coupling_ent->elts,
445 dim_shift,
446 false,
447 false,
448 1,
449 writer_ids);
450
451 /* Register post processing function */
452
453 cs_post_add_time_dep_output(_cs_syr_coupling_post_function,
454 (void *)syr_coupling);
455
456 /* Update start and end (negative) numbers associated with
457 dedicated post processing meshes */
458
459 if (_syr_post_mesh_ext[0] == 0)
460 _syr_post_mesh_ext[0] = coupling_ent->post_mesh_id;
461
462 _syr_post_mesh_ext[1] = coupling_ent->post_mesh_id;
463
464 cs_timer_stats_switch(t_top_id);
465 }
466
467 /*----------------------------------------------------------------------------
468 * Check if coupling location is complete
469 *
470 * parameters:
471 * syr_coupling <-- partially initialized SYRTHES coupling structure
472 * coupling_ent <-- coupling entity
473 * n_ext --> number of unlocated points
474 * ext_syr --> 1 if SYRTHES has some unlocted elements, 0 otherwise
475 *
476 * returns:
477 * true if location is complete, false otherwise
478 *----------------------------------------------------------------------------*/
479
480 static bool
_is_location_complete(cs_syr_coupling_t * syr_coupling,cs_syr_coupling_ent_t * coupling_ent,cs_gnum_t * n_ext,bool * ext_syr)481 _is_location_complete(cs_syr_coupling_t *syr_coupling,
482 cs_syr_coupling_ent_t *coupling_ent,
483 cs_gnum_t *n_ext,
484 bool *ext_syr)
485 {
486 bool location_complete = true;
487
488 /* Check that all points are effectively located */
489
490 ple_lnum_t n_exterior = ple_locator_get_n_exterior(coupling_ent->locator);
491 *n_ext = n_exterior;
492
493 char op_name_send[32 + 1];
494 char op_name_recv[32 + 1];
495
496 cs_parall_counter(n_ext, 1);
497
498 if (*n_ext > 0) {
499 strcpy(op_name_send, "coupling:location:incomplete");
500 location_complete = false;
501 }
502 else
503 strcpy(op_name_send, "coupling:location:ok");
504
505 _exchange_sync(syr_coupling, op_name_send, op_name_recv);
506 if (!strcmp(op_name_recv, "coupling:location:incomplete")) {
507 location_complete = false;
508 *ext_syr = true;
509 }
510 else
511 *ext_syr = false;
512
513 return location_complete;
514 }
515
516 /*----------------------------------------------------------------------------
517 * Define nodal mesh for Syrthes coupling from selection criteria.
518 *
519 * parameters:
520 * syr_coupling <-- partially initialized SYRTHES coupling structure
521 * n_locations <-- number of associated locations
522 * location_ids <-- associated location ids
523 * elt_dim <-- element dimension
524 *
525 * returns:
526 * pointer to created Syrthes coupling entity helper structure
527 *----------------------------------------------------------------------------*/
528
529 static cs_syr_coupling_ent_t *
_create_coupled_ent(cs_syr_coupling_t * syr_coupling,int n_locations,int location_ids[],int elt_dim)530 _create_coupled_ent(cs_syr_coupling_t *syr_coupling,
531 int n_locations,
532 int location_ids[],
533 int elt_dim)
534 {
535 char *coupled_mesh_name = NULL;
536 bool ext_syr = false;
537 cs_lnum_t n_exterior = 0;
538 cs_gnum_t n_ext = 0;
539 cs_coord_t *elt_centers = NULL;
540 fvm_nodal_t *location_elts = NULL;
541 float *cs_to_syr_dist = NULL;
542 float *syr_to_cs_dist = NULL;
543
544 bool location_complete = false;
545 cs_syr_coupling_ent_t *coupling_ent = NULL;
546
547 int locator_options[PLE_LOCATOR_N_OPTIONS];
548 locator_options[PLE_LOCATOR_NUMBERING] = 1;
549
550 assert(syr_coupling != NULL);
551
552 /* Initialization */
553
554 BFT_MALLOC(coupling_ent, 1, cs_syr_coupling_ent_t);
555
556 coupling_ent->locator = NULL;
557 coupling_ent->elt_dim = elt_dim;
558
559 coupling_ent->n_elts = 0;
560 coupling_ent->elts = NULL;
561
562 coupling_ent->post_mesh_id = 0;
563 coupling_ent->solid_temp = NULL;
564 coupling_ent->flux = NULL;
565 coupling_ent->tfluid_tmp = NULL;
566
567 coupling_ent->hvol = NULL;
568
569 if (syr_coupling->verbosity > 0) {
570 bft_printf(_("\nExtracting coupled mesh ..."));
571 bft_printf_flush();
572 }
573
574 /* Select elements */
575
576 cs_lnum_t n_elts = 0;
577 cs_lnum_t *elt_list = NULL;
578
579 for (int l_i = 0; l_i < n_locations; l_i++)
580 n_elts += cs_mesh_location_get_n_elts(location_ids[l_i])[0];
581
582 BFT_MALLOC(elt_list, n_elts, cs_lnum_t);
583
584 n_elts = 0;
585 for (int l_i = 0; l_i < n_locations; l_i++) {
586 int loc_id = location_ids[l_i];
587 const cs_lnum_t n = cs_mesh_location_get_n_elts(loc_id)[0];
588 const cs_lnum_t *ids = cs_mesh_location_get_elt_list(loc_id);
589 if (ids != NULL) {
590 for (cs_lnum_t i = 0; i < n; i++)
591 elt_list[n_elts++] = ids[i] + 1;
592 }
593 else {
594 for (cs_lnum_t i = 0; i < n; i++)
595 elt_list[n_elts++] = i + 1;
596 }
597 }
598
599 /* Creation of a new nodal mesh from selected cells */
600
601 if (elt_dim == syr_coupling->dim) {
602
603 BFT_MALLOC(coupled_mesh_name,
604 strlen(_("SYRTHES %s cells"))
605 + strlen(syr_coupling->syr_name) + 1, char);
606 sprintf(coupled_mesh_name, _("SYRTHES %s cells"), syr_coupling->syr_name);
607
608 coupling_ent->n_elts = n_elts;
609
610 coupling_ent->elts
611 = cs_mesh_connect_cells_to_nodal(cs_glob_mesh,
612 coupled_mesh_name,
613 false,
614 coupling_ent->n_elts,
615 elt_list);
616
617 /* Allocate additional buffers */
618
619 BFT_MALLOC(coupling_ent->hvol, coupling_ent->n_elts, double);
620 BFT_MALLOC(coupling_ent->solid_temp, coupling_ent->n_elts, cs_real_t);
621
622 }
623
624 /* Creation of a new nodal mesh from selected border faces */
625
626 else if (elt_dim == syr_coupling->dim - 1) {
627
628 BFT_MALLOC(coupled_mesh_name,
629 strlen("SYRTHES faces") + strlen(syr_coupling->syr_name) + 1,
630 char);
631 sprintf(coupled_mesh_name, _("SYRTHES %s faces"), syr_coupling->syr_name);
632
633 coupling_ent->n_elts = n_elts;
634
635 coupling_ent->elts
636 = cs_mesh_connect_faces_to_nodal(cs_glob_mesh,
637 coupled_mesh_name,
638 false,
639 0,
640 coupling_ent->n_elts,
641 NULL,
642 elt_list);
643
644 }
645
646 BFT_FREE(elt_list);
647
648 BFT_FREE(coupled_mesh_name);
649
650 if (syr_coupling->verbosity > 0) {
651 bft_printf(" [ok]\n");
652 bft_printf_flush();
653 }
654
655 if (fvm_nodal_get_n_g_vertices(coupling_ent->elts) == 0)
656 bft_error(__FILE__, __LINE__, 0,
657 _(" Selected mesh locations:\n"
658 " lead to an empty mesh for SYRTHES coupling .\n"
659 " \"%s\"\n"),
660 syr_coupling->syr_name);
661
662 /* In case of 2D coupling, project coupled elements to 2D */
663
664 location_elts = coupling_ent->elts;
665
666 if (syr_coupling->dim == 2) {
667
668 double a[6];
669 cs_lnum_t n_errors = 0;
670
671 if (syr_coupling->verbosity > 0) {
672 bft_printf(_("Projecting the extracted mesh to 2D ..."));
673 bft_printf_flush();
674 }
675
676 fvm_nodal_project(coupling_ent->elts, syr_coupling->ref_axis, &n_errors);
677
678 if (n_errors > 0)
679 bft_error(__FILE__, __LINE__, 0,
680 _("Error projecting the extracted mesh."));
681
682 if (syr_coupling->verbosity > 0) {
683 bft_printf(_(" [ok]\n"));
684 bft_printf_flush();
685 }
686
687 location_elts = fvm_nodal_copy(coupling_ent->elts);
688
689 if (syr_coupling->ref_axis == 0) {
690 a[0] = 0.; a[1] = 1.; a[2] = 0.; a[3] = 0.; a[4] = 0.; a[5] = 1.;
691 }
692 else if (syr_coupling->ref_axis == 1) {
693 a[0] = 1.; a[1] = 0.; a[2] = 0.; a[3] = 0.; a[4] = 0.; a[5] = 1.;
694 }
695 else if (syr_coupling->ref_axis == 2) {
696 a[0] = 1.; a[1] = 0.; a[2] = 0.; a[3] = 0.; a[4] = 1.; a[5] = 0.;
697 }
698
699 fvm_nodal_project_coords(location_elts, a);
700 }
701
702 /* Element information */
703
704 if (syr_coupling->verbosity > 0) {
705 cs_gnum_t n_g_elts = coupling_ent->n_elts;
706 cs_parall_counter(&n_g_elts, 1);
707 bft_printf(_("\nExtracted mesh built of %llu elements.\n"),
708 (unsigned long long)n_g_elts);
709 bft_printf_flush();
710 }
711
712 /* Initialize post-processing */
713
714 /* Precaution: deactivate visualization for time-dependent meshes,
715 as this would require regenerating visualization at each time step */
716
717 if (cs_post_get_writer_time_dep(-1) != FVM_WRITER_FIXED_MESH)
718 syr_coupling->visualization = 0;
719
720 if (syr_coupling->visualization != 0)
721 _post_init(syr_coupling, coupling_ent);
722
723 /* Build and initialize associated locator */
724
725 if (syr_coupling->verbosity > 0) {
726 bft_printf(_("\nLocator structure and mesh creation ..."));
727 bft_printf_flush();
728 }
729
730 /* Retrieve coordinates using FVM functions rather than previous list and
731 coordinates, in case the extracted nodal mesh contains elements in a
732 different order (multiple element types) or coordinates are projected
733 in 2D. */
734
735 if (coupling_ent->n_elts > 0) {
736
737 if (syr_coupling->visualization != 0)
738 BFT_MALLOC(cs_to_syr_dist, coupling_ent->n_elts, float);
739
740 BFT_MALLOC(elt_centers,
741 coupling_ent->n_elts*syr_coupling->dim,
742 cs_coord_t);
743 fvm_nodal_get_element_centers(location_elts,
744 CS_INTERLACE,
745 coupling_ent->elt_dim,
746 elt_centers);
747 }
748
749 /* Locate entities */
750
751 #if defined(PLE_HAVE_MPI)
752 coupling_ent->locator = ple_locator_create(syr_coupling->comm,
753 syr_coupling->n_syr_ranks,
754 syr_coupling->syr_root_rank);
755 #else
756 coupling_ent->locator = ple_locator_create();
757 #endif
758
759 ple_locator_set_mesh(coupling_ent->locator,
760 location_elts,
761 locator_options,
762 0.,
763 syr_coupling->tolerance,
764 syr_coupling->dim,
765 coupling_ent->n_elts,
766 NULL,
767 NULL,
768 elt_centers,
769 cs_to_syr_dist,
770 cs_coupling_mesh_extents,
771 cs_coupling_point_in_mesh);
772
773 /* Check that all points are effectively located */
774
775 location_complete = _is_location_complete(syr_coupling,
776 coupling_ent,
777 &n_ext,
778 &ext_syr);
779
780 if (syr_coupling->allow_nearest) {
781
782 float tolerance = syr_coupling->tolerance;
783
784 while (location_complete == false) {
785
786 tolerance *= 4;
787
788 if (syr_coupling->verbosity > 0) {
789 bft_printf(_(" [failed]\n"));
790 if (n_ext > 0)
791 bft_printf(_(" %llu fluid mesh elements not located on solid mesh\n"),
792 (unsigned long long) n_ext);
793 if (ext_syr)
794 bft_printf(_(" Some solid mesh elements not located on fluid mesh\n"));
795 bft_printf(_("\n Extending search with tolerance factor %f..."),
796 tolerance);
797 bft_printf_flush();
798 }
799
800 ple_locator_extend_search(coupling_ent->locator,
801 location_elts,
802 locator_options,
803 0.,
804 tolerance,
805 coupling_ent->n_elts,
806 NULL,
807 NULL,
808 elt_centers,
809 cs_to_syr_dist,
810 cs_coupling_mesh_extents,
811 cs_coupling_point_in_mesh);
812
813 location_complete = _is_location_complete(syr_coupling,
814 coupling_ent,
815 &n_ext,
816 &ext_syr);
817
818 }
819
820 }
821
822 if (syr_coupling->verbosity > 0) {
823 bft_printf(_(" [ok]\n"));
824 bft_printf_flush();
825 }
826
827 if (location_elts != coupling_ent->elts)
828 fvm_nodal_destroy(location_elts);
829
830 if (elt_centers != NULL)
831 BFT_FREE(elt_centers);
832
833 if (syr_coupling->visualization != 0) {
834
835 cs_post_activate_writer(-1, 1);
836 cs_post_write_meshes(cs_glob_time_step);
837
838 const float *b_dist = NULL, *v_dist = NULL;
839
840 if (coupling_ent->elt_dim == syr_coupling->dim - 1)
841 b_dist = cs_to_syr_dist;
842 else if (coupling_ent->elt_dim == syr_coupling->dim)
843 v_dist = cs_to_syr_dist;
844
845 cs_post_write_var(coupling_ent->post_mesh_id,
846 CS_POST_WRITER_ALL_ASSOCIATED,
847 _("distance_to_solid"),
848 1,
849 false,
850 false, /* use_parent, */
851 CS_POST_TYPE_float,
852 v_dist,
853 NULL,
854 b_dist,
855 NULL); /* time-independent variable */
856
857 BFT_FREE(cs_to_syr_dist);
858
859 }
860
861 /* Post-process distances from SYRTHES points to Code_Saturne faces */
862
863 if (elt_dim == syr_coupling->dim - 1) {
864
865 cs_lnum_t n_dist_elts = ple_locator_get_n_dist_points(coupling_ent->locator);
866
867 BFT_MALLOC(syr_to_cs_dist, n_dist_elts, float);
868
869 ple_locator_exchange_point_var(coupling_ent->locator,
870 syr_to_cs_dist,
871 NULL,
872 NULL,
873 sizeof(float),
874 1,
875 1);
876
877 if ( syr_coupling->visualization != 0
878 && syr_coupling->allow_nearest == false) {
879
880 cs_lnum_t i;
881 int writer_ids[] = {-1};
882 int mesh_id = coupling_ent->post_mesh_id - 1;
883 cs_lnum_t *p_vtx_num = NULL;
884 fvm_io_num_t *vtx_io_num = NULL;
885 fvm_nodal_t *syr_points = fvm_nodal_create("SYRTHES face centers",
886 syr_coupling->dim);
887
888 BFT_MALLOC(p_vtx_num, n_dist_elts, cs_lnum_t);
889
890 for (i = 0; i < (cs_lnum_t)n_dist_elts; i++)
891 p_vtx_num[i] = i+1;
892
893 fvm_nodal_define_vertex_list(syr_points, n_dist_elts, p_vtx_num);
894 fvm_nodal_set_shared_vertices
895 (syr_points,
896 ple_locator_get_dist_coords(coupling_ent->locator));
897
898 if (cs_glob_n_ranks > 1) {
899
900 vtx_io_num = fvm_io_num_create_from_scan(n_dist_elts);
901
902 fvm_nodal_init_io_num(syr_points,
903 fvm_io_num_get_global_num(vtx_io_num),
904 0);
905
906 }
907
908 cs_post_define_existing_mesh(mesh_id,
909 syr_points,
910 0,
911 true,
912 false,
913 1,
914 writer_ids);
915
916 cs_post_activate_writer(-1, 1);
917 cs_post_write_meshes(cs_glob_time_step);
918
919 cs_post_write_vertex_var(mesh_id,
920 CS_POST_WRITER_ALL_ASSOCIATED,
921 _("distance_to_fluid"),
922 1,
923 false,
924 false, /* use parent */
925 CS_POST_TYPE_float,
926 syr_to_cs_dist,
927 NULL); /* time-independent variable */
928
929 cs_post_free_mesh(mesh_id);
930
931 if (cs_glob_n_ranks > 1)
932 fvm_io_num_destroy(vtx_io_num);
933
934 } /* Do post-processing */
935
936 BFT_FREE(syr_to_cs_dist);
937
938 }
939
940 if (n_ext) {
941
942 int i;
943 int writer_ids[] = {-1};
944 int mesh_id = cs_post_get_free_mesh_id();
945 cs_lnum_t *post_vtx_num = NULL;
946 cs_coord_t *exterior_coords = NULL;
947 cs_coord_t *el_list = NULL;
948 fvm_io_num_t *vtx_io_num = NULL;
949 fvm_nodal_t *ulck_points = fvm_nodal_create("unlocated elements (centers)",
950 3);
951 n_exterior = ple_locator_get_n_exterior(coupling_ent->locator);
952 const ple_lnum_t *exterior_list
953 = ple_locator_get_exterior_list(coupling_ent->locator);
954
955 BFT_MALLOC(post_vtx_num, n_exterior, cs_lnum_t);
956 BFT_MALLOC(exterior_coords, 3*n_exterior, cs_coord_t);
957 BFT_MALLOC(el_list,
958 coupling_ent->n_elts*3,
959 cs_coord_t);
960
961 fvm_nodal_get_element_centers(coupling_ent->elts,
962 CS_INTERLACE,
963 coupling_ent->elt_dim,
964 el_list);
965
966 for (i = 0; i < (cs_lnum_t)n_exterior; i++) {
967 post_vtx_num[i] = i+1;
968 if (exterior_list[i] > coupling_ent->n_elts)
969 bft_error(__FILE__, __LINE__, 0,
970 _("Error: invalid exterior elements selection."));
971 exterior_coords[3*i ] = el_list[3*(exterior_list[i]-1) ];
972 exterior_coords[3*i +1] = el_list[3*(exterior_list[i]-1) +1];
973 exterior_coords[3*i +2] = el_list[3*(exterior_list[i]-1) +2];
974 }
975
976 fvm_nodal_define_vertex_list(ulck_points,
977 (cs_lnum_t)n_exterior,
978 post_vtx_num);
979 fvm_nodal_set_shared_vertices
980 (ulck_points,
981 exterior_coords);
982
983 if (cs_glob_n_ranks > 1) {
984 vtx_io_num = fvm_io_num_create_from_scan(n_exterior);
985 fvm_nodal_init_io_num(ulck_points,
986 fvm_io_num_get_global_num(vtx_io_num),
987 0);
988 }
989
990 cs_post_define_existing_mesh(mesh_id,
991 ulck_points,
992 0,
993 true,
994 false,
995 1,
996 writer_ids);
997
998 cs_post_activate_writer(writer_ids[0], 1);
999 cs_post_write_meshes(cs_glob_time_step);
1000 cs_post_free_mesh(mesh_id);
1001
1002 if (cs_glob_n_ranks > 1)
1003 fvm_io_num_destroy(vtx_io_num);
1004
1005 BFT_FREE(el_list);
1006 BFT_FREE(exterior_coords);
1007
1008 cs_base_warn(__FILE__, __LINE__);
1009 bft_printf(_("Coupling with SYRTHES impossible:\n"
1010 "%llu element centers from mesh \"%s\"\n"
1011 "not located on SYRTHES mesh."),
1012 (unsigned long long)n_ext,
1013 fvm_nodal_get_name(coupling_ent->elts));
1014
1015 }
1016
1017 /* Ensure clean stop */
1018
1019 if (location_complete == false)
1020 cs_coupling_set_sync_flag(PLE_COUPLING_STOP);
1021
1022 return coupling_ent;
1023 }
1024
1025 /*----------------------------------------------------------------------------
1026 * Log timing info
1027 *----------------------------------------------------------------------------*/
1028
1029 static void
_all_comm_times(void)1030 _all_comm_times(void)
1031 {
1032 int coupl_id, ent_id;
1033 cs_syr_coupling_t *syr_coupling = NULL;
1034
1035 if (_syr_n_couplings == 0)
1036 return;
1037
1038 cs_log_printf(CS_LOG_PERFORMANCE, "\n");
1039 cs_log_separator(CS_LOG_PERFORMANCE);
1040
1041 cs_log_printf(CS_LOG_PERFORMANCE,
1042 _("\nSYRTHES coupling overheads\n"));
1043
1044 for (coupl_id = 0; coupl_id < _syr_n_couplings; coupl_id++) {
1045
1046 syr_coupling = _syr_couplings[coupl_id];
1047
1048 for (ent_id = 0; ent_id < 2; ent_id++) {
1049
1050 cs_syr_coupling_ent_t
1051 *ce = (ent_id == 0) ? syr_coupling->faces : syr_coupling->cells;
1052 const char *ent_type[] = {N_("surface"), N_("volume")};
1053
1054 if (ce != NULL) {
1055
1056 double location_wtime, exchange_wtime;
1057 double location_comm_wtime, exchange_comm_wtime;
1058
1059 if (syr_coupling->syr_name != NULL)
1060 cs_log_printf(CS_LOG_PERFORMANCE,
1061 _("\n %s (%s):\n\n"),
1062 syr_coupling->syr_name, _(ent_type[ent_id]));
1063 else
1064 cs_log_printf(CS_LOG_PERFORMANCE,
1065 _("\n coupling %d (%s):\n\n"),
1066 coupl_id, _(ent_type[ent_id]));
1067
1068 ple_locator_get_times(ce->locator,
1069 &location_wtime,
1070 NULL,
1071 &exchange_wtime,
1072 NULL);
1073
1074 ple_locator_get_comm_times(ce->locator,
1075 &location_comm_wtime,
1076 NULL,
1077 &exchange_comm_wtime,
1078 NULL);
1079
1080 cs_log_printf(CS_LOG_PERFORMANCE,
1081 _(" location time: %12.3f\n"
1082 " communication and wait: %12.3f\n"
1083 " variable exchange time: %12.3f\n"
1084 " communication and wait: %12.3f\n"),
1085 location_wtime, location_comm_wtime,
1086 exchange_wtime, exchange_comm_wtime);
1087
1088 }
1089
1090 }
1091
1092 }
1093 }
1094
1095 /*----------------------------------------------------------------------------
1096 * Destroy coupled entity helper structure.
1097 *
1098 * parameters:
1099 * coupling ent <-> pointer to structure pointer
1100 *----------------------------------------------------------------------------*/
1101
1102 static void
_destroy_coupled_ent(cs_syr_coupling_ent_t ** coupling_ent)1103 _destroy_coupled_ent(cs_syr_coupling_ent_t **coupling_ent)
1104 {
1105 cs_syr_coupling_ent_t *ce = *coupling_ent;
1106
1107 if (ce == NULL)
1108 return;
1109
1110 if (ce->locator != NULL)
1111 ce->locator = ple_locator_destroy(ce->locator);
1112
1113 if (ce->solid_temp != NULL)
1114 BFT_FREE(ce->solid_temp);
1115 if (ce->flux != NULL)
1116 BFT_FREE(ce->flux);
1117
1118 if (ce->hvol != NULL)
1119 BFT_FREE(ce->hvol);
1120
1121 if (ce->elts != NULL)
1122 ce->elts = fvm_nodal_destroy(ce->elts);
1123
1124 BFT_FREE(*coupling_ent);
1125 }
1126
1127 /*----------------------------------------------------------------------------
1128 * Update post-processing variables of a Syrthes coupling
1129 *
1130 * Note that only the solid temperature is output for surface coupling,
1131 * while the heat flux is also output for volume coupling.
1132 *
1133 * parameters:
1134 * coupling_ent <-- Syrthes coupling structure
1135 * step <-- 0: var = wall temperature
1136 * 1: var = fluid temperature
1137 * 2: var = exchange coefficient
1138 * var <-- Pointer to variable values
1139 *----------------------------------------------------------------------------*/
1140
1141 static void
_post_var_update(cs_syr_coupling_ent_t * coupling_ent,int step,const cs_real_t * var)1142 _post_var_update(cs_syr_coupling_ent_t *coupling_ent,
1143 int step,
1144 const cs_real_t *var)
1145 {
1146 cs_lnum_t n_elts, ii;
1147
1148 assert(coupling_ent != NULL);
1149
1150 if (coupling_ent->post_mesh_id == 0)
1151 return;
1152
1153 assert(coupling_ent->solid_temp != NULL);
1154 assert(step == 0 || coupling_ent->flux != NULL);
1155
1156 n_elts = coupling_ent->n_elts;
1157
1158 /* Allocate arrays */
1159
1160 switch(step) {
1161
1162 case 0:
1163 for (ii = 0; ii < n_elts; ii++)
1164 coupling_ent->solid_temp[ii] = var[ii];
1165 break;
1166
1167 case 1:
1168 coupling_ent->tfluid_tmp = coupling_ent->flux;
1169 for (ii = 0; ii < n_elts; ii++)
1170 coupling_ent->tfluid_tmp[ii] = var[ii];
1171 break;
1172
1173 case 2:
1174 assert(coupling_ent->tfluid_tmp == coupling_ent->flux);
1175 for (ii = 0; ii < n_elts; ii++)
1176 coupling_ent->flux[ii] = var[ii] * ( coupling_ent->solid_temp[ii]
1177 - coupling_ent->flux[ii]);
1178 coupling_ent->tfluid_tmp = NULL;
1179 break;
1180
1181 default:
1182 assert(0);
1183 }
1184
1185 }
1186
1187 /*----------------------------------------------------------------------------
1188 * Ensure conservativity thanks to a corrector coefficient computed by SYRTHES
1189 * SYRTHES computes a global flux for a given tfluid and hfluid field.
1190 * Code_Saturne sent before its computed global flux for this time step.
1191 *
1192 * parameters:
1193 * syr_coupling <-- SYRTHES coupling structure
1194 * coupl_face_ids <-- ids of coupled boundary faces
1195 *----------------------------------------------------------------------------*/
1196
1197 static void
_ensure_conservativity(cs_syr_coupling_t * syr_coupling,const cs_lnum_t coupl_face_ids[])1198 _ensure_conservativity(cs_syr_coupling_t *syr_coupling,
1199 const cs_lnum_t coupl_face_ids[])
1200 {
1201 cs_lnum_t ii, face_id;
1202
1203 double g_flux = 0.0, _flux = 0.0, coef = 0.0;
1204
1205 double *surf = cs_glob_mesh_quantities->b_face_surf;
1206 cs_syr_coupling_ent_t *coupling_ent = NULL;
1207
1208 /* Sanity checks */
1209
1210 assert(surf != NULL);
1211 assert(coupl_face_ids != NULL);
1212 assert(syr_coupling != NULL);
1213 coupling_ent = syr_coupling->faces;
1214 assert(coupling_ent != NULL);
1215
1216 /* Compute Code_Saturne's global flux */
1217
1218 for (ii = 0; ii < coupling_ent->n_elts; ii++) {
1219 face_id = coupl_face_ids[ii];
1220 _flux += coupling_ent->flux[ii] * surf[face_id];
1221 }
1222
1223 #if defined(HAVE_MPI)
1224 if (cs_glob_n_ranks > 1)
1225 MPI_Reduce(&_flux, &g_flux, 1, MPI_DOUBLE, MPI_SUM, 0, cs_glob_mpi_comm);
1226 #endif
1227
1228 if (cs_glob_n_ranks == 1)
1229 g_flux = _flux;
1230
1231 /* Send the computed global flux to SYRTHES and receive the corrector
1232 coefficient */
1233
1234 #if defined(HAVE_MPI)
1235 if (cs_glob_rank_id < 1) {
1236
1237 MPI_Status status;
1238
1239 /* Send global flux */
1240
1241 MPI_Send(&g_flux, 1, MPI_DOUBLE,
1242 syr_coupling->syr_root_rank,
1243 cs_syr_coupling_tag,
1244 syr_coupling->comm);
1245
1246 if (syr_coupling->verbosity > 1)
1247 bft_printf(_(" Global heat flux exchanged with SYRTHES in W: %5.3e\n"),
1248 g_flux);
1249
1250 /* Receive corrector coefficient */
1251
1252 MPI_Recv(&coef, 1, MPI_DOUBLE,
1253 syr_coupling->syr_root_rank,
1254 cs_syr_coupling_tag,
1255 syr_coupling->comm,
1256 &status);
1257
1258 }
1259 #endif
1260
1261 /* Print message */
1262
1263 if (syr_coupling->verbosity > 1)
1264 bft_printf(_(" Correction coefficient used to force conservativity during"
1265 " coupling with SYRTHES: %5.3e\n"), coef);
1266 }
1267
1268 /*----------------------------------------------------------------------------
1269 * Exchange location synchronization status
1270 *
1271 * parameters:
1272 * syr_coupling <-- SYRTHES coupling structure
1273 *
1274 * returns:
1275 * 0 in case of success, 1 otherwise
1276 *----------------------------------------------------------------------------*/
1277
1278 static int
_sync_after_location(cs_syr_coupling_t * syr_coupling)1279 _sync_after_location(cs_syr_coupling_t *syr_coupling)
1280 {
1281 int retval = 1;
1282
1283 char op_name_send[32 + 1];
1284 char op_name_recv[32 + 1];
1285
1286 /* Communication with SYRTHES */
1287 /*----------------------------*/
1288
1289 /* Ready to start time iterations */
1290
1291 strcpy(op_name_send, "coupling:start");
1292
1293 _exchange_sync(syr_coupling, op_name_send, op_name_recv);
1294
1295 if (!strcmp(op_name_recv, "coupling:error:location")) {
1296
1297 cs_coupling_set_sync_flag(PLE_COUPLING_STOP);
1298
1299 cs_base_warn(__FILE__, __LINE__);
1300
1301 bft_printf(_(" Message received from SYRTHES: \"%s\"\n"
1302 " indicates meshes have not been matched correctly.\n\n"
1303 " The calculation will not run.\n\n"),
1304 op_name_recv);
1305
1306 }
1307 else if (strcmp(op_name_recv, "coupling:start"))
1308 bft_error(__FILE__, __LINE__, 0,
1309 _(" Message received from SYRTHES: \"%s\"\n"
1310 " indicates an error or is unexpected."),
1311 op_name_recv);
1312
1313 else
1314 retval = 0;
1315
1316 return retval;
1317 }
1318
1319 /*----------------------------------------------------------------------------
1320 * Get pointer to SYRTHES coupling.
1321 *
1322 * parameters:
1323 * coupling_id <-- Id (0 to n-1) of SYRTHES coupling
1324 *
1325 * returns:
1326 * pointer to SYRTHES coupling structure
1327 *----------------------------------------------------------------------------*/
1328
1329 static cs_syr_coupling_t *
_syr_coupling_by_id(int coupling_id)1330 _syr_coupling_by_id(int coupling_id)
1331 {
1332 cs_syr_coupling_t *retval = NULL;
1333
1334 if ( coupling_id > -1
1335 && coupling_id < _syr_n_couplings)
1336 retval = _syr_couplings[coupling_id];
1337
1338 return retval;
1339 }
1340
1341 /*----------------------------------------------------------------------------
1342 * Create or redefine a syr_coupling_t structure.
1343 *
1344 * If a structure is redefined, associated locations are reset.
1345 *
1346 * parameters:
1347 * dim <-- spatial mesh dimension
1348 * ref_axis <-- reference axis
1349 * syr_name <-- SYRTHES application name
1350 * allow_nonmatching <-- nearest-neighbor search for non-matching faces flag
1351 * tolerance <-- addition to local extents of each element
1352 * extent = base_extent * (1 + tolerance)
1353 * verbosity <-- verbosity level
1354 * visualization <-- visualization output flag
1355 *----------------------------------------------------------------------------*/
1356
1357 static cs_syr_coupling_t *
_syr_coupling_define(int dim,int ref_axis,const char * syr_name,bool allow_nonmatching,float tolerance,int verbosity,int visualization)1358 _syr_coupling_define(int dim,
1359 int ref_axis,
1360 const char *syr_name,
1361 bool allow_nonmatching,
1362 float tolerance,
1363 int verbosity,
1364 int visualization)
1365 {
1366 cs_syr_coupling_t *syr_coupling = NULL;
1367
1368 /* Search in existing couplings */
1369
1370 for (int i = 0;
1371 i < _syr_n_couplings && syr_coupling == NULL;
1372 i++) {
1373
1374 if (strcmp(_syr_couplings[i]->syr_name, syr_name) == 0) {
1375 syr_coupling = _syr_couplings[i];
1376
1377 BFT_FREE(syr_coupling->syr_name);
1378 BFT_FREE(syr_coupling->b_location_ids);
1379 BFT_FREE(syr_coupling->v_location_ids);
1380
1381 assert(syr_coupling->faces == NULL); /* not built yet at this stage */
1382 assert(syr_coupling->cells == NULL);
1383 }
1384 }
1385
1386 /* Allocate _cs_syr_coupling_t structure */
1387
1388 if (syr_coupling == NULL) {
1389 BFT_REALLOC(_syr_couplings,
1390 _syr_n_couplings + 1, cs_syr_coupling_t *);
1391 BFT_MALLOC(syr_coupling, 1, cs_syr_coupling_t);
1392
1393 _syr_couplings[_syr_n_couplings] = syr_coupling;
1394
1395 _syr_n_couplings++;
1396 }
1397
1398 syr_coupling->dim = dim;
1399 syr_coupling->ref_axis = ref_axis;
1400
1401 syr_coupling->syr_name = NULL;
1402
1403 if (syr_name != NULL) {
1404 BFT_MALLOC(syr_coupling->syr_name, strlen(syr_name) + 1, char);
1405 strcpy(syr_coupling->syr_name, syr_name);
1406 }
1407 else {
1408 BFT_MALLOC(syr_coupling->syr_name, 1, char);
1409 syr_coupling->syr_name[0] = '\0';
1410 }
1411
1412 /* Selection criteria */
1413
1414 syr_coupling->n_b_locations = 0;
1415 syr_coupling->n_v_locations = 0;
1416 syr_coupling->b_location_ids = NULL;
1417 syr_coupling->v_location_ids = NULL;
1418
1419 syr_coupling->faces = NULL;
1420 syr_coupling->cells = NULL;
1421
1422 syr_coupling->allow_nearest = allow_nonmatching;
1423 syr_coupling->tolerance = tolerance;
1424 syr_coupling->verbosity = verbosity;
1425 syr_coupling->visualization = visualization;
1426
1427 /* Initialize communicators */
1428
1429 #if defined(HAVE_MPI)
1430
1431 syr_coupling->comm = MPI_COMM_NULL;
1432 syr_coupling->n_syr_ranks = 0;
1433 syr_coupling->syr_root_rank = -1;
1434
1435 #endif
1436
1437 return syr_coupling;
1438 }
1439
1440 /*----------------------------------------------------------------------------
1441 * Add a mesh location to a syr_coupling_t structure.
1442 *
1443 * parameters:
1444 * syr_coupling <-- SYRTHES coupling structure
1445 * location_id <-- id of mesh location to add (boundary faces or cells)
1446 *----------------------------------------------------------------------------*/
1447
1448 static void
_add_mesh_location(cs_syr_coupling_t * syr_coupling,int location_id)1449 _add_mesh_location(cs_syr_coupling_t *syr_coupling,
1450 int location_id)
1451 {
1452 cs_mesh_location_type_t l_type = cs_mesh_location_get_type(location_id);
1453
1454 if (l_type == CS_MESH_LOCATION_BOUNDARY_FACES) {
1455 int i = syr_coupling->n_b_locations;
1456 syr_coupling->n_b_locations += 1;
1457 BFT_REALLOC(syr_coupling->b_location_ids, syr_coupling->n_b_locations, int);
1458
1459 syr_coupling->b_location_ids[i] = location_id;
1460 }
1461
1462 else if (l_type == CS_MESH_LOCATION_CELLS) {
1463 int i = syr_coupling->n_v_locations;
1464 syr_coupling->n_v_locations += 1;
1465 BFT_REALLOC(syr_coupling->v_location_ids, syr_coupling->n_v_locations, int);
1466
1467 syr_coupling->v_location_ids[i] = location_id;
1468 }
1469 }
1470
1471 /*----------------------------------------------------------------------------
1472 * Destroy cs_syr_coupling_t structures
1473 *----------------------------------------------------------------------------*/
1474
1475 static void
_syr_coupling_all_destroy(void)1476 _syr_coupling_all_destroy(void)
1477 {
1478 cs_lnum_t i_coupl;
1479 cs_syr_coupling_t *syr_coupling = NULL;
1480
1481 if (_syr_n_couplings == 0)
1482 return;
1483
1484 _all_comm_times();
1485
1486 for (i_coupl = 0; i_coupl < _syr_n_couplings; i_coupl++) {
1487
1488 syr_coupling = _syr_couplings[i_coupl];
1489
1490 /* Free _cs_syr_coupling structure */
1491
1492 BFT_FREE(syr_coupling->syr_name);
1493 BFT_FREE(syr_coupling->b_location_ids);
1494 BFT_FREE(syr_coupling->v_location_ids);
1495
1496 if (syr_coupling->faces != NULL)
1497 _destroy_coupled_ent(&(syr_coupling->faces));
1498 if (syr_coupling->cells != NULL)
1499 _destroy_coupled_ent(&(syr_coupling->cells));
1500
1501 /* Close communication */
1502
1503 _finalize_comm(syr_coupling);
1504
1505 BFT_FREE(syr_coupling);
1506
1507 } /* End of loop on _syr_couplings */
1508
1509 _syr_n_couplings = 0;
1510 BFT_FREE(_syr_couplings);
1511
1512 bft_printf(_("\nStructures associated with SYRTHES coupling freed.\n"));
1513 bft_printf_flush();
1514 }
1515
1516 /*----------------------------------------------------------------------------
1517 * Get name of SYRTHES coupling.
1518 *
1519 * parameters:
1520 * syr_coupling <-- SYRTHES coupling structure
1521 *
1522 * returns:
1523 * pointer to SYRTHES coupling name
1524 *----------------------------------------------------------------------------*/
1525
1526 static const char *
_syr_coupling_get_name(cs_syr_coupling_t * syr_coupling)1527 _syr_coupling_get_name(cs_syr_coupling_t *syr_coupling)
1528 {
1529 const char *retval = cs_empty_string;
1530
1531 if (syr_coupling->syr_name != NULL)
1532 retval = syr_coupling->syr_name;
1533
1534 return retval;
1535 }
1536
1537 /*----------------------------------------------------------------------------
1538 * Initialize communicator for SYRTHES coupling
1539 *
1540 * parameters:
1541 * syr_coupling <-> SYRTHES coupling structure
1542 * coupling_id <-- id of this coupling (for log file message)
1543 * syr_root_rank <-- SYRTHES root rank
1544 * n_syr_ranks <-- Number of ranks associated with SYRTHES
1545 *----------------------------------------------------------------------------*/
1546
1547 static void
_syr_coupling_init_comm(cs_syr_coupling_t * syr_coupling,int coupling_id,int syr_root_rank,int n_syr_ranks)1548 _syr_coupling_init_comm(cs_syr_coupling_t *syr_coupling,
1549 int coupling_id,
1550 int syr_root_rank,
1551 int n_syr_ranks)
1552 {
1553 #if defined(HAVE_MPI)
1554
1555 char volume_flag = ' ';
1556 char boundary_flag = ' ';
1557 char conservativity_flag = '1';
1558 char allow_nearest_flag = '1';
1559 char op_name_send[32 + 1];
1560 char op_name_recv[32 + 1];
1561
1562 syr_coupling->n_syr_ranks = n_syr_ranks;
1563 syr_coupling->syr_root_rank = syr_root_rank;
1564
1565 _init_comm(syr_coupling, coupling_id);
1566
1567 /* Exchange coupling options */
1568
1569 if (syr_coupling->n_b_locations > 0)
1570 boundary_flag = 'b';
1571 if (syr_coupling->n_v_locations > 0)
1572 volume_flag = 'v';
1573 if (_syr_coupling_conservativity == 0)
1574 conservativity_flag = '0';
1575 if (syr_coupling->allow_nearest == false)
1576 allow_nearest_flag = '0';
1577
1578 snprintf(op_name_send, 32, "coupling:type:%c%c%c \2\2%c(%6.2g)",
1579 boundary_flag, volume_flag, conservativity_flag,
1580 allow_nearest_flag, (double)syr_coupling->tolerance);
1581
1582 _exchange_sync(syr_coupling, op_name_send, op_name_recv);
1583
1584 if (strncmp(op_name_recv, op_name_send, 16))
1585 bft_error
1586 (__FILE__, __LINE__, 0,
1587 _("========================================================\n"
1588 " ** Incompatible SYRTHES coupling options:\n"
1589 " ------------------------------------------------\n"
1590 " Code_Saturne options: \"%s\"\n"
1591 " SYRTHES options: \"%s\"\n"
1592 "========================================================\n"),
1593 op_name_send, op_name_recv);
1594
1595 #endif
1596 }
1597
1598 /*----------------------------------------------------------------------------
1599 * Define coupled mesh and send it to SYRTHES
1600 *
1601 * Optional post-processing output is also built at this stage.
1602 *
1603 * parameters:
1604 * syr_coupling <-- SYRTHES coupling structure
1605 *----------------------------------------------------------------------------*/
1606
1607 static void
_syr_coupling_init_mesh(cs_syr_coupling_t * syr_coupling)1608 _syr_coupling_init_mesh(cs_syr_coupling_t *syr_coupling)
1609 {
1610 const cs_lnum_t verbosity = syr_coupling->verbosity;
1611
1612 if (verbosity > 0)
1613 bft_printf(_("\n ** Processing the mesh for SYRTHES coupling "
1614 "\"%s\"\n\n"),
1615 syr_coupling->syr_name);
1616
1617 /* Define coupled mesh */
1618
1619 assert(syr_coupling->dim == 3 || syr_coupling->dim == 2);
1620
1621 int match_flag = 0;
1622
1623 if (syr_coupling->n_b_locations > 0) {
1624 syr_coupling->faces = _create_coupled_ent(syr_coupling,
1625 syr_coupling->n_b_locations,
1626 syr_coupling->b_location_ids,
1627 syr_coupling->dim - 1);
1628 match_flag += _sync_after_location(syr_coupling);
1629 }
1630
1631 if (syr_coupling->n_v_locations > 0) {
1632 syr_coupling->cells = _create_coupled_ent(syr_coupling,
1633 syr_coupling->n_v_locations,
1634 syr_coupling->v_location_ids,
1635 syr_coupling->dim);
1636 match_flag += _sync_after_location(syr_coupling);
1637 }
1638
1639 /* Communication with SYRTHES */
1640 /*----------------------------*/
1641
1642 if (match_flag == 0 && verbosity > 0) {
1643 bft_printf(_("\n ** Mesh located for SYRTHES coupling \"%s\".\n\n"),
1644 syr_coupling->syr_name);
1645 bft_printf_flush();
1646 }
1647 }
1648
1649 /*----------------------------------------------------------------------------
1650 * Return 1 if this coupling is a surface coupling else 0
1651 *
1652 * parameters:
1653 * syr_coupling <-- SYRTHES coupling structure
1654 *
1655 * returns:
1656 * 1 or 0
1657 *----------------------------------------------------------------------------*/
1658
1659 static inline int
_syr_coupling_is_surf(const cs_syr_coupling_t * syr_coupling)1660 _syr_coupling_is_surf(const cs_syr_coupling_t *syr_coupling)
1661 {
1662 int retval = 0;
1663
1664 if (syr_coupling->n_b_locations > 0)
1665 retval = 1;
1666
1667 return retval;
1668 }
1669
1670 /*----------------------------------------------------------------------------
1671 * Return 1 if this coupling is a volume coupling else 0
1672 *
1673 * parameters:
1674 * syr_coupling <-- SYRTHES coupling structure
1675 *
1676 * returns:
1677 * 1 or 0
1678 *----------------------------------------------------------------------------*/
1679
1680 static inline int
_syr_coupling_is_vol(const cs_syr_coupling_t * syr_coupling)1681 _syr_coupling_is_vol(const cs_syr_coupling_t *syr_coupling)
1682 {
1683 int retval = 0;
1684
1685 if (syr_coupling->n_v_locations > 0)
1686 retval = 1;
1687
1688 return retval;
1689 }
1690
1691 /*----------------------------------------------------------------------------
1692 * Receive coupling variables from SYRTHES
1693 *
1694 * parameters:
1695 * syr_coupling <-- SYRTHES coupling structure
1696 * tsolid --> solid temperature
1697 * mode <-- 0: surface coupling; 1: volume coupling
1698 *----------------------------------------------------------------------------*/
1699
1700 static void
_syr_coupling_recv_tsolid(cs_syr_coupling_t * syr_coupling,cs_real_t tsolid[],int mode)1701 _syr_coupling_recv_tsolid(cs_syr_coupling_t *syr_coupling,
1702 cs_real_t tsolid[],
1703 int mode)
1704 {
1705 cs_syr_coupling_ent_t *coupling_ent = NULL;
1706
1707 assert(mode == 0 || mode == 1);
1708
1709 if (mode == 0)
1710 coupling_ent = syr_coupling->faces;
1711 else
1712 coupling_ent = syr_coupling->cells;
1713
1714 if (coupling_ent == NULL)
1715 return;
1716
1717 /* Receive data */
1718
1719 ple_locator_exchange_point_var(coupling_ent->locator,
1720 NULL,
1721 tsolid,
1722 NULL,
1723 sizeof(cs_real_t),
1724 1,
1725 0);
1726
1727 if (coupling_ent->n_elts > 0) {
1728 if (mode == 1) { /* Save tsolid for a future used
1729 in source term definition */
1730 cs_lnum_t i;
1731 assert(coupling_ent->solid_temp != NULL);
1732 for (i = 0; i < coupling_ent->n_elts; i++)
1733 coupling_ent->solid_temp[i] = tsolid[i];
1734 }
1735 else
1736 _post_var_update(coupling_ent, 0, tsolid);
1737 }
1738 }
1739
1740 /*----------------------------------------------------------------------------
1741 * Send coupling variables to SYRTHES
1742 *
1743 * parameters:
1744 * syr_coupling <-- SYRTHES coupling structure
1745 * cpl_elt_ids <-- ids of coupled elements
1746 * tf <-- fluid temperature
1747 * hf <-- fluid heat exchange coef. (numerical or user-defined)
1748 * mode <-- 0: surface coupling; 1: volume coupling
1749 *----------------------------------------------------------------------------*/
1750
1751 static void
_syr_coupling_send_tf_hf(cs_syr_coupling_t * syr_coupling,const cs_lnum_t cpl_elt_ids[],cs_real_t tf[],cs_real_t hf[],int mode)1752 _syr_coupling_send_tf_hf(cs_syr_coupling_t *syr_coupling,
1753 const cs_lnum_t cpl_elt_ids[],
1754 cs_real_t tf[],
1755 cs_real_t hf[],
1756 int mode)
1757 {
1758 cs_lnum_t ii;
1759
1760 cs_lnum_t n_dist = 0;
1761
1762 double *send_var = NULL;
1763 cs_syr_coupling_ent_t *coupling_ent = NULL;
1764
1765 const cs_lnum_t *dist_loc = NULL;
1766
1767 assert(mode == 0 || mode == 1);
1768
1769 if (mode == 0)
1770 coupling_ent = syr_coupling->faces;
1771 else
1772 coupling_ent = syr_coupling->cells;
1773
1774 if (coupling_ent == NULL)
1775 return;
1776
1777 n_dist = ple_locator_get_n_dist_points(coupling_ent->locator);
1778 dist_loc = ple_locator_get_dist_locations(coupling_ent->locator);
1779
1780 /* Prepare and send data */
1781
1782 BFT_MALLOC(send_var, n_dist*2, double);
1783
1784 for (ii = 0; ii < n_dist; ii++) {
1785 send_var[ii*2] = tf[dist_loc[ii] - 1];
1786 send_var[ii*2 + 1] = hf[dist_loc[ii] - 1];
1787 }
1788
1789 ple_locator_exchange_point_var(coupling_ent->locator,
1790 send_var,
1791 NULL,
1792 NULL,
1793 sizeof(double),
1794 2,
1795 0);
1796
1797 BFT_FREE(send_var);
1798
1799 if (mode == 1 && coupling_ent->n_elts > 0) {
1800
1801 _post_var_update(coupling_ent, 1, tf);
1802 _post_var_update(coupling_ent, 2, hf);
1803
1804 /* Saved hf for a future used in source term definition */
1805
1806 assert(coupling_ent->hvol != NULL);
1807 for (ii = 0; ii < coupling_ent->n_elts; ii++)
1808 coupling_ent->hvol[ii] = hf[ii];
1809
1810 }
1811
1812 /* Exchange flux and corrector coefficient to ensure conservativity */
1813
1814 if (_syr_coupling_conservativity > 0 && mode == 0)
1815 _ensure_conservativity(syr_coupling, cpl_elt_ids);
1816 }
1817
1818 /*----------------------------------------------------------------------------
1819 * Compute the explicit/implicit contribution to source terms in case of
1820 * volume coupling with SYRTHES
1821 *
1822 * parameters:
1823 * syr_coupling <-- SYRTHES coupling structure
1824 * tf <-- fluid temperature
1825 * ctbimp <-> implicit contribution
1826 * ctbexp <-> explicit contribution
1827 *----------------------------------------------------------------------------*/
1828
1829 static void
_syr_coupling_ts_contrib(const cs_syr_coupling_t * syr_coupling,const cs_real_t tf[],cs_real_t ctbimp[],cs_real_t ctbexp[])1830 _syr_coupling_ts_contrib(const cs_syr_coupling_t *syr_coupling,
1831 const cs_real_t tf[],
1832 cs_real_t ctbimp[],
1833 cs_real_t ctbexp[])
1834 {
1835 int i;
1836
1837 const double *hvol = NULL;
1838 const cs_real_t *solid_temp = NULL;
1839 cs_syr_coupling_ent_t *ent = NULL;
1840
1841 /* sanity checks */
1842
1843 assert(syr_coupling != NULL);
1844 assert(syr_coupling->cells != NULL);
1845
1846 ent = syr_coupling->cells;
1847 hvol = ent->hvol;
1848 solid_temp = ent->solid_temp;
1849
1850 assert(hvol != NULL || ent->n_elts == 0);
1851 assert(solid_temp != NULL || ent->n_elts == 0);
1852
1853 /* Compute contribution */
1854
1855 if (_syr_coupling_implicit == 0) { /* Explicit treatment */
1856
1857 for (i = 0; i < ent->n_elts; i++) {
1858 ctbexp[i] = -hvol[i] * (tf[i] - solid_temp[i]);
1859 ctbimp[i] = 0.0;
1860 }
1861
1862 }
1863 else { /* Implicit treatment */
1864
1865 for (i = 0; i < ent->n_elts; i++) {
1866 ctbexp[i] = hvol[i] * solid_temp[i];
1867 ctbimp[i] = hvol[i];
1868 }
1869
1870 } /* Test if implicit */
1871
1872 }
1873
1874 /*----------------------------------------------------------------------------
1875 * Print information on yet unmatched SYRTHES couplings.
1876 *
1877 * parameters:
1878 * n_unmatched <-- number of unmatched couplings
1879 * unmatched_ids <-- array of unmatched couplings
1880 *----------------------------------------------------------------------------*/
1881
1882 static void
_print_all_unmatched_syr(int n_unmatched,const int unmatched_ids[])1883 _print_all_unmatched_syr(int n_unmatched,
1884 const int unmatched_ids[])
1885 {
1886 /* Loop on defined SYRTHES instances */
1887
1888 for (int i = 0; i < n_unmatched; i++) {
1889
1890 cs_syr_coupling_t *syr_coupling
1891 = _syr_coupling_by_id(unmatched_ids[i]);
1892 const char *local_name = _syr_coupling_get_name(syr_coupling);
1893
1894 bft_printf(_(" SYRTHES coupling:\n"
1895 " coupling id: %d\n"
1896 " local name: \"%s\"\n\n"),
1897 i, local_name);
1898
1899 }
1900
1901 bft_printf_flush();
1902 }
1903
1904 #if defined(HAVE_MPI)
1905
1906 /*----------------------------------------------------------------------------
1907 * Initialize MPI SYRTHES couplings using MPI.
1908 *
1909 * This function may be called once all couplings have been defined,
1910 * and it will match defined couplings with available applications.
1911 *
1912 * parameters:
1913 * n_unmatched <-> pointer to number of unmatched couplings
1914 * unmatched_ids <-> pointer to array of unmatched couplings
1915 *----------------------------------------------------------------------------*/
1916
1917 static void
_init_all_mpi_syr(int * n_unmatched,int ** unmatched_ids)1918 _init_all_mpi_syr(int *n_unmatched,
1919 int **unmatched_ids)
1920 {
1921 int _n_unmatched = *n_unmatched;
1922 int *_unmatched_ids = *unmatched_ids;
1923
1924 const int n_couplings = _syr_n_couplings;
1925
1926 const ple_coupling_mpi_set_t *mpi_apps = cs_coupling_get_mpi_apps();
1927
1928 if (mpi_apps == NULL)
1929 return;
1930
1931 const int n_apps = ple_coupling_mpi_set_n_apps(mpi_apps);
1932
1933 /* Loop on applications */
1934
1935 for (int i = 0; i < n_apps; i++) {
1936
1937 ple_coupling_mpi_set_info_t ai = ple_coupling_mpi_set_get_info(mpi_apps, i);
1938
1939 if (strncmp(ai.app_type, "SYRTHES 4", 9) == 0) {
1940
1941 int match_queue_id = -1;
1942 int coupling_id = -1;
1943
1944 if (n_apps == 2 && n_couplings == 1 && _n_unmatched == 1) {
1945 match_queue_id = 0;
1946 coupling_id = 0;
1947 }
1948 else if (ai.app_name != NULL) {
1949 for (int j = 0; j < _n_unmatched; j++) {
1950 int k = _unmatched_ids[j];
1951 cs_syr_coupling_t *scpl = _syr_coupling_by_id(k);
1952 if (strcmp(ai.app_name, _syr_coupling_get_name(scpl)) == 0) {
1953 coupling_id = k;
1954 match_queue_id = j;
1955 break;
1956 }
1957 }
1958 }
1959
1960 if (coupling_id > -1) {
1961
1962 /* Remove from unmatched queue */
1963 _n_unmatched -= 1;
1964 for (int l = match_queue_id; l < _n_unmatched; l++)
1965 _unmatched_ids[l] = _unmatched_ids[l+1];
1966 if (_n_unmatched == 0)
1967 BFT_FREE(_unmatched_ids);
1968
1969 /* Set communicator */
1970 _syr_coupling_init_comm(_syr_coupling_by_id(coupling_id),
1971 coupling_id,
1972 ai.root_rank,
1973 ai.n_ranks);
1974
1975 /* Print matching info */
1976
1977 const char *syr_version = cs_empty_string;
1978 const char *local_name = cs_empty_string;
1979 const char *distant_name = cs_empty_string;
1980
1981 if (ai.app_name != NULL)
1982 local_name = ai.app_name;
1983 if (ai.app_type != NULL)
1984 syr_version = ai.app_type;
1985 if (ai.app_name != NULL)
1986 distant_name = ai.app_name;
1987
1988 bft_printf(_(" SYRTHES coupling:\n"
1989 " coupling id: %d\n"
1990 " version: \"%s\"\n"
1991 " local name: \"%s\"\n"
1992 " distant application name: \"%s\"\n"
1993 " MPI application id: %d\n"
1994 " MPI root rank: %d\n"
1995 " number of MPI ranks: %d\n\n"),
1996 coupling_id, syr_version, local_name, distant_name,
1997 i, ai.root_rank, ai.n_ranks);
1998 }
1999
2000 /* Note that a SYRTHES app may be present in the coupling set, but
2001 not coupled to the current code_saturne instance, so
2002 coupling_id < 0 here should not be reported as an error or
2003 complained about here. In case if missing matches, only the
2004 codes having defined and missing couplings should complain. */
2005
2006 }
2007
2008 } /* End of loop on applications */
2009
2010 bft_printf_flush();
2011
2012 /* Set return values */
2013
2014 *n_unmatched = _n_unmatched;
2015 *unmatched_ids = _unmatched_ids;
2016 }
2017
2018 /*----------------------------------------------------------------------------
2019 * Find name of single SYRTHES coupling using MPI.
2020 *
2021 * If no coupling or multiple couplings are present, the default cannot be
2022 * determined, so NULL is returned.
2023 *----------------------------------------------------------------------------*/
2024
2025 static const char *
_mpi_syr_default_name(void)2026 _mpi_syr_default_name(void)
2027 {
2028 const char *retval = NULL;
2029
2030 int n_syr_apps = 0;
2031
2032 const ple_coupling_mpi_set_t *mpi_apps = cs_coupling_get_mpi_apps();
2033
2034 if (mpi_apps == NULL)
2035 return NULL;
2036
2037 int n_apps = ple_coupling_mpi_set_n_apps(mpi_apps);
2038
2039 /* First pass to count available SYRTHES couplings */
2040
2041 for (int i = 0; i < n_apps; i++) {
2042 const ple_coupling_mpi_set_info_t
2043 ai = ple_coupling_mpi_set_get_info(mpi_apps, i);
2044 if (strncmp(ai.app_type, "SYRTHES 4", 9) == 0) {
2045 if (n_syr_apps == 0)
2046 retval = ai.app_name;
2047 else
2048 retval = NULL;
2049 n_syr_apps += 1;
2050 }
2051 }
2052
2053 return retval;
2054 }
2055
2056 #endif /* defined(HAVE_MPI) */
2057
2058 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
2059
2060 /*============================================================================
2061 * Public function definitions
2062 *============================================================================*/
2063
2064 /*----------------------------------------------------------------------------*/
2065 /*!
2066 * \brief Define new SYRTHES coupling.
2067 *
2068 * \param[in] syrthes_name matching SYRTHES application name
2069 * \param[in] boundary_criteria surface selection criteria, or NULL
2070 * \param[in] volume_criteria volume selection criteria, or NULL
2071 * \param[in] projection_axis x', 'y', or 'y' for 2D projection axis (case
2072 * independent), or ' ' for standard 3D coupling
2073 * \param[in] allow_nonmatching allow nearest-neighbor mapping where matching
2074 * within tolerance is not available (useful
2075 * when meshes have a different level of detail)
2076 * \param[in] tolerance addition to local extents of each element
2077 * extent = base_extent * (1 + tolerance)
2078 * \param[in] verbosity verbosity level
2079 * \param[in] visualization visualization output level (0 or 1)
2080 *
2081 * In the case of a single Code_Saturne and single SYRTHES instance, the
2082 * 'syrthes_name' argument is ignored, as there is only one matching
2083 * possibility.
2084 *
2085 * In case of multiple couplings, a coupling will be matched with available
2086 * SYRTHES instances based on the 'syrthes_name' argument.
2087 */
2088 /*----------------------------------------------------------------------------*/
2089
2090 void
cs_syr_coupling_define(const char * syrthes_name,const char * boundary_criteria,const char * volume_criteria,char projection_axis,bool allow_nonmatching,float tolerance,int verbosity,int visualization)2091 cs_syr_coupling_define(const char *syrthes_name,
2092 const char *boundary_criteria,
2093 const char *volume_criteria,
2094 char projection_axis,
2095 bool allow_nonmatching,
2096 float tolerance,
2097 int verbosity,
2098 int visualization)
2099 {
2100 int dim = 3;
2101 int ref_axis = -1;
2102
2103 switch (projection_axis) {
2104 case 'x':
2105 case 'X':
2106 dim = 2;
2107 ref_axis = 0;
2108 break;
2109 case 'y':
2110 case 'Y':
2111 dim = 2;
2112 ref_axis = 1;
2113 break;
2114 case 'z':
2115 case 'Z':
2116 dim = 2;
2117 ref_axis = 2;
2118 break;
2119 default:
2120 break;
2121 }
2122
2123 /* Ensure name is available */
2124
2125 #if defined(HAVE_MPI)
2126 if (syrthes_name == NULL)
2127 syrthes_name = _mpi_syr_default_name();
2128 #endif
2129
2130 if (syrthes_name == NULL)
2131 syrthes_name = cs_empty_string;
2132
2133 /* Define additional coupling */
2134
2135 cs_syr_coupling_t *syr_coupling = _syr_coupling_define(dim,
2136 ref_axis,
2137 syrthes_name,
2138 allow_nonmatching,
2139 tolerance,
2140 verbosity,
2141 visualization);
2142
2143 /* Add locations if done at that stage (deprecated) */
2144
2145 int n_locations = cs_mesh_location_n_locations();
2146
2147 const char *sel_criteria[2] = {boundary_criteria, volume_criteria};
2148 const char *type_name[2] = {"faces", "cells"};
2149 cs_mesh_location_type_t type_filter[2] = {CS_MESH_LOCATION_BOUNDARY_FACES,
2150 CS_MESH_LOCATION_CELLS};
2151
2152 for (int i = 0; i < 2; i++) {
2153
2154 if (sel_criteria[i] != NULL) {
2155 for (int j = 0; j < n_locations && sel_criteria[i] != NULL; j++) {
2156 cs_mesh_location_type_t l_type = cs_mesh_location_get_type(j);
2157 if (l_type & type_filter[i]) {
2158 const char *c = cs_mesh_location_get_selection_string(j);
2159 if (c != NULL) {
2160 if (strcmp(c, sel_criteria[i]) == 0) {
2161 _add_mesh_location(syr_coupling, j);
2162 sel_criteria[i] = NULL;
2163 }
2164 }
2165 }
2166 }
2167 }
2168
2169 if (sel_criteria[i] != NULL) {
2170
2171 char *name;
2172 size_t l = strlen(syrthes_name) + strlen(type_name[i]) + 2;
2173 BFT_MALLOC(name, l, char);
2174 snprintf(name, l, "%s_%s", syrthes_name, type_name[i]);
2175
2176 int j = cs_mesh_location_add(name, type_filter[i], sel_criteria[i]);
2177
2178 BFT_FREE(name);
2179
2180 _add_mesh_location(syr_coupling, j);
2181
2182 }
2183
2184 }
2185
2186 }
2187
2188 /*----------------------------------------------------------------------------*/
2189 /*!
2190 * \brief Associated a zone to a defined SYRTHES coupling.
2191 *
2192 * \param[in] syrthes_name matching SYRTHES application name
2193 * \param[in] z pointer to matching zone
2194 * (boundary or volume)
2195 */
2196 /*----------------------------------------------------------------------------*/
2197
2198 void
cs_syr_coupling_add_zone(const char * syrthes_name,const cs_zone_t * z)2199 cs_syr_coupling_add_zone(const char *syrthes_name,
2200 const cs_zone_t *z)
2201 {
2202 /* Ensure name is available */
2203
2204 #if defined(HAVE_MPI)
2205 if (syrthes_name == NULL)
2206 syrthes_name = _mpi_syr_default_name();
2207 #endif
2208
2209 if (syrthes_name == NULL)
2210 syrthes_name = cs_empty_string;
2211
2212 /* Search for matching name in existing couplings */
2213
2214 int n_couplings = _syr_n_couplings;
2215 bool match = false;
2216
2217 for (int i = 0; i < n_couplings; i++) {
2218
2219 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(i);
2220 const char *cmp_name = _syr_coupling_get_name(syr_coupling);
2221
2222 if (strcmp(syrthes_name, cmp_name) == 0) {
2223 _add_mesh_location(syr_coupling, z->location_id);
2224 match = true;
2225 break;
2226 }
2227
2228 }
2229
2230 if (match == false)
2231 bft_error(__FILE__, __LINE__, 0,
2232 _("%s: no defined SYRTHES coupling named \"%s\"."),
2233 __func__, syrthes_name);
2234 }
2235
2236 /*----------------------------------------------------------------------------
2237 * Initialize SYRTHES couplings.
2238 *
2239 * This function may be called once all couplings have been defined,
2240 * and it will match defined couplings with available applications.
2241 *----------------------------------------------------------------------------*/
2242
2243 void
cs_syr_coupling_all_init(void)2244 cs_syr_coupling_all_init(void)
2245 {
2246 int n_unmatched = _syr_n_couplings;
2247
2248 int *unmatched_ids;
2249 BFT_MALLOC(unmatched_ids, n_unmatched, int);
2250
2251 for (int i = 0; i < n_unmatched; i++)
2252 unmatched_ids[i] = i;
2253
2254 /* First try using MPI */
2255
2256 #if defined(HAVE_MPI)
2257
2258 if (n_unmatched > 0)
2259 _init_all_mpi_syr(&n_unmatched, &unmatched_ids);
2260
2261 #endif
2262
2263 if (n_unmatched > 0) {
2264
2265 bft_printf("Unmatched SYRTHES couplings:\n"
2266 "----------------------------\n\n");
2267
2268 _print_all_unmatched_syr(n_unmatched, unmatched_ids);
2269
2270 BFT_FREE(unmatched_ids);
2271
2272 bft_error(__FILE__, __LINE__, 0,
2273 _("At least 1 SYRTHES coupling was defined for which\n"
2274 "no communication with a SYRTHES instance is possible."));
2275 }
2276 }
2277
2278 /*----------------------------------------------------------------------------
2279 * Finalize all SYRTHES couplings.
2280 *----------------------------------------------------------------------------*/
2281
2282 void
cs_syr_coupling_all_finalize(void)2283 cs_syr_coupling_all_finalize(void)
2284 {
2285 _syr_coupling_all_destroy();
2286 }
2287
2288 /*----------------------------------------------------------------------------
2289 * Return number of SYRTHES couplings.
2290 *
2291 * return:
2292 * number of SYRTHES couplings defined
2293 *----------------------------------------------------------------------------*/
2294
2295 int
cs_syr_coupling_n_couplings(void)2296 cs_syr_coupling_n_couplings(void)
2297 {
2298 return _syr_n_couplings;
2299 }
2300
2301 /*----------------------------------------------------------------------------
2302 * Set conservativity forcing flag to True (1) or False (0) for all defined
2303 * SYRTHES couplings
2304 *
2305 * parameter:
2306 * flag <-- Conservativity forcing flag to set
2307 *----------------------------------------------------------------------------*/
2308
2309 void
cs_syr_coupling_set_conservativity(int flag)2310 cs_syr_coupling_set_conservativity(int flag)
2311 {
2312 assert(flag == 0 || flag == 1);
2313 _syr_coupling_conservativity = flag;
2314 }
2315
2316 /*----------------------------------------------------------------------------
2317 * Set explicit treatment for the source terms in SYRTHES volume couplings
2318 *----------------------------------------------------------------------------*/
2319
2320 void
cs_syr_coupling_set_explicit_treatment(void)2321 cs_syr_coupling_set_explicit_treatment(void)
2322 {
2323 _syr_coupling_implicit = 0;
2324 }
2325
2326 /*----------------------------------------------------------------------------*/
2327 /*!
2328 * \brief Log SYRTHES coupling setup information.
2329 */
2330 /*----------------------------------------------------------------------------*/
2331
2332 void
cs_syr_coupling_log_setup(void)2333 cs_syr_coupling_log_setup(void)
2334 {
2335 /* Get the number of SYRTHES couplings */
2336 int n_coupl = _syr_n_couplings;
2337 const int keysca = cs_field_key_id("scalar_id");
2338 const int kcpsyr = cs_field_key_id("syrthes_coupling");
2339 int icpsyr;
2340
2341 if (n_coupl >= 1) {
2342
2343 cs_log_printf
2344 (CS_LOG_SETUP,
2345 _("SYRTHES coupling\n"
2346 "----------------\n\n"
2347 " number of couplings: %d\n"),
2348 n_coupl);
2349
2350 int n_surf_coupl = 0, n_vol_coupl = 0, issurf, isvol;
2351
2352 for (int coupl_id = 0; coupl_id < n_coupl; coupl_id++) {
2353 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(coupl_id);
2354
2355 /* Add a new surface coupling if detected */
2356 issurf = _syr_coupling_is_surf(syr_coupling);
2357 n_surf_coupl += issurf;
2358
2359 /* Add a new volume coupling if detected */
2360 isvol = _syr_coupling_is_vol(syr_coupling);
2361 n_vol_coupl += isvol;
2362 }
2363
2364 cs_log_printf
2365 (CS_LOG_SETUP,
2366 _(" with %d surface coupling(s)\n"
2367 " with %d volume coupling(s)\n"),
2368 n_surf_coupl, n_vol_coupl);
2369
2370 cs_log_printf
2371 (CS_LOG_SETUP,
2372 _("\n"
2373 " Coupled scalars\n"
2374 "------------------------\n"
2375 " Scalar Number icpsyr\n"
2376 "------------------------\n"));
2377
2378
2379 for (int f_id = 0 ; f_id < cs_field_n_fields() ; f_id++) {
2380 cs_field_t *f = cs_field_by_id(f_id);
2381 if ((f->type & CS_FIELD_VARIABLE) || (f->type & CS_FIELD_USER)) {
2382 int ii = cs_field_get_key_int(f, keysca);
2383 if (ii > 0) {
2384 icpsyr = cs_field_get_key_int(f, kcpsyr);
2385 cs_log_printf
2386 (CS_LOG_SETUP,
2387 _(" %s %7d %7d\n"),cs_field_get_label(f),ii, icpsyr);
2388 }
2389 }
2390 }
2391 cs_log_printf
2392 (CS_LOG_SETUP,
2393 _("------------------------\n\n"
2394 " icpsyr = 0 or 1 (1: scalar coupled to SYRTHES)\n"));
2395 }
2396 }
2397
2398 /*----------------------------------------------------------------------------*/
2399 /*!
2400 * \brief Create coupled meshes and setup PLE locator for Syrthes couplings.
2401 */
2402 /*----------------------------------------------------------------------------*/
2403
2404 void
cs_syr_coupling_init_meshes(void)2405 cs_syr_coupling_init_meshes(void)
2406 {
2407 int n_coupl = _syr_n_couplings;
2408
2409 for (int coupl_id = 0; coupl_id < n_coupl; coupl_id++) {
2410 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(coupl_id);
2411 _syr_coupling_init_mesh(syr_coupling);
2412 }
2413 }
2414
2415 /*----------------------------------------------------------------------------*/
2416 /*!
2417 * \brief Check if the given SYRTHES coupling number is a surface couplings.
2418 *
2419 * \param[in] cpl_id matching SYRTHES coupling id
2420 *
2421 * \return 1 if the coupling includes the surface, 0 otherwise.
2422 */
2423 /*----------------------------------------------------------------------------*/
2424
2425 int
cs_syr_coupling_is_surf(int cpl_id)2426 cs_syr_coupling_is_surf(int cpl_id)
2427 {
2428 int retval = 0; /* Default initialization */
2429
2430 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
2431
2432 if (syr_coupling == NULL)
2433 bft_error(__FILE__, __LINE__, 0,
2434 _("SYRTHES coupling id %d impossible; "
2435 "there are %d couplings"),
2436 cpl_id, _syr_n_couplings);
2437
2438 retval = _syr_coupling_is_surf(syr_coupling);
2439
2440 return retval;
2441 }
2442
2443 /*----------------------------------------------------------------------------*/
2444 /*!
2445 * \brief Read boundary field/variable values relative to a SYRTHES coupling.
2446 *
2447 * \param[in] nvar number of variables
2448 * \param[in] bc_type boundary condition type
2449 * \param[in, out] icodcl boundary condition codes
2450 * \param[in, out] rcodcl boundary condition values
2451 */
2452 /*----------------------------------------------------------------------------*/
2453
2454 void
cs_syr_coupling_recv_boundary(int nvar,int bc_type[],int icodcl[],cs_real_t rcodcl[])2455 cs_syr_coupling_recv_boundary(int nvar,
2456 int bc_type[],
2457 int icodcl[],
2458 cs_real_t rcodcl[])
2459 {
2460 /* SYRTHES coupling: get wall temperature
2461 ====================================== */
2462
2463 const int kcpsyr = cs_field_key_id("syrthes_coupling");
2464
2465 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
2466 const cs_lnum_t n_vars = nvar; /* cast to cs_lnum_t because this
2467 is used in address computations here */
2468
2469 /* Get number of coupling cases */
2470
2471 int n_cpl = cs_syr_coupling_n_couplings();
2472
2473 /* Loop on fields, handling only those coupled with Syrthes */
2474
2475 int n_fields = cs_field_n_fields();
2476 for (int field_id = 0 ; field_id < n_fields; field_id++) {
2477
2478 cs_field_t *f = cs_field_by_id(field_id);
2479
2480 int icpsyr = 0;
2481 if (f->type & CS_FIELD_VARIABLE)
2482 icpsyr = cs_field_get_key_int(f, kcpsyr);
2483
2484 if (icpsyr < 1)
2485 continue;
2486
2487 /* Loop on couplings: get wall temperature array for each coupling
2488 and apply matching boundary condition. */
2489
2490 for (int cpl_id = 0; cpl_id < n_cpl; cpl_id++) {
2491
2492 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
2493 cs_syr_coupling_ent_t *coupling_ent = syr_coupling->faces;
2494
2495 if (coupling_ent == NULL) /* ignore if volume-only */
2496 continue;
2497
2498 cs_lnum_t n_cpl_faces = coupling_ent->n_elts;
2499
2500 /* Get list of coupled faces */
2501
2502 cs_lnum_t *f_ids;
2503 BFT_MALLOC(f_ids, n_cpl_faces, cs_lnum_t);
2504 fvm_nodal_get_parent_id(coupling_ent->elts,
2505 coupling_ent->elt_dim,
2506 f_ids);
2507
2508 /* Read wall temperature and interpolate if necessary */
2509
2510 cs_real_t *t_solid;
2511 BFT_MALLOC(t_solid, n_cpl_faces, cs_real_t);
2512 _syr_coupling_recv_tsolid(syr_coupling, t_solid, 0);
2513
2514 /* For scalars coupled with SYRTHES, prescribe a Dirichlet
2515 condition at coupled faces.
2516 For the time being, pass here only once, as only one scalar is
2517 coupled with SYRTHES.
2518 For the compressible module, solve in energy, but save the
2519 temperature separately, for BC's to be clearer. */
2520
2521 const int k_var_id = cs_field_key_id("variable_id");
2522 int var_id = cs_field_get_key_int(f, k_var_id) - 1;
2523
2524 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] >= 0) {
2525 if (f == CS_F_(e_tot)) {
2526 const cs_field_t *f_t_kelvin = CS_F_(t_kelvin);
2527 var_id = cs_field_get_key_int(f_t_kelvin, k_var_id);
2528 }
2529 else
2530 bft_error
2531 (__FILE__, __LINE__, 0,
2532 _("With the compressible module, only the \"total energy\"\n"
2533 "scalar field may be coupled with SYRTHES.\n"
2534 "Here, one tries to couple with the field \"%s\"."),
2535 f->name);
2536 }
2537
2538 int *_icodcl = icodcl + (var_id*n_b_faces);
2539 cs_real_t *_rcodcl1 = rcodcl + (var_id*n_b_faces);
2540 cs_real_t *_rcodcl2 = rcodcl + (n_b_faces*n_vars + var_id*n_b_faces);
2541 cs_real_t *_rcodcl3 = rcodcl + (2*n_b_faces*n_vars + var_id*n_b_faces);
2542
2543 for (cs_lnum_t i = 0; i < n_cpl_faces; i++) {
2544
2545 cs_lnum_t face_id = f_ids[i];
2546
2547 if ( _icodcl[face_id] != CS_INDEF
2548 && _icodcl[face_id] != CS_SMOOTHWALL
2549 && _icodcl[face_id] != CS_ROUGHWALL) {
2550 if (bc_type[face_id] == CS_SMOOTHWALL)
2551 _icodcl[face_id] = CS_SMOOTHWALL;
2552 else if (bc_type[face_id] == CS_ROUGHWALL)
2553 _icodcl[face_id] = CS_ROUGHWALL;
2554 }
2555
2556 _rcodcl1[face_id] = t_solid[i];
2557 _rcodcl2[face_id] = cs_math_infinite_r;
2558 _rcodcl3[face_id] = 0.;
2559
2560 }
2561
2562 /* Require temperature -> enthalpy conversion */
2563
2564 if (cs_glob_thermal_model->itherm == CS_THERMAL_MODEL_ENTHALPY) {
2565
2566 if (f == cs_thermal_model_field()) {
2567 for (cs_lnum_t i = 0; i < n_cpl_faces; i++) {
2568 cs_lnum_t face_id = f_ids[i];
2569 _icodcl[face_id] *= -1;
2570 }
2571 }
2572
2573 } /* End case for enthalpy */
2574
2575 BFT_FREE(f_ids);
2576 BFT_FREE(t_solid);
2577
2578 } /* End loop on couplings */
2579
2580 } /* End loop on fields */
2581 }
2582
2583 /*----------------------------------------------------------------------------*/
2584 /*!
2585 * \brief Send field/variable values relative to a SYRTHES coupling.
2586 *
2587 * \param[in] h_wall wall thermal exchange coefficient
2588 * \param[in] v_fluid near-wall fluid thermal variable
2589 */
2590 /*----------------------------------------------------------------------------*/
2591
2592 void
cs_syr_coupling_send_boundary(const cs_real_t h_wall[],cs_real_t v_fluid[])2593 cs_syr_coupling_send_boundary(const cs_real_t h_wall[],
2594 cs_real_t v_fluid[])
2595 {
2596 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
2597 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
2598
2599 /* Get number of coupling cases */
2600
2601 int n_cpl = cs_syr_coupling_n_couplings();
2602
2603 /* Check if we have a boundary coupling. */
2604
2605 bool have_boundary_cpl = false;
2606 for (int cpl_id = 0; cpl_id < n_cpl; cpl_id++) {
2607 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
2608 if (_syr_coupling_is_surf(syr_coupling)) {
2609 have_boundary_cpl = true;
2610 break;
2611 }
2612 }
2613
2614 if (! have_boundary_cpl)
2615 return;
2616
2617 /* Build arrays large enough for all cases */
2618
2619 cs_lnum_t *f_ids;
2620 cs_real_t *t_fluid, *h_cpl;
2621 BFT_MALLOC(f_ids, n_b_faces, cs_lnum_t);
2622 BFT_MALLOC(t_fluid, n_b_faces, cs_real_t);
2623 BFT_MALLOC(h_cpl, n_b_faces, cs_real_t);
2624
2625 /* Prepare conversion to temperature for enthalpy or energy
2626 (check for surface couplings to make sure it is needed,
2627 exit earlier otherwise) */
2628
2629 cs_real_t *wa = NULL;
2630 if (cs_glob_thermal_model->itherm == CS_THERMAL_MODEL_ENTHALPY) {
2631 BFT_MALLOC(wa, n_b_faces, cs_real_t);
2632 cs_ht_convert_h_to_t_faces(v_fluid, wa);
2633 }
2634 else if (cs_glob_thermal_model->itherm == CS_THERMAL_MODEL_TOTAL_ENERGY) {
2635 /* Epsilon sup for perfect gas at cells */
2636 BFT_MALLOC(wa, n_cells, cs_real_t);
2637 cs_cf_thermo_eps_sup(CS_F_(rho)->val, wa, n_cells);
2638 }
2639
2640 /* Loop on couplings */
2641
2642 for (int cpl_id = 0; cpl_id < n_cpl; cpl_id++) {
2643
2644 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
2645 cs_syr_coupling_ent_t *coupling_ent = syr_coupling->faces;
2646
2647 if (coupling_ent == NULL) /* ignore if volume-only */
2648 continue;
2649
2650 cs_lnum_t n_cpl_faces = coupling_ent->n_elts;
2651
2652 /* Get list of coupled faces */
2653
2654 fvm_nodal_get_parent_id(coupling_ent->elts,
2655 coupling_ent->elt_dim,
2656 f_ids);
2657
2658 switch (cs_glob_thermal_model->itherm) {
2659
2660 case CS_THERMAL_MODEL_TEMPERATURE:
2661 {
2662 for (cs_lnum_t i = 0; i < n_cpl_faces; i++) {
2663 cs_lnum_t face_id = f_ids[i];
2664
2665 /* Saved fluid temperatures and exchange coefficients */
2666 t_fluid[i] = v_fluid[face_id];
2667 h_cpl[i] = h_wall[face_id];
2668 }
2669 }
2670 break;
2671
2672 case CS_THERMAL_MODEL_ENTHALPY:
2673 {
2674 /* In enthalpy formulation, transform to temperatures for SYRTHES
2675 * To conserve flux Phi = (lambda/d ) Delta T
2676 * or Phi = (lambda/(d Cp)) Delta H
2677 * recall hbord = lambda/d.
2678 * Conservation is not guaranteed, so we add a warning. */
2679
2680 for (cs_lnum_t i = 0; i < n_cpl_faces; i++) {
2681 cs_lnum_t face_id = f_ids[i];
2682
2683 t_fluid[i] = wa[face_id];
2684 h_cpl[i] = h_wall[face_id];
2685 }
2686 }
2687 break;
2688
2689 case CS_THERMAL_MODEL_TOTAL_ENERGY:
2690 {
2691 /* In energy formulation, transform to temperatures for SYRTHES
2692 * To conserve flux Phi = (lambda/d ) Delta T
2693 * or Phi = (lambda/(d Cp)) Delta H
2694 * Recall hbord = lambda/ d
2695 * Note that Ei = Cv Ti + 1/2 Ui*Ui + Epsilon_sup_i
2696 * and that Ep = Cv Tp + 1/2 Ui*Ui + Epsilon_sup_i
2697 * (the difference is thus Cv Delta T)
2698
2699 * Modify temperature and exchange coefficient
2700
2701 * Compute e - CvT */
2702
2703 const cs_lnum_t *b_face_cells = cs_glob_mesh->b_face_cells;
2704
2705 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
2706 const cs_real_t *cv = NULL;
2707 cs_lnum_t cv_step = 0;
2708
2709 const cs_real_3_t *cvar_vel = (const cs_real_3_t *)CS_F_(vel)->val;
2710
2711 if (CS_F_(cv) != NULL) {
2712 cv = (const cs_real_t *)CS_F_(cv)->val;
2713 cv_step = 1;
2714 }
2715 else
2716 cv = &cv0;
2717
2718 const cs_real_t *cvar_e_tot = (const cs_real_t *)CS_F_(e_tot)->val;
2719
2720 for (cs_lnum_t i = 0; i < n_cpl_faces; i++) {
2721 cs_lnum_t face_id = f_ids[i];
2722 cs_lnum_t cell_id = b_face_cells[face_id];
2723
2724 cs_real_t cvt = cvar_e_tot[face_id]
2725 - 0.5*cs_math_3_square_norm(cvar_vel[cell_id])
2726 + wa[cell_id];
2727
2728 t_fluid[i] = cvt / cv[cell_id * cv_step];
2729 h_cpl[i] = h_wall[face_id];
2730 }
2731 }
2732 break;
2733
2734 default:
2735 break;
2736 }
2737
2738 /* Fluxes are multiplied by porosity if present.
2739 * Here as the flux is expressed as h.(Tw-Tf), the exchange coefficient
2740 * is multipled by the porosity. */
2741
2742 const cs_field_t *f_poro = CS_F_(poro);
2743
2744 if (f_poro != NULL) {
2745
2746 const cs_lnum_t *b_face_cells = cs_glob_mesh->b_face_cells;
2747
2748 if (f_poro->dim == 1) {
2749 const cs_real_t * cpro_poro = (const cs_real_t *)f_poro->val;
2750 for (cs_lnum_t i = 0; i < n_cpl_faces; i++) {
2751 cs_lnum_t face_id = f_ids[i];
2752 cs_lnum_t cell_id = b_face_cells[face_id];
2753 h_cpl[i] *= cpro_poro[cell_id];
2754 }
2755 }
2756 else if (f_poro->dim == 6) {
2757 const cs_real_6_t * cpro_poro = (const cs_real_6_t *)f_poro->val;
2758 for (cs_lnum_t i = 0; i < n_cpl_faces; i++) {
2759 cs_lnum_t face_id = f_ids[i];
2760 cs_lnum_t cell_id = b_face_cells[face_id];
2761 /* TODO: using a product between the porosity
2762 and the boundary face normal would be more precise. */
2763 h_cpl[i] *= 1./3. * ( cpro_poro[cell_id][0]
2764 + cpro_poro[cell_id][1]
2765 + cpro_poro[cell_id][2]);
2766 }
2767 }
2768
2769 }
2770
2771 _syr_coupling_send_tf_hf(syr_coupling, f_ids, t_fluid, h_cpl, 0);
2772
2773 } /* End loop on couplings */
2774
2775 BFT_FREE(wa);
2776
2777 BFT_FREE(f_ids);
2778 BFT_FREE(t_fluid);
2779 BFT_FREE(h_cpl);
2780 }
2781
2782 /*----------------------------------------------------------------------------*/
2783 /*!
2784 * \brief Exchange volume values relative to a SYRTHES coupling.
2785 */
2786 /*----------------------------------------------------------------------------*/
2787
2788 void
cs_syr_coupling_exchange_volume(void)2789 cs_syr_coupling_exchange_volume(void)
2790 {
2791 const int kcpsyr = cs_field_key_id("syrthes_coupling");
2792
2793 /* Get number of coupling cases */
2794
2795 int n_cpl = cs_syr_coupling_n_couplings();
2796
2797 /* Loop on fields, handling only those coupled with Syrthes */
2798
2799 int n_fields = cs_field_n_fields();
2800 for (int field_id = 0 ; field_id < n_fields; field_id++) {
2801
2802 cs_field_t *f = cs_field_by_id(field_id);
2803
2804 int icpsyr = 0;
2805 if (f->type & CS_FIELD_VARIABLE)
2806 icpsyr = cs_field_get_key_int(f, kcpsyr);
2807
2808 if (icpsyr < 1)
2809 continue;
2810
2811 /* Sanity check : only temperature is possible when doing a
2812 volume coupling with SYRTHES */
2813 if (f != cs_thermal_model_field())
2814 bft_error
2815 (__FILE__, __LINE__, 0,
2816 _("SYRTHES volume coupling possible only with temperature variable,\n"
2817 "not \"%s\"."),
2818 f->name);
2819
2820 /* Loop on couplings: get wall temperature array for each coupling
2821 and apply matching boundary condition. */
2822
2823 for (int cpl_id = 0; cpl_id < n_cpl; cpl_id++) {
2824
2825 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
2826 cs_syr_coupling_ent_t *coupling_ent = syr_coupling->cells;
2827
2828 if (coupling_ent == NULL) /* ignore if surface-only */
2829 continue;
2830
2831 cs_lnum_t n_cpl_cells = coupling_ent->n_elts;
2832
2833 /* Get list of coupled cells */
2834
2835 cs_lnum_t *c_ids;
2836 cs_real_t *t_fluid, *h_vol;
2837 BFT_MALLOC(c_ids, n_cpl_cells, cs_lnum_t);
2838 BFT_MALLOC(t_fluid, n_cpl_cells, cs_real_t);
2839 BFT_MALLOC(h_vol, n_cpl_cells, cs_real_t);
2840
2841 fvm_nodal_get_parent_id(coupling_ent->elts,
2842 coupling_ent->elt_dim,
2843 c_ids);
2844
2845 for (cs_lnum_t i = 0; i < n_cpl_cells; i++) {
2846 h_vol[i] = 0.;
2847 }
2848
2849 /* Receive solid temperature.
2850 * This temperature is stored in a C structure for a future
2851 * use in source term definition. */
2852
2853 _syr_coupling_recv_tsolid(syr_coupling, t_fluid, 1);
2854
2855 const cs_real_t *cvar_t = (const cs_real_t *)f->val;
2856
2857 const char *syrthes_name = _syr_coupling_get_name(syr_coupling);
2858
2859
2860 cs_user_syrthes_coupling_volume_h(cpl_id,
2861 syrthes_name,
2862 n_cpl_cells,
2863 c_ids,
2864 h_vol);
2865
2866 for (cs_lnum_t i = 0; i < n_cpl_cells; i++)
2867 t_fluid[i] = cvar_t[c_ids[i]];
2868
2869 _syr_coupling_send_tf_hf(syr_coupling, c_ids, t_fluid, h_vol, 1);
2870
2871 BFT_FREE(c_ids);
2872 BFT_FREE(t_fluid);
2873 BFT_FREE(h_vol);
2874
2875 } /* End loop on couplings */
2876
2877 } /* End loop on fields */
2878 }
2879
2880 /*----------------------------------------------------------------------------*/
2881 /*!
2882 * \brief Compute the source term (implicit and/or explicit part) for a
2883 * volume coupling with SYRTHES.
2884 *
2885 * \param[in] field_id field id
2886 * \param[in, out] st_exp explicit source term
2887 * \param[in, out] st_imp implicit source term
2888 */
2889 /*----------------------------------------------------------------------------*/
2890
2891 void
cs_syr_coupling_volume_source_terms(int field_id,cs_real_t st_exp[],cs_real_t st_imp[])2892 cs_syr_coupling_volume_source_terms(int field_id,
2893 cs_real_t st_exp[],
2894 cs_real_t st_imp[])
2895 {
2896 cs_field_t *f = cs_field_by_id(field_id);
2897
2898 const cs_real_t *cell_f_vol = cs_glob_mesh_quantities->cell_f_vol;
2899
2900 /* Get number of coupling cases */
2901
2902 int n_cpl = cs_syr_coupling_n_couplings();
2903
2904 /* Sanity check : only temperature is possible when doing a
2905 volume coupling with SYRTHES */
2906 if (f != cs_thermal_model_field())
2907 bft_error
2908 (__FILE__, __LINE__, 0,
2909 _("SYRTHES volume coupling possible only with temperature variable,\n"
2910 "not \"%s\"."),
2911 f->name);
2912
2913 /* Loop on couplings: get wall temperature array for each coupling
2914 and apply matching boundary condition. */
2915
2916 for (int cpl_id = 0; cpl_id < n_cpl; cpl_id++) {
2917
2918 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
2919 cs_syr_coupling_ent_t *coupling_ent = syr_coupling->cells;
2920
2921 if (coupling_ent == NULL) /* ignore if surface-only */
2922 continue;
2923
2924 /* Get list of coupled cells */
2925
2926 cs_lnum_t n_cpl_cells = coupling_ent->n_elts;
2927
2928 cs_lnum_t *c_ids;
2929 cs_real_t *t_fluid, *ctbimp, *ctbexp;
2930 BFT_MALLOC(c_ids, n_cpl_cells, cs_lnum_t);
2931 BFT_MALLOC(t_fluid, n_cpl_cells, cs_real_t);
2932 BFT_MALLOC(ctbimp, n_cpl_cells, cs_real_t);
2933 BFT_MALLOC(ctbexp, n_cpl_cells, cs_real_t);
2934
2935 fvm_nodal_get_parent_id(coupling_ent->elts,
2936 coupling_ent->elt_dim,
2937 c_ids);
2938
2939 /* Compute implicit and explicit contribution to source terms */
2940
2941 const cs_real_t *cvara_vart = (const cs_real_t *)f->vals[1];
2942
2943 for (cs_lnum_t i = 0; i < n_cpl_cells; i++) {
2944 t_fluid[i] = cvara_vart[c_ids[i]];
2945 }
2946
2947 _syr_coupling_ts_contrib(syr_coupling, t_fluid, ctbimp, ctbexp);
2948
2949 /* Loop on coupled cells to compute crvexp and crvimp */
2950
2951 for (cs_lnum_t i = 0; i < n_cpl_cells; i++) {
2952
2953 cs_lnum_t c_id = c_ids[i];
2954
2955 cs_real_t tsexp = (ctbexp[i] - ctbimp[i]*t_fluid[i]) * cell_f_vol[c_id];
2956 cs_real_t tsimp = ctbimp[i] * cell_f_vol[c_id];
2957
2958 st_exp[c_id] += tsexp;
2959 st_imp[c_id] += tsimp;
2960
2961 }
2962
2963 BFT_FREE(c_ids);
2964 BFT_FREE(t_fluid);
2965 BFT_FREE(ctbimp);
2966 BFT_FREE(ctbexp);
2967
2968 } /* End loop on couplings */
2969 }
2970
2971 /*----------------------------------------------------------------------------*/
2972 /*!
2973 * \brief Get number of coupled elements with SYRTHES.
2974 *
2975 * \param[in] cpl_id coupling id
2976 * \param[in] mode 0 for boundary, 1 for volume
2977 *
2978 * \return number of coupled elements for this coupling
2979 */
2980 /*----------------------------------------------------------------------------*/
2981
2982 cs_lnum_t
cs_syr_coupling_n_elts(int cpl_id,int mode)2983 cs_syr_coupling_n_elts(int cpl_id,
2984 int mode)
2985 {
2986 cs_lnum_t retval = 0;
2987
2988 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
2989
2990 if (syr_coupling == NULL)
2991 bft_error(__FILE__, __LINE__, 0,
2992 _("SYRTHES coupling id %d impossible; "
2993 "there are %d couplings"),
2994 cpl_id, _syr_n_couplings);
2995
2996 else {
2997
2998 cs_syr_coupling_ent_t *coupling_ent = NULL;
2999
3000 assert(mode == 0 || mode == 1);
3001
3002 if (mode == 0)
3003 coupling_ent = syr_coupling->faces;
3004 else
3005 coupling_ent = syr_coupling->cells;
3006
3007 if (coupling_ent != NULL)
3008 retval = coupling_ent->n_elts;
3009
3010 }
3011
3012 return retval;
3013 }
3014
3015 /*----------------------------------------------------------------------------*/
3016 /*!
3017 * \brief Get local ids of elements coupled with SYRTHES
3018 *
3019 * \param[in] cpl_id coupling id
3020 * \param[in] mode 0 for boundary, 1 for volume
3021 * \param[out] elt_ids ids of coupled elements (preallocated)
3022 */
3023 /*----------------------------------------------------------------------------*/
3024
3025 void
cs_syr_coupling_elt_ids(int cpl_id,int mode,cs_lnum_t elt_ids[])3026 cs_syr_coupling_elt_ids(int cpl_id,
3027 int mode,
3028 cs_lnum_t elt_ids[])
3029 {
3030 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
3031
3032 if (syr_coupling == NULL)
3033 bft_error(__FILE__, __LINE__, 0,
3034 _("SYRTHES coupling id %d impossible; "
3035 "there are %d couplings"),
3036 cpl_id, _syr_n_couplings);
3037
3038 else {
3039
3040 cs_syr_coupling_ent_t *coupling_ent = NULL;
3041
3042 assert(mode == 0 || mode == 1);
3043
3044 if (mode == 0)
3045 coupling_ent = syr_coupling->faces;
3046 else
3047 coupling_ent = syr_coupling->cells;
3048
3049 if (coupling_ent != NULL)
3050 fvm_nodal_get_parent_id(coupling_ent->elts,
3051 coupling_ent->elt_dim,
3052 elt_ids);
3053
3054 }
3055 }
3056
3057 /*----------------------------------------------------------------------------*/
3058 /*!
3059 * \brief Receive coupling variables from SYRTHES.
3060 *
3061 * \param[in] cpl_id coupling id
3062 * \param[in] mode 0 for boundary, 1 for volume
3063 * \param[out] t_solid solid temperature
3064 */
3065 /*----------------------------------------------------------------------------*/
3066
3067 void
cs_syr_coupling_recv_tsolid(int cpl_id,int mode,cs_real_t t_solid[])3068 cs_syr_coupling_recv_tsolid(int cpl_id,
3069 int mode,
3070 cs_real_t t_solid[])
3071 {
3072 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
3073
3074 if (syr_coupling == NULL)
3075 bft_error(__FILE__, __LINE__, 0,
3076 _("SYRTHES coupling id %d impossible; "
3077 "there are %d couplings"),
3078 cpl_id, _syr_n_couplings);
3079
3080 else
3081 _syr_coupling_recv_tsolid(syr_coupling, t_solid, mode);
3082 }
3083
3084 /*----------------------------------------------------------------------------*/
3085 /*!
3086 * \brief Send coupling variables to SYRTHES.
3087 *
3088 * \param[in] cpl_id coupling id
3089 * \param[in] mode 0 for boundary, 1 for volume
3090 * \param[in] elt_ids ids of coupled elements
3091 * \param[in] t_fluid fluid temperature
3092 * \param[in] h_fluid fluid exchange coefficient
3093 */
3094 /*----------------------------------------------------------------------------*/
3095
3096 void
cs_syr_coupling_send_tf_hf(int cpl_id,int mode,const cs_lnum_t elt_ids[],cs_real_t t_fluid[],cs_real_t h_fluid[])3097 cs_syr_coupling_send_tf_hf(int cpl_id,
3098 int mode,
3099 const cs_lnum_t elt_ids[],
3100 cs_real_t t_fluid[],
3101 cs_real_t h_fluid[])
3102 {
3103 cs_syr_coupling_t *syr_coupling = _syr_coupling_by_id(cpl_id);
3104
3105 if (syr_coupling == NULL)
3106 bft_error(__FILE__, __LINE__, 0,
3107 _("SYRTHES coupling id %d impossible; "
3108 "there are %d couplings"),
3109 cpl_id, _syr_n_couplings);
3110
3111 else
3112 _syr_coupling_send_tf_hf(syr_coupling, elt_ids,
3113 t_fluid, h_fluid, mode);
3114 }
3115
3116 /*----------------------------------------------------------------------------*/
3117
3118 END_C_DECLS
3119