1 /*============================================================================
2 * Functions to handle common features for building algebraic system in CDO
3 * schemes
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <assert.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <float.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include <bft_mem.h>
44
45 #include "cs_blas.h"
46 #include "cs_boundary_zone.h"
47 #include "cs_cdo_local.h"
48 #include "cs_log.h"
49 #include "cs_math.h"
50 #include "cs_parall.h"
51 #include "cs_parameters.h"
52 #include "cs_xdef_eval.h"
53
54 #if defined(DEBUG) && !defined(NDEBUG)
55 #include "cs_dbg.h"
56 #endif
57
58 /*----------------------------------------------------------------------------*/
59
60 BEGIN_C_DECLS
61
62 /*----------------------------------------------------------------------------
63 * Header for the current file
64 *----------------------------------------------------------------------------*/
65
66 #include "cs_equation_common.h"
67
68 /*----------------------------------------------------------------------------*/
69
70 BEGIN_C_DECLS
71
72 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
73
74 /*============================================================================
75 * Type definitions and macros
76 *============================================================================*/
77
78 #define CS_EQUATION_COMMON_DBG 0 /* Debug level */
79
80 /*============================================================================
81 * Local private variables
82 *============================================================================*/
83
84 /* Temporary buffers useful during the building of all algebraic systems */
85
86 static size_t cs_equation_common_work_buffer_size = 0;
87 static cs_real_t *cs_equation_common_work_buffer = NULL;
88
89 /* Pointer to shared structures (owned by a cs_domain_t structure) */
90
91 static const cs_cdo_quantities_t *cs_shared_quant;
92 static const cs_cdo_connect_t *cs_shared_connect;
93 static const cs_time_step_t *cs_shared_time_step;
94
95 /*============================================================================
96 * Private function prototypes
97 *============================================================================*/
98
99 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
100
101 /*============================================================================
102 * Public function prototypes
103 *============================================================================*/
104
105 /*----------------------------------------------------------------------------*/
106 /*!
107 * \brief Allocate a pointer to a buffer of size at least the n_cells for
108 * managing temporary usage of memory when dealing with equations
109 * The size of the temporary buffer can be bigger according to the
110 * numerical settings
111 * Set also shared pointers from the main domain members
112 *
113 * \param[in] connect pointer to a cs_cdo_connect_t structure
114 * \param[in] quant pointer to additional mesh quantities struct.
115 * \param[in] time_step pointer to a time step structure
116 * \param[in] eb_flag metadata for Edge-based schemes
117 * \param[in] fb_flag metadata for Face-based schemes
118 * \param[in] vb_flag metadata for Vertex-based schemes
119 * \param[in] vcb_flag metadata for Vertex+Cell-basde schemes
120 * \param[in] hho_flag metadata for HHO schemes
121 */
122 /*----------------------------------------------------------------------------*/
123
124 void
cs_equation_common_init(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step,cs_flag_t eb_flag,cs_flag_t fb_flag,cs_flag_t vb_flag,cs_flag_t vcb_flag,cs_flag_t hho_flag)125 cs_equation_common_init(const cs_cdo_connect_t *connect,
126 const cs_cdo_quantities_t *quant,
127 const cs_time_step_t *time_step,
128 cs_flag_t eb_flag,
129 cs_flag_t fb_flag,
130 cs_flag_t vb_flag,
131 cs_flag_t vcb_flag,
132 cs_flag_t hho_flag)
133 {
134 assert(connect != NULL && quant != NULL); /* Sanity check */
135
136 /* Allocate cell-wise and face-wise view of a mesh */
137
138 cs_cdo_local_initialize(connect);
139
140 const cs_lnum_t n_cells = connect->n_cells;
141 const cs_lnum_t n_faces = connect->n_faces[CS_ALL_FACES];
142 const cs_lnum_t n_vertices = connect->n_vertices;
143 const cs_lnum_t n_edges = connect->n_edges;
144
145 /* Allocate shared buffer and initialize shared structures */
146
147 size_t cwb_size = n_cells; /* initial cell-wise buffer size */
148
149 /* Allocate and initialize matrix assembler and matrix structures */
150
151 if (vb_flag > 0 || vcb_flag > 0) {
152
153 if (vb_flag & CS_FLAG_SCHEME_SCALAR || vcb_flag & CS_FLAG_SCHEME_SCALAR) {
154
155 if (vb_flag & CS_FLAG_SCHEME_SCALAR)
156 cwb_size = CS_MAX(cwb_size, (size_t)n_vertices);
157
158 if (vcb_flag & CS_FLAG_SCHEME_SCALAR)
159 cwb_size = CS_MAX(cwb_size, (size_t)(n_vertices + n_cells));
160
161 } /* scalar-valued equations */
162
163 if (vb_flag & CS_FLAG_SCHEME_VECTOR || vcb_flag & CS_FLAG_SCHEME_VECTOR) {
164
165 cwb_size = CS_MAX(cwb_size, (size_t)3*n_cells);
166 if (vb_flag & CS_FLAG_SCHEME_VECTOR)
167 cwb_size = CS_MAX(cwb_size, (size_t)3*n_vertices);
168
169 if (vcb_flag & CS_FLAG_SCHEME_VECTOR)
170 cwb_size = CS_MAX(cwb_size, (size_t)3*(n_vertices + n_cells));
171
172 } /* vector-valued equations */
173
174 } /* Vertex-based schemes and related ones */
175
176 if (eb_flag > 0) {
177
178 if (eb_flag & CS_FLAG_SCHEME_SCALAR) {
179
180 /* This is a vector-valued equation but the DoF is scalar-valued since
181 * it is a circulation associated to each edge */
182
183 cwb_size = CS_MAX(cwb_size, (size_t)3*n_cells);
184 cwb_size = CS_MAX(cwb_size, (size_t)n_edges);
185
186 } /* vector-valued equations with scalar-valued DoFs */
187
188 } /* Edge-based schemes */
189
190 if (fb_flag > 0 || hho_flag > 0) {
191
192 if (cs_flag_test(fb_flag, CS_FLAG_SCHEME_POLY0 | CS_FLAG_SCHEME_SCALAR) ||
193 cs_flag_test(hho_flag, CS_FLAG_SCHEME_POLY0 | CS_FLAG_SCHEME_SCALAR)) {
194
195 assert(n_faces > n_cells);
196 if (fb_flag & CS_FLAG_SCHEME_SCALAR)
197 cwb_size = CS_MAX(cwb_size, (size_t)n_faces);
198
199 if (hho_flag & CS_FLAG_SCHEME_SCALAR)
200 cwb_size = CS_MAX(cwb_size, (size_t)n_faces);
201
202 } /* Scalar-valued CDO-Fb or HHO-P0 */
203
204 if (cs_flag_test(fb_flag, CS_FLAG_SCHEME_POLY0 | CS_FLAG_SCHEME_VECTOR) ||
205 cs_flag_test(hho_flag, CS_FLAG_SCHEME_POLY1 | CS_FLAG_SCHEME_SCALAR) ||
206 cs_flag_test(hho_flag, CS_FLAG_SCHEME_POLY0 | CS_FLAG_SCHEME_VECTOR)) {
207
208 assert((CS_CDO_CONNECT_FACE_SP1 == CS_CDO_CONNECT_FACE_VP0) &&
209 (CS_CDO_CONNECT_FACE_SP1 == CS_CDO_CONNECT_FACE_VHP0));
210
211 cwb_size = CS_MAX(cwb_size, (size_t)CS_N_FACE_DOFS_1ST * n_faces);
212
213 } /* Vector CDO-Fb or HHO-P1 or vector HHO-P0 */
214
215 if (cs_flag_test(hho_flag,
216 CS_FLAG_SCHEME_POLY2 | CS_FLAG_SCHEME_SCALAR))
217 cwb_size = CS_MAX(cwb_size, (size_t)CS_N_FACE_DOFS_2ND * n_faces);
218
219 /* For vector equations and HHO */
220
221 if (cs_flag_test(hho_flag, CS_FLAG_SCHEME_VECTOR | CS_FLAG_SCHEME_POLY1) ||
222 cs_flag_test(hho_flag, CS_FLAG_SCHEME_VECTOR | CS_FLAG_SCHEME_POLY2)) {
223
224 if (hho_flag & CS_FLAG_SCHEME_POLY1)
225 cwb_size = CS_MAX(cwb_size, (size_t)3*CS_N_FACE_DOFS_1ST*n_faces);
226
227 else if (hho_flag & CS_FLAG_SCHEME_POLY2)
228 cwb_size = CS_MAX(cwb_size, (size_t)3*CS_N_FACE_DOFS_2ND*n_faces);
229
230 }
231
232 } /* Face-based schemes (CDO or HHO) */
233
234 /* Assign static const pointers: shared pointers with a cs_domain_t */
235
236 cs_shared_quant = quant;
237 cs_shared_connect = connect;
238 cs_shared_time_step = time_step;
239
240 /* Common buffer for temporary usage */
241
242 cs_equation_common_work_buffer_size = cwb_size;
243 BFT_MALLOC(cs_equation_common_work_buffer, cwb_size, double);
244 }
245
246 /*----------------------------------------------------------------------------*/
247 /*!
248 * \brief Free buffers shared among the equations solved with CDO schemes
249 */
250 /*----------------------------------------------------------------------------*/
251
252 void
cs_equation_common_finalize(void)253 cs_equation_common_finalize(void)
254 {
255 /* Free cell-wise and face-wise view of a mesh */
256
257 cs_cdo_local_finalize();
258
259 /* Free common buffer */
260
261 BFT_FREE(cs_equation_common_work_buffer);
262 }
263
264 /*----------------------------------------------------------------------------*/
265 /*!
266 * \brief Allocate a new structure to handle the building of algebraic system
267 * related to a cs_equation_t structure
268 *
269 * \param[in] eqp pointer to a cs_equation_param_t structure
270 * \param[in] mesh pointer to a cs_mesh_t structure
271 *
272 * \return a pointer to a new allocated cs_equation_builder_t structure
273 */
274 /*----------------------------------------------------------------------------*/
275
276 cs_equation_builder_t *
cs_equation_builder_init(const cs_equation_param_t * eqp,const cs_mesh_t * mesh)277 cs_equation_builder_init(const cs_equation_param_t *eqp,
278 const cs_mesh_t *mesh)
279 {
280 cs_equation_builder_t *eqb = NULL;
281
282 BFT_MALLOC(eqb, 1, cs_equation_builder_t);
283
284 eqb->init_step = true;
285
286 /* Initialize flags used to knows what kind of cell quantities to build */
287
288 eqb->msh_flag = 0;
289 eqb->bd_msh_flag = 0;
290 eqb->st_msh_flag = 0;
291 if (eqp->dim > 1)
292 eqb->sys_flag = CS_FLAG_SYS_VECTOR;
293 else
294 eqb->sys_flag = 0;
295
296 /* Handle properties */
297
298 eqb->diff_pty_uniform = true;
299 if (cs_equation_param_has_diffusion(eqp))
300 eqb->diff_pty_uniform = cs_property_is_uniform(eqp->diffusion_property);
301
302 eqb->curlcurl_pty_uniform = true;
303 if (cs_equation_param_has_curlcurl(eqp))
304 eqb->curlcurl_pty_uniform = cs_property_is_uniform(eqp->curlcurl_property);
305
306 eqb->graddiv_pty_uniform = true;
307 if (cs_equation_param_has_graddiv(eqp))
308 eqb->graddiv_pty_uniform = cs_property_is_uniform(eqp->graddiv_property);
309
310 eqb->time_pty_uniform = true;
311 if (cs_equation_param_has_time(eqp))
312 eqb->time_pty_uniform = cs_property_is_uniform(eqp->time_property);
313
314 if (eqp->n_reaction_terms > CS_CDO_N_MAX_REACTIONS)
315 bft_error(__FILE__, __LINE__, 0,
316 " %s: Number of reaction terms for an equation is too high.\n"
317 " Current value: %d (max: %d)\n"
318 " Change the value of CS_CDO_N_MAX_REACTIONS in the code or\n"
319 " modify your settings or contact the developpement team.",
320 __func__, eqp->n_reaction_terms, CS_CDO_N_MAX_REACTIONS);
321
322 for (int i = 0; i < eqp->n_reaction_terms; i++)
323 eqb->reac_pty_uniform[i]
324 = cs_property_is_uniform(eqp->reaction_properties[i]);
325
326 /* Enforcement of DoFs */
327
328 eqb->enforced_values = NULL;
329
330 /* Handle source terms */
331
332 eqb->source_mask = NULL;
333 if (cs_equation_param_has_sourceterm(eqp)) {
334
335 /* Default initialization */
336
337 eqb->st_msh_flag = cs_source_term_init(eqp->space_scheme,
338 eqp->n_source_terms,
339 (cs_xdef_t *const *)eqp->source_terms,
340 eqb->compute_source,
341 &(eqb->sys_flag),
342 &(eqb->source_mask));
343
344 } /* There is at least one source term */
345
346 /* Set members and structures related to the management of the BCs
347 Translate user-defined information about BC into a structure well-suited
348 for computation. We make the distinction between homogeneous and
349 non-homogeneous BCs. */
350
351 eqb->dir_values = NULL;
352 eqb->face_bc = cs_cdo_bc_face_define(eqp->default_bc,
353 true, /* Steady BC up to now */
354 eqp->dim,
355 eqp->n_bc_defs,
356 eqp->bc_defs,
357 mesh->n_b_faces);
358
359 /* Monitoring */
360
361 CS_TIMER_COUNTER_INIT(eqb->tcb); /* build system */
362 CS_TIMER_COUNTER_INIT(eqb->tcs); /* solve system */
363 CS_TIMER_COUNTER_INIT(eqb->tce); /* extra operations */
364
365 return eqb;
366 }
367
368 /*----------------------------------------------------------------------------*/
369 /*!
370 * \brief Free a cs_equation_builder_t structure
371 *
372 * \param[in, out] p_builder pointer of pointer to the cs_equation_builder_t
373 * structure to free
374 */
375 /*----------------------------------------------------------------------------*/
376
377 void
cs_equation_builder_free(cs_equation_builder_t ** p_builder)378 cs_equation_builder_free(cs_equation_builder_t **p_builder)
379 {
380 if (p_builder == NULL)
381 return;
382 if (*p_builder == NULL)
383 return;
384
385 cs_equation_builder_t *eqb = *p_builder;
386
387 cs_equation_builder_reset(eqb);
388
389 if (eqb->source_mask != NULL)
390 BFT_FREE(eqb->source_mask);
391
392 /* Free BC structure */
393
394 eqb->face_bc = cs_cdo_bc_free(eqb->face_bc);
395
396 BFT_FREE(eqb);
397
398 *p_builder = NULL;
399 }
400
401 /*----------------------------------------------------------------------------*/
402 /*!
403 * \brief Free some members of a cs_equation_builder_t structure
404 *
405 * \param[in, out] eqb pointer to the cs_equation_builder_t structure
406 */
407 /*----------------------------------------------------------------------------*/
408
409 void
cs_equation_builder_reset(cs_equation_builder_t * eqb)410 cs_equation_builder_reset(cs_equation_builder_t *eqb)
411 {
412 if (eqb == NULL)
413 return;
414
415 BFT_FREE(eqb->enforced_values);
416 BFT_FREE(eqb->dir_values);
417 }
418
419 /*----------------------------------------------------------------------------*/
420 /*!
421 * \brief Compute the value of the renormalization coefficient for the
422 * the residual norm of the linear system
423 * A pre-computation arising from cellwise contribution may have
424 * been done before this call according to the requested type of
425 * renormalization.
426 *
427 * \param[in] type type of renormalization
428 * \param[in] rhs_size size of the rhs array
429 * \param[in] rhs array related to the right-hand side
430 * \param[in, out] normalization value of the residual normalization
431 */
432 /*----------------------------------------------------------------------------*/
433
434 void
cs_equation_sync_rhs_normalization(cs_param_resnorm_type_t type,cs_lnum_t rhs_size,const cs_real_t rhs[],double * normalization)435 cs_equation_sync_rhs_normalization(cs_param_resnorm_type_t type,
436 cs_lnum_t rhs_size,
437 const cs_real_t rhs[],
438 double *normalization)
439 {
440 switch (type) {
441
442 case CS_PARAM_RESNORM_NORM2_RHS:
443 *normalization = cs_dot_xx(rhs_size, rhs);
444 cs_parall_sum(1, CS_REAL_TYPE, normalization);
445 if (*normalization < 100*DBL_MIN)
446 *normalization = 1.0;
447 else
448 *normalization = sqrt((*normalization));
449 break;
450
451 case CS_PARAM_RESNORM_FILTERED_RHS:
452 cs_parall_sum(1, CS_REAL_TYPE, normalization);
453 if (*normalization < 100*DBL_MIN)
454 *normalization = 1.0;
455 else
456 *normalization = sqrt((*normalization));
457 break;
458
459 case CS_PARAM_RESNORM_WEIGHTED_RHS:
460 cs_parall_sum(1, CS_REAL_TYPE, normalization);
461 if (*normalization < 100*DBL_MIN)
462 *normalization = 1.0;
463 else
464 *normalization = sqrt((*normalization)/cs_shared_quant->vol_tot);
465 break;
466
467 default:
468 *normalization = 1.0;
469 break;
470
471 } /* Type of normalization */
472 }
473
474 /*----------------------------------------------------------------------------*/
475 /*!
476 * \brief Prepare a linear system and synchronize buffers to handle parallelism.
477 * Transfer a mesh-based description of arrays x0 and rhs into an
478 * algebraic description for the linear system in x and b.
479 *
480 * \param[in] stride stride to apply to the range set operations
481 * \param[in] x_size size of the vector unknowns (scatter view)
482 * \param[in] matrix pointer to a cs_matrix_t structure
483 * \param[in] rset pointer to a range set structure
484 * \param[in] rhs_redux do or not a parallel sum reduction on the RHS
485 * \param[in, out] x array of unknowns (in: initial guess)
486 * \param[in, out] b right-hand side
487 */
488 /*----------------------------------------------------------------------------*/
489
490 void
cs_equation_prepare_system(int stride,cs_lnum_t x_size,const cs_matrix_t * matrix,const cs_range_set_t * rset,bool rhs_redux,cs_real_t * x,cs_real_t * b)491 cs_equation_prepare_system(int stride,
492 cs_lnum_t x_size,
493 const cs_matrix_t *matrix,
494 const cs_range_set_t *rset,
495 bool rhs_redux,
496 cs_real_t *x,
497 cs_real_t *b)
498 {
499 const cs_lnum_t n_scatter_elts = x_size; /* size of x and rhs */
500
501 #if defined(DEBUG) && !defined(NDEBUG) && CS_EQUATION_COMMON_DBG > 0
502 const cs_lnum_t n_gather_elts = cs_matrix_get_n_rows(matrix);
503 assert(n_gather_elts <= n_scatter_elts);
504
505 cs_log_printf(CS_LOG_DEFAULT,
506 " n_gather_elts: %d\n"
507 " n_scatter_elts: %d\n"
508 " n_matrix_rows: %d\n"
509 " n_matrix_columns: %d\n",
510 n_gather_elts, n_scatter_elts, cs_matrix_get_n_rows(matrix),
511 cs_matrix_get_n_columns(matrix));
512 #else
513 CS_UNUSED(matrix);
514 #endif
515
516 if (rset != NULL) { /* Parallel or periodic mode
517 ========================= */
518
519 /* x and b should be changed to have a "gathered" view through the range set
520 operation. Their size is equal to n_sles_gather_elts <=
521 n_sles_scatter_elts */
522
523 /* Compact numbering to fit the algebraic decomposition */
524
525 cs_range_set_gather(rset,
526 CS_REAL_TYPE, /* type */
527 stride, /* stride */
528 x, /* in: size = n_sles_scatter_elts */
529 x); /* out: size = n_sles_gather_elts */
530
531 /* The right-hand side stems from a cellwise building on this rank.
532 Other contributions from distant ranks may contribute to an element
533 owned by the local rank */
534
535 if (rhs_redux && rset->ifs != NULL)
536 cs_interface_set_sum(rset->ifs,
537 n_scatter_elts, stride, false, CS_REAL_TYPE,
538 b);
539
540 cs_range_set_gather(rset,
541 CS_REAL_TYPE,/* type */
542 stride, /* stride */
543 b, /* in: size = n_sles_scatter_elts */
544 b); /* out: size = n_sles_gather_elts */
545 }
546
547 #if defined(DEBUG) && !defined(NDEBUG) && CS_EQUATION_COMMON_DBG > 2
548 const cs_lnum_t *row_index, *col_id;
549 const cs_real_t *d_val, *x_val;
550
551 cs_matrix_get_msr_arrays(matrix, &row_index, &col_id, &d_val, &x_val);
552
553 cs_dbg_dump_linear_system("Dump linear system",
554 n_gather_elts, CS_EQUATION_COMMON_DBG,
555 x, b,
556 row_index, col_id, x_val, d_val);
557 #endif
558 }
559
560 /*----------------------------------------------------------------------------*/
561 /*!
562 * \brief Solve a linear system arising with scalar-valued cell-based DoFs*
563 * No rotation is taken into account when synchronizing halo.
564 *
565 * \param[in] n_dofs local number of DoFs
566 * \param[in] slesp pointer to a cs_param_sles_t structure
567 * \param[in] matrix pointer to a cs_matrix_t structure
568 * \param[in] normalization value used for the residual normalization
569 * \param[in, out] sles pointer to a cs_sles_t structure
570 * \param[in, out] x solution of the linear system (in: initial guess)
571 * \param[in, out] b right-hand side (scatter/gather if needed)
572 *
573 * \return the number of iterations of the linear solver
574 */
575 /*----------------------------------------------------------------------------*/
576
577 int
cs_equation_solve_scalar_cell_system(cs_lnum_t n_dofs,const cs_param_sles_t * slesp,const cs_matrix_t * matrix,cs_real_t normalization,cs_sles_t * sles,cs_real_t * x,cs_real_t * b)578 cs_equation_solve_scalar_cell_system(cs_lnum_t n_dofs,
579 const cs_param_sles_t *slesp,
580 const cs_matrix_t *matrix,
581 cs_real_t normalization,
582 cs_sles_t *sles,
583 cs_real_t *x,
584 cs_real_t *b)
585 {
586 /* Retrieve the solving info structure stored in the cs_field_t structure */
587
588 cs_solving_info_t sinfo;
589 cs_field_t *fld = NULL;
590 if (slesp->field_id > -1) {
591 fld = cs_field_by_id(slesp->field_id);
592 cs_field_get_key_struct(fld, cs_field_key_id("solving_info"), &sinfo);
593 }
594
595 sinfo.n_it = 0;
596 sinfo.res_norm = DBL_MAX;
597 sinfo.rhs_norm = normalization;
598
599 const cs_halo_t *halo = cs_matrix_get_halo(matrix);
600 const cs_lnum_t n_rows = cs_matrix_get_n_rows(matrix);
601 const cs_lnum_t n_cols_ext = cs_matrix_get_n_columns(matrix);
602
603 assert(n_dofs == n_rows);
604
605 /* Handle parallelism */
606
607 cs_real_t *_x = x, *_b = b;
608 if (n_cols_ext > n_rows) {
609
610 BFT_MALLOC(_b, n_cols_ext, cs_real_t);
611 BFT_MALLOC(_x, n_cols_ext, cs_real_t);
612
613 memcpy(_x, x, n_dofs*sizeof(cs_real_t));
614 memcpy(_b, b, n_dofs*sizeof(cs_real_t));
615
616 cs_matrix_pre_vector_multiply_sync(matrix, _b);
617 cs_halo_sync_var(halo, CS_HALO_STANDARD, _x);
618 }
619
620 /* Solve the linear solver */
621
622 cs_sles_convergence_state_t code = cs_sles_solve(sles,
623 matrix,
624 slesp->eps,
625 sinfo.rhs_norm,
626 &(sinfo.n_it),
627 &(sinfo.res_norm),
628 _b,
629 _x,
630 0, /* aux. size */
631 NULL); /* aux. buffers */
632
633 if (n_cols_ext > n_rows) {
634 BFT_FREE(_b);
635 memcpy(x, _x, n_dofs*sizeof(cs_real_t));
636 BFT_FREE(_x);
637 }
638
639 /* Output information about the convergence of the resolution */
640
641 if (slesp->verbosity > 0)
642 cs_log_printf(CS_LOG_DEFAULT, " <%20s/sles_cvg_code=%-d>"
643 " n_iter %3d | res.norm % -8.4e | rhs.norm % -8.4e\n",
644 slesp->name, code,
645 sinfo.n_it, sinfo.res_norm, sinfo.rhs_norm);
646
647 if (slesp->field_id > -1)
648 cs_field_set_key_struct(fld, cs_field_key_id("solving_info"), &sinfo);
649
650 return (sinfo.n_it);
651 }
652
653 /*----------------------------------------------------------------------------*/
654 /*!
655 * \brief Solve a linear system arising from CDO schemes with scalar-valued
656 * degrees of freedom
657 *
658 * \param[in] n_scatter_dofs local number of DoFs (may be != n_gather_elts)
659 * \param[in] slesp pointer to a cs_param_sles_t structure
660 * \param[in] matrix pointer to a cs_matrix_t structure
661 * \param[in] rs pointer to a cs_range_set_t structure
662 * \param[in] normalization value used for the residual normalization
663 * \param[in] rhs_redux do or not a parallel sum reduction on the RHS
664 * \param[in, out] sles pointer to a cs_sles_t structure
665 * \param[in, out] x solution of the linear system (in: initial guess)
666 * \param[in, out] b right-hand side (scatter/gather if needed)
667 *
668 * \return the number of iterations of the linear solver
669 */
670 /*----------------------------------------------------------------------------*/
671
672 int
cs_equation_solve_scalar_system(cs_lnum_t n_scatter_dofs,const cs_param_sles_t * slesp,const cs_matrix_t * matrix,const cs_range_set_t * rset,cs_real_t normalization,bool rhs_redux,cs_sles_t * sles,cs_real_t * x,cs_real_t * b)673 cs_equation_solve_scalar_system(cs_lnum_t n_scatter_dofs,
674 const cs_param_sles_t *slesp,
675 const cs_matrix_t *matrix,
676 const cs_range_set_t *rset,
677 cs_real_t normalization,
678 bool rhs_redux,
679 cs_sles_t *sles,
680 cs_real_t *x,
681 cs_real_t *b)
682 {
683 const cs_lnum_t n_cols = cs_matrix_get_n_columns(matrix);
684
685 /* Set xsol */
686
687 cs_real_t *xsol = NULL;
688 if (n_cols > n_scatter_dofs) {
689 assert(cs_glob_n_ranks > 1);
690 BFT_MALLOC(xsol, n_cols, cs_real_t);
691 memcpy(xsol, x, n_scatter_dofs*sizeof(cs_real_t));
692 }
693 else
694 xsol = x;
695
696 /* Retrieve the solving info structure stored in the cs_field_t structure */
697
698 cs_field_t *fld = cs_field_by_id(slesp->field_id);
699 cs_solving_info_t sinfo;
700 cs_field_get_key_struct(fld, cs_field_key_id("solving_info"), &sinfo);
701
702 sinfo.n_it = 0;
703 sinfo.res_norm = DBL_MAX;
704 sinfo.rhs_norm = normalization;
705
706 /* Prepare solving (handle parallelism)
707 * stride = 1 for scalar-valued */
708
709 cs_equation_prepare_system(1, n_scatter_dofs, matrix, rset, rhs_redux,
710 xsol, b);
711
712 /* Solve the linear solver */
713
714 cs_sles_convergence_state_t code = cs_sles_solve(sles,
715 matrix,
716 slesp->eps,
717 sinfo.rhs_norm,
718 &(sinfo.n_it),
719 &(sinfo.res_norm),
720 b,
721 xsol,
722 0, /* aux. size */
723 NULL); /* aux. buffers */
724
725 /* Output information about the convergence of the resolution */
726
727 if (slesp->verbosity > 0)
728 cs_log_printf(CS_LOG_DEFAULT, " <%20s/sles_cvg_code=%-d>"
729 " n_iter %3d | res.norm % -8.4e | rhs.norm % -8.4e\n",
730 slesp->name, code,
731 sinfo.n_it, sinfo.res_norm, sinfo.rhs_norm);
732
733 cs_range_set_scatter(rset,
734 CS_REAL_TYPE, 1, /* type and stride */
735 xsol, x);
736 cs_range_set_scatter(rset,
737 CS_REAL_TYPE, 1, /* type and stride */
738 b, b);
739
740 #if defined(DEBUG) && !defined(NDEBUG) && CS_EQUATION_COMMON_DBG > 1
741 cs_dbg_fprintf_system(slesp->name, cs_shared_time_step->nt_cur,
742 slesp->verbosity,
743 x, b, n_scatter_dofs);
744 #endif
745
746 if (n_cols > n_scatter_dofs)
747 BFT_FREE(xsol);
748
749 cs_field_set_key_struct(fld, cs_field_key_id("solving_info"), &sinfo);
750
751 return (sinfo.n_it);
752 }
753
754 /*----------------------------------------------------------------------------*/
755 /*!
756 * \brief Print a message in the performance output file related to the
757 * monitoring of equation
758 *
759 * \param[in] eqname pointer to the name of the current equation
760 * \param[in] eqb pointer to a cs_equation_builder_t structure
761 */
762 /*----------------------------------------------------------------------------*/
763
764 void
cs_equation_write_monitoring(const char * eqname,const cs_equation_builder_t * eqb)765 cs_equation_write_monitoring(const char *eqname,
766 const cs_equation_builder_t *eqb)
767 {
768 double t[3] = {eqb->tcb.nsec, eqb->tcs.nsec, eqb->tce.nsec};
769 for (int i = 0; i < 3; i++) t[i] *= 1e-9;
770
771 if (eqname == NULL)
772 cs_log_printf(CS_LOG_PERFORMANCE, " %-35s %10.4f %10.4f %10.4f (seconds)\n",
773 "<CDO/Equation> Monitoring", t[0], t[1], t[2]);
774 else {
775
776 char *msg = NULL;
777 int len = 1 + strlen("<CDO/> Monitoring") + strlen(eqname);
778
779 BFT_MALLOC(msg, len, char);
780 sprintf(msg, "<CDO/%s> Monitoring", eqname);
781 cs_log_printf(CS_LOG_PERFORMANCE, " %-35s %10.4f %10.4f %10.4f (seconds)\n",
782 msg, t[0], t[1], t[2]);
783 BFT_FREE(msg);
784
785 }
786 }
787
788 /*----------------------------------------------------------------------------*/
789 /*!
790 * \brief Initialize all reaction properties. This function is shared across
791 * all CDO schemes. The \ref cs_cell_builder_t structure stores the
792 * computed property values.
793 *
794 * \param[in] eqp pointer to a cs_equation_param_t structure
795 * \param[in] eqb pointer to a cs_equation_builder_t structure
796 * \param[in] t_eval time at which one performs the evaluation
797 * \param[in, out] cb pointer to a cs_cell_builder_t structure
798 */
799 /*----------------------------------------------------------------------------*/
800
801 void
cs_equation_init_reaction_properties(const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,cs_real_t t_eval,cs_cell_builder_t * cb)802 cs_equation_init_reaction_properties(const cs_equation_param_t *eqp,
803 const cs_equation_builder_t *eqb,
804 cs_real_t t_eval,
805 cs_cell_builder_t *cb)
806 {
807 assert(cs_equation_param_has_reaction(eqp));
808
809 /* Preparatory step for the reaction term(s) */
810
811 for (int i = 0; i < CS_CDO_N_MAX_REACTIONS; i++) cb->rpty_vals[i] = 1.0;
812
813 for (int r = 0; r < eqp->n_reaction_terms; r++)
814 if (eqb->reac_pty_uniform[r])
815 cb->rpty_vals[r] = cs_property_get_cell_value(0, t_eval,
816 eqp->reaction_properties[r]);
817 }
818
819 /*----------------------------------------------------------------------------*/
820 /*!
821 * \brief Initialize all reaction properties. This function is shared across
822 * all CDO schemes. The \ref cs_cell_builder_t structure stores the
823 * computed property values.
824 * If the property is uniform, a first call to the function
825 * \ref cs_equation_init_reaction_properties or to the function
826 * \ref cs_equation_init_properties has to be done before the
827 * loop on cells
828 *
829 * \param[in] eqp pointer to a cs_equation_param_t structure
830 * \param[in] eqb pointer to a cs_equation_builder_t structure
831 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
832 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
833 */
834 /*----------------------------------------------------------------------------*/
835
836 void
cs_equation_set_reaction_properties_cw(const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_cell_mesh_t * cm,cs_cell_builder_t * cb)837 cs_equation_set_reaction_properties_cw(const cs_equation_param_t *eqp,
838 const cs_equation_builder_t *eqb,
839 const cs_cell_mesh_t *cm,
840 cs_cell_builder_t *cb)
841 {
842 assert(cs_equation_param_has_reaction(eqp));
843
844 /* Set the (linear) reaction property */
845
846 cb->rpty_val = 0;
847 for (int r = 0; r < eqp->n_reaction_terms; r++)
848 if (eqb->reac_pty_uniform[r])
849 cb->rpty_val += cb->rpty_vals[r];
850 else
851 cb->rpty_val += cs_property_value_in_cell(cm,
852 eqp->reaction_properties[r],
853 cb->t_pty_eval);
854 }
855
856 /*----------------------------------------------------------------------------*/
857 /*!
858 * \brief Initialize all properties potentially useful to build the algebraic
859 * system. This function is shared across all CDO schemes.
860 * The \ref cs_cell_builder_t structure stores property values related
861 * to the reaction term, unsteady term and grad-div term.
862 *
863 * \param[in] eqp pointer to a cs_equation_param_t structure
864 * \param[in] eqb pointer to a cs_equation_builder_t structure
865 * \param[in, out] diffusion_hodge pointer to the diffusion hodge structure
866 * \param[in, out] cb pointer to a cs_cell_builder_t structure
867 */
868 /*----------------------------------------------------------------------------*/
869
870 void
cs_equation_init_properties(const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,cs_hodge_t * diffusion_hodge,cs_cell_builder_t * cb)871 cs_equation_init_properties(const cs_equation_param_t *eqp,
872 const cs_equation_builder_t *eqb,
873 cs_hodge_t *diffusion_hodge,
874 cs_cell_builder_t *cb)
875 {
876 /* Preparatory step for diffusion term
877 * One calls this function with the boundary tag to examine all tests */
878
879 if (diffusion_hodge != NULL && eqb->diff_pty_uniform)
880 cs_hodge_set_property_value(0, /* cell_id */
881 cb->t_pty_eval,
882 CS_FLAG_BOUNDARY_CELL_BY_FACE,
883 diffusion_hodge);
884
885 /* Preparatory step for the grad-div term */
886
887 if (cs_equation_param_has_graddiv(eqp) && eqb->graddiv_pty_uniform)
888 cb->gpty_val = cs_property_get_cell_value(0, cb->t_pty_eval,
889 eqp->graddiv_property);
890
891 /* Preparatory step for the unsteady term */
892
893 if (cs_equation_param_has_time(eqp) && eqb->time_pty_uniform)
894 cb->tpty_val = cs_property_get_cell_value(0, cb->t_pty_eval,
895 eqp->time_property);
896
897 /* Preparatory step for the reaction term(s) */
898
899 if (cs_equation_param_has_reaction(eqp)) {
900
901 for (int i = 0; i < CS_CDO_N_MAX_REACTIONS; i++) cb->rpty_vals[i] = 1.0;
902
903 for (int r = 0; r < eqp->n_reaction_terms; r++) {
904 if (eqb->reac_pty_uniform[r]) {
905 cb->rpty_vals[r] =
906 cs_property_get_cell_value(0, cb->t_pty_eval,
907 eqp->reaction_properties[r]);
908 }
909 } /* Loop on reaction properties */
910
911 }
912 }
913
914 /*----------------------------------------------------------------------------*/
915 /*!
916 * \brief Take into account the enforcement of internal DoFs. Apply an
917 * algebraic manipulation. Update members of the cs_cell_sys_t
918 * structure related to the internal enforcement.
919 *
920 * | | | | | | | | | |
921 * | Aii | Aie | | Aii | 0 | |bi| |bi -Aid.x_enf|
922 * |------------| --> |------------| and |--| --> |-------------|
923 * | | | | | | | | | |
924 * | Aei | Aee | | 0 | Id | |be| | x_enf |
925 *
926 * where x_enf is the value of the enforcement for the selected internal DoFs
927 *
928 * \param[in] eqb pointer to a cs_equation_builder_t structure
929 * \param[in, out] cb pointer to a cs_cell_builder_t structure
930 * \param[in, out] csys structure storing the cell-wise system
931 */
932 /*----------------------------------------------------------------------------*/
933
934 void
cs_equation_enforced_internal_dofs(const cs_equation_builder_t * eqb,cs_cell_builder_t * cb,cs_cell_sys_t * csys)935 cs_equation_enforced_internal_dofs(const cs_equation_builder_t *eqb,
936 cs_cell_builder_t *cb,
937 cs_cell_sys_t *csys)
938 {
939 /* Enforcement of internal DoFs */
940
941 double *x_vals = cb->values; /* define with cs_enforcement_dofs_cw() */
942 double *ax = cb->values + csys->n_dofs;
943
944 memset(cb->values, 0, 2*csys->n_dofs*sizeof(double));
945
946 bool do_enforcement = cs_enforcement_dofs_cw(eqb->enforced_values,
947 csys,
948 cb->values);
949
950 csys->has_internal_enforcement = do_enforcement;
951
952 if (!do_enforcement)
953 return;
954
955 /* Contribution of the DoFs which are enforced */
956
957 cs_sdm_matvec(csys->mat, x_vals, ax);
958
959 /* Second pass: Replace the block of enforced DoFs by a diagonal block */
960
961 for (int i = 0; i < csys->n_dofs; i++) {
962
963 if (csys->dof_is_forced[i]) {
964
965 /* Reset row i */
966
967 memset(csys->mat->val + csys->n_dofs*i, 0, csys->n_dofs*sizeof(double));
968
969 /* Reset column i */
970
971 for (int j = 0; j < csys->n_dofs; j++)
972 csys->mat->val[i + csys->n_dofs*j] = 0;
973 csys->mat->val[i*(1 + csys->n_dofs)] = 1;
974
975 /* Set the RHS */
976
977 csys->rhs[i] = x_vals[i];
978
979 } /* DoF associated to a Dirichlet BC */
980 else
981 csys->rhs[i] -= ax[i]; /* Update RHS */
982
983 } /* Loop on degrees of freedom */
984 }
985
986 /*----------------------------------------------------------------------------*/
987 /*!
988 * \brief Take into account the enforcement of internal DoFs. Case of matrices
989 * defined by blocks. Apply an algebraic manipulation. Update members
990 * of the cs_cell_sys_t structure related to the internal enforcement.
991 *
992 * | | | | | | | | | |
993 * | Aii | Aie | | Aii | 0 | |bi| |bi -Aid.x_enf|
994 * |------------| --> |------------| and |--| --> |-------------|
995 * | | | | | | | | | |
996 * | Aei | Aee | | 0 | Id | |be| | x_enf |
997 *
998 * where x_enf is the value of the enforcement for the selected internal DoFs
999 *
1000 * \param[in] eqb pointer to a cs_equation_builder_t structure
1001 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1002 * \param[in, out] csys structure storing the cell-wise system
1003 */
1004 /*----------------------------------------------------------------------------*/
1005
1006 void
cs_equation_enforced_internal_block_dofs(const cs_equation_builder_t * eqb,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1007 cs_equation_enforced_internal_block_dofs(const cs_equation_builder_t *eqb,
1008 cs_cell_builder_t *cb,
1009 cs_cell_sys_t *csys)
1010 {
1011 /* Enforcement of internal DoFs */
1012
1013 double *x_vals = cb->values; /* define with cs_enforcement_dofs_cw() */
1014 double *ax = cb->values + csys->n_dofs;
1015
1016 memset(cb->values, 0, 2*csys->n_dofs*sizeof(double));
1017
1018 bool do_enforcement = cs_enforcement_dofs_cw(eqb->enforced_values,
1019 csys,
1020 cb->values);
1021
1022 csys->has_internal_enforcement = do_enforcement;
1023
1024 if (!do_enforcement)
1025 return;
1026
1027 /* Contribution of the DoFs which are enforced */
1028
1029 cs_sdm_block_matvec(csys->mat, x_vals, ax);
1030
1031 /* Define the new right-hand side (rhs) */
1032
1033 for (int i = 0; i < csys->n_dofs; i++) {
1034 if (csys->dof_is_forced[i])
1035 csys->rhs[i] = x_vals[i];
1036 else
1037 csys->rhs[i] -= ax[i]; /* Update RHS */
1038 }
1039
1040 const cs_sdm_block_t *bd = csys->mat->block_desc;
1041
1042 /* Second pass: Replace the block of enforced DoFs by a diagonal block */
1043
1044 int s = 0;
1045 for (int ii = 0; ii < bd->n_row_blocks; ii++) {
1046
1047 cs_sdm_t *db = cs_sdm_get_block(csys->mat, ii, ii);
1048 const int bsize = db->n_rows*db->n_cols;
1049
1050 if (csys->dof_is_forced[s]) {
1051
1052 /* Identity for the diagonal block */
1053
1054 memset(db->val, 0, sizeof(cs_real_t)*bsize);
1055 for (int i = 0; i < db->n_rows; i++) {
1056 db->val[i*(1+db->n_rows)] = 1;
1057 assert(csys->dof_is_forced[s+i]);
1058 }
1059
1060 /* Reset column and row block jj < ii */
1061
1062 for (int jj = 0; jj < ii; jj++) {
1063
1064 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, ii, jj);
1065 memset(bij->val, 0, sizeof(cs_real_t)*bsize);
1066
1067 cs_sdm_t *bji = cs_sdm_get_block(csys->mat, jj, ii);
1068 memset(bji->val, 0, sizeof(cs_real_t)*bsize);
1069
1070 }
1071
1072 /* Reset column and row block jj < ii */
1073
1074 for (int jj = ii+1; jj < db->n_rows; jj++) {
1075
1076 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, ii, jj);
1077 memset(bij->val, 0, sizeof(cs_real_t)*bsize);
1078
1079 cs_sdm_t *bji = cs_sdm_get_block(csys->mat, jj, ii);
1080 memset(bji->val, 0, sizeof(cs_real_t)*bsize);
1081
1082 }
1083
1084 } /* DoF associated to an enforcement of their values*/
1085
1086 s += db->n_rows;
1087
1088 } /* Loop on degrees of freedom */
1089 }
1090
1091 /*----------------------------------------------------------------------------*/
1092 /*!
1093 * \brief Retrieve a pointer to a buffer of size at least the 2*n_cells
1094 * The size of the temporary buffer can be bigger according to the
1095 * numerical settings
1096 *
1097 * \return a pointer to an array of double
1098 */
1099 /*----------------------------------------------------------------------------*/
1100
1101 cs_real_t *
cs_equation_get_tmpbuf(void)1102 cs_equation_get_tmpbuf(void)
1103 {
1104 return cs_equation_common_work_buffer;
1105 }
1106
1107 /*----------------------------------------------------------------------------*/
1108 /*!
1109 * \brief Get the allocation size of the temporary buffer
1110 *
1111 * \return the size of the temporary buffer
1112 */
1113 /*----------------------------------------------------------------------------*/
1114
1115 size_t
cs_equation_get_tmpbuf_size(void)1116 cs_equation_get_tmpbuf_size(void)
1117 {
1118 return cs_equation_common_work_buffer_size;
1119 }
1120
1121 /*----------------------------------------------------------------------------*/
1122 /*!
1123 * \brief Allocate a cs_equation_balance_t structure
1124 *
1125 * \param[in] location where the balance is performed
1126 * \param[in] size size of arrays in the structure
1127 *
1128 * \return a pointer to the new allocated structure
1129 */
1130 /*----------------------------------------------------------------------------*/
1131
1132 cs_equation_balance_t *
cs_equation_balance_create(cs_flag_t location,cs_lnum_t size)1133 cs_equation_balance_create(cs_flag_t location,
1134 cs_lnum_t size)
1135 {
1136 cs_equation_balance_t *b = NULL;
1137
1138 BFT_MALLOC(b, 1, cs_equation_balance_t);
1139
1140 b->size = size;
1141 b->location = location;
1142 if (cs_flag_test(location, cs_flag_primal_cell) == false &&
1143 cs_flag_test(location, cs_flag_primal_vtx) == false)
1144 bft_error(__FILE__, __LINE__, 0, " %s: Invalid location", __func__);
1145
1146 BFT_MALLOC(b->balance, 7*size, cs_real_t);
1147 b->unsteady_term = b->balance + size;
1148 b->reaction_term = b->balance + 2*size;
1149 b->diffusion_term = b->balance + 3*size;
1150 b->advection_term = b->balance + 4*size;
1151 b->source_term = b->balance + 5*size;
1152 b->boundary_term = b->balance + 6*size;
1153
1154 /* Set to zero all members */
1155
1156 cs_equation_balance_reset(b);
1157
1158 return b;
1159 }
1160
1161 /*----------------------------------------------------------------------------*/
1162 /*!
1163 * \brief Reset a cs_equation_balance_t structure
1164 *
1165 * \param[in, out] b pointer to a cs_equation_balance_t to reset
1166 */
1167 /*----------------------------------------------------------------------------*/
1168
1169 void
cs_equation_balance_reset(cs_equation_balance_t * b)1170 cs_equation_balance_reset(cs_equation_balance_t *b)
1171 {
1172 if (b == NULL)
1173 return;
1174 if (b->size < 1)
1175 return;
1176
1177 if (b->balance == NULL)
1178 bft_error(__FILE__, __LINE__, 0, " %s: array is not allocated.", __func__);
1179
1180 size_t bufsize = b->size *sizeof(cs_real_t);
1181
1182 memset(b->balance, 0, 7*bufsize);
1183 }
1184
1185 /*----------------------------------------------------------------------------*/
1186 /*!
1187 * \brief Synchronize balance terms if this is a parallel computation
1188 *
1189 * \param[in] connect pointer to a cs_cdo_connect_t structure
1190 * \param[in, out] b pointer to a cs_equation_balance_t to rsync
1191 */
1192 /*----------------------------------------------------------------------------*/
1193
1194 void
cs_equation_balance_sync(const cs_cdo_connect_t * connect,cs_equation_balance_t * b)1195 cs_equation_balance_sync(const cs_cdo_connect_t *connect,
1196 cs_equation_balance_t *b)
1197 {
1198 if (b == NULL)
1199 bft_error(__FILE__, __LINE__, 0, " %s: structure not allocated", __func__);
1200
1201 if (cs_flag_test(b->location, cs_flag_primal_vtx)) {
1202
1203 assert(b->size == connect->n_vertices);
1204
1205 if (connect->interfaces[CS_CDO_CONNECT_VTX_SCAL] != NULL)
1206 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
1207 b->size,
1208 7, /* stride: 1 for each kind of balance */
1209 false,
1210 CS_REAL_TYPE,
1211 b->balance);
1212
1213 }
1214 }
1215
1216 /*----------------------------------------------------------------------------*/
1217 /*!
1218 * \brief Free a cs_equation_balance_t structure
1219 *
1220 * \param[in, out] p_balance pointer to the pointer to free
1221 */
1222 /*----------------------------------------------------------------------------*/
1223
1224 void
cs_equation_balance_destroy(cs_equation_balance_t ** p_balance)1225 cs_equation_balance_destroy(cs_equation_balance_t **p_balance)
1226 {
1227 cs_equation_balance_t *b = *p_balance;
1228
1229 if (b == NULL)
1230 return;
1231
1232 BFT_FREE(b->balance);
1233
1234 BFT_FREE(b);
1235 *p_balance = NULL;
1236 }
1237
1238 /*----------------------------------------------------------------------------*/
1239 /*!
1240 * \brief Synchronize the volumetric definitions to consider at each vertex
1241 *
1242 * \param[in] connect pointer to a cs_cdo_connect_t structure
1243 * \param[in] n_defs number of definitions
1244 * \param[in] defs number of times the values has been updated
1245 * \param[in, out] def2v_idx index array to define
1246 * \param[in, out] def2v_ids array of ids to define
1247 */
1248 /*----------------------------------------------------------------------------*/
1249
1250 void
cs_equation_sync_vol_def_at_vertices(const cs_cdo_connect_t * connect,int n_defs,cs_xdef_t ** defs,cs_lnum_t def2v_idx[],cs_lnum_t def2v_ids[])1251 cs_equation_sync_vol_def_at_vertices(const cs_cdo_connect_t *connect,
1252 int n_defs,
1253 cs_xdef_t **defs,
1254 cs_lnum_t def2v_idx[],
1255 cs_lnum_t def2v_ids[])
1256 {
1257 if (n_defs == 0)
1258 return;
1259
1260 const cs_lnum_t n_vertices = connect->n_vertices;
1261 const cs_adjacency_t *c2v = connect->c2v;
1262
1263 int *v2def_ids = NULL;
1264 BFT_MALLOC(v2def_ids, n_vertices, int);
1265 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
1266 for (cs_lnum_t v = 0; v < n_vertices; v++)
1267 v2def_ids[v] = -1; /* default */
1268
1269 for (int def_id = 0; def_id < n_defs; def_id++) {
1270
1271 /* Get and then set the definition of the initial condition */
1272
1273 const cs_xdef_t *def = defs[def_id];
1274 assert(def->support == CS_XDEF_SUPPORT_VOLUME);
1275
1276 if (def->meta & CS_FLAG_FULL_LOC) {
1277
1278 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
1279 for (cs_lnum_t v = 0; v < n_vertices; v++)
1280 v2def_ids[v] = def_id;
1281
1282 }
1283 else {
1284
1285 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1286
1287 for (cs_lnum_t i = 0; i < z->n_elts; i++) { /* Loop on selected cells */
1288 const cs_lnum_t c_id = z->elt_ids[i];
1289 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++)
1290 v2def_ids[c2v->ids[j]] = def_id;
1291 }
1292
1293 }
1294
1295 } /* Loop on definitions */
1296
1297 if (connect->interfaces[CS_CDO_CONNECT_VTX_SCAL] != NULL) {
1298
1299 /* Last definition is used in case of conflict */
1300
1301 cs_interface_set_max(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
1302 n_vertices,
1303 1, /* stride */
1304 false, /* interlace (not useful here) */
1305 CS_INT_TYPE, /* int */
1306 v2def_ids);
1307
1308 }
1309
1310 /* 0. Initialization */
1311
1312 cs_lnum_t *count = NULL;
1313 BFT_MALLOC(count, n_defs, cs_lnum_t);
1314 memset(count, 0, n_defs*sizeof(cs_lnum_t));
1315 memset(def2v_idx, 0, (n_defs+1)*sizeof(cs_lnum_t));
1316
1317 /* 1. Count number of vertices related to each definition */
1318
1319 for (cs_lnum_t v = 0; v < n_vertices; v++)
1320 if (v2def_ids[v] > -1)
1321 def2v_idx[v2def_ids[v]+1] += 1;
1322
1323 /* 2. Build index */
1324
1325 for (int def_id = 0; def_id < n_defs; def_id++)
1326 def2v_idx[def_id+1] += def2v_idx[def_id];
1327
1328 /* 3. Build list */
1329
1330 for (cs_lnum_t v = 0; v < n_vertices; v++) {
1331 const int def_id = v2def_ids[v];
1332 if (def_id > -1) {
1333 def2v_ids[def2v_idx[def_id] + count[def_id]] = v;
1334 count[def_id] += 1;
1335 }
1336 }
1337
1338 BFT_FREE(v2def_ids);
1339 BFT_FREE(count);
1340 }
1341
1342 /*----------------------------------------------------------------------------*/
1343 /*!
1344 * \brief Synchronize the volumetric definitions to consider at each edge
1345 *
1346 * \param[in] connect pointer to a cs_cdo_connect_t structure
1347 * \param[in] n_defs number of definitions
1348 * \param[in] defs number of times the values has been updated
1349 * \param[in, out] def2e_idx index array to define
1350 * \param[in, out] def2e_ids array of ids to define
1351 */
1352 /*----------------------------------------------------------------------------*/
1353
1354 void
cs_equation_sync_vol_def_at_edges(const cs_cdo_connect_t * connect,int n_defs,cs_xdef_t ** defs,cs_lnum_t def2e_idx[],cs_lnum_t def2e_ids[])1355 cs_equation_sync_vol_def_at_edges(const cs_cdo_connect_t *connect,
1356 int n_defs,
1357 cs_xdef_t **defs,
1358 cs_lnum_t def2e_idx[],
1359 cs_lnum_t def2e_ids[])
1360 {
1361 if (n_defs == 0)
1362 return;
1363
1364 const cs_lnum_t n_edges = connect->n_edges;
1365 const cs_adjacency_t *c2e = connect->c2e;
1366
1367 int *e2def_ids = NULL;
1368 BFT_MALLOC(e2def_ids, n_edges, int);
1369 # pragma omp parallel for if (n_edges > CS_THR_MIN)
1370 for (cs_lnum_t e = 0; e < n_edges; e++)
1371 e2def_ids[e] = -1; /* default: not associated to a definition */
1372
1373 for (int def_id = 0; def_id < n_defs; def_id++) {
1374
1375 /* Get and then set the definition of the initial condition */
1376
1377 const cs_xdef_t *def = defs[def_id];
1378 assert(def->support == CS_XDEF_SUPPORT_VOLUME);
1379
1380 if (def->meta & CS_FLAG_FULL_LOC) {
1381
1382 # pragma omp parallel for if (n_edges > CS_THR_MIN)
1383 for (cs_lnum_t e = 0; e < n_edges; e++)
1384 e2def_ids[e] = def_id;
1385
1386 }
1387 else {
1388
1389 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1390
1391 for (cs_lnum_t i = 0; i < z->n_elts; i++) { /* Loop on selected cells */
1392 const cs_lnum_t c_id = z->elt_ids[i];
1393 for (cs_lnum_t j = c2e->idx[c_id]; j < c2e->idx[c_id+1]; j++)
1394 e2def_ids[c2e->ids[j]] = def_id;
1395 }
1396
1397 }
1398
1399 } /* Loop on definitions */
1400
1401 if (connect->interfaces[CS_CDO_CONNECT_EDGE_SCAL] != NULL) {
1402
1403 /* Last definition is used in case of conflict */
1404
1405 cs_interface_set_max(connect->interfaces[CS_CDO_CONNECT_EDGE_SCAL],
1406 n_edges,
1407 1, /* stride */
1408 false, /* interlace (not useful here) */
1409 CS_INT_TYPE, /* int */
1410 e2def_ids);
1411
1412 }
1413
1414 /* 0. Initialization */
1415
1416 cs_lnum_t *count = NULL;
1417 BFT_MALLOC(count, n_defs, cs_lnum_t);
1418 memset(count, 0, n_defs*sizeof(cs_lnum_t));
1419 memset(def2e_idx, 0, (n_defs+1)*sizeof(cs_lnum_t));
1420
1421 /* 1. Count the number of edges related to each definition */
1422
1423 for (cs_lnum_t e = 0; e < n_edges; e++)
1424 if (e2def_ids[e] > -1)
1425 def2e_idx[e2def_ids[e]+1] += 1;
1426
1427 /* 2. Build the index */
1428
1429 for (int def_id = 0; def_id < n_defs; def_id++)
1430 def2e_idx[def_id+1] += def2e_idx[def_id];
1431
1432 /* 3. Build the list */
1433
1434 for (cs_lnum_t e = 0; e < n_edges; e++) {
1435 const int def_id = e2def_ids[e];
1436 if (def_id > -1) {
1437 def2e_ids[def2e_idx[def_id] + count[def_id]] = e;
1438 count[def_id] += 1;
1439 }
1440 }
1441
1442 BFT_FREE(e2def_ids);
1443 BFT_FREE(count);
1444 }
1445
1446 /*----------------------------------------------------------------------------*/
1447 /*!
1448 * \brief Synchronize the volumetric definitions to consider at each face
1449 *
1450 * \param[in] connect pointer to a cs_cdo_connect_t structure
1451 * \param[in] n_defs number of definitions
1452 * \param[in] defs number of times the values has been updated
1453 * \param[in, out] def2f_idx index array to define
1454 * \param[in, out] def2f_ids array of ids to define
1455 */
1456 /*----------------------------------------------------------------------------*/
1457
1458 void
cs_equation_sync_vol_def_at_faces(const cs_cdo_connect_t * connect,int n_defs,cs_xdef_t ** defs,cs_lnum_t def2f_idx[],cs_lnum_t def2f_ids[])1459 cs_equation_sync_vol_def_at_faces(const cs_cdo_connect_t *connect,
1460 int n_defs,
1461 cs_xdef_t **defs,
1462 cs_lnum_t def2f_idx[],
1463 cs_lnum_t def2f_ids[])
1464 {
1465 if (n_defs == 0)
1466 return;
1467
1468 const cs_lnum_t n_faces = connect->n_faces[CS_ALL_FACES];
1469 const cs_adjacency_t *c2f = connect->c2f;
1470
1471 int *f2def_ids = NULL;
1472 BFT_MALLOC(f2def_ids, n_faces, int);
1473 # pragma omp parallel for if (n_faces > CS_THR_MIN)
1474 for (cs_lnum_t f = 0; f < n_faces; f++)
1475 f2def_ids[f] = -1; /* default */
1476
1477 for (int def_id = 0; def_id < n_defs; def_id++) {
1478
1479 /* Get and then set the definition of the initial condition */
1480
1481 const cs_xdef_t *def = defs[def_id];
1482 assert(def->support == CS_XDEF_SUPPORT_VOLUME);
1483
1484 if (def->meta & CS_FLAG_FULL_LOC) {
1485
1486 # pragma omp parallel for if (n_faces > CS_THR_MIN)
1487 for (cs_lnum_t f = 0; f < n_faces; f++)
1488 f2def_ids[f] = def_id;
1489
1490 }
1491 else {
1492
1493 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1494
1495 for (cs_lnum_t i = 0; i < z->n_elts; i++) { /* Loop on selected cells */
1496 const cs_lnum_t c_id = z->elt_ids[i];
1497 for (cs_lnum_t j = c2f->idx[c_id]; j < c2f->idx[c_id+1]; j++)
1498 f2def_ids[c2f->ids[j]] = def_id;
1499 }
1500
1501 }
1502
1503 } /* Loop on definitions */
1504
1505 if (connect->interfaces[CS_CDO_CONNECT_FACE_SP0] != NULL) {
1506
1507 /* Last definition is used in case of conflict */
1508
1509 cs_interface_set_max(connect->interfaces[CS_CDO_CONNECT_FACE_SP0],
1510 n_faces,
1511 1, /* stride */
1512 false, /* interlace (not useful here) */
1513 CS_INT_TYPE, /* int */
1514 f2def_ids);
1515
1516 }
1517
1518 /* 0. Initialization */
1519
1520 cs_lnum_t *count = NULL;
1521 BFT_MALLOC(count, n_defs, cs_lnum_t);
1522 memset(count, 0, n_defs*sizeof(cs_lnum_t));
1523 memset(def2f_idx, 0, (n_defs+1)*sizeof(cs_lnum_t));
1524
1525 /* 1. Count number of faces related to each definition */
1526
1527 for (cs_lnum_t f = 0; f < n_faces; f++)
1528 if (f2def_ids[f] > -1)
1529 def2f_idx[f2def_ids[f]+1] += 1;
1530
1531 /* 2. Build index */
1532
1533 for (int def_id = 0; def_id < n_defs; def_id++)
1534 def2f_idx[def_id+1] += def2f_idx[def_id];
1535
1536 /* 3. Build list */
1537
1538 for (cs_lnum_t f = 0; f < n_faces; f++) {
1539 const int def_id = f2def_ids[f];
1540 if (def_id > -1) {
1541 def2f_ids[def2f_idx[def_id] + count[def_id]] = f;
1542 count[def_id] += 1;
1543 }
1544 }
1545
1546 BFT_FREE(f2def_ids);
1547 BFT_FREE(count);
1548 }
1549
1550 /*----------------------------------------------------------------------------*/
1551 /*!
1552 * \brief Compute the mean-value across ranks at each vertex
1553 *
1554 * \param[in] connect pointer to a cs_cdo_connect_t structure
1555 * \param[in] dim number of entries for each vertex
1556 * \param[in] counter number of occurences on this rank
1557 * \param[in, out] values array to update
1558 */
1559 /*----------------------------------------------------------------------------*/
1560
1561 void
cs_equation_sync_vertex_mean_values(const cs_cdo_connect_t * connect,int dim,int * counter,cs_real_t * values)1562 cs_equation_sync_vertex_mean_values(const cs_cdo_connect_t *connect,
1563 int dim,
1564 int *counter,
1565 cs_real_t *values)
1566 {
1567 const cs_lnum_t n_vertices = connect->n_vertices;
1568
1569 if (connect->interfaces[CS_CDO_CONNECT_VTX_SCAL] != NULL) {
1570
1571 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
1572 n_vertices,
1573 1, /* stride */
1574 false, /* interlace (not useful here) */
1575 CS_INT_TYPE, /* int */
1576 counter);
1577
1578 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
1579 n_vertices,
1580 dim, /* stride */
1581 true, /* interlace */
1582 CS_REAL_TYPE,
1583 values);
1584
1585 }
1586
1587 if (dim == 1) {
1588
1589 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
1590 for (cs_lnum_t v_id = 0; v_id < n_vertices; v_id++)
1591 if (counter[v_id] > 1)
1592 values[v_id] /= counter[v_id];
1593
1594 }
1595 else { /* dim > 1 */
1596
1597 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
1598 for (cs_lnum_t v_id = 0; v_id < n_vertices; v_id++) {
1599 if (counter[v_id] > 1) {
1600 const cs_real_t inv_count = 1./counter[v_id];
1601 for (int k = 0; k < dim; k++)
1602 values[dim*v_id + k] *= inv_count;
1603 }
1604 }
1605
1606 }
1607 }
1608
1609 /*----------------------------------------------------------------------------*/
1610
1611 END_C_DECLS
1612