1 /*============================================================================
2 * Functions to handle the cs_equation_t structure and its related structures
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 <assert.h>
34 #include <ctype.h>
35 #include <float.h>
36 #include <math.h>
37 #include <stdlib.h>
38 #include <string.h>
39
40 #if defined(HAVE_MPI)
41 #include <mpi.h>
42 #endif
43
44 /*----------------------------------------------------------------------------
45 * Local headers
46 *----------------------------------------------------------------------------*/
47
48 #include <bft_mem.h>
49
50 #include "cs_cdovb_scaleq.h"
51 #include "cs_cdovb_vecteq.h"
52 #include "cs_cdovcb_scaleq.h"
53 #include "cs_cdoeb_vecteq.h"
54 #include "cs_cdofb_scaleq.h"
55 #include "cs_cdofb_vecteq.h"
56 #include "cs_equation_bc.h"
57 #include "cs_equation_common.h"
58 #include "cs_equation_priv.h"
59 #include "cs_evaluate.h"
60 #include "cs_field_default.h"
61 #include "cs_hho_scaleq.h"
62 #include "cs_hho_vecteq.h"
63 #include "cs_log.h"
64 #include "cs_parall.h"
65 #include "cs_post.h"
66 #include "cs_prototypes.h"
67 #include "cs_range_set.h"
68 #include "cs_sles.h"
69 #include "cs_timer_stats.h"
70
71 #if defined(DEBUG) && !defined(NDEBUG)
72 #include "cs_dbg.h"
73 #endif
74
75 /*----------------------------------------------------------------------------
76 * Header for the current file
77 *----------------------------------------------------------------------------*/
78
79 #include "cs_equation.h"
80
81 /*----------------------------------------------------------------------------*/
82
83 BEGIN_C_DECLS
84
85 /*=============================================================================
86 * Local macro definitions
87 *============================================================================*/
88
89 #define CS_EQUATION_DBG 0
90
91 /*============================================================================
92 * Type definitions
93 *============================================================================*/
94
95 /*============================================================================
96 * Local variables
97 *============================================================================*/
98
99 static int _n_equations = 0;
100 static int _n_predef_equations = 0;
101 static int _n_user_equations = 0;
102 static cs_equation_t **_equations = NULL;
103
104 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
105
106 /*============================================================================
107 * Private variables
108 *============================================================================*/
109
110 static const char _err_empty_eq[] =
111 N_(" %s: Stop setting an empty cs_equation_t structure.\n"
112 " Please check your settings.\n");
113
114 /*============================================================================
115 * Prototypes for functions intended for use only by Fortran wrappers.
116 * (descriptions follow, with function bodies).
117 *============================================================================*/
118
119 void
120 solve_steady_state_cdo_equation(const char *eqname);
121
122 void
123 solve_cdo_equation(bool cur2prev,
124 const char *eqname);
125
126 /*============================================================================
127 * Private function prototypes
128 *============================================================================*/
129
130 /*----------------------------------------------------------------------------*/
131 /*!
132 * \brief Set the initial values for the variable related to an equation
133 *
134 * \param[in] eq pointer to a cs_equation_t structure
135 * \param[in] ts pointer to cs_time_step_t structure
136 * \param[in] tag tag to add to the equation name to build the label
137 * \param[in, out] label label for the postprocessing
138 * \param[in] values pointer to the array of values to post
139 */
140 /*----------------------------------------------------------------------------*/
141
142 static inline void
_post_balance_at_vertices(const cs_equation_t * eq,const cs_time_step_t * ts,const char * tag,char * label,const cs_real_t * values)143 _post_balance_at_vertices(const cs_equation_t *eq,
144 const cs_time_step_t *ts,
145 const char *tag,
146 char *label,
147 const cs_real_t *values)
148
149 {
150 sprintf(label, "%s.Balance.%s", eq->param->name, tag);
151
152 cs_post_write_vertex_var(CS_POST_MESH_VOLUME,
153 CS_POST_WRITER_DEFAULT,
154 label,
155 eq->param->dim,
156 false,
157 false,
158 CS_POST_TYPE_cs_real_t,
159 values,
160 ts);
161 }
162
163 /*----------------------------------------------------------------------------*/
164 /*!
165 * \brief Carry out operations for allocating and/or initializing the solution
166 * array and the right hand side of the linear system to solve.
167 * Handle parallelism thanks to cs_range_set_t structure.
168 *
169 * \param[in, out] eq_to_cast pointer to generic builder structure
170 * \param[in, out] p_x pointer of pointer to the solution array
171 * \param[in, out] p_rhs pointer of pointer to the RHS array
172 */
173 /*----------------------------------------------------------------------------*/
174
175 static void
_prepare_fb_solving(void * eq_to_cast,cs_real_t * p_x[],cs_real_t * p_rhs[])176 _prepare_fb_solving(void *eq_to_cast,
177 cs_real_t *p_x[],
178 cs_real_t *p_rhs[])
179 {
180 cs_equation_t *eq = (cs_equation_t *)eq_to_cast;
181
182 const cs_real_t *f_values = eq->get_face_values(eq->scheme_context, false);
183 const int stride = 1; /* Since the global numbering is adapted in each
184 case (scalar-, vector-valued equations) */
185
186 /* Sanity check */
187 assert(f_values != NULL);
188
189 cs_real_t *x = NULL, *b = NULL;
190 BFT_MALLOC(x, CS_MAX(eq->n_sles_scatter_elts,
191 cs_matrix_get_n_columns(eq->matrix)), cs_real_t);
192
193 /* x and b are a "gathered" view of field->val and eq->rhs respectively
194 through the range set operation.
195 Their size is equal to n_sles_gather_elts <= n_sles_scatter_elts */
196
197 if (cs_glob_n_ranks > 1) { /* Parallel mode */
198
199 /* Compact numbering to fit the algebraic decomposition */
200 cs_range_set_gather(eq->rset,
201 CS_REAL_TYPE, /* type */
202 stride, /* stride */
203 f_values, /* in: size = n_sles_scatter_elts */
204 x); /* out: size = n_sles_gather_elts */
205
206 /* The right-hand side stems from a cellwise building on this rank.
207 Other contributions from distant ranks may contribute to an element
208 owned by the local rank */
209 BFT_MALLOC(b, eq->n_sles_scatter_elts, cs_real_t);
210
211 #if defined(HAVE_OPENMP)
212 # pragma omp parallel for if (eq->n_sles_scatter_elts > CS_THR_MIN)
213 for (cs_lnum_t i = 0; i < eq->n_sles_scatter_elts; i++)
214 b[i] = eq->rhs[i];
215 #else
216 memcpy(b, eq->rhs, eq->n_sles_scatter_elts * sizeof(cs_real_t));
217 #endif
218
219 cs_interface_set_sum(eq->rset->ifs,
220 eq->n_sles_scatter_elts, stride, false, CS_REAL_TYPE,
221 b);
222
223 cs_range_set_gather(eq->rset,
224 CS_REAL_TYPE, /* type */
225 stride, /* stride */
226 b, /* in: size = n_sles_scatter_elts */
227 b); /* out: size = n_sles_gather_elts */
228
229 }
230 else { /* Serial mode *** without periodicity *** */
231
232 assert(eq->n_sles_gather_elts == eq->n_sles_scatter_elts);
233
234 #if defined(HAVE_OPENMP)
235 # pragma omp parallel for if (eq->n_sles_scatter_elts > CS_THR_MIN)
236 for (cs_lnum_t i = 0; i < eq->n_sles_scatter_elts; i++)
237 x[i] = f_values[i];
238 #else
239 memcpy(x, f_values, eq->n_sles_scatter_elts * sizeof(cs_real_t));
240 #endif
241
242 /* Nothing to do for the right-hand side */
243 b = eq->rhs;
244
245 }
246
247 /* Return pointers */
248 *p_x = x;
249 *p_rhs = b;
250 }
251
252 /*----------------------------------------------------------------------------*/
253 /*!
254 * \brief Set the pointers of function for the given equation.
255 * Case of scalar-valued HHO schemes
256 *
257 * \param[in, out] eq pointer to a cs_equation_t structure
258 */
259 /*----------------------------------------------------------------------------*/
260
261 static void
_set_scal_hho_function_pointers(cs_equation_t * eq)262 _set_scal_hho_function_pointers(cs_equation_t *eq)
263 {
264 if (eq == NULL)
265 return;
266
267 eq->init_context = cs_hho_scaleq_init_context;
268 eq->free_context = cs_hho_scaleq_free_context;
269
270 /* New functions */
271 eq->init_field_values = cs_hho_scaleq_init_values;
272 eq->solve = NULL;
273 eq->solve_steady_state = NULL;
274
275 eq->postprocess = cs_hho_scaleq_extra_post;
276 eq->read_restart = cs_hho_scaleq_read_restart;
277 eq->write_restart = cs_hho_scaleq_write_restart;
278
279 /* Function pointers to retrieve values at mesh locations */
280 eq->get_vertex_values = NULL;
281 eq->get_edge_values = NULL;
282 eq->get_face_values = cs_hho_scaleq_get_face_values;
283 eq->get_cell_values = cs_hho_scaleq_get_cell_values;
284
285 /* Deprecated functions */
286 eq->initialize_system = cs_hho_scaleq_initialize_system;
287 eq->set_dir_bc = NULL;
288 eq->build_system = cs_hho_scaleq_build_system;
289 eq->prepare_solving = _prepare_fb_solving;
290 eq->update_field = cs_hho_scaleq_update_field;
291 }
292
293 /*----------------------------------------------------------------------------*/
294 /*!
295 * \brief Set the pointers of function for the given equation.
296 * Case of vector-valued HHO schemes
297 *
298 * \param[in, out] eq pointer to a cs_equation_t structure
299 */
300 /*----------------------------------------------------------------------------*/
301
302 static void
_set_vect_hho_function_pointers(cs_equation_t * eq)303 _set_vect_hho_function_pointers(cs_equation_t *eq)
304 {
305 if (eq == NULL)
306 return;
307
308 eq->init_context = cs_hho_vecteq_init_context;
309 eq->free_context = cs_hho_vecteq_free_context;
310
311 /* New functions */
312 eq->init_field_values = cs_hho_vecteq_init_values;
313 eq->solve = NULL;
314 eq->solve_steady_state = NULL;
315
316 eq->postprocess = cs_hho_vecteq_extra_post;
317 eq->read_restart = cs_hho_vecteq_read_restart;
318 eq->write_restart = cs_hho_vecteq_write_restart;
319
320 /* Function pointers to retrieve values at mesh locations */
321 eq->get_vertex_values = NULL;
322 eq->get_edge_values = NULL;
323 eq->get_face_values = cs_hho_vecteq_get_face_values;
324 eq->get_cell_values = cs_hho_vecteq_get_cell_values;
325
326 /* Deprecated functions */
327 eq->initialize_system = cs_hho_vecteq_initialize_system;
328 eq->build_system = cs_hho_vecteq_build_system;
329 eq->prepare_solving = _prepare_fb_solving;
330 eq->update_field = cs_hho_vecteq_update_field;
331 }
332
333 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
334
335 /*============================================================================
336 * Public function prototypes
337 *============================================================================*/
338
339 /*----------------------------------------------------------------------------*/
340 /*!
341 * \brief Retrieve the number of equations
342 *
343 * \return the current number of cs_equation_t structure allocated
344 */
345 /*----------------------------------------------------------------------------*/
346
347 int
cs_equation_get_n_equations(void)348 cs_equation_get_n_equations(void)
349 {
350 return _n_equations;
351 }
352
353 /*----------------------------------------------------------------------------*/
354 /*!
355 * \brief Find the cs_equation_t structure with name eqname
356 * Return NULL if not find
357 *
358 * \param[in] eqname name of the equation to find
359 *
360 * \return a pointer to a cs_equation_t structure or NULL if not found
361 */
362 /*----------------------------------------------------------------------------*/
363
364 cs_equation_t *
cs_equation_by_name(const char * eqname)365 cs_equation_by_name(const char *eqname)
366 {
367 cs_equation_t *eq = NULL;
368 if (eqname == NULL)
369 return eq;
370
371 size_t len_in = strlen(eqname);
372 for (int i = 0; i < _n_equations; i++) {
373
374 cs_equation_t *_eq = _equations[i];
375 assert(_eq != NULL);
376 cs_equation_param_t *eqp = _eq->param;
377 assert(eqp != NULL && eqp->name != NULL);
378 if (strlen(eqp->name) == len_in)
379 if (strcmp(eqname, eqp->name) == 0)
380 return _eq;
381
382 }
383
384 return eq;
385 }
386
387 /*----------------------------------------------------------------------------*/
388 /*!
389 * \brief Return the pointer to a cs_equation_t structure thanks to the field
390 * name of the variable field associated to a cs_equation_t structure
391 *
392 * \param[in] field_name name of the field
393 *
394 * \return a pointer to a cs_equation_t structure or NULL if not found
395 */
396 /*----------------------------------------------------------------------------*/
397
398 cs_equation_t *
cs_equation_by_field_name(const char * field_name)399 cs_equation_by_field_name(const char *field_name)
400 {
401 if (field_name == NULL)
402 return NULL;
403
404 for (int i = 0; i < _n_equations; i++) {
405
406 cs_equation_t *eq = _equations[i];
407 assert(eq != NULL);
408
409 if (cs_equation_has_field_name(eq, field_name))
410 return eq;
411
412 } /* Loop on equations */
413
414 return NULL;
415 }
416
417 /*----------------------------------------------------------------------------*/
418 /*!
419 * \brief Check if the asociated field to a \ref cs_equation_t structure
420 * has name equal to fld_name
421 *
422 * \param[in] eq pointer to a \ref cs_equation_t structure to test
423 * \param[in] fld_name name of the field
424 *
425 * \return true if the \ref cs_equation_t structure has an associated field
426 * named fld_name, otherwise false
427 */
428 /*----------------------------------------------------------------------------*/
429
430 bool
cs_equation_has_field_name(const cs_equation_t * eq,const char * fld_name)431 cs_equation_has_field_name(const cs_equation_t *eq,
432 const char *fld_name)
433 {
434 if (eq == NULL)
435 return false;
436
437 cs_field_t *fld = cs_field_by_id(eq->field_id);
438 if (fld == NULL)
439 return false;
440
441 if (strcmp(fld->name, fld_name) == 0)
442 return true;
443 else
444 return false;
445 }
446
447 /*----------------------------------------------------------------------------*/
448 /*!
449 * \brief Return the cs_equation_param_t structure associated to a
450 * cs_equation_t structure based on the equation name
451 *
452 * If no equation matches the given name but a field does, equation
453 * parameter structure associated to the field will be returned instead.
454 * This allows using this function with non-CDO (legacy) fields.
455 *
456 * \param[in] eqname name of the equation
457 *
458 * \return a cs_equation_param_t structure or NULL if not found
459 */
460 /*----------------------------------------------------------------------------*/
461
462 cs_equation_param_t *
cs_equation_param_by_name(const char * eqname)463 cs_equation_param_by_name(const char *eqname)
464 {
465 cs_equation_param_t *eq_param = NULL;
466
467 if (eqname != NULL) {
468
469 cs_equation_t *eq = cs_equation_by_name(eqname);
470 if (eq != NULL)
471 eq_param = eq->param;
472
473 else {
474 cs_field_t *f = cs_field_by_name_try(eqname);
475 if (f != NULL)
476 eq_param = cs_field_get_equation_param(f);
477 }
478
479 }
480
481 return eq_param;
482 }
483
484 /*----------------------------------------------------------------------------*/
485 /*!
486 * \brief Return the cs_equation_param_t structure related to a
487 * cs_equation_t structure thanks to the field name of the variable
488 * field associated to a cs_equation_t structure
489 *
490 * \param[in] field_name name of the field
491 *
492 * \return a cs_equation_param_t structure or NULL if not found
493 */
494 /*----------------------------------------------------------------------------*/
495
496 cs_equation_param_t *
cs_equation_param_by_field_name(const char * field_name)497 cs_equation_param_by_field_name(const char *field_name)
498 {
499 if (field_name == NULL)
500 return NULL;
501
502 cs_equation_t *eq = cs_equation_by_field_name(field_name);
503
504 if (eq == NULL)
505 return NULL;
506 else
507 return eq->param;
508 }
509
510 /*----------------------------------------------------------------------------*/
511 /*!
512 * \brief Return the cs_equation_param_t structure associated to a
513 * cs_equation_t structure
514 *
515 * \param[in] eq pointer to a cs_equation_t structure
516 *
517 * \return a cs_equation_param_t structure or NULL if not found
518 */
519 /*----------------------------------------------------------------------------*/
520
521 cs_equation_param_t *
cs_equation_get_param(const cs_equation_t * eq)522 cs_equation_get_param(const cs_equation_t *eq)
523 {
524 if (eq == NULL)
525 return NULL;
526 else
527 return eq->param;
528 }
529
530 /*----------------------------------------------------------------------------*/
531 /*!
532 * \brief Find the \ref cs_equation_t structure with id eq_id
533 * Return NULL if not find
534 *
535 * \param[in] eq_id id of the equation to find
536 *
537 * \return a pointer to a \ref cs_equation_t structure or NULL if not found
538 */
539 /*----------------------------------------------------------------------------*/
540
541 cs_equation_t *
cs_equation_by_id(int eq_id)542 cs_equation_by_id(int eq_id)
543 {
544 if (eq_id < 0 || eq_id > _n_equations - 1)
545 return NULL;
546
547 return _equations[eq_id];
548 }
549
550
551 /*----------------------------------------------------------------------------*/
552 /*!
553 * \brief Return the name related to the given cs_equation_t structure
554 *
555 * \param[in] eq pointer to a cs_equation_t structure
556 *
557 * \return a name or NULL if not found
558 */
559 /*----------------------------------------------------------------------------*/
560
561 const char *
cs_equation_get_name(const cs_equation_t * eq)562 cs_equation_get_name(const cs_equation_t *eq)
563 {
564 cs_equation_param_t *eqp = cs_equation_get_param(eq);
565 if (eqp == NULL)
566 return NULL;
567 else
568 return eqp->name;
569 }
570
571 /*----------------------------------------------------------------------------*/
572 /*!
573 * \brief Return the id number related to the given cs_equation_t structure
574 *
575 * \param[in] eq pointer to a cs_equation_t structure
576 *
577 * \return an id (0 ... n-1) or -1 if not found
578 */
579 /*----------------------------------------------------------------------------*/
580
581 int
cs_equation_get_id(const cs_equation_t * eq)582 cs_equation_get_id(const cs_equation_t *eq)
583 {
584 if (eq == NULL)
585 return -1;
586 else
587 return eq->id;
588 }
589
590
591 /*----------------------------------------------------------------------------*/
592 /*!
593 * \brief Return the field structure associated to a cs_equation_t structure
594 *
595 * \param[in] eq pointer to a cs_equation_t structure
596 *
597 * \return a cs_field_t structure or NULL if not found
598 */
599 /*----------------------------------------------------------------------------*/
600
601 cs_field_t *
cs_equation_get_field(const cs_equation_t * eq)602 cs_equation_get_field(const cs_equation_t *eq)
603 {
604 if (eq == NULL)
605 return NULL;
606 else
607 return cs_field_by_id(eq->field_id);
608 }
609
610 /*----------------------------------------------------------------------------*/
611 /*!
612 * \brief Return the id related to the variable field structure associated to
613 * the cs_equation_t structure
614 *
615 * \param[in] eq pointer to a cs_equation_t structure
616 *
617 * \return an integer (-1 if the field is not defined)
618 */
619 /*----------------------------------------------------------------------------*/
620
621 int
cs_equation_get_field_id(const cs_equation_t * eq)622 cs_equation_get_field_id(const cs_equation_t *eq)
623 {
624 if (eq == NULL)
625 return -1;
626 else
627 return eq->field_id;
628 }
629
630 /*----------------------------------------------------------------------------*/
631 /*!
632 * \brief Return the range set structure associated to a cs_equation_t
633 * structure
634 *
635 * \param[in] eq pointer to a cs_equation_t structure
636 *
637 * \return a cs_range_set_t structure or NULL if not found
638 */
639 /*----------------------------------------------------------------------------*/
640
641 const cs_range_set_t *
cs_equation_get_range_set(const cs_equation_t * eq)642 cs_equation_get_range_set(const cs_equation_t *eq)
643 {
644 if (eq == NULL)
645 return NULL;
646 else
647 return eq->rset;
648 }
649
650 /*----------------------------------------------------------------------------*/
651 /*!
652 * \brief Return the global number of degrees of freedom associated to this
653 * cs_equation_t structure
654 *
655 * \param[in] eq pointer to a cs_equation_t structure
656 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
657 *
658 * \return a global number of degrees of freedom (DoFs)
659 */
660 /*----------------------------------------------------------------------------*/
661
662 cs_gnum_t
cs_equation_get_global_n_dofs(const cs_equation_t * eq,const cs_cdo_quantities_t * cdoq)663 cs_equation_get_global_n_dofs(const cs_equation_t *eq,
664 const cs_cdo_quantities_t *cdoq)
665 {
666 if (eq == NULL || cdoq == NULL)
667 return 0;
668
669 switch (cs_equation_get_space_scheme(eq)) {
670
671 case CS_SPACE_SCHEME_CDOVB:
672 if (cs_glob_n_ranks > 1)
673 return cdoq->n_g_vertices;
674 else
675 return (cs_gnum_t)cdoq->n_vertices;
676 break;
677
678 case CS_SPACE_SCHEME_CDOVCB:
679 if (cs_glob_n_ranks > 1)
680 return cdoq->n_g_vertices + cdoq->n_g_cells;
681 else
682 return (cs_gnum_t)(cdoq->n_vertices + cdoq->n_cells);
683 break;
684
685 case CS_SPACE_SCHEME_CDOEB:
686 if (cs_glob_n_ranks > 1)
687 return cdoq->n_g_edges;
688 else
689 return (cs_gnum_t)(cdoq->n_edges);
690 break;
691
692 case CS_SPACE_SCHEME_CDOFB:
693 case CS_SPACE_SCHEME_HHO_P0:
694 if (cs_glob_n_ranks > 1)
695 return cdoq->n_g_faces + cdoq->n_g_cells;
696 else
697 return (cs_gnum_t)(cdoq->n_faces + cdoq->n_cells);
698 break;
699
700 case CS_SPACE_SCHEME_HHO_P1:
701 if (cs_glob_n_ranks > 1)
702 return CS_N_FACE_DOFS_1ST*cdoq->n_g_faces
703 + CS_N_CELL_DOFS_1ST*cdoq->n_g_cells;
704 else
705 return (cs_gnum_t)(CS_N_FACE_DOFS_1ST*cdoq->n_faces
706 + CS_N_CELL_DOFS_1ST*cdoq->n_cells);
707 break;
708
709 case CS_SPACE_SCHEME_HHO_P2:
710 if (cs_glob_n_ranks > 1)
711 return CS_N_FACE_DOFS_2ND*cdoq->n_g_faces
712 + CS_N_CELL_DOFS_2ND*cdoq->n_g_cells;
713 else
714 return (cs_gnum_t)(CS_N_FACE_DOFS_2ND*cdoq->n_faces
715 + CS_N_CELL_DOFS_2ND*cdoq->n_cells);
716 break;
717
718 default:
719 return 0;
720 }
721
722 }
723
724 /*----------------------------------------------------------------------------*/
725 /*!
726 * \brief Return the field structure for the (normal) boundary flux associated
727 * to a cs_equation_t structure
728 *
729 * \param[in] eq pointer to a cs_equation_t structure
730 *
731 * \return a cs_field_t structure or NULL
732 */
733 /*----------------------------------------------------------------------------*/
734
735 cs_field_t *
cs_equation_get_boundary_flux(const cs_equation_t * eq)736 cs_equation_get_boundary_flux(const cs_equation_t *eq)
737 {
738 if (eq == NULL)
739 return NULL;
740 else if (eq->boundary_flux_id < 0)
741 return NULL;
742 else
743 return cs_field_by_id(eq->boundary_flux_id);
744 }
745
746 /*----------------------------------------------------------------------------*/
747 /*!
748 * \brief Return the flag associated to an equation
749 *
750 * \param[in] eq pointer to a cs_equation_t structure
751 *
752 * \return a flag (cs_flag_t type)
753 */
754 /*----------------------------------------------------------------------------*/
755
756 cs_flag_t
cs_equation_get_flag(const cs_equation_t * eq)757 cs_equation_get_flag(const cs_equation_t *eq)
758 {
759 cs_flag_t ret_flag = 0;
760
761 if (eq == NULL)
762 return ret_flag;
763
764 ret_flag = eq->param->flag;
765
766 return ret_flag;
767 }
768
769 /*----------------------------------------------------------------------------*/
770 /*!
771 * \brief Redefine the flag associated to an equation
772 *
773 * \param[in, out] eq pointer to a cs_equation_t structure
774 * \param[in] flag new flag to set
775 */
776 /*----------------------------------------------------------------------------*/
777
778 void
cs_equation_set_flag(cs_equation_t * eq,cs_flag_t flag)779 cs_equation_set_flag(cs_equation_t *eq,
780 cs_flag_t flag)
781 {
782 if (eq == NULL)
783 bft_error(__FILE__, __LINE__, 0, _err_empty_eq, __func__);
784 assert(eq->param != NULL);
785
786 eq->param->flag = flag;;
787 }
788
789 /*----------------------------------------------------------------------------*/
790 /*!
791 * \brief Add a user hook to enable an advanced user to get a fine control of
792 * the cellwise system building.
793 * Only for an advanced usage. The context may be set to NULL if there
794 * is no need to get additional information.
795 *
796 * \param[in, out] eq pointer to the cs_equation_t stucture to update
797 * \param[in] context pointer to a structure for additional information
798 * \param[in] func pointer to the user function
799 */
800 /*----------------------------------------------------------------------------*/
801
802 void
cs_equation_add_user_hook(cs_equation_t * eq,void * context,cs_equation_user_hook_t * func)803 cs_equation_add_user_hook(cs_equation_t *eq,
804 void *context,
805 cs_equation_user_hook_t *func)
806 {
807 if (eq == NULL)
808 return;
809
810 cs_equation_param_t *eqp = eq->param;
811 assert(eqp != NULL);
812
813 if (eq->builder == NULL)
814 bft_error(__FILE__, __LINE__, 0,
815 " %s: Initialization of equation %s has not been done yet.\n"
816 " Please call this operation later in"
817 " cs_user_extra_operations_initialize() for instance.",
818 __func__, eqp->name);
819
820 cs_equation_builder_t *eqb = eq->builder;
821
822 eqb->user_hook_context = context;
823 eqb->user_hook_function = func;
824 eqp->flag |= CS_EQUATION_USER_HOOK;
825
826 /* Add an entry in the setup log file (this is done after the main setup
827 * log but one needs to initialize equations before calling this function) */
828 cs_log_printf(CS_LOG_SETUP, " Equation %s: Add a user hook function\n",
829 eqp->name);
830 }
831
832 /*----------------------------------------------------------------------------*/
833 /*!
834 * \brief Return the cs_equation_builder_t structure associated to a
835 * cs_equation_t structure. Only for an advanced usage.
836 *
837 * \param[in] eq pointer to a cs_equation_t structure
838 *
839 * \return a cs_equation_builder_t structure or NULL if not found
840 */
841 /*----------------------------------------------------------------------------*/
842
843 cs_equation_builder_t *
cs_equation_get_builder(const cs_equation_t * eq)844 cs_equation_get_builder(const cs_equation_t *eq)
845 {
846 if (eq == NULL)
847 return NULL;
848 else
849 return eq->builder;
850 }
851
852 /*----------------------------------------------------------------------------*/
853 /*!
854 * \brief Return a pointer to a structure useful to handle low-level
855 * operations for the given equation
856 *
857 * \param[in] eq pointer to a cs_equation_t structure
858 *
859 * \return a pointer to a structure to cast on-the-fly or NULL if not found
860 */
861 /*----------------------------------------------------------------------------*/
862
863 void *
cs_equation_get_scheme_context(const cs_equation_t * eq)864 cs_equation_get_scheme_context(const cs_equation_t *eq)
865 {
866 if (eq == NULL)
867 return NULL;
868 else
869 return eq->scheme_context;
870 }
871
872 /*----------------------------------------------------------------------------*/
873 /*!
874 * \brief Return a pointer to the cs_property_t structure associated to the
875 * diffusion term for this equation (NULL if not activated).
876 *
877 * \param[in] eq pointer to a cs_equation_t structure
878 *
879 * \return a pointer to a cs_property_t structure
880 */
881 /*----------------------------------------------------------------------------*/
882
883 cs_property_t *
cs_equation_get_diffusion_property(const cs_equation_t * eq)884 cs_equation_get_diffusion_property(const cs_equation_t *eq)
885 {
886 if (eq == NULL)
887 return NULL;
888 else if (eq->param == NULL)
889 return NULL;
890 else
891 return eq->param->diffusion_property;
892 }
893
894 /*----------------------------------------------------------------------------*/
895 /*!
896 * \brief Return a pointer to the cs_property_t structure associated to the
897 * unsteady term for this equation (NULL if not activated).
898 *
899 * \param[in] eq pointer to a cs_equation_t structure
900 *
901 * \return a pointer to a cs_property_t structure
902 */
903 /*----------------------------------------------------------------------------*/
904
905 cs_property_t *
cs_equation_get_time_property(const cs_equation_t * eq)906 cs_equation_get_time_property(const cs_equation_t *eq)
907 {
908 if (eq == NULL)
909 return NULL;
910 else if (eq->param == NULL)
911 return NULL;
912 else
913 return eq->param->time_property;
914 }
915
916 /*----------------------------------------------------------------------------*/
917 /*!
918 * \brief Return a pointer to the cs_property_t structure associated to the
919 * reaction term with id equal to reaction_id and related to this
920 * equation
921 *
922 *
923 * \param[in] eq pointer to a cs_equation_t structure
924 * \param[in] reaction_id id related to this reaction term
925 *
926 * \return a pointer to a cs_property_t structure or NULL if not found
927 */
928 /*----------------------------------------------------------------------------*/
929
930 cs_property_t *
cs_equation_get_reaction_property(const cs_equation_t * eq,const int reaction_id)931 cs_equation_get_reaction_property(const cs_equation_t *eq,
932 const int reaction_id)
933 {
934 if (eq == NULL)
935 return NULL;
936
937 const cs_equation_param_t *eqp = eq->param;
938 if (reaction_id < 0 || reaction_id > eqp->n_reaction_terms - 1)
939 return NULL;
940
941 return eqp->reaction_properties[reaction_id];
942 }
943
944 /*----------------------------------------------------------------------------*/
945 /*!
946 * \brief Return the type of numerical scheme used for the discretization in
947 * time
948 *
949 * \param[in] eq pointer to a cs_equation_t structure
950 *
951 * \return a cs_param_time_scheme_t variable
952 */
953 /*----------------------------------------------------------------------------*/
954
955 cs_param_time_scheme_t
cs_equation_get_time_scheme(const cs_equation_t * eq)956 cs_equation_get_time_scheme(const cs_equation_t *eq)
957 {
958 if (eq == NULL)
959 return CS_TIME_N_SCHEMES;
960 else if (eq->param == NULL)
961 return CS_TIME_N_SCHEMES;
962 else
963 return eq->param->time_scheme;
964 }
965
966 /*----------------------------------------------------------------------------*/
967 /*!
968 * \brief Return the value of the theta parameter in theta time scheme
969 * discretization
970 *
971 * \param[in] eq pointer to a cs_equation_t structure
972 *
973 * \return the value of the theta coefficient. -1 if the time scheme is not
974 * a theta time scheme
975 */
976 /*----------------------------------------------------------------------------*/
977
978 cs_real_t
cs_equation_get_theta_time_val(const cs_equation_t * eq)979 cs_equation_get_theta_time_val(const cs_equation_t *eq)
980 {
981 cs_real_t theta = -1;
982
983 if (eq == NULL)
984 return theta;
985 else if (eq->param == NULL)
986 return theta;
987
988 else {
989
990 switch (eq->param->time_scheme) {
991
992 case CS_TIME_SCHEME_THETA:
993 theta = eq->param->theta;
994 break;
995 case CS_TIME_SCHEME_CRANKNICO:
996 theta = 0.5;
997 break;
998 case CS_TIME_SCHEME_BDF2:
999 case CS_TIME_SCHEME_EULER_IMPLICIT:
1000 theta = 1;
1001 break;
1002 case CS_TIME_SCHEME_EULER_EXPLICIT:
1003 theta = 0.;
1004 break;
1005
1006 default:
1007 break;
1008 }
1009
1010 }
1011
1012 return theta;
1013 }
1014
1015 /*----------------------------------------------------------------------------*/
1016 /*!
1017 * \brief Return the type of numerical scheme used for the discretization in
1018 * space
1019 *
1020 * \param[in] eq pointer to a cs_equation_t structure
1021 *
1022 * \return a cs_param_space_scheme_t variable
1023 */
1024 /*----------------------------------------------------------------------------*/
1025
1026 cs_param_space_scheme_t
cs_equation_get_space_scheme(const cs_equation_t * eq)1027 cs_equation_get_space_scheme(const cs_equation_t *eq)
1028 {
1029 if (eq == NULL)
1030 return CS_SPACE_N_SCHEMES;
1031 else if (eq->param == NULL)
1032 return CS_SPACE_N_SCHEMES;
1033 else
1034 return eq->param->space_scheme;
1035 }
1036
1037 /*----------------------------------------------------------------------------*/
1038 /*!
1039 * \brief Return the max. degree used in the polynomial basis for the space
1040 * discretization
1041 *
1042 * \param[in] eq pointer to a cs_equation_t structure
1043 *
1044 * \return the polynomial order
1045 */
1046 /*----------------------------------------------------------------------------*/
1047
1048 int
cs_equation_get_space_poly_degree(const cs_equation_t * eq)1049 cs_equation_get_space_poly_degree(const cs_equation_t *eq)
1050 {
1051 if (eq == NULL)
1052 return -1;
1053 else if (eq->param == NULL)
1054 return -1;
1055 else
1056 return eq->param->space_poly_degree;
1057 }
1058
1059 /*----------------------------------------------------------------------------*/
1060 /*!
1061 * \brief Return the dimension of the variable solved by this equation
1062 *
1063 * \param[in] eq pointer to a cs_equation_t structure
1064 *
1065 * \return an integer corresponding to the dimension of the variable
1066 */
1067 /*----------------------------------------------------------------------------*/
1068
1069 int
cs_equation_get_var_dim(const cs_equation_t * eq)1070 cs_equation_get_var_dim(const cs_equation_t *eq)
1071 {
1072 if (eq == NULL)
1073 return 0;
1074 else if (eq->param == NULL)
1075 return 0;
1076 else
1077 return eq->param->dim;
1078 }
1079
1080 /*----------------------------------------------------------------------------*/
1081 /*!
1082 * \brief Return the type of equation for the given equation structure
1083 *
1084 * \param[in] eq pointer to a cs_equation_t structure
1085 *
1086 * \return the type of the given equation
1087 */
1088 /*----------------------------------------------------------------------------*/
1089
1090 cs_equation_type_t
cs_equation_get_type(const cs_equation_t * eq)1091 cs_equation_get_type(const cs_equation_t *eq)
1092 {
1093 if (eq == NULL)
1094 return CS_EQUATION_N_TYPES;
1095 else if (eq->param == NULL)
1096 return CS_EQUATION_N_TYPES;
1097 else
1098 return eq->param->type;
1099 }
1100
1101 /*----------------------------------------------------------------------------*/
1102 /*!
1103 * \brief Return true is the given equation is steady otherwise false
1104 *
1105 * \param[in] eq pointer to a cs_equation_t structure
1106 *
1107 * \return true or false
1108 */
1109 /*----------------------------------------------------------------------------*/
1110
1111 bool
cs_equation_is_steady(const cs_equation_t * eq)1112 cs_equation_is_steady(const cs_equation_t *eq)
1113 {
1114 if (eq == NULL)
1115 return true;
1116 if (eq->param == NULL)
1117 return true;
1118
1119 if (eq->param->flag & CS_EQUATION_UNSTEADY)
1120 return false;
1121 else
1122 return true;
1123 }
1124
1125 /*----------------------------------------------------------------------------*/
1126 /*!
1127 * \brief Return true is the given equation is steady otherwise false
1128 *
1129 * \param[in] eq pointer to a cs_equation_t structure
1130 *
1131 * \return true or false
1132 */
1133 /*----------------------------------------------------------------------------*/
1134
1135 bool
cs_equation_uses_new_mechanism(const cs_equation_t * eq)1136 cs_equation_uses_new_mechanism(const cs_equation_t *eq)
1137 {
1138 if (eq == NULL)
1139 return false;
1140 assert(eq->param != NULL);
1141
1142 if (eq->param->dim == 1) {
1143 if ((eq->param->space_scheme == CS_SPACE_SCHEME_CDOVB) ||
1144 (eq->param->space_scheme == CS_SPACE_SCHEME_CDOVCB) ||
1145 (eq->param->space_scheme == CS_SPACE_SCHEME_CDOFB))
1146 return true;
1147 }
1148 else if (eq->param->dim == 3) {
1149 if ((eq->param->space_scheme == CS_SPACE_SCHEME_CDOVB) ||
1150 (eq->param->space_scheme == CS_SPACE_SCHEME_CDOFB) ||
1151 (eq->param->space_scheme == CS_SPACE_SCHEME_CDOEB))
1152 return true;
1153 }
1154
1155 return false;
1156 }
1157
1158 /*----------------------------------------------------------------------------*/
1159 /*!
1160 * \brief Add a new equation structure and set a first set of parameters
1161 *
1162 * \param[in] eqname name of the equation
1163 * \param[in] varname name of the variable associated to this equation
1164 * \param[in] eqtype type of equation (user, predefined...)
1165 * \param[in] dim dimension of the unknow attached to this equation
1166 * \param[in] default_bc type of boundary condition set by default
1167 *
1168 * \return a pointer to the new allocated cs_equation_t structure
1169 */
1170 /*----------------------------------------------------------------------------*/
1171
1172 cs_equation_t *
cs_equation_add(const char * eqname,const char * varname,cs_equation_type_t eqtype,int dim,cs_param_bc_type_t default_bc)1173 cs_equation_add(const char *eqname,
1174 const char *varname,
1175 cs_equation_type_t eqtype,
1176 int dim,
1177 cs_param_bc_type_t default_bc)
1178 {
1179 /* Sanity checks */
1180 if (varname == NULL)
1181 bft_error(__FILE__, __LINE__, 0,
1182 _(" %s: No variable name associated to an equation structure.\n"
1183 " Check your initialization."), __func__);
1184 if (eqname == NULL)
1185 bft_error(__FILE__, __LINE__, 0,
1186 _(" %s No equation name associated to an equation structure.\n"
1187 " Check your initialization."), __func__);
1188 if (cs_equation_by_name(eqname) != NULL)
1189 bft_error(__FILE__, __LINE__, 0,
1190 _(" %s: Stop adding a new equation.\n"
1191 " Equation name %s is already defined."), __func__, eqname);
1192
1193 cs_equation_t *eq = NULL;
1194
1195 BFT_MALLOC(eq, 1, cs_equation_t);
1196
1197 int eq_id = _n_equations;
1198 _n_equations++;
1199 BFT_REALLOC(_equations, _n_equations, cs_equation_t *);
1200 _equations[eq_id] = eq;
1201
1202 switch (eqtype) {
1203
1204 case CS_EQUATION_TYPE_USER:
1205 _n_user_equations++;
1206 break;
1207
1208 case CS_EQUATION_TYPE_GROUNDWATER:
1209 case CS_EQUATION_TYPE_MAXWELL:
1210 case CS_EQUATION_TYPE_NAVSTO:
1211 case CS_EQUATION_TYPE_PREDEFINED:
1212 case CS_EQUATION_TYPE_SOLIDIFICATION:
1213 case CS_EQUATION_TYPE_THERMAL:
1214 _n_predef_equations++;
1215 break;
1216
1217 default:
1218 bft_error(__FILE__, __LINE__, 0,
1219 " %s: This type of equation is not handled.\n"
1220 " Stop adding a new equation.", __func__);
1221 break;
1222
1223 }
1224
1225 eq->id = eq_id;
1226
1227 /* Store varname */
1228
1229 size_t len = strlen(varname);
1230 BFT_MALLOC(eq->varname, len + 1, char);
1231 strncpy(eq->varname, varname, len);
1232 eq->varname[len] = '\0';
1233
1234 eq->param = cs_equation_param_create(eqname, eqtype, dim, default_bc);
1235
1236 eq->field_id = -1; /* This field is created in a second step */
1237 eq->boundary_flux_id = -1; /* Not always defined (done in a second step) */
1238
1239 /* Algebraic system: allocated later */
1240
1241 eq->rset = NULL;
1242 eq->n_sles_gather_elts = eq->n_sles_scatter_elts = 0;
1243
1244 /* Builder structure for this equation */
1245
1246 eq->builder = NULL;
1247 eq->scheme_context = NULL;
1248
1249 /* Pointers of function */
1250
1251 eq->init_context = NULL;
1252 eq->free_context = NULL;
1253
1254 /* Extra-operations */
1255
1256 eq->compute_balance = NULL;
1257 eq->postprocess = NULL;
1258 eq->current_to_previous = NULL;
1259
1260 /* Restart */
1261
1262 eq->read_restart = NULL;
1263 eq->write_restart = NULL;
1264
1265 /* Function pointers to retrieve values at mesh locations */
1266
1267 eq->get_vertex_values = NULL;
1268 eq->get_edge_values = NULL;
1269 eq->get_face_values = NULL;
1270 eq->get_cell_values = NULL;
1271
1272 /* New functions */
1273
1274 eq->init_field_values = NULL;
1275 eq->solve = NULL;
1276 eq->solve_steady_state = NULL;
1277
1278 /* Deprecated functions */
1279
1280 eq->initialize_system = NULL;
1281 eq->set_dir_bc = NULL;
1282 eq->build_system = NULL;
1283 eq->update_field = NULL;
1284
1285 /* Deprecated members related to deprecated functions */
1286
1287 eq->matrix = NULL;
1288 eq->rhs = NULL;
1289
1290 /* Set timer statistic structure to a default value */
1291
1292 eq->main_ts_id = cs_timer_stats_id_by_name(eqname);
1293 if (eq->main_ts_id < 0)
1294 eq->main_ts_id = cs_timer_stats_create(NULL, /* new root */
1295 eqname,
1296 eqname);
1297
1298 return eq;
1299 }
1300
1301 /*----------------------------------------------------------------------------*/
1302 /*!
1303 * \brief Add a new user equation structure and set a first set of parameters
1304 *
1305 * \param[in] eqname name of the equation
1306 * \param[in] varname name of the variable associated to this equation
1307 * \param[in] dim dimension of the unknow attached to this equation
1308 * \param[in] default_bc type of boundary condition set by default
1309 *
1310 * \return a pointer to the new allocated cs_equation_t structure
1311 */
1312 /*----------------------------------------------------------------------------*/
1313
1314 cs_equation_t *
cs_equation_add_user(const char * eqname,const char * varname,int dim,cs_param_bc_type_t default_bc)1315 cs_equation_add_user(const char *eqname,
1316 const char *varname,
1317 int dim,
1318 cs_param_bc_type_t default_bc)
1319 {
1320 if (eqname == NULL)
1321 bft_error(__FILE__, __LINE__, 0, " %s: Empty equation name.", __func__);
1322 if (varname == NULL)
1323 bft_error(__FILE__, __LINE__, 0, " %s: Empty variable name.", __func__);
1324
1325 if ((default_bc != CS_PARAM_BC_HMG_DIRICHLET) &&
1326 (default_bc != CS_PARAM_BC_HMG_NEUMANN))
1327 bft_error(__FILE__, __LINE__, 0,
1328 _(" %s: Invalid type of boundary condition by default.\n"
1329 " Valid choices are CS_PARAM_BC_HMG_DIRICHLET or"
1330 " CS_PARAM_BC_HMG_NEUMANN"), __func__);
1331
1332 /* Add a new user equation */
1333 cs_equation_t *eq =
1334 cs_equation_add(eqname, /* equation name */
1335 varname, /* variable name */
1336 CS_EQUATION_TYPE_USER, /* type of equation */
1337 dim, /* dimension of the variable */
1338 default_bc); /* default BC */
1339
1340 return eq;
1341 }
1342
1343 /*----------------------------------------------------------------------------*/
1344 /*!
1345 * \brief Add a new user transport equation and set a first set of parameters
1346 * If time_pty is NULL, then no unsteady term is added.
1347 * If adv is NULL, then no advection term is added.
1348 * If diff_pty is NULL, then no diffusion term is added.
1349 *
1350 * \param[in] eqname name of the equation
1351 * \param[in] varname name of the variable associated to this equation
1352 * \param[in] dim dimension of the unknow attached to this equation
1353 * \param[in] default_bc type of boundary condition set by default
1354 * \param[in] time_pty property related to the unsteady term
1355 * \param[in] adv advection field
1356 * \param[in] diff_pty property related to the diffusion term
1357 *
1358 * \return a pointer to the new allocated cs_equation_t structure
1359 */
1360 /*----------------------------------------------------------------------------*/
1361
1362 cs_equation_t *
cs_equation_add_user_tracer(const char * eqname,const char * varname,int dim,cs_param_bc_type_t default_bc,cs_property_t * time_pty,cs_adv_field_t * adv,cs_property_t * diff_pty)1363 cs_equation_add_user_tracer(const char *eqname,
1364 const char *varname,
1365 int dim,
1366 cs_param_bc_type_t default_bc,
1367 cs_property_t *time_pty,
1368 cs_adv_field_t *adv,
1369 cs_property_t *diff_pty)
1370 {
1371 cs_equation_t *eq = cs_equation_add_user(eqname, varname, dim, default_bc);
1372 assert(eq != NULL);
1373
1374 /* Add an advection term */
1375 if (adv != NULL)
1376 cs_equation_add_advection(eq->param, adv);
1377
1378 /* Add the unsteady term */
1379 if (time_pty != NULL)
1380 cs_equation_add_time(eq->param, time_pty);
1381
1382 /* Add the diffusion term */
1383 if (diff_pty != NULL)
1384 cs_equation_add_diffusion(eq->param, diff_pty);
1385
1386 return eq;
1387 }
1388
1389 /*----------------------------------------------------------------------------*/
1390 /*!
1391 * \brief Destroy all cs_equation_t structures
1392 */
1393 /*----------------------------------------------------------------------------*/
1394
1395 void
cs_equation_destroy_all(void)1396 cs_equation_destroy_all(void)
1397 {
1398 if (_n_equations == 0)
1399 return;
1400
1401 for (int i = 0; i < _n_equations; i++) {
1402
1403 cs_equation_t *eq = _equations[i];
1404
1405 if (eq->main_ts_id > -1)
1406 cs_timer_stats_start(eq->main_ts_id);
1407
1408 eq->param = cs_equation_free_param(eq->param);
1409
1410 /* Sanity check */
1411 assert(eq->matrix == NULL && eq->rhs == NULL);
1412 /* Since eq->rset is only shared, no free is done at this stage */
1413
1414 /* Free the associated builder structure */
1415 cs_equation_builder_free(&(eq->builder));
1416 eq->scheme_context = eq->free_context(eq->scheme_context);
1417
1418 if (eq->main_ts_id > -1)
1419 cs_timer_stats_stop(eq->main_ts_id);
1420
1421 BFT_FREE(eq->varname);
1422 BFT_FREE(eq);
1423
1424 } /* Loop on equations */
1425
1426 BFT_FREE(_equations);
1427
1428 _n_equations = 0;
1429 _n_user_equations = 0;
1430 _n_predef_equations = 0;
1431 }
1432
1433 /*----------------------------------------------------------------------------*/
1434 /*!
1435 * \brief Check if a steady-state computation is requested according to the
1436 * setting
1437 *
1438 * \return true or false
1439 */
1440 /*----------------------------------------------------------------------------*/
1441
1442 bool
cs_equation_needs_steady_state_solve(void)1443 cs_equation_needs_steady_state_solve(void)
1444 {
1445 for (int eq_id = 0; eq_id < _n_equations; eq_id++) {
1446
1447 cs_equation_t *eq = _equations[eq_id];
1448
1449 if (cs_equation_is_steady(eq))
1450 return true;
1451
1452 } /* Loop on equations */
1453
1454 return false;
1455 }
1456
1457 /*----------------------------------------------------------------------------*/
1458 /*!
1459 * \brief Print a synthesis of the monitoring information in the performance
1460 * file
1461 */
1462 /*----------------------------------------------------------------------------*/
1463
1464 void
cs_equation_log_monitoring(void)1465 cs_equation_log_monitoring(void)
1466 {
1467 cs_log_printf(CS_LOG_PERFORMANCE, "%-36s %9s %9s %9s\n",
1468 " ", "Build", "Solve", "Extra");
1469
1470 for (int i = 0; i < _n_equations; i++) {
1471
1472 cs_equation_t *eq = _equations[i];
1473
1474 /* Display high-level timer counter related to the current equation
1475 before deleting the structure */
1476 cs_equation_write_monitoring(eq->param->name, eq->builder);
1477
1478 } /* Loop on equations */
1479 }
1480
1481 /*----------------------------------------------------------------------------*/
1482 /*!
1483 * \brief Get the count of equations of each macro type
1484 *
1485 * \param[out] n_equations total number of equations
1486 * \param[out] n_predef_equations number of predefined equations
1487 * \param[out] n_user_equations number of user equations
1488 */
1489 /*----------------------------------------------------------------------------*/
1490
1491 void
cs_equation_get_count(int * n_equations,int * n_predef_equations,int * n_user_equations)1492 cs_equation_get_count(int *n_equations,
1493 int *n_predef_equations,
1494 int *n_user_equations)
1495 {
1496 *n_equations = _n_equations;
1497 *n_predef_equations = _n_predef_equations;
1498 *n_user_equations = _n_user_equations;
1499 }
1500
1501 /*----------------------------------------------------------------------------*/
1502 /*!
1503 * \brief Summarize all cs_equation_t structures
1504 */
1505 /*----------------------------------------------------------------------------*/
1506
1507 void
cs_equation_log_setup(void)1508 cs_equation_log_setup(void)
1509 {
1510 cs_log_printf(CS_LOG_SETUP, "\nSettings for equations\n");
1511 cs_log_printf(CS_LOG_SETUP, "%s\n", cs_sep_h1);
1512
1513 for (int eq_id = 0; eq_id < _n_equations; eq_id++) {
1514
1515 cs_equation_t *eq = _equations[eq_id];
1516 assert(eq->param != NULL);
1517
1518 if (eq->main_ts_id > -1)
1519 cs_timer_stats_start(eq->main_ts_id);
1520
1521 cs_log_printf(CS_LOG_SETUP,
1522 "Summary of settings for %s eq. (variable %s)\n",
1523 eq->param->name, eq->varname);
1524 cs_log_printf(CS_LOG_SETUP, "%s", cs_sep_h2);
1525
1526 cs_equation_param_log(eq->param);
1527
1528 if (eq->main_ts_id > -1)
1529 cs_timer_stats_stop(eq->main_ts_id);
1530
1531 } /* Loop on equations */
1532
1533 }
1534
1535 /*----------------------------------------------------------------------------*/
1536 /*!
1537 * \brief Set a parameter attached to a keyname for the default settigns
1538 *
1539 * \param[in] key key related to the member of eq to set
1540 * \param[in] keyval accessor to the value to set
1541 */
1542 /*----------------------------------------------------------------------------*/
1543
1544 void
cs_equation_set_default_param(cs_equation_key_t key,const char * keyval)1545 cs_equation_set_default_param(cs_equation_key_t key,
1546 const char *keyval)
1547 {
1548 if (_n_equations == 0)
1549 return;
1550
1551 for (int eq_id = 0; eq_id < _n_equations; eq_id++) {
1552
1553 cs_equation_t *eq = _equations[eq_id];
1554 if (eq == NULL)
1555 continue;
1556
1557 if (eq->main_ts_id > -1)
1558 cs_timer_stats_start(eq->main_ts_id);
1559
1560 cs_equation_param_set(eq->param, key, keyval);
1561
1562 if (eq->main_ts_id > -1)
1563 cs_timer_stats_stop(eq->main_ts_id);
1564
1565 } /* Loop on equations */
1566 }
1567
1568 /*----------------------------------------------------------------------------*/
1569 /*!
1570 * \brief Setup the linear algebra requirements
1571 */
1572 /*----------------------------------------------------------------------------*/
1573
1574 void
cs_equation_set_sles(void)1575 cs_equation_set_sles(void)
1576 {
1577 if (_n_equations == 0)
1578 return;
1579
1580 for (int eq_id = 0; eq_id < _n_equations; eq_id++) {
1581
1582 cs_equation_t *eq = _equations[eq_id];
1583 cs_equation_param_t *eqp = eq->param;
1584
1585 if (eq->main_ts_id > -1)
1586 cs_timer_stats_start(eq->main_ts_id);
1587
1588 /* Initialize cs_sles_t structure */
1589 if (eqp->type != CS_EQUATION_TYPE_NAVSTO)
1590 cs_equation_param_set_sles(eqp);
1591
1592 if (eq->main_ts_id > -1)
1593 cs_timer_stats_stop(eq->main_ts_id);
1594
1595 } /* Loop on equations */
1596 }
1597
1598 /*----------------------------------------------------------------------------*/
1599 /*!
1600 * \brief Set shared structures among the activated class of discretization
1601 * schemes
1602 *
1603 * \param[in] connect pointer to a cs_cdo_connect_t structure
1604 * \param[in] quant pointer to additional mesh quantities struct.
1605 * \param[in] time_step pointer to a time step structure
1606 * \param[in] eb_scheme_flag metadata for Eb schemes
1607 * \param[in] fb_scheme_flag metadata for Fb schemes
1608 * \param[in] vb_scheme_flag metadata for Vb schemes
1609 * \param[in] vcb_scheme_flag metadata for V+C schemes
1610 * \param[in] hho_scheme_flag metadata for HHO schemes
1611 */
1612 /*----------------------------------------------------------------------------*/
1613
1614 void
cs_equation_set_shared_structures(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step,cs_flag_t eb_scheme_flag,cs_flag_t fb_scheme_flag,cs_flag_t vb_scheme_flag,cs_flag_t vcb_scheme_flag,cs_flag_t hho_scheme_flag)1615 cs_equation_set_shared_structures(const cs_cdo_connect_t *connect,
1616 const cs_cdo_quantities_t *quant,
1617 const cs_time_step_t *time_step,
1618 cs_flag_t eb_scheme_flag,
1619 cs_flag_t fb_scheme_flag,
1620 cs_flag_t vb_scheme_flag,
1621 cs_flag_t vcb_scheme_flag,
1622 cs_flag_t hho_scheme_flag)
1623 {
1624 if (vb_scheme_flag > 0 || vcb_scheme_flag > 0) {
1625
1626 if (vb_scheme_flag & CS_FLAG_SCHEME_SCALAR ||
1627 vcb_scheme_flag & CS_FLAG_SCHEME_SCALAR) {
1628
1629 cs_matrix_structure_t *ms
1630 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_VTX_SCAL);
1631
1632 if (vb_scheme_flag & CS_FLAG_SCHEME_SCALAR)
1633 cs_cdovb_scaleq_init_common(quant, connect, time_step, ms);
1634
1635 if (vcb_scheme_flag & CS_FLAG_SCHEME_SCALAR)
1636 cs_cdovcb_scaleq_init_common(quant, connect, time_step, ms);
1637
1638 } /* scalar-valued DoFs */
1639
1640 if (vb_scheme_flag & CS_FLAG_SCHEME_VECTOR ||
1641 vcb_scheme_flag & CS_FLAG_SCHEME_VECTOR) {
1642
1643 cs_matrix_structure_t *ms
1644 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_VTX_VECT);
1645
1646 if (vb_scheme_flag & CS_FLAG_SCHEME_VECTOR)
1647 cs_cdovb_vecteq_init_common(quant, connect, time_step, ms);
1648
1649 /* TODO: VCb schemes DoFs */
1650
1651 } /* vector-valued DoFs */
1652
1653 } /* Vertex-based class of discretization schemes */
1654
1655 if (eb_scheme_flag > 0) {
1656
1657 if (eb_scheme_flag & CS_FLAG_SCHEME_SCALAR) {
1658 /* This is a vector-valued equation but the DoF is scalar-valued since
1659 * it is a circulation associated to each edge */
1660
1661 cs_matrix_structure_t *ms
1662 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_EDGE_SCAL);
1663
1664 cs_cdoeb_vecteq_init_common(quant, connect, time_step, ms);
1665
1666 } /* scalar-valued DoFs */
1667
1668 } /* Edge-based class of discretization schemes */
1669
1670 if (fb_scheme_flag > 0 || hho_scheme_flag > 0) {
1671
1672 if (cs_flag_test(fb_scheme_flag,
1673 CS_FLAG_SCHEME_POLY0 | CS_FLAG_SCHEME_SCALAR)) {
1674
1675 cs_matrix_structure_t *ms
1676 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_FACE_SP0);
1677
1678 cs_cdofb_scaleq_init_common(quant, connect, time_step, ms);
1679
1680 } /* Scalar-valued CDO-Fb DoFs */
1681
1682 if (cs_flag_test(fb_scheme_flag,
1683 CS_FLAG_SCHEME_POLY0 | CS_FLAG_SCHEME_VECTOR)) {
1684
1685 cs_matrix_structure_t *ms
1686 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_FACE_SP1);
1687
1688 cs_cdofb_vecteq_init_common(quant, connect, time_step, ms);
1689
1690 } /* Vector-valued CDO-Fb DoFs */
1691
1692 if (hho_scheme_flag & CS_FLAG_SCHEME_SCALAR) {
1693
1694 cs_matrix_structure_t *ms0
1695 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_FACE_SP0);
1696 cs_matrix_structure_t *ms1
1697 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_FACE_SP1);
1698 cs_matrix_structure_t *ms2
1699 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_FACE_SP2);
1700
1701 cs_hho_scaleq_init_common(hho_scheme_flag,
1702 quant, connect, time_step,
1703 ms0, ms1, ms2);
1704
1705 } /* Scalar-valued HHO schemes DoFs */
1706
1707 if (hho_scheme_flag & CS_FLAG_SCHEME_VECTOR) {
1708
1709 cs_matrix_structure_t *ms0
1710 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_FACE_VHP0);
1711 cs_matrix_structure_t *ms1
1712 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_FACE_VHP1);
1713 cs_matrix_structure_t *ms2
1714 = cs_equation_get_matrix_structure(CS_CDO_CONNECT_FACE_VHP2);
1715
1716 cs_hho_vecteq_init_common(hho_scheme_flag,
1717 quant, connect, time_step,
1718 ms0, ms1, ms2);
1719
1720 } /* Vector-valued HHO schemes DoFs */
1721
1722 } /* Face-based class of discretization schemes */
1723
1724 }
1725
1726 /*----------------------------------------------------------------------------*/
1727 /*!
1728 * \brief Release shared structures among the activated class of discretization
1729 * schemes
1730 *
1731 * \param[in] vb_scheme_flag metadata for Vb schemes
1732 * \param[in] vcb_scheme_flag metadata for V+C schemes
1733 * \param[in] eb_scheme_flag metadata for Eb schemes
1734 * \param[in] fb_scheme_flag metadata for Fb schemes
1735 * \param[in] hho_scheme_flag metadata for HHO schemes
1736 */
1737 /*----------------------------------------------------------------------------*/
1738
1739 void
cs_equation_unset_shared_structures(cs_flag_t vb_scheme_flag,cs_flag_t vcb_scheme_flag,cs_flag_t eb_scheme_flag,cs_flag_t fb_scheme_flag,cs_flag_t hho_scheme_flag)1740 cs_equation_unset_shared_structures(cs_flag_t vb_scheme_flag,
1741 cs_flag_t vcb_scheme_flag,
1742 cs_flag_t eb_scheme_flag,
1743 cs_flag_t fb_scheme_flag,
1744 cs_flag_t hho_scheme_flag)
1745 {
1746 /* Free common structures specific to a numerical scheme */
1747 if (vb_scheme_flag & CS_FLAG_SCHEME_SCALAR)
1748 cs_cdovb_scaleq_finalize_common();
1749
1750 if (vb_scheme_flag & CS_FLAG_SCHEME_VECTOR)
1751 cs_cdovb_vecteq_finalize_common();
1752
1753 if (vcb_scheme_flag & CS_FLAG_SCHEME_SCALAR)
1754 cs_cdovcb_scaleq_finalize_common();
1755
1756 if (eb_scheme_flag & CS_FLAG_SCHEME_SCALAR)
1757 cs_cdoeb_vecteq_finalize_common();
1758
1759 if (fb_scheme_flag & CS_FLAG_SCHEME_SCALAR)
1760 cs_cdofb_scaleq_finalize_common();
1761
1762 if (fb_scheme_flag & CS_FLAG_SCHEME_VECTOR)
1763 cs_cdofb_vecteq_finalize_common();
1764
1765 if (hho_scheme_flag & CS_FLAG_SCHEME_SCALAR)
1766 cs_hho_scaleq_finalize_common();
1767
1768 if (hho_scheme_flag & CS_FLAG_SCHEME_VECTOR)
1769 cs_hho_vecteq_finalize_common();
1770 }
1771
1772 /*----------------------------------------------------------------------------*/
1773 /*!
1774 * \brief Assign a \ref cs_range_set_t structures for synchronization when
1775 * computing in parallel mode.
1776 *
1777 * \param[in] connect pointer to a cs_cdo_connect_t structure
1778 */
1779 /*----------------------------------------------------------------------------*/
1780
1781 void
cs_equation_set_range_set(const cs_cdo_connect_t * connect)1782 cs_equation_set_range_set(const cs_cdo_connect_t *connect)
1783
1784 {
1785 if (_n_equations == 0)
1786 return;
1787
1788 const char s_err_msg[] =
1789 "%s: Only the scalar-valued case is handled for this scheme.\n";
1790 const char sv_err_msg[] =
1791 "%s: Only the scalar-valued and vector-valued case are handled"
1792 "for this scheme.\n";
1793
1794 const cs_lnum_t n_faces = connect->n_faces[CS_ALL_FACES];
1795
1796 for (int eq_id = 0; eq_id < _n_equations; eq_id++) {
1797
1798 cs_equation_t *eq = _equations[eq_id];
1799 cs_equation_param_t *eqp = eq->param;
1800
1801 if (eq->main_ts_id > -1)
1802 cs_timer_stats_start(eq->main_ts_id);
1803
1804 /* Set function pointers */
1805 switch(eqp->space_scheme) {
1806
1807 case CS_SPACE_SCHEME_CDOVB:
1808 if (eqp->dim == 1) {
1809
1810 /* Set the cs_range_set_t structure */
1811 eq->rset = connect->range_sets[CS_CDO_CONNECT_VTX_SCAL];
1812
1813 /* Set the size of the algebraic system arising from the cellwise
1814 process */
1815 eq->n_sles_gather_elts = eq->n_sles_scatter_elts = connect->n_vertices;
1816
1817 }
1818 else if (eqp->dim == 3) {
1819
1820 /* Set the cs_range_set_t structure */
1821 eq->rset = connect->range_sets[CS_CDO_CONNECT_VTX_VECT];
1822
1823 /* Set the size of the algebraic system arising from the cellwise
1824 process */
1825 eq->n_sles_gather_elts = eqp->dim * connect->n_vertices;
1826 eq->n_sles_scatter_elts = eqp->dim * connect->n_vertices;
1827
1828 }
1829 else
1830 bft_error(__FILE__, __LINE__, 0, sv_err_msg, __func__);
1831
1832 break;
1833
1834 case CS_SPACE_SCHEME_CDOVCB:
1835 if (eqp->dim == 1) {
1836
1837 /* Set the cs_range_set_t structure */
1838 eq->rset = connect->range_sets[CS_CDO_CONNECT_VTX_SCAL];
1839
1840 /* Set the size of the algebraic system arising from the cellwise
1841 process */
1842 eq->n_sles_gather_elts = eq->n_sles_scatter_elts = connect->n_vertices;
1843
1844 }
1845 else
1846 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
1847
1848 break;
1849
1850 case CS_SPACE_SCHEME_CDOEB:
1851 if (eqp->dim == 3) {
1852 /* This is a vector-valued equation but the DoF is scalar-valued since
1853 * it is a circulation associated to each edge */
1854
1855 /* Set the cs_range_set_t structure */
1856 eq->rset = connect->range_sets[CS_CDO_CONNECT_EDGE_SCAL];
1857
1858 /* Set the size of the algebraic system arising from the cellwise
1859 process */
1860 eq->n_sles_gather_elts = eq->n_sles_scatter_elts = connect->n_vertices;
1861
1862 }
1863 else
1864 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
1865
1866 break;
1867
1868 case CS_SPACE_SCHEME_CDOFB:
1869 if (eqp->dim == 1) {
1870
1871 /* Set the cs_range_set_t structure */
1872 eq->rset = connect->range_sets[CS_CDO_CONNECT_FACE_SP0];
1873
1874 /* Set the size of the algebraic system arising from the cellwise
1875 process */
1876 eq->n_sles_gather_elts = eq->n_sles_scatter_elts = n_faces;
1877
1878 }
1879 else if (eqp->dim == 3) {
1880
1881 /* Set the cs_range_set_t structure */
1882 eq->rset = connect->range_sets[CS_CDO_CONNECT_FACE_VP0];
1883
1884 /* Set the size of the algebraic system arising from the cellwise
1885 process (OK for a sequential run) */
1886 eq->n_sles_gather_elts = eqp->dim * n_faces;
1887 eq->n_sles_scatter_elts = eqp->dim * n_faces;
1888
1889 }
1890 else
1891 bft_error(__FILE__, __LINE__, 0, sv_err_msg, __func__);
1892
1893 break;
1894
1895 case CS_SPACE_SCHEME_HHO_P0:
1896 if (eqp->dim == 1) {
1897
1898 /* Set the cs_range_set_t structure */
1899 eq->rset = connect->range_sets[CS_CDO_CONNECT_FACE_SP0];
1900
1901 /* Set the size of the algebraic system arising from the cellwise
1902 process */
1903 eq->n_sles_gather_elts = eq->n_sles_scatter_elts = n_faces;
1904
1905 }
1906 else
1907 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
1908
1909 break;
1910
1911 case CS_SPACE_SCHEME_HHO_P1:
1912 if (eqp->dim == 1) {
1913
1914 /* Set the cs_range_set_t structure */
1915 eq->rset = connect->range_sets[CS_CDO_CONNECT_FACE_SP1];
1916
1917 /* Set the size of the algebraic system arising from the cellwise
1918 process (OK for a sequential run) */
1919 eq->n_sles_gather_elts = CS_N_FACE_DOFS_1ST * n_faces;
1920 eq->n_sles_scatter_elts = CS_N_FACE_DOFS_1ST * n_faces;
1921
1922 }
1923 else if (eqp->dim == 3) {
1924
1925 /* Set the cs_range_set_t structure */
1926 eq->rset = connect->range_sets[CS_CDO_CONNECT_FACE_VHP1];
1927
1928 /* Set the size of the algebraic system arising from the cellwise
1929 process (OK for a sequential run) */
1930 eq->n_sles_gather_elts = 3*CS_N_FACE_DOFS_1ST * n_faces;
1931 eq->n_sles_scatter_elts = 3*CS_N_FACE_DOFS_1ST * n_faces;
1932
1933 }
1934 else
1935 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
1936
1937 break;
1938
1939 case CS_SPACE_SCHEME_HHO_P2:
1940 if (eqp->dim == 1) {
1941
1942 /* Set the cs_range_set_t structure */
1943 eq->rset = connect->range_sets[CS_CDO_CONNECT_FACE_SP2];
1944
1945 /* Set the size of the algebraic system arising from the cellwise
1946 process (OK for a sequential run) */
1947 eq->n_sles_gather_elts = CS_N_FACE_DOFS_2ND * n_faces;
1948 eq->n_sles_scatter_elts = CS_N_FACE_DOFS_2ND * n_faces;
1949
1950 }
1951 else if (eqp->dim == 3) {
1952
1953 /* Set the cs_range_set_t structure */
1954 eq->rset = connect->range_sets[CS_CDO_CONNECT_FACE_VHP2];
1955
1956 /* Set the size of the algebraic system arising from the cellwise
1957 process (OK for a sequential run) */
1958 eq->n_sles_gather_elts = 3*CS_N_FACE_DOFS_2ND * n_faces;
1959 eq->n_sles_scatter_elts = 3*CS_N_FACE_DOFS_2ND * n_faces;
1960
1961 }
1962 else
1963 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
1964
1965 break;
1966
1967 default:
1968 bft_error(__FILE__, __LINE__, 0,
1969 _(" %s: Invalid scheme for the space discretization.\n"
1970 " Please check your settings."), __func__);
1971 break;
1972 }
1973
1974 if (cs_glob_n_ranks > 1)
1975 eq->n_sles_gather_elts = eq->rset->n_elts[0];
1976
1977 if (eq->main_ts_id > -1)
1978 cs_timer_stats_stop(eq->main_ts_id);
1979
1980 } /* Loop on equations */
1981 }
1982
1983 /*----------------------------------------------------------------------------*/
1984 /*!
1985 * \brief Assign a set of pointer functions for managing the cs_equation_t
1986 * structure during the computation
1987 * After this call, parameters related to an equation are set once for
1988 * all
1989 *
1990 * \return true if all equations are steady-state otherwise false
1991 */
1992 /*----------------------------------------------------------------------------*/
1993
1994 bool
cs_equation_set_functions(void)1995 cs_equation_set_functions(void)
1996 {
1997 if (_n_equations == 0)
1998 return true;
1999
2000 const char s_err_msg[] =
2001 "%s: Only the scalar-valued case is handled for this scheme.\n";
2002 const char sv_err_msg[] =
2003 "%s: Only the scalar-valued and vector-valued case are handled"
2004 "for this scheme.\n";
2005
2006 bool all_are_steady = true;
2007
2008 for (int eq_id = 0; eq_id < _n_equations; eq_id++) {
2009
2010 cs_equation_t *eq = _equations[eq_id];
2011 cs_equation_param_t *eqp = eq->param;
2012
2013 if (eq->main_ts_id > -1)
2014 cs_timer_stats_start(eq->main_ts_id);
2015
2016 if (eqp->flag & CS_EQUATION_UNSTEADY)
2017 all_are_steady = false;
2018 else
2019 cs_equation_param_set(eqp, CS_EQKEY_TIME_SCHEME, "steady");
2020
2021 /* Apply the last modifications to the cs_equation_param_t structure */
2022 cs_equation_param_last_stage(eqp);
2023
2024 /* Set function pointers */
2025 switch(eqp->space_scheme) {
2026
2027 case CS_SPACE_SCHEME_CDOVB:
2028 if (eqp->dim == 1) {
2029
2030 eq->init_context = cs_cdovb_scaleq_init_context;
2031 eq->free_context = cs_cdovb_scaleq_free_context;
2032 eq->init_field_values = cs_cdovb_scaleq_init_values;
2033
2034 /* deprecated */
2035 eq->set_dir_bc = NULL;
2036 eq->initialize_system = NULL;
2037 eq->build_system = NULL;
2038 eq->prepare_solving = NULL;
2039 eq->update_field = NULL;
2040
2041 /* New mechanism */
2042 eq->solve_steady_state = cs_cdovb_scaleq_solve_steady_state;
2043 switch (eqp->time_scheme) {
2044 case CS_TIME_SCHEME_STEADY:
2045 eq->solve = eq->solve_steady_state;
2046 break;
2047 case CS_TIME_SCHEME_EULER_IMPLICIT:
2048 eq->solve = cs_cdovb_scaleq_solve_implicit;
2049 break;
2050 case CS_TIME_SCHEME_THETA:
2051 case CS_TIME_SCHEME_CRANKNICO:
2052 eq->solve = cs_cdovb_scaleq_solve_theta;
2053 break;
2054
2055 case CS_TIME_SCHEME_BDF2:
2056 default:
2057 bft_error(__FILE__, __LINE__, 0,
2058 "%s: Eq. %s. This time scheme is not yet implemented",
2059 __func__, eqp->name);
2060 }
2061
2062 eq->compute_balance = cs_cdovb_scaleq_balance;
2063 eq->postprocess = cs_cdovb_scaleq_extra_post;
2064 eq->current_to_previous = cs_cdovb_scaleq_current_to_previous;
2065
2066 eq->read_restart = NULL;
2067 eq->write_restart = NULL;
2068
2069 /* Function pointers to retrieve values at mesh locations */
2070 eq->get_vertex_values = cs_cdovb_scaleq_get_vertex_values;
2071 eq->get_edge_values = NULL;
2072 eq->get_face_values = NULL;
2073 eq->get_cell_values = cs_cdovb_scaleq_get_cell_values;
2074
2075 eq->get_cw_build_structures = cs_cdovb_scaleq_get;
2076
2077 }
2078 else if (eqp->dim == 3) {
2079
2080 eq->init_context = cs_cdovb_vecteq_init_context;
2081 eq->free_context = cs_cdovb_vecteq_free_context;
2082 eq->init_field_values = cs_cdovb_vecteq_init_values;
2083
2084 /* Deprecated */
2085 eq->set_dir_bc = NULL;
2086 eq->initialize_system = NULL;
2087 eq->build_system = NULL;
2088 eq->prepare_solving = NULL;
2089 eq->update_field = NULL;
2090
2091 /* New mechanism */
2092 eq->solve_steady_state = cs_cdovb_vecteq_solve_steady_state;
2093 switch (eqp->time_scheme) {
2094 case CS_TIME_SCHEME_STEADY:
2095 eq->solve = eq->solve_steady_state;
2096 break;
2097 case CS_TIME_SCHEME_BDF2:
2098 case CS_TIME_SCHEME_EULER_IMPLICIT:
2099 case CS_TIME_SCHEME_THETA:
2100 case CS_TIME_SCHEME_CRANKNICO:
2101 default:
2102 bft_error(__FILE__, __LINE__, 0,
2103 "%s: Eq. %s. This time scheme is not yet implemented",
2104 __func__, eqp->name);
2105 }
2106
2107 eq->postprocess = cs_cdovb_vecteq_extra_post;
2108 eq->current_to_previous = cs_cdovb_vecteq_current_to_previous;
2109
2110 eq->read_restart = NULL;
2111 eq->write_restart = NULL;
2112
2113 /* Function pointers to retrieve values at mesh locations */
2114 eq->get_vertex_values = cs_cdovb_vecteq_get_vertex_values;
2115 eq->get_edge_values = NULL;
2116 eq->get_face_values = NULL;
2117 eq->get_cell_values = cs_cdovb_vecteq_get_cell_values;
2118
2119 eq->get_cw_build_structures = cs_cdovb_vecteq_get;
2120
2121 }
2122 else
2123 bft_error(__FILE__, __LINE__, 0, sv_err_msg, __func__);
2124
2125 break;
2126
2127 case CS_SPACE_SCHEME_CDOVCB:
2128 if (eqp->dim == 1) {
2129
2130 eq->init_context = cs_cdovcb_scaleq_init_context;
2131 eq->free_context = cs_cdovcb_scaleq_free_context;
2132 eq->init_field_values = cs_cdovcb_scaleq_init_values;
2133
2134 /* Deprecated */
2135 eq->set_dir_bc = NULL;
2136 eq->initialize_system = NULL;
2137 eq->build_system = NULL;
2138 eq->prepare_solving = NULL;
2139 eq->update_field = NULL;
2140
2141 /* New mechanism */
2142 eq->solve_steady_state = cs_cdovcb_scaleq_solve_steady_state;
2143 switch (eqp->time_scheme) {
2144 case CS_TIME_SCHEME_STEADY:
2145 eq->solve = eq->solve_steady_state;
2146 break;
2147
2148 case CS_TIME_SCHEME_EULER_IMPLICIT:
2149 eq->solve = cs_cdovcb_scaleq_solve_implicit;
2150 break;
2151
2152 case CS_TIME_SCHEME_THETA:
2153 case CS_TIME_SCHEME_CRANKNICO:
2154 eq->solve = cs_cdovcb_scaleq_solve_theta;
2155 break;
2156
2157 case CS_TIME_SCHEME_BDF2:
2158 default:
2159 bft_error(__FILE__, __LINE__, 0,
2160 "%s: Eq. %s. This time scheme is not yet implemented",
2161 __func__, eqp->name);
2162 }
2163
2164 eq->postprocess = cs_cdovcb_scaleq_extra_post;
2165 eq->current_to_previous = cs_cdovcb_scaleq_current_to_previous;
2166
2167 eq->read_restart = cs_cdovcb_scaleq_read_restart;
2168 eq->write_restart = cs_cdovcb_scaleq_write_restart;
2169
2170 /* Function pointers to retrieve values at mesh locations */
2171 eq->get_vertex_values = cs_cdovcb_scaleq_get_vertex_values;
2172 eq->get_edge_values = NULL;
2173 eq->get_face_values = NULL;
2174 eq->get_cell_values = cs_cdovcb_scaleq_get_cell_values;
2175
2176 eq->get_cw_build_structures = cs_cdovcb_scaleq_get;
2177
2178 }
2179 else
2180 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
2181
2182 break;
2183
2184 case CS_SPACE_SCHEME_CDOFB:
2185 if (eqp->dim == 1) {
2186
2187 eq->init_context = cs_cdofb_scaleq_init_context;
2188 eq->free_context = cs_cdofb_scaleq_free_context;
2189 eq->init_field_values = cs_cdofb_scaleq_init_values;
2190
2191 /* Deprecated */
2192 eq->set_dir_bc = NULL;
2193 eq->initialize_system = NULL;
2194 eq->build_system = NULL;
2195 eq->prepare_solving = NULL;
2196 eq->update_field = NULL;
2197
2198 /* New mechanism */
2199 eq->solve_steady_state = cs_cdofb_scaleq_solve_steady_state;
2200 switch (eqp->time_scheme) {
2201 case CS_TIME_SCHEME_STEADY:
2202 eq->solve = eq->solve_steady_state;
2203 break;
2204
2205 case CS_TIME_SCHEME_EULER_IMPLICIT:
2206 eq->solve = cs_cdofb_scaleq_solve_implicit;
2207 break;
2208
2209 case CS_TIME_SCHEME_THETA:
2210 case CS_TIME_SCHEME_CRANKNICO:
2211 eq->solve = cs_cdofb_scaleq_solve_theta;
2212 break;
2213
2214 case CS_TIME_SCHEME_BDF2:
2215 default:
2216 bft_error(__FILE__, __LINE__, 0,
2217 "%s: Eq. %s. This time scheme is not yet implemented",
2218 __func__, eqp->name);
2219 }
2220
2221 eq->compute_balance = cs_cdofb_scaleq_balance;
2222 eq->postprocess = cs_cdofb_scaleq_extra_post;
2223 eq->current_to_previous = cs_cdofb_scaleq_current_to_previous;
2224
2225 eq->read_restart = cs_cdofb_scaleq_read_restart;
2226 eq->write_restart = cs_cdofb_scaleq_write_restart;
2227
2228 /* Function pointers to retrieve values at mesh locations */
2229 eq->get_vertex_values = NULL;
2230 eq->get_edge_values = NULL;
2231 eq->get_face_values = cs_cdofb_scaleq_get_face_values;
2232 eq->get_cell_values = cs_cdofb_scaleq_get_cell_values;
2233
2234 eq->get_cw_build_structures = cs_cdofb_scaleq_get;
2235
2236 }
2237 else if (eqp->dim == 3) {
2238
2239 eq->init_context = cs_cdofb_vecteq_init_context;
2240 eq->free_context = cs_cdofb_vecteq_free_context;
2241 eq->init_field_values = cs_cdofb_vecteq_init_values;
2242
2243 /* Deprecated */
2244 eq->initialize_system = NULL;
2245 eq->set_dir_bc = NULL;
2246 eq->build_system = NULL;
2247 eq->prepare_solving = NULL;
2248 eq->update_field = NULL;
2249
2250 /* New mechanism */
2251 eq->solve_steady_state = cs_cdofb_vecteq_solve_steady_state;
2252 switch (eqp->time_scheme) {
2253 case CS_TIME_SCHEME_STEADY:
2254 eq->solve = eq->solve_steady_state;
2255 break;
2256
2257 case CS_TIME_SCHEME_EULER_IMPLICIT:
2258 eq->solve = cs_cdofb_vecteq_solve_implicit;
2259 break;
2260
2261 case CS_TIME_SCHEME_THETA:
2262 case CS_TIME_SCHEME_CRANKNICO:
2263 eq->solve = cs_cdofb_vecteq_solve_theta;
2264 break;
2265
2266 case CS_TIME_SCHEME_BDF2:
2267 eq->solve = NULL; /* cs_cdofb_vecteq_solve_bdf2 */
2268 bft_error(__FILE__, __LINE__, 0,
2269 "%s: Eq. %s. This time scheme is not yet implemented",
2270 __func__, eqp->name);
2271 break;
2272
2273 default:
2274 bft_error(__FILE__, __LINE__, 0,
2275 "%s: Eq. %s. This time scheme is not yet implemented",
2276 __func__, eqp->name);
2277 }
2278
2279 eq->postprocess = cs_cdofb_vecteq_extra_post;
2280 eq->current_to_previous = cs_cdofb_vecteq_current_to_previous;
2281
2282 eq->read_restart = cs_cdofb_vecteq_read_restart;
2283 eq->write_restart = cs_cdofb_vecteq_write_restart;
2284
2285 /* Function pointers to retrieve values at mesh locations */
2286 eq->get_vertex_values = NULL;
2287 eq->get_edge_values = NULL;
2288 eq->get_face_values = cs_cdofb_vecteq_get_face_values;
2289 eq->get_cell_values = cs_cdofb_vecteq_get_cell_values;
2290
2291 eq->get_cw_build_structures = cs_cdofb_vecteq_get;
2292 }
2293 else
2294 bft_error(__FILE__, __LINE__, 0, sv_err_msg, __func__);
2295
2296 break;
2297
2298 case CS_SPACE_SCHEME_CDOEB:
2299 if (eqp->dim == 3) {
2300 /* This is a vector-valued equation but the DoF is scalar-valued since
2301 * it is a circulation associated to each edge */
2302
2303 eq->init_context = cs_cdoeb_vecteq_init_context;
2304 eq->free_context = cs_cdoeb_vecteq_free_context;
2305 eq->init_field_values = cs_cdoeb_vecteq_init_values;
2306
2307 /* Deprecated */
2308 eq->set_dir_bc = NULL;
2309 eq->initialize_system = NULL;
2310 eq->build_system = NULL;
2311 eq->prepare_solving = NULL;
2312 eq->update_field = NULL;
2313
2314 /* New mechanism */
2315 eq->solve_steady_state = cs_cdoeb_vecteq_solve_steady_state;
2316 switch (eqp->time_scheme) {
2317 case CS_TIME_SCHEME_STEADY:
2318 eq->solve = eq->solve_steady_state;
2319 break;
2320
2321 default:
2322 bft_error(__FILE__, __LINE__, 0,
2323 "%s: Eq. %s. This time scheme is not yet implemented",
2324 __func__, eqp->name);
2325 }
2326
2327 eq->postprocess = cs_cdoeb_vecteq_extra_post;
2328 eq->current_to_previous = cs_cdoeb_vecteq_current_to_previous;
2329
2330 eq->read_restart = cs_cdoeb_vecteq_read_restart;
2331 eq->write_restart = cs_cdoeb_vecteq_write_restart;
2332
2333 /* Function pointers to retrieve values at mesh locations */
2334 eq->get_vertex_values = NULL;
2335 eq->get_edge_values = cs_cdoeb_vecteq_get_edge_values;
2336 eq->get_face_values = NULL;
2337 eq->get_cell_values = cs_cdoeb_vecteq_get_cell_values;
2338
2339 eq->get_cw_build_structures = cs_cdoeb_vecteq_get;
2340
2341 }
2342 else
2343 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
2344
2345 break;
2346
2347 case CS_SPACE_SCHEME_HHO_P0:
2348 if (eqp->dim == 1) /* Set pointers of function */
2349 _set_scal_hho_function_pointers(eq);
2350 else
2351 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
2352
2353 break;
2354
2355 case CS_SPACE_SCHEME_HHO_P1:
2356 if (eqp->dim == 1) /* Set pointers of function */
2357 _set_scal_hho_function_pointers(eq);
2358
2359 else if (eqp->dim == 3) /* Set pointers of function */
2360 _set_vect_hho_function_pointers(eq);
2361
2362 else
2363 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
2364
2365 break;
2366
2367 case CS_SPACE_SCHEME_HHO_P2:
2368 if (eqp->dim == 1) /* Set pointers of function */
2369 _set_scal_hho_function_pointers(eq);
2370
2371 else if (eqp->dim == 3) /* Set pointers of function */
2372 _set_vect_hho_function_pointers(eq);
2373
2374 else
2375 bft_error(__FILE__, __LINE__, 0, s_err_msg, __func__);
2376
2377 break;
2378
2379 default:
2380 bft_error(__FILE__, __LINE__, 0,
2381 _(" %s: Eq. %s. Invalid scheme for the space discretization.\n"
2382 " Please check your settings."), __func__, eqp->name);
2383 break;
2384 }
2385
2386 /* Flag this equation such that parametrization is not modifiable anymore */
2387 eqp->flag |= CS_EQUATION_LOCKED;
2388
2389 if (eq->main_ts_id > -1)
2390 cs_timer_stats_stop(eq->main_ts_id);
2391
2392 } /* Loop on equations */
2393
2394 return all_are_steady;
2395 }
2396
2397 /*----------------------------------------------------------------------------*/
2398 /*!
2399 * \brief Create a field structure related to all cs_equation_t structures
2400 */
2401 /*----------------------------------------------------------------------------*/
2402
2403 void
cs_equation_create_fields(void)2404 cs_equation_create_fields(void)
2405 {
2406 for (int eq_id = 0; eq_id < _n_equations; eq_id++) {
2407
2408 int location_id = -1; /* initialize values to avoid a warning */
2409 cs_equation_t *eq = _equations[eq_id];
2410
2411 /* Sanity check */
2412 assert(eq != NULL);
2413
2414 cs_equation_param_t *eqp = eq->param;
2415
2416 /* Redundant definition to handle C/FORTRAN */
2417 int has_previous = (eqp->flag & CS_EQUATION_UNSTEADY) ? 1 : 0;
2418 bool previous = (eqp->flag & CS_EQUATION_UNSTEADY) ? true : false;
2419
2420 if (eq->main_ts_id > -1)
2421 cs_timer_stats_start(eq->main_ts_id);
2422
2423 /* Associate a predefined mesh_location_id to this field */
2424 switch (eqp->space_scheme) {
2425 case CS_SPACE_SCHEME_CDOVB:
2426 case CS_SPACE_SCHEME_CDOVCB:
2427 location_id = cs_mesh_location_get_id_by_name("vertices");
2428 break;
2429 case CS_SPACE_SCHEME_CDOEB:
2430 case CS_SPACE_SCHEME_CDOFB:
2431 case CS_SPACE_SCHEME_HHO_P0:
2432 case CS_SPACE_SCHEME_HHO_P1:
2433 case CS_SPACE_SCHEME_HHO_P2:
2434 location_id = cs_mesh_location_get_id_by_name("cells");
2435 break;
2436 default:
2437 bft_error(__FILE__, __LINE__, 0,
2438 _(" Space scheme for eq. %s is incompatible with a field.\n"
2439 " Stop adding a cs_field_t structure.\n"), eqp->name);
2440 break;
2441 }
2442
2443 if (location_id == -1)
2444 bft_error(__FILE__, __LINE__, 0,
2445 _(" Invalid mesh location id (= -1) for the current field\n"));
2446
2447 /* Store the related field id */
2448 eq->field_id = cs_variable_cdo_field_create(eq->varname,
2449 NULL, /* label */
2450 location_id,
2451 eqp->dim,
2452 has_previous);
2453
2454 /* SLES is associated to a field_id */
2455 eqp->sles_param->field_id = eq->field_id;
2456
2457 if (eqp->process_flag & CS_EQUATION_POST_NORMAL_FLUX) {
2458
2459 /* Add a field for the normal boundary flux */
2460 location_id = cs_mesh_location_get_id_by_name("boundary_faces");
2461
2462 char *bdy_flux_name = NULL;
2463 int len = strlen(eq->varname) + strlen("_normal_boundary_flux") + 2;
2464
2465 BFT_MALLOC(bdy_flux_name, len, char);
2466 sprintf(bdy_flux_name, "%s_normal_boundary_flux", eq->varname);
2467
2468 /* If a scalar: the scalar diffusive flux across the boundary
2469 * If a vector: the vector dot the normal of the boundary face */
2470 int flx_dim = (eqp->dim > 5) ? 3 : 1;
2471 cs_field_t *bdy_flux_fld = cs_field_find_or_create(bdy_flux_name,
2472 0, /* field_mask */
2473 location_id,
2474 flx_dim,
2475 previous);
2476
2477 eq->boundary_flux_id = cs_field_id_by_name(bdy_flux_name);
2478
2479 const int post_flag = CS_POST_ON_LOCATION | CS_POST_MONITOR;
2480 cs_field_set_key_int(bdy_flux_fld, cs_field_key_id("log"), 1);
2481 cs_field_set_key_int(bdy_flux_fld, cs_field_key_id("post_vis"),
2482 post_flag);
2483
2484 BFT_FREE(bdy_flux_name);
2485
2486 } /* Create a boundary flux field */
2487
2488 if (eq->main_ts_id > -1)
2489 cs_timer_stats_stop(eq->main_ts_id);
2490
2491 } /* Loop on equations */
2492
2493 }
2494
2495 /*----------------------------------------------------------------------------*/
2496 /*!
2497 * \brief Allocate and initialize the builder of the algebraic system.
2498 * Set the initialize condition to all variable fields associated to
2499 * each cs_equation_t structure.
2500 *
2501 * \param[in] mesh pointer to a cs_mesh_t structure
2502 * \param[in] ts pointer to a cs_time_step_t structure
2503 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2504 * \param[in, out] connect pointer to a cs_cdo_connect_t structure
2505 */
2506 /*----------------------------------------------------------------------------*/
2507
2508 void
cs_equation_initialize(const cs_mesh_t * mesh,const cs_time_step_t * ts,const cs_cdo_quantities_t * quant,cs_cdo_connect_t * connect)2509 cs_equation_initialize(const cs_mesh_t *mesh,
2510 const cs_time_step_t *ts,
2511 const cs_cdo_quantities_t *quant,
2512 cs_cdo_connect_t *connect)
2513 {
2514 CS_UNUSED(quant);
2515
2516 for (int i = 0; i < _n_equations; i++) {
2517
2518 cs_equation_t *eq = _equations[i];
2519 assert(eq != NULL); /* Sanity check */
2520
2521 if (eq->main_ts_id > -1)
2522 cs_timer_stats_start(eq->main_ts_id);
2523
2524 const cs_equation_param_t *eqp = eq->param;
2525
2526 /* Allocate and initialize a system builder */
2527 /* Not initialized here if it is a restart */
2528 if (eq->builder == NULL)
2529 eq->builder = cs_equation_builder_init(eqp, mesh);
2530
2531 if (eq->scheme_context == NULL)
2532 eq->scheme_context = eq->init_context(eqp,
2533 eq->field_id,
2534 eq->boundary_flux_id,
2535 eq->builder);
2536
2537 /* Define a face interface if not already allocated in this specific case */
2538 if (eqp->n_ic_defs > 0) {
2539 if (cs_param_space_scheme_is_face_based(eqp->space_scheme)) {
2540 if (connect->interfaces[CS_CDO_CONNECT_FACE_SP0] == NULL)
2541 connect->interfaces[CS_CDO_CONNECT_FACE_SP0] =
2542 cs_cdo_connect_define_face_interface(mesh);
2543 }
2544 }
2545
2546 /* Assign an initial value for the variable fields */
2547 if (ts->nt_cur < 1)
2548 eq->init_field_values(ts->t_cur, eq->field_id,
2549 mesh, eqp, eq->builder, eq->scheme_context);
2550
2551 if (eq->main_ts_id > -1)
2552 cs_timer_stats_stop(eq->main_ts_id);
2553
2554 } /* Loop on equations */
2555
2556 }
2557
2558 /*----------------------------------------------------------------------------*/
2559 /*!
2560 * \brief Build the linear system for this equation
2561 *
2562 * \param[in] mesh pointer to a cs_mesh_t structure
2563 * \param[in, out] eq pointer to a cs_equation_t structure
2564 */
2565 /*----------------------------------------------------------------------------*/
2566
2567 void
cs_equation_build_system(const cs_mesh_t * mesh,cs_equation_t * eq)2568 cs_equation_build_system(const cs_mesh_t *mesh,
2569 cs_equation_t *eq)
2570 {
2571 assert(eq != NULL);
2572
2573 const cs_field_t *fld = cs_field_by_id(eq->field_id);
2574
2575 if (eq->main_ts_id > -1)
2576 cs_timer_stats_start(eq->main_ts_id);
2577
2578 /* Sanity checks */
2579 assert(eq->matrix == NULL && eq->rhs == NULL);
2580
2581 /* Initialize the algebraic system to build */
2582 eq->initialize_system(eq->param,
2583 eq->builder,
2584 eq->scheme_context,
2585 &(eq->matrix),
2586 &(eq->rhs));
2587
2588 /* Build the algebraic system to solve */
2589 eq->build_system(mesh, fld->val,
2590 eq->param,
2591 eq->builder,
2592 eq->scheme_context,
2593 eq->rhs,
2594 eq->matrix);
2595
2596 if (eq->main_ts_id > -1)
2597 cs_timer_stats_stop(eq->main_ts_id);
2598 }
2599
2600 /*----------------------------------------------------------------------------*/
2601 /*!
2602 * \brief Solve the linear system for this equation
2603 *
2604 * \param[in, out] eq pointer to a cs_equation_t structure
2605 */
2606 /*----------------------------------------------------------------------------*/
2607
2608 void
cs_equation_solve_deprecated(cs_equation_t * eq)2609 cs_equation_solve_deprecated(cs_equation_t *eq)
2610 {
2611 int n_iters = 0;
2612 double residual = DBL_MAX;
2613 cs_sles_t *sles = cs_sles_find_or_add(eq->field_id, NULL);
2614 cs_field_t *fld = cs_field_by_id(eq->field_id);
2615 cs_real_t *x = NULL, *b = NULL;
2616
2617 if (eq->main_ts_id > -1)
2618 cs_timer_stats_start(eq->main_ts_id);
2619
2620 const cs_equation_param_t *eqp = eq->param;
2621 const double r_norm = 1.0; /* No renormalization by default (TODO) */
2622 const cs_param_sles_t *sles_param = eqp->sles_param;
2623
2624 /* Sanity checks (up to now, only scalar field are handled) */
2625 assert(eq->n_sles_gather_elts <= eq->n_sles_scatter_elts);
2626 assert(eq->n_sles_gather_elts == cs_matrix_get_n_rows(eq->matrix));
2627
2628 #if defined(DEBUG) && !defined(NDEBUG) && CS_EQUATION_DBG > 0
2629 cs_log_printf(CS_LOG_DEFAULT,
2630 " n_sles_gather_elts: %d\n"
2631 " n_sles_scatter_elts: %d\n"
2632 " n_matrix_rows: %d\n"
2633 " n_matrix_columns: %d\n",
2634 eq->n_sles_gather_elts, eq->n_sles_scatter_elts,
2635 cs_matrix_get_n_rows(eq->matrix),
2636 cs_matrix_get_n_columns(eq->matrix));
2637 #endif
2638
2639 /* Handle parallelism */
2640 eq->prepare_solving(eq, &x, &b);
2641
2642 cs_sles_convergence_state_t code = cs_sles_solve(sles,
2643 eq->matrix,
2644 sles_param->eps,
2645 r_norm,
2646 &n_iters,
2647 &residual,
2648 b,
2649 x,
2650 0, /* aux. size */
2651 NULL); /* aux. buffers */
2652
2653 if (sles_param->verbosity > 0) {
2654
2655 const cs_lnum_t size = eq->n_sles_gather_elts;
2656 const cs_lnum_t *row_index, *col_id;
2657 const cs_real_t *d_val, *x_val;
2658
2659 cs_matrix_get_msr_arrays(eq->matrix, &row_index, &col_id, &d_val, &x_val);
2660
2661 #if defined(DEBUG) && !defined(NDEBUG) && CS_EQUATION_DBG > 1
2662 cs_dbg_dump_linear_system(eqp->name, size, CS_EQUATION_DBG,
2663 x, b,
2664 row_index, col_id, x_val, d_val);
2665 #endif
2666
2667 cs_gnum_t nnz = row_index[size];
2668 cs_parall_counter(&nnz, 1);
2669 cs_log_printf(CS_LOG_DEFAULT, " <%s/sles_cvg> code %-d n_iters %d"
2670 " residual % -8.4e nnz %lu\n",
2671 eqp->name, code, n_iters, residual, nnz);
2672
2673 }
2674
2675 if (cs_glob_n_ranks > 1) { /* Parallel mode */
2676
2677 cs_range_set_scatter(eq->rset,
2678 CS_REAL_TYPE, 1, // type and stride
2679 x,
2680 x);
2681
2682 cs_range_set_scatter(eq->rset,
2683 CS_REAL_TYPE, 1, // type and stride
2684 b,
2685 eq->rhs);
2686
2687 }
2688
2689 /* Copy current field values to previous values */
2690 cs_field_current_to_previous(fld);
2691
2692 /* Define the new field value for the current time */
2693 eq->update_field(x, eq->rhs, eq->param,
2694 eq->builder, eq->scheme_context, fld->val);
2695
2696 if (eq->main_ts_id > -1)
2697 cs_timer_stats_stop(eq->main_ts_id);
2698
2699 /* Free memory */
2700 BFT_FREE(x);
2701 if (b != eq->rhs)
2702 BFT_FREE(b);
2703 BFT_FREE(eq->rhs);
2704 cs_sles_free(sles);
2705 cs_matrix_destroy(&(eq->matrix));
2706 }
2707
2708 /*----------------------------------------------------------------------------*/
2709 /*!
2710 * \brief Build and then solve the linear system for this equation when the
2711 * goal is to find the steady state
2712 *
2713 * \param[in] mesh pointer to a cs_mesh_t structure
2714 * \param[in, out] eq pointer to a cs_equation_t structure
2715 */
2716 /*----------------------------------------------------------------------------*/
2717
2718 void
cs_equation_solve_steady_state(const cs_mesh_t * mesh,cs_equation_t * eq)2719 cs_equation_solve_steady_state(const cs_mesh_t *mesh,
2720 cs_equation_t *eq)
2721 {
2722 if (eq == NULL)
2723 bft_error(__FILE__, __LINE__, 0, "%s: Empty equation structure", __func__);
2724 assert(cs_equation_uses_new_mechanism(eq));
2725
2726 if (eq->main_ts_id > -1)
2727 cs_timer_stats_start(eq->main_ts_id);
2728
2729 /* Allocate, build and solve the algebraic system:
2730 * The linear solver is called inside and the field value is updated inside
2731 */
2732
2733 eq->solve_steady_state(false, /* current to previous */
2734 mesh,
2735 eq->field_id,
2736 eq->param,
2737 eq->builder,
2738 eq->scheme_context);
2739
2740 if (eq->main_ts_id > -1)
2741 cs_timer_stats_stop(eq->main_ts_id);
2742 }
2743
2744 /*----------------------------------------------------------------------------*/
2745 /*!
2746 * \brief Build and then solve the linear system for an equation with an
2747 * unsteady term
2748 *
2749 * \param[in] cur2prev true="current to previous" operation is performed
2750 * \param[in] mesh pointer to a cs_mesh_t structure
2751 * \param[in, out] eq pointer to a cs_equation_t structure
2752 */
2753 /*----------------------------------------------------------------------------*/
2754
2755 void
cs_equation_solve(bool cur2prev,const cs_mesh_t * mesh,cs_equation_t * eq)2756 cs_equation_solve(bool cur2prev,
2757 const cs_mesh_t *mesh,
2758 cs_equation_t *eq)
2759 {
2760 if (eq == NULL)
2761 bft_error(__FILE__, __LINE__, 0, "%s: Empty equation structure", __func__);
2762 assert(cs_equation_uses_new_mechanism(eq));
2763
2764 if (eq->main_ts_id > -1)
2765 cs_timer_stats_start(eq->main_ts_id);
2766
2767 /* Allocate, build and solve the algebraic system:
2768 * The linear solver is called inside and the field value is updated inside
2769 */
2770
2771 eq->solve(cur2prev,
2772 mesh,
2773 eq->field_id,
2774 eq->param,
2775 eq->builder,
2776 eq->scheme_context);
2777
2778 if (eq->main_ts_id > -1)
2779 cs_timer_stats_stop(eq->main_ts_id);
2780 }
2781
2782 /*----------------------------------------------------------------------------*/
2783 /*!
2784 * \brief Build and then solve the linear system for a steady-state equation.
2785 * This is wrapper for the FORTRAN interface (limitation of the
2786 * parameters to simple types).
2787 *
2788 * \param[in] eqname name of the equation to solve
2789 */
2790 /*----------------------------------------------------------------------------*/
2791
2792 void
cs_equation_solve_steady_state_wrapper(const char * eqname)2793 cs_equation_solve_steady_state_wrapper(const char *eqname)
2794 {
2795 cs_equation_t *eq = cs_equation_by_name(eqname);
2796
2797 if (eq == NULL)
2798 bft_error(__FILE__, __LINE__, 0, "%s: Empty equation structure", __func__);
2799 assert(cs_equation_uses_new_mechanism(eq));
2800
2801 if (eq->main_ts_id > -1)
2802 cs_timer_stats_start(eq->main_ts_id);
2803
2804 /* Allocate, build and solve the algebraic system:
2805 The linear solver is called inside and the field value is updated inside
2806 */
2807
2808 eq->solve_steady_state(false, /* current to previous */
2809 cs_glob_mesh,
2810 eq->field_id,
2811 eq->param,
2812 eq->builder,
2813 eq->scheme_context);
2814
2815 if (eq->main_ts_id > -1)
2816 cs_timer_stats_stop(eq->main_ts_id);
2817 }
2818
2819 /*----------------------------------------------------------------------------*/
2820 /*!
2821 * \brief Build and then solve the linear system for an equation with an
2822 * unsteady term. This is wrapper for the FORTRAN interface (limitation
2823 * of the parameters to simple types)
2824 *
2825 * \param[in] cur2prev true="current to previous" operation is performed
2826 * \param[in] eqname name of the equation to solve
2827 */
2828 /*----------------------------------------------------------------------------*/
2829
2830 void
cs_equation_solve_wrapper(bool cur2prev,const char * eqname)2831 cs_equation_solve_wrapper(bool cur2prev,
2832 const char *eqname)
2833 {
2834 cs_equation_t *eq = cs_equation_by_name(eqname);
2835
2836 if (eq == NULL)
2837 bft_error(__FILE__, __LINE__, 0, "%s: Empty equation structure", __func__);
2838 assert(cs_equation_uses_new_mechanism(eq));
2839
2840 if (eq->main_ts_id > -1)
2841 cs_timer_stats_start(eq->main_ts_id);
2842
2843 /* Allocate, build and solve the algebraic system:
2844 The linear solver is called inside and the field value is updated inside
2845 */
2846 eq->solve(cur2prev,
2847 cs_glob_mesh,
2848 eq->field_id,
2849 eq->param,
2850 eq->builder,
2851 eq->scheme_context);
2852
2853 if (eq->main_ts_id > -1)
2854 cs_timer_stats_stop(eq->main_ts_id);
2855 }
2856
2857 /*----------------------------------------------------------------------------*/
2858 /*!
2859 * \brief Apply the current to previous to all fields (and potentially arrays)
2860 * related to an equation. This function fas to be called when a solve
2861 * step is called with the parameter: cur2prev = false
2862 *
2863 * \param[in] eq pointer to a \ref cs_equation_t structure
2864 */
2865 /*----------------------------------------------------------------------------*/
2866
2867 void
cs_equation_current_to_previous(const cs_equation_t * eq)2868 cs_equation_current_to_previous(const cs_equation_t *eq)
2869 {
2870 if (eq == NULL)
2871 bft_error(__FILE__, __LINE__, 0, "%s: Empty equation structure", __func__);
2872
2873 if (eq->current_to_previous == NULL)
2874 return;
2875
2876 if (eq->main_ts_id > -1)
2877 cs_timer_stats_start(eq->main_ts_id);
2878
2879 eq->current_to_previous(eq->param, eq->builder, eq->scheme_context);
2880
2881 if (eq->main_ts_id > -1)
2882 cs_timer_stats_stop(eq->main_ts_id);
2883 }
2884
2885 /*----------------------------------------------------------------------------*/
2886 /*!
2887 * \brief For a given equation, retrieve the related cellwise builder
2888 * structures: cs_cell_builder_t and cs_cell_system_t structures
2889 *
2890 * \param[in] eq pointer to a \ref cs_equation_t structure
2891 * \param[out] cb pointer to a pointer on a cs_cell_sys_t structure
2892 * \param[out] csys pointer to a pointer on a cs_cell_builder_t structure
2893 */
2894 /*----------------------------------------------------------------------------*/
2895
2896 void
cs_equation_get_cellwise_builders(const cs_equation_t * eq,cs_cell_sys_t ** csys,cs_cell_builder_t ** cb)2897 cs_equation_get_cellwise_builders(const cs_equation_t *eq,
2898 cs_cell_sys_t **csys,
2899 cs_cell_builder_t **cb)
2900 {
2901 *csys = NULL;
2902 *cb = NULL;
2903
2904 if (eq == NULL)
2905 return;
2906
2907 if (eq->get_cw_build_structures != NULL)
2908 eq->get_cw_build_structures(csys, cb);
2909 }
2910
2911 /*----------------------------------------------------------------------------*/
2912 /*!
2913 * \brief For a given equation, retrieve an array of values related to each
2914 * cell of the mesh for the unknowns
2915 *
2916 * \param[in] eq pointer to a \ref cs_equation_t structure
2917 * \param[in] previous retrieve the previous state (true/false)
2918 *
2919 * \return a pointer to an array of cell values
2920 */
2921 /*----------------------------------------------------------------------------*/
2922
2923 cs_real_t *
cs_equation_get_cell_values(const cs_equation_t * eq,bool previous)2924 cs_equation_get_cell_values(const cs_equation_t *eq,
2925 bool previous)
2926 {
2927 if (eq == NULL)
2928 return NULL;
2929
2930 cs_real_t *c_values = NULL;
2931 if (eq->get_cell_values != NULL)
2932 c_values = eq->get_cell_values(eq->scheme_context, previous);
2933
2934 return c_values;
2935 }
2936
2937 /*----------------------------------------------------------------------------*/
2938 /*!
2939 * \brief For a given equation, retrieve an array of values related to each
2940 * face of the mesh for the unknowns
2941 *
2942 * \param[in] eq pointer to a \ref cs_equation_t structure
2943 * \param[in] previous retrieve the previous state (true/false)
2944 *
2945 * \return a pointer to an array of face values
2946 */
2947 /*----------------------------------------------------------------------------*/
2948
2949 cs_real_t *
cs_equation_get_face_values(const cs_equation_t * eq,bool previous)2950 cs_equation_get_face_values(const cs_equation_t *eq,
2951 bool previous)
2952 {
2953 if (eq == NULL)
2954 return NULL;
2955
2956 cs_real_t *f_values = NULL;
2957 if (eq->get_face_values != NULL)
2958 f_values = eq->get_face_values(eq->scheme_context, previous);
2959
2960 return f_values;
2961 }
2962
2963 /*----------------------------------------------------------------------------*/
2964 /*!
2965 * \brief For a given equation, retrieve an array of values related to each
2966 * edge of the mesh for the unknowns
2967 *
2968 * \param[in] eq pointer to a \ref cs_equation_t structure
2969 * \param[in] previous retrieve the previous state (true/false)
2970 *
2971 * \return a pointer to an array of edge values
2972 */
2973 /*----------------------------------------------------------------------------*/
2974
2975 cs_real_t *
cs_equation_get_edge_values(const cs_equation_t * eq,bool previous)2976 cs_equation_get_edge_values(const cs_equation_t *eq,
2977 bool previous)
2978 {
2979 if (eq == NULL)
2980 return NULL;
2981
2982 cs_real_t *e_values = NULL;
2983 if (eq->get_edge_values != NULL)
2984 e_values = eq->get_edge_values(eq->scheme_context, previous);
2985
2986 return e_values;
2987 }
2988
2989 /*----------------------------------------------------------------------------*/
2990 /*!
2991 * \brief For a given equation, retrieve an array of values related to each
2992 * vertex of the mesh for the unknowns
2993 *
2994 * \param[in] eq pointer to a \ref cs_equation_t structure
2995 * \param[in] previous retrieve the previous state (true/false)
2996 *
2997 * \return a pointer to an array of vertex values
2998 */
2999 /*----------------------------------------------------------------------------*/
3000
3001 cs_real_t *
cs_equation_get_vertex_values(const cs_equation_t * eq,bool previous)3002 cs_equation_get_vertex_values(const cs_equation_t *eq,
3003 bool previous)
3004 {
3005 if (eq == NULL)
3006 return NULL;
3007
3008 cs_real_t *v_values = NULL;
3009 if (eq->get_vertex_values != NULL)
3010 v_values = eq->get_vertex_values(eq->scheme_context, previous);
3011
3012 return v_values;
3013 }
3014
3015 /*----------------------------------------------------------------------------*/
3016 /*!
3017 * \brief Compute the integral over the domain of the current variable field
3018 * associated to the given equation.
3019 *
3020 * \param[in] connect pointer to a \ref cs_cdo_connect_t structure
3021 * \param[in] cdoq pointer to a \ref cs_cdo_quantities_t structure
3022 * \param[in] eq pointer to a \ref cs_equation_t structure
3023 * \param[in, out] integral result of the computation
3024 */
3025 /*----------------------------------------------------------------------------*/
3026
3027 void
cs_equation_integrate_variable(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,const cs_equation_t * eq,cs_real_t * result)3028 cs_equation_integrate_variable(const cs_cdo_connect_t *connect,
3029 const cs_cdo_quantities_t *cdoq,
3030 const cs_equation_t *eq,
3031 cs_real_t *result)
3032 {
3033 /* Initialize the value to return */
3034 *result = 0.;
3035
3036 if (eq == NULL)
3037 return;
3038
3039 const cs_equation_param_t *eqp = eq->param;
3040 assert(eqp != NULL);
3041 if (eqp->dim > 1)
3042 bft_error(__FILE__, __LINE__, 0, "%s: (Eq. %s) Not implemented",
3043 __func__, eqp->name);
3044
3045 /* Scalar-valued equation */
3046 switch (eqp->space_scheme) {
3047
3048 case CS_SPACE_SCHEME_CDOVB:
3049 {
3050 const cs_real_t *p_v = cs_equation_get_vertex_values(eq, false);
3051 const cs_adjacency_t *c2v = connect->c2v;
3052
3053 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
3054
3055 cs_real_t int_cell = 0.;
3056 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++)
3057 int_cell += cdoq->dcell_vol[j] * p_v[c2v->ids[j]];
3058
3059 *result += int_cell;
3060
3061 } /* Loop on cells */
3062
3063 }
3064 break;
3065
3066 case CS_SPACE_SCHEME_CDOVCB:
3067 {
3068 const cs_real_t *p_v = cs_equation_get_vertex_values(eq, false);
3069 const cs_real_t *p_c = cs_equation_get_cell_values(eq, false);
3070 const cs_adjacency_t *c2v = connect->c2v;
3071
3072 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
3073
3074 /* Shares between cell and vertex unknowns:
3075 - the cell unknown stands for 1/4 of the cell volume
3076 - the vertex unknown stands for 3/4 of the dual cell volume
3077 */
3078 cs_real_t int_cell = 0.25*cdoq->cell_vol[c_id]*p_c[c_id];
3079 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++)
3080 int_cell += 0.75*cdoq->dcell_vol[j] * p_v[c2v->ids[j]];
3081
3082 *result += int_cell;
3083
3084 } /* Loop on cells */
3085
3086 }
3087 break;
3088
3089 case CS_SPACE_SCHEME_CDOFB:
3090 {
3091 const cs_real_t *p_f = cs_equation_get_face_values(eq, false);
3092 const cs_real_t *p_c = cs_equation_get_cell_values(eq, false);
3093 const cs_adjacency_t *c2f = connect->c2f;
3094
3095 assert(cdoq->pvol_fc != NULL); /* Sanity check */
3096
3097 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
3098
3099 /* Shares between cell and face unknowns (assuming stab.coeff = 1/3):
3100 - the cell unknown stands for 1/4 of the cell volume
3101 - the face unknown stands for 3/4 of the dual cell volume
3102 */
3103 cs_real_t int_cell = 0.25*cdoq->cell_vol[c_id]*p_c[c_id];
3104 for (cs_lnum_t j = c2f->idx[c_id]; j < c2f->idx[c_id+1]; j++)
3105 int_cell += 0.75*cdoq->pvol_fc[j] * p_f[c2f->ids[j]];
3106
3107 *result += int_cell;
3108
3109 } /* Loop on cells */
3110
3111 }
3112 break;
3113
3114 default:
3115 bft_error(__FILE__, __LINE__, 0,
3116 "%s: (Eq. %s). Not implemented.", __func__, eqp->name);
3117
3118 } /* End of switch */
3119
3120 /* Parallel synchronization */
3121 cs_parall_sum(1, CS_REAL_TYPE, result);
3122 }
3123
3124 /*----------------------------------------------------------------------------*/
3125 /*!
3126 * \brief Compute the diffusive flux across all boundary faces
3127 * According to the space discretization scheme, the size of the
3128 * resulting array differs.
3129 * For Vb and VCb schemes, this array relies on the bf2v adjacency.
3130 *
3131 * \param[in] t_eval time at which one performs the property evaluation
3132 * \param[in] eq pointer to a cs_equation_t structure
3133 * \param[in, out] diff_flux value of the diffusive part of the flux
3134 */
3135 /*----------------------------------------------------------------------------*/
3136
3137 void
cs_equation_compute_boundary_diff_flux(cs_real_t t_eval,const cs_equation_t * eq,cs_real_t * diff_flux)3138 cs_equation_compute_boundary_diff_flux(cs_real_t t_eval,
3139 const cs_equation_t *eq,
3140 cs_real_t *diff_flux)
3141 {
3142 if (diff_flux == NULL)
3143 return;
3144
3145 if (eq == NULL)
3146 bft_error(__FILE__, __LINE__, 0, _err_empty_eq, __func__);
3147
3148 const cs_equation_param_t *eqp = eq->param;
3149
3150 assert(eqp != NULL);
3151 if (eqp->dim > 1)
3152 bft_error(__FILE__, __LINE__, 0, "%s: (Eq. %s) Not implemented",
3153 __func__, eqp->name);
3154
3155 /* Scalar-valued equation */
3156 switch (eqp->space_scheme) {
3157
3158 case CS_SPACE_SCHEME_CDOVB:
3159 {
3160 const cs_real_t *p_v = cs_equation_get_vertex_values(eq, false);
3161
3162 cs_cdovb_scaleq_boundary_diff_flux(t_eval,
3163 eqp,
3164 p_v,
3165 eq->builder,
3166 eq->scheme_context,
3167 diff_flux);
3168 }
3169 break;
3170
3171 case CS_SPACE_SCHEME_CDOVCB:
3172 {
3173 const cs_real_t *p_v = cs_equation_get_vertex_values(eq, false);
3174 const cs_real_t *p_c = cs_equation_get_cell_values(eq, false);
3175
3176 cs_cdovcb_scaleq_boundary_diff_flux(t_eval,
3177 eqp,
3178 p_v,
3179 p_c,
3180 eq->builder,
3181 eq->scheme_context,
3182 diff_flux);
3183 }
3184 break;
3185
3186 case CS_SPACE_SCHEME_CDOFB:
3187 {
3188 const cs_real_t *p_f = cs_equation_get_face_values(eq, false);
3189 const cs_real_t *p_c = cs_equation_get_cell_values(eq, false);
3190
3191 cs_cdofb_scaleq_boundary_diff_flux(t_eval,
3192 eqp,
3193 p_f,
3194 p_c,
3195 eq->builder,
3196 eq->scheme_context,
3197 diff_flux);
3198 }
3199 break;
3200
3201 default:
3202 bft_error(__FILE__, __LINE__, 0,
3203 "%s: (Eq. %s). Not implemented.", __func__, eqp->name);
3204
3205 } /* End of switch */
3206
3207 }
3208
3209 /*----------------------------------------------------------------------------*/
3210 /*!
3211 * \brief Compute the diffusive and convective flux across a plane defined
3212 * by a mesh location structure attached to the name ml_name.
3213 *
3214 * \param[in] eq pointer to a cs_equation_t structure
3215 * \param[in] ml_name name of the related mesh location
3216 * \param[in] direction vector indicating in which direction flux is > 0
3217 * \param[in, out] diff_flux value of the diffusive part of the flux
3218 * \param[in, out] conv_flux value of the convective part of the flux
3219 */
3220 /*----------------------------------------------------------------------------*/
3221
3222 void
cs_equation_compute_flux_across_plane(const cs_equation_t * eq,const char * ml_name,const cs_real_3_t direction,cs_real_t * diff_flux,cs_real_t * conv_flux)3223 cs_equation_compute_flux_across_plane(const cs_equation_t *eq,
3224 const char *ml_name,
3225 const cs_real_3_t direction,
3226 cs_real_t *diff_flux,
3227 cs_real_t *conv_flux)
3228 {
3229 if (eq == NULL)
3230 bft_error(__FILE__, __LINE__, 0, _err_empty_eq, __func__);
3231
3232 /* Get the mesh location id from its name */
3233 const int ml_id = cs_mesh_location_get_id_by_name(ml_name);
3234 if (ml_id == -1)
3235 bft_error(__FILE__, __LINE__, 0,
3236 " %s: Invalid mesh location name %s.\n"
3237 " This mesh location is not already defined.\n",
3238 __func__, ml_name);
3239
3240 const char emsg[] = " %s: Computation of the diffusive and convective flux"
3241 " across a plane\n is not available for equation %s\n";
3242
3243 /* Retrieve the field from its id */
3244 cs_field_t *fld = cs_field_by_id(eq->field_id);
3245
3246 cs_equation_param_t *eqp = eq->param;
3247 assert(eqp != NULL && fld != NULL);
3248
3249 if (eqp->dim > 1)
3250 bft_error(__FILE__, __LINE__, 0, emsg, __func__, eqp->name);
3251
3252
3253 /* Scalar-valued equation */
3254 switch (eqp->space_scheme) {
3255
3256 case CS_SPACE_SCHEME_CDOVB:
3257 cs_cdovb_scaleq_flux_across_plane(direction,
3258 fld->val,
3259 eqp,
3260 ml_id,
3261 eq->builder,
3262 eq->scheme_context,
3263 diff_flux, conv_flux);
3264 break;
3265
3266 case CS_SPACE_SCHEME_CDOVCB:
3267 cs_cdovcb_scaleq_flux_across_plane(direction,
3268 fld->val,
3269 eqp,
3270 ml_id,
3271 eq->builder,
3272 eq->scheme_context,
3273 diff_flux, conv_flux);
3274 break;
3275
3276 default:
3277 bft_error(__FILE__, __LINE__, 0, emsg, __func__, eqp->name);
3278
3279 } /* End of switch */
3280
3281 }
3282
3283 /*----------------------------------------------------------------------------*/
3284 /*!
3285 * \brief Cellwise computation of the diffusive flux across all cell faces.
3286 * Primal or dual faces are considered according to the space scheme.
3287 *
3288 * \param[in] eq pointer to a cs_equation_t structure
3289 * \param[in] location indicate where the flux has to be computed
3290 * \param[in] t_eval time at which one performs the evaluation
3291 * \param[in, out] diff_flux value of the diffusive flux
3292 */
3293 /*----------------------------------------------------------------------------*/
3294
3295 void
cs_equation_compute_diff_flux_cellwise(const cs_equation_t * eq,cs_flag_t location,cs_real_t t_eval,cs_real_t * diff_flux)3296 cs_equation_compute_diff_flux_cellwise(const cs_equation_t *eq,
3297 cs_flag_t location,
3298 cs_real_t t_eval,
3299 cs_real_t *diff_flux)
3300 {
3301 if (diff_flux == NULL)
3302 return;
3303 if (eq == NULL)
3304 bft_error(__FILE__, __LINE__, 0, _err_empty_eq, __func__);
3305
3306 cs_equation_param_t *eqp = eq->param;
3307 assert(eqp != NULL);
3308
3309 /* Retrieve the field from its id */
3310 cs_field_t *fld = cs_field_by_id(eq->field_id);
3311
3312 const char fmsg[] = " %s: (Eq. %s) Stop computing the diffusive flux.\n"
3313 " This functionality is not available for this scheme.";
3314 const char lmsg[] = " %s: (Eq. %s) Stop computing the diffusive flux.\n"
3315 " This mesh location is not available for this scheme.";
3316
3317 if (eqp->dim > 1)
3318 bft_error(__FILE__, __LINE__, 0, fmsg, __func__, eqp->name);
3319
3320 switch (eqp->space_scheme) {
3321
3322 case CS_SPACE_SCHEME_CDOVB:
3323 if (cs_flag_test(location, cs_flag_primal_cell))
3324 cs_cdovb_scaleq_diff_flux_in_cells(fld->val,
3325 eqp,
3326 t_eval,
3327 eq->builder,
3328 eq->scheme_context,
3329 diff_flux);
3330 else if (cs_flag_test(location, cs_flag_dual_face_byc))
3331 cs_cdovb_scaleq_diff_flux_dfaces(fld->val,
3332 eqp,
3333 t_eval,
3334 eq->builder,
3335 eq->scheme_context,
3336 diff_flux);
3337 else
3338 bft_error(__FILE__, __LINE__, 0, lmsg, __func__, eqp->name);
3339 break;
3340
3341 case CS_SPACE_SCHEME_CDOVCB:
3342 if (cs_flag_test(location, cs_flag_primal_cell))
3343 cs_cdovcb_scaleq_diff_flux_in_cells(fld->val,
3344 eqp,
3345 t_eval,
3346 eq->builder,
3347 eq->scheme_context,
3348 diff_flux);
3349 else if (cs_flag_test(location, cs_flag_dual_face_byc))
3350 cs_cdovcb_scaleq_diff_flux_dfaces(fld->val,
3351 eqp,
3352 t_eval,
3353 eq->builder,
3354 eq->scheme_context,
3355 diff_flux);
3356 else
3357 bft_error(__FILE__, __LINE__, 0, lmsg, __func__, eqp->name);
3358 break;
3359
3360 default:
3361 bft_error(__FILE__, __LINE__, 0, fmsg, __func__, eqp->name);
3362
3363 } /* End of switch on the space scheme */
3364
3365 }
3366
3367 /*----------------------------------------------------------------------------*/
3368 /*!
3369 * \brief Cellwise computation of the discrete gradient at vertices
3370 *
3371 * \param[in] eq pointer to a cs_equation_t structure
3372 * \param[in, out] v_gradient gradient at vertices
3373 */
3374 /*----------------------------------------------------------------------------*/
3375
3376 void
cs_equation_compute_vtx_field_gradient(const cs_equation_t * eq,cs_real_t * v_gradient)3377 cs_equation_compute_vtx_field_gradient(const cs_equation_t *eq,
3378 cs_real_t *v_gradient)
3379 {
3380 if (eq == NULL)
3381 bft_error(__FILE__, __LINE__, 0, _err_empty_eq, __func__);
3382 assert(v_gradient != NULL);
3383
3384 const cs_equation_param_t *eqp = eq->param;
3385
3386 /* Retrieve the field from its id */
3387 cs_field_t *fld = cs_field_by_id(eq->field_id);
3388
3389 switch (eqp->space_scheme) {
3390
3391 case CS_SPACE_SCHEME_CDOVCB:
3392 cs_cdovcb_scaleq_vtx_gradient(fld->val, eq->builder, eq->scheme_context,
3393 v_gradient);
3394 break;
3395
3396 case CS_SPACE_SCHEME_CDOVB: /* --> wall distance */
3397 default:
3398 bft_error(__FILE__, __LINE__, 0,
3399 " %s: Invalid type of scheme for equation %s when computing"
3400 " the gradient at vertices", __func__, eqp->name);
3401 break;
3402
3403 }
3404
3405 }
3406
3407 /*----------------------------------------------------------------------------*/
3408 /*!
3409 * \brief Compute and post-process Peclet number if requested
3410 *
3411 * \param[in] eq pointer to a cs_equation_t structure
3412 * \param[in] ts pointer to a cs_time_step_t struct.
3413 * \param[in, out] peclet pointer to an array storing the resulting Peclet
3414 * number in each cell
3415 */
3416 /*----------------------------------------------------------------------------*/
3417
3418 void
cs_equation_compute_peclet(const cs_equation_t * eq,const cs_time_step_t * ts,cs_real_t peclet[])3419 cs_equation_compute_peclet(const cs_equation_t *eq,
3420 const cs_time_step_t *ts,
3421 cs_real_t peclet[])
3422 {
3423 if (eq == NULL)
3424 bft_error(__FILE__, __LINE__, 0, _err_empty_eq, __func__);
3425 assert(peclet != NULL);
3426
3427 const cs_equation_param_t *eqp = eq->param;
3428
3429 /* Check if the computation of the Peclet number is requested */
3430 if (!(eqp->process_flag & CS_EQUATION_POST_PECLET))
3431 return;
3432
3433 if (eqp->diffusion_property == NULL)
3434 bft_error(__FILE__, __LINE__, 0,
3435 "%s: Computation of the Peclet number is requested for\n"
3436 " equation %s but no diffusion property is set.\n",
3437 __func__, eqp->name);
3438 if (eqp->adv_field == NULL)
3439 bft_error(__FILE__, __LINE__, 0,
3440 "%s: Computation of the Peclet number is requested for\n"
3441 " equation %s but no advection field is set.\n",
3442 __func__, eqp->name);
3443
3444 if (eq->main_ts_id > -1) /* Activate timer statistics */
3445 cs_timer_stats_start(eq->main_ts_id);
3446
3447 /* Compute the Peclet number in each cell */
3448 cs_advection_get_peclet(eqp->adv_field,
3449 eqp->diffusion_property,
3450 ts->t_cur,
3451 peclet);
3452
3453 if (eq->main_ts_id > -1)
3454 cs_timer_stats_stop(eq->main_ts_id);
3455 }
3456
3457 /*----------------------------------------------------------------------------*/
3458 /*!
3459 * \brief Write into the restart file additionnal arrays (not defined as
3460 * fields) but useful for the checkpoint/restart process
3461 *
3462 * \param[in, out] restart pointer to a \ref cs_restart_t structure
3463 */
3464 /*----------------------------------------------------------------------------*/
3465
3466 void
cs_equation_read_extra_restart(cs_restart_t * restart)3467 cs_equation_read_extra_restart(cs_restart_t *restart)
3468 {
3469 for (int i = 0; i < _n_equations; i++) {
3470
3471 cs_equation_t *eq = _equations[i];
3472 assert(eq != NULL); /* Sanity check */
3473
3474 if (eq->read_restart != NULL)
3475 eq->read_restart(restart, eq->param->name, eq->scheme_context);
3476
3477 } /* Loop on equations */
3478 }
3479
3480 /*----------------------------------------------------------------------------*/
3481 /*!
3482 * \brief Write into the restart file additionnal arrays (not defined as
3483 * fields) but useful for the checkpoint/restart process
3484 *
3485 * \param[in, out] restart pointer to a \ref cs_restart_t structure
3486 */
3487 /*----------------------------------------------------------------------------*/
3488
3489 void
cs_equation_write_extra_restart(cs_restart_t * restart)3490 cs_equation_write_extra_restart(cs_restart_t *restart)
3491 {
3492 for (int i = 0; i < _n_equations; i++) {
3493
3494 cs_equation_t *eq = _equations[i];
3495 assert(eq != NULL); /* Sanity check */
3496
3497 if (eq->write_restart != NULL)
3498 eq->write_restart(restart, eq->param->name, eq->scheme_context);
3499
3500 } /* Loop on equations */
3501 }
3502
3503 /*----------------------------------------------------------------------------*/
3504 /*!
3505 * \brief Predefined extra-operations related to all equations
3506 *
3507 * \param[in] mesh pointer to a cs_mesh_t structure
3508 * \param[in] connect pointer to a cs_cdo_connect_t structure
3509 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
3510 * \param[in] ts pointer to a cs_time_step_t struct.
3511 */
3512 /*----------------------------------------------------------------------------*/
3513
3514 void
cs_equation_post_balance(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,const cs_time_step_t * ts)3515 cs_equation_post_balance(const cs_mesh_t *mesh,
3516 const cs_cdo_connect_t *connect,
3517 const cs_cdo_quantities_t *cdoq,
3518 const cs_time_step_t *ts)
3519 {
3520 CS_UNUSED(connect);
3521 CS_UNUSED(cdoq);
3522
3523 for (int i = 0; i < _n_equations; i++) {
3524
3525 cs_equation_t *eq = _equations[i];
3526 assert(eq != NULL); /* Sanity check */
3527 const cs_equation_param_t *eqp = eq->param;
3528
3529 /* Check if the computation of the balance is requested */
3530 if (!(eqp->process_flag & CS_EQUATION_POST_BALANCE))
3531 continue;
3532
3533 if (eq->compute_balance != NULL)
3534 bft_error(__FILE__, __LINE__, 0,
3535 "%s: Balance for equation %s is requested but\n"
3536 " this functionality is not available yet.\n",
3537 __func__, eqp->name);
3538
3539 if (eq->main_ts_id > -1) /* Activate timer statistics */
3540 cs_timer_stats_start(eq->main_ts_id);
3541
3542 cs_equation_balance_t *b = eq->compute_balance(eqp,
3543 eq->builder,
3544 eq->scheme_context);
3545
3546 char *postlabel = NULL;
3547 int len = strlen(eqp->name) + 13 + 1;
3548 BFT_MALLOC(postlabel, len, char);
3549
3550 switch (eqp->space_scheme) {
3551
3552 case CS_SPACE_SCHEME_CDOVB:
3553 {
3554 sprintf(postlabel, "%s.Balance", eqp->name);
3555
3556 cs_post_write_vertex_var(CS_POST_MESH_VOLUME,
3557 CS_POST_WRITER_DEFAULT,
3558 postlabel,
3559 eqp->dim,
3560 false,
3561 false,
3562 CS_POST_TYPE_cs_real_t,
3563 b->balance,
3564 ts);
3565
3566 if (cs_equation_param_has_diffusion(eqp))
3567 _post_balance_at_vertices(eq, ts, "Diff", postlabel,
3568 b->diffusion_term);
3569
3570 if (cs_equation_param_has_convection(eqp))
3571 _post_balance_at_vertices(eq, ts, "Adv", postlabel,
3572 b->advection_term);
3573
3574 if (cs_equation_param_has_time(eqp))
3575 _post_balance_at_vertices(eq, ts, "Time", postlabel,
3576 b->unsteady_term);
3577
3578 if (cs_equation_param_has_reaction(eqp))
3579 _post_balance_at_vertices(eq, ts, "Reac", postlabel,
3580 b->reaction_term);
3581
3582 if (cs_equation_param_has_sourceterm(eqp))
3583 _post_balance_at_vertices(eq, ts, "Src", postlabel,
3584 b->source_term);
3585
3586 }
3587 break;
3588
3589 default:
3590 break;
3591 }
3592
3593 /* In case of postprocessing of the border faces, one has to check if there
3594 is a mesh modification. In particular, a removal of 2D extruded border
3595 faces*/
3596
3597 bool use_parent = (mesh->n_g_b_faces_all > mesh->n_g_b_faces) ?
3598 false : true;
3599
3600 sprintf(postlabel, "%s.BdyFlux", eqp->name);
3601
3602 /* Post-process the boundary fluxes (diffusive and convective) */
3603
3604 cs_post_write_var(CS_POST_MESH_BOUNDARY,
3605 CS_POST_WRITER_DEFAULT,
3606 postlabel,
3607 1,
3608 true, /* interlace */
3609 use_parent,
3610 CS_POST_TYPE_cs_real_t,
3611 NULL, /* values on cells */
3612 NULL, /* values at internal faces */
3613 b->boundary_term, /* values at border faces */
3614 ts); /* time step structure */
3615
3616 /* Free buffers */
3617
3618 BFT_FREE(postlabel);
3619 cs_equation_balance_destroy(&b);
3620
3621 if (eq->main_ts_id > -1)
3622 cs_timer_stats_stop(eq->main_ts_id);
3623
3624 } /* Loop on equations */
3625
3626 }
3627
3628 /*----------------------------------------------------------------------------*/
3629 /*!
3630 * \brief Predefined extra-operations related to equations according to the
3631 * type of numerical scheme (for the space discretization)
3632 */
3633 /*----------------------------------------------------------------------------*/
3634
3635 void
cs_equation_extra_post(void)3636 cs_equation_extra_post(void)
3637 {
3638 for (int i = 0; i < _n_equations; i++) {
3639
3640 cs_equation_t *eq = _equations[i];
3641 assert(eq != NULL); /* Sanity check */
3642 const cs_equation_param_t *eqp = eq->param;
3643
3644 if (eq->main_ts_id > -1) /* Activate timer statistics */
3645 cs_timer_stats_start(eq->main_ts_id);
3646 assert(eq->postprocess != NULL);
3647
3648 /* Perform post-processing specific to a numerical scheme */
3649 eq->postprocess(eqp,
3650 eq->builder,
3651 eq->scheme_context);
3652
3653 if (eq->main_ts_id > -1)
3654 cs_timer_stats_stop(eq->main_ts_id);
3655
3656 } /* Loop on equations */
3657 }
3658
3659 /*----------------------------------------------------------------------------*/
3660
3661 END_C_DECLS
3662