1 /*============================================================================
2 * Functions related to the structure cs_equation_param_t storing the settings
3 * related to an equation.
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <assert.h>
35 #include <ctype.h>
36 #include <stdlib.h>
37 #include <string.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include <bft_error.h>
44 #include <bft_mem.h>
45 #include <bft_printf.h>
46
47 #include "cs_boundary_zone.h"
48 #include "cs_cdo_advection.h"
49 #include "cs_cdo_bc.h"
50 #include "cs_fp_exception.h"
51 #include "cs_hodge.h"
52 #include "cs_log.h"
53 #include "cs_mesh_location.h"
54 #include "cs_multigrid.h"
55 #include "cs_sles.h"
56 #include "cs_source_term.h"
57 #include "cs_volume_zone.h"
58
59 /*----------------------------------------------------------------------------
60 * Header for the current file
61 *----------------------------------------------------------------------------*/
62
63 #include "cs_equation_param.h"
64
65 /*----------------------------------------------------------------------------*/
66
67 BEGIN_C_DECLS
68
69 /*=============================================================================
70 * Additional doxygen documentation
71 *============================================================================*/
72
73 /*!
74 \file cs_equation_param.c
75
76 \brief Structure and functions handling the specific settings related
77 to a cs_equation_t structure
78
79 */
80
81 /*============================================================================
82 * Type definitions
83 *============================================================================*/
84
85 /*============================================================================
86 * Local private variables
87 *============================================================================*/
88
89 static const cs_real_t _weak_pena_bc_coef_by_default = 100.;
90 static const cs_real_t _strong_pena_bc_coef_by_default = 1e12;
91
92 static const char _err_empty_eqp[] =
93 N_(" Stop setting an empty cs_equation_param_t structure.\n"
94 " Please check your settings.\n");
95
96 /*============================================================================
97 * Private function prototypes
98 *============================================================================*/
99
100 /*----------------------------------------------------------------------------*/
101 /*!
102 * \brief Check if PETSc or HYPRE is available and return the possible
103 * solver class.
104 *
105 * \param[in] slesp pointer to a \ref cs_param_sles_t structure
106 * \param[in] keyname name of the key to handle
107 *
108 * \return a solver class
109 */
110 /*----------------------------------------------------------------------------*/
111
112 static cs_param_sles_class_t
_get_petsc_or_hypre(const cs_param_sles_t * slesp,const char * keyname)113 _get_petsc_or_hypre(const cs_param_sles_t *slesp,
114 const char *keyname)
115 {
116 assert(slesp != NULL);
117
118 /* Either with PETSc or with PETSc/HYPRE using Euclid */
119 cs_param_sles_class_t ret_class =
120 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_PETSC);
121
122 if (ret_class != CS_PARAM_SLES_CLASS_PETSC)
123 bft_error(__FILE__, __LINE__, 0,
124 " %s(): Eq. %s Error detected while setting \"%s\" key.\n"
125 " PETSc is not available with your installation.\n"
126 " Please check your installation settings.\n",
127 __func__, slesp->name, keyname);
128
129 if (slesp->solver_class == CS_PARAM_SLES_CLASS_HYPRE)
130 ret_class = cs_param_sles_check_class(CS_PARAM_SLES_CLASS_HYPRE);
131
132 return ret_class;
133 }
134
135 /*----------------------------------------------------------------------------*/
136 /*!
137 * \brief Set a parameter attached to a keyname in a \ref cs_equation_param_t
138 * structure
139 *
140 * \param[in, out] eqp pointer to a \ref cs_equation_param_t structure
141 * \param[in] key key related to the member of eq to set
142 * \param[in] keyval accessor to the value to set
143 */
144 /*----------------------------------------------------------------------------*/
145
146 static void
_set_key(cs_equation_param_t * eqp,cs_equation_key_t key,const char * keyval)147 _set_key(cs_equation_param_t *eqp,
148 cs_equation_key_t key,
149 const char *keyval)
150 {
151 const char *eqname = eqp->name;
152 const char emsg[] = " %s: %s equation --> Invalid key value %s for"
153 " keyword %s.\n";
154
155 switch(key) {
156
157 case CS_EQKEY_ADV_EXTRAPOL:
158 if (strcmp(keyval, "none") == 0)
159 eqp->adv_extrapol = CS_PARAM_ADVECTION_EXTRAPOL_NONE;
160 else if (strcmp(keyval, "taylor") == 0)
161 eqp->adv_formulation = CS_PARAM_ADVECTION_EXTRAPOL_TAYLOR_2;
162 else if (strcmp(keyval, "adams_bashforth") == 0)
163 eqp->adv_formulation = CS_PARAM_ADVECTION_EXTRAPOL_ADAMS_BASHFORTH_2;
164 else {
165 const char *_val = keyval;
166 bft_error(__FILE__, __LINE__, 0,
167 emsg, __func__, eqname, _val, "CS_EQKEY_ADV_EXTRAPOL");
168 }
169 break;
170
171 case CS_EQKEY_ADV_FORMULATION:
172 if (strcmp(keyval, "conservative") == 0)
173 eqp->adv_formulation = CS_PARAM_ADVECTION_FORM_CONSERV;
174 else if (strcmp(keyval, "non_conservative") == 0)
175 eqp->adv_formulation = CS_PARAM_ADVECTION_FORM_NONCONS;
176 else if (strcmp(keyval, "skew_symmetric") == 0)
177 eqp->adv_formulation = CS_PARAM_ADVECTION_FORM_SKEWSYM;
178 else {
179 const char *_val = keyval;
180 bft_error(__FILE__, __LINE__, 0,
181 emsg, __func__, eqname, _val, "CS_EQKEY_ADV_FORMULATION");
182 }
183 break;
184
185 case CS_EQKEY_ADV_SCHEME:
186 if (strcmp(keyval, "upwind") == 0)
187 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_UPWIND;
188 else if (strcmp(keyval, "samarskii") == 0)
189 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_SAMARSKII;
190 else if (strcmp(keyval, "sg") == 0)
191 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_SG;
192 else if (strcmp(keyval, "centered") == 0)
193 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_CENTERED;
194 else if (strcmp(keyval, "mix_centered_upwind") == 0 ||
195 strcmp(keyval, "hybrid_centered_upwind") == 0)
196 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_HYBRID_CENTERED_UPWIND;
197 else if (strcmp(keyval, "cip") == 0) {
198 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_CIP;
199 /* Automatically switch to a non-conservative formulation */
200 eqp->adv_formulation = CS_PARAM_ADVECTION_FORM_NONCONS;
201 }
202 else if (strcmp(keyval, "cip_cw") == 0) {
203 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_CIP_CW;
204 /* Automatically switch to a non-conservative formulation */
205 eqp->adv_formulation = CS_PARAM_ADVECTION_FORM_NONCONS;
206 }
207 else {
208 const char *_val = keyval;
209 bft_error(__FILE__, __LINE__, 0,
210 emsg, __func__, eqname, _val, "CS_EQKEY_ADV_SCHEME");
211 }
212 break;
213
214 case CS_EQKEY_ADV_STRATEGY:
215 if (strcmp(keyval, "fully_implicit") == 0 ||
216 strcmp(keyval, "implicit") == 0)
217 eqp->adv_strategy = CS_PARAM_ADVECTION_IMPLICIT_FULL;
218 else if (strcmp(keyval, "implicit_linear") == 0 ||
219 strcmp(keyval, "linearized") == 0)
220 eqp->adv_strategy = CS_PARAM_ADVECTION_IMPLICIT_LINEARIZED;
221 else if (strcmp(keyval, "explicit") == 0)
222 eqp->adv_strategy = CS_PARAM_ADVECTION_EXPLICIT;
223 else {
224 const char *_val = keyval;
225 bft_error(__FILE__, __LINE__, 0,
226 emsg, __func__, eqname, _val, "CS_EQKEY_ADV_STRATEGY");
227 }
228 break;
229
230 case CS_EQKEY_ADV_UPWIND_PORTION:
231 eqp->upwind_portion = atof(keyval);
232 /* Automatic witch to a hybrid upwind/centered scheme for advection */
233 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_HYBRID_CENTERED_UPWIND;
234 break;
235
236 case CS_EQKEY_AMG_TYPE:
237 if (strcmp(keyval, "none") == 0 || strcmp(keyval, "") == 0)
238 eqp->sles_param->amg_type = CS_PARAM_AMG_NONE;
239 else if (strcmp(keyval, "v_cycle") == 0) {
240 eqp->sles_param->amg_type = CS_PARAM_AMG_HOUSE_V;
241 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_CS;
242 }
243 else if (strcmp(keyval, "k_cycle") == 0) {
244 eqp->sles_param->amg_type = CS_PARAM_AMG_HOUSE_K;
245 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_CS;
246 }
247 else if (strcmp(keyval, "boomer") == 0 || strcmp(keyval, "boomer_v") == 0) {
248
249 cs_param_sles_class_t ret_class =
250 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_HYPRE);
251
252 if (ret_class == CS_PARAM_SLES_CLASS_HYPRE) {
253 eqp->sles_param->amg_type = CS_PARAM_AMG_HYPRE_BOOMER_V;
254 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_HYPRE;
255 }
256 else if (ret_class == CS_PARAM_SLES_CLASS_PETSC) {
257 eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_GAMG_V;
258 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
259 }
260 else
261 bft_error(__FILE__, __LINE__, 0,
262 "%s(): Eq. %s\n Invalid choice of AMG type.\n"
263 " HYPRE/PETSc are not available."
264 " Please check your settings.", __func__, eqname);
265
266 }
267 else if (strcmp(keyval, "boomer_w") == 0) {
268
269 cs_param_sles_class_t ret_class =
270 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_HYPRE);
271
272 if (ret_class == CS_PARAM_SLES_CLASS_HYPRE) {
273 eqp->sles_param->amg_type = CS_PARAM_AMG_HYPRE_BOOMER_W;
274 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_HYPRE;
275 }
276 else if (ret_class == CS_PARAM_SLES_CLASS_PETSC) {
277 eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_GAMG_W;
278 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
279 }
280 else
281 bft_error(__FILE__, __LINE__, 0,
282 "%s(): Eq. %s\n Invalid choice of AMG type.\n"
283 " HYPRE/PETSc are not available."
284 " Please check your settings.", __func__, eqname);
285
286 }
287 else if (strcmp(keyval, "gamg") == 0 || strcmp(keyval, "gamg_v") == 0) {
288
289 cs_param_sles_class_t ret_class =
290 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_PETSC);
291
292 if (ret_class != CS_PARAM_SLES_CLASS_PETSC)
293 bft_error(__FILE__, __LINE__, 0,
294 "%s(): Eq. %s\n Invalid choice of AMG type.\n"
295 " PETSc is not available."
296 " Please check your settings.", __func__, eqname);
297
298 eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_GAMG_V;
299 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
300
301 }
302 else if (strcmp(keyval, "gamg_w") == 0) {
303
304 cs_param_sles_class_t ret_class =
305 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_PETSC);
306
307 if (ret_class != CS_PARAM_SLES_CLASS_PETSC)
308 bft_error(__FILE__, __LINE__, 0,
309 "%s(): Eq. %s\n Invalid choice of AMG type.\n"
310 " PETSc is not available."
311 " Please check your settings.", __func__, eqname);
312
313 eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_GAMG_W;
314 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
315
316 }
317 else if (strcmp(keyval, "pcmg") == 0) {
318
319 cs_param_sles_class_t ret_class =
320 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_PETSC);
321
322 if (ret_class != CS_PARAM_SLES_CLASS_PETSC)
323 bft_error(__FILE__, __LINE__, 0,
324 "%s(): Eq. %s\n Invalid choice of AMG type.\n"
325 " PETSc is not available."
326 " Please check your settings.", __func__, eqname);
327
328 eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_PCMG;
329 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
330
331 }
332 else {
333 const char *_val = keyval;
334 bft_error(__FILE__, __LINE__, 0,
335 emsg, __func__, eqname, _val, "CS_EQKEY_AMG_TYPE");
336 }
337 break;
338
339 case CS_EQKEY_BC_ENFORCEMENT:
340 if (strcmp(keyval, "algebraic") == 0)
341 eqp->default_enforcement = CS_PARAM_BC_ENFORCE_ALGEBRAIC;
342 else if (strcmp(keyval, "penalization") == 0)
343 eqp->default_enforcement = CS_PARAM_BC_ENFORCE_PENALIZED;
344 else if (strcmp(keyval, "weak_sym") == 0)
345 eqp->default_enforcement = CS_PARAM_BC_ENFORCE_WEAK_SYM;
346 else if (strcmp(keyval, "weak") == 0)
347 eqp->default_enforcement = CS_PARAM_BC_ENFORCE_WEAK_NITSCHE;
348 else {
349 const char *_val = keyval;
350 bft_error(__FILE__, __LINE__, 0,
351 emsg, __func__, eqname, _val, "CS_EQKEY_BC_ENFORCEMENT");
352 }
353 break;
354
355 case CS_EQKEY_BC_QUADRATURE:
356 {
357 cs_quadrature_type_t qtype = CS_QUADRATURE_NONE;
358
359 if (strcmp(keyval, "bary") == 0)
360 qtype = CS_QUADRATURE_BARY;
361 else if (strcmp(keyval, "bary_subdiv") == 0)
362 qtype = CS_QUADRATURE_BARY_SUBDIV;
363 else if (strcmp(keyval, "higher") == 0)
364 qtype = CS_QUADRATURE_HIGHER;
365 else if (strcmp(keyval, "highest") == 0)
366 qtype = CS_QUADRATURE_HIGHEST;
367 else {
368 const char *_val = keyval;
369 bft_error(__FILE__, __LINE__, 0,
370 emsg, __func__, eqname, _val, "CS_EQKEY_BC_QUADRATURE");
371 }
372
373 for (int i = 0; i < eqp->n_bc_defs; i++)
374 cs_xdef_set_quadrature(eqp->bc_defs[i], qtype);
375
376 }
377 break;
378
379 case CS_EQKEY_BC_STRONG_PENA_COEFF:
380 eqp->strong_pena_bc_coeff = atof(keyval);
381 if (eqp->strong_pena_bc_coeff < 1.)
382 bft_error(__FILE__, __LINE__, 0,
383 " %s: Invalid value of the penalization coefficient %5.3e\n"
384 " This should be positive and large.\n"
385 " Equation: %s\n",
386 __func__, eqp->strong_pena_bc_coeff, eqname);
387 break;
388
389 case CS_EQKEY_BC_WEAK_PENA_COEFF:
390 eqp->weak_pena_bc_coeff = atof(keyval);
391 if (eqp->weak_pena_bc_coeff < 0.)
392 bft_error(__FILE__, __LINE__, 0,
393 " %s: Invalid value of the penalization coefficient %5.3e\n"
394 " This should be positive.\n"
395 " Equation: %s\n",
396 __func__, eqp->weak_pena_bc_coeff, eqname);
397 break;
398
399 case CS_EQKEY_DO_LUMPING:
400 if (strcmp(keyval, "true") == 0 || strcmp(keyval, "1") == 0)
401 eqp->do_lumping = true;
402 else
403 eqp->do_lumping = false; /* Should be the default behavior */
404 break;
405
406 case CS_EQKEY_DOF_REDUCTION:
407 if (strcmp(keyval, "derham") == 0)
408 eqp->dof_reduction = CS_PARAM_REDUCTION_DERHAM;
409 else if (strcmp(keyval, "average") == 0)
410 eqp->dof_reduction = CS_PARAM_REDUCTION_AVERAGE;
411 else {
412 const char *_val = keyval;
413 bft_error(__FILE__, __LINE__, 0,
414 emsg, __func__, eqname, _val, "CS_EQKEY_DOF_REDUCTION");
415 }
416 break;
417
418 case CS_EQKEY_EXTRA_OP:
419 if (strcmp(keyval, "balance") == 0)
420 eqp->process_flag |= CS_EQUATION_POST_BALANCE;
421 else if (strcmp(keyval, "peclet") == 0)
422 eqp->process_flag |= CS_EQUATION_POST_PECLET;
423 else if (strcmp(keyval, "upwind_coef") == 0)
424 eqp->process_flag |= CS_EQUATION_POST_UPWIND_COEF;
425 else if (strcmp(keyval, "normal_flux") == 0)
426 eqp->process_flag |= CS_EQUATION_POST_NORMAL_FLUX;
427 else {
428 const char *_val = keyval;
429 bft_error(__FILE__, __LINE__, 0,
430 emsg, __func__, eqname, _val, "CS_EQKEY_EXTRA_OP");
431 }
432 break;
433
434 case CS_EQKEY_HODGE_DIFF_ALGO:
435 if (strcmp(keyval,"cost") == 0 || strcmp(keyval,"ocs") == 0)
436 eqp->diffusion_hodgep.algo = CS_HODGE_ALGO_COST;
437 else if (strcmp(keyval, "ocs2") == 0)
438 eqp->diffusion_hodgep.algo = CS_HODGE_ALGO_OCS2;
439 else if (strcmp(keyval, "bubble") == 0)
440 eqp->diffusion_hodgep.algo = CS_HODGE_ALGO_BUBBLE;
441 else if (strcmp(keyval, "voronoi") == 0)
442 eqp->diffusion_hodgep.algo = CS_HODGE_ALGO_VORONOI;
443 else if (strcmp(keyval, "wbs") == 0)
444 eqp->diffusion_hodgep.algo = CS_HODGE_ALGO_WBS;
445 else if (strcmp(keyval, "auto") == 0)
446 eqp->diffusion_hodgep.algo = CS_HODGE_ALGO_AUTO;
447 else {
448 const char *_val = keyval;
449 bft_error(__FILE__, __LINE__, 0,
450 emsg, __func__, eqname, _val, "CS_EQKEY_HODGE_DIFF_ALGO");
451 }
452 break;
453
454 case CS_EQKEY_HODGE_DIFF_COEF:
455 if (strcmp(keyval, "dga") == 0)
456 eqp->diffusion_hodgep.coef = 1./3.;
457 else if (strcmp(keyval, "sushi") == 0)
458 eqp->diffusion_hodgep.coef = 1./sqrt(3.);
459 else if (strcmp(keyval, "gcr") == 0)
460 eqp->diffusion_hodgep.coef = 1.0;
461 else if (strcmp(keyval, "frac23") == 0 || strcmp(keyval, "2/3") == 0)
462 eqp->diffusion_hodgep.coef = 2./3.;
463 else
464 eqp->diffusion_hodgep.coef = atof(keyval);
465 break;
466
467 case CS_EQKEY_HODGE_TIME_ALGO:
468 if (strcmp(keyval, "voronoi") == 0)
469 eqp->time_hodgep.algo = CS_HODGE_ALGO_VORONOI;
470 else if (strcmp(keyval,"cost") == 0 || strcmp(keyval,"ocs") == 0)
471 eqp->time_hodgep.algo = CS_HODGE_ALGO_COST;
472 else if (strcmp(keyval, "wbs") == 0)
473 eqp->time_hodgep.algo = CS_HODGE_ALGO_WBS;
474 else {
475 const char *_val = keyval;
476 bft_error(__FILE__, __LINE__, 0,
477 emsg, __func__, eqname, _val, "CS_EQKEY_HODGE_TIME_ALGO");
478 }
479 break;
480
481 case CS_EQKEY_HODGE_REAC_ALGO:
482 if (strcmp(keyval, "voronoi") == 0)
483 eqp->reaction_hodgep.algo = CS_HODGE_ALGO_VORONOI;
484 else if (strcmp(keyval, "wbs") == 0)
485 eqp->reaction_hodgep.algo = CS_HODGE_ALGO_WBS;
486 else {
487 const char *_val = keyval;
488 bft_error(__FILE__, __LINE__, 0,
489 emsg, __func__, eqname, _val, "CS_EQKEY_HODGE_REAC_ALGO");
490 }
491 break;
492
493 case CS_EQKEY_ITSOL:
494 if (strcmp(keyval, "amg") == 0)
495 eqp->sles_param->solver = CS_PARAM_ITSOL_AMG;
496 else if (strcmp(keyval, "bicg") == 0)
497 eqp->sles_param->solver = CS_PARAM_ITSOL_BICG;
498 else if (strcmp(keyval, "bicgstab2") == 0)
499 eqp->sles_param->solver = CS_PARAM_ITSOL_BICGSTAB2;
500 else if (strcmp(keyval, "cg") == 0)
501 eqp->sles_param->solver = CS_PARAM_ITSOL_CG;
502 else if (strcmp(keyval, "cr3") == 0)
503 eqp->sles_param->solver = CS_PARAM_ITSOL_CR3;
504 else if (strcmp(keyval, "fcg") == 0)
505 eqp->sles_param->solver = CS_PARAM_ITSOL_FCG;
506 else if (strcmp(keyval, "gauss_seidel") == 0 ||
507 strcmp(keyval, "gs") == 0) {
508 eqp->sles_param->solver = CS_PARAM_ITSOL_GAUSS_SEIDEL;
509 eqp->sles_param->precond = CS_PARAM_PRECOND_NONE;
510 }
511 else if (strcmp(keyval, "gcr") == 0)
512 eqp->sles_param->solver = CS_PARAM_ITSOL_GCR;
513 else if (strcmp(keyval, "gmres") == 0)
514 eqp->sles_param->solver = CS_PARAM_ITSOL_GMRES;
515 else if (strcmp(keyval, "fgmres") == 0) {
516
517 cs_param_sles_class_t ret_class =
518 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_PETSC);
519
520 if (ret_class != CS_PARAM_SLES_CLASS_PETSC) {
521
522 cs_base_warn(__FILE__, __LINE__);
523 bft_printf(" Switch to the GCR implemantation of code_saturne\n");
524 eqp->sles_param->solver = CS_PARAM_ITSOL_GCR;
525 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_CS;
526
527 }
528 else {
529 eqp->sles_param->solver = CS_PARAM_ITSOL_FGMRES;
530 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
531 }
532
533 }
534 else if (strcmp(keyval, "jacobi") == 0 || strcmp(keyval, "diag") == 0 ||
535 strcmp(keyval, "diagonal") == 0) {
536 eqp->sles_param->solver = CS_PARAM_ITSOL_JACOBI;
537 eqp->sles_param->precond = CS_PARAM_PRECOND_NONE;
538 }
539 else if (strcmp(keyval, "minres") == 0)
540 eqp->sles_param->solver = CS_PARAM_ITSOL_MINRES;
541
542 else if (strcmp(keyval, "mumps") == 0 ||
543 strcmp(keyval, "mumps_float") == 0 ||
544 strcmp(keyval, "mumps_float_ldlt") == 0 ||
545 strcmp(keyval, "mumps_ldlt") == 0) {
546
547 eqp->sles_param->precond = CS_PARAM_PRECOND_NONE;
548
549 /* Modify the default and check availability of MUMPS solvers */
550 if (eqp->sles_param->solver_class == CS_PARAM_SLES_CLASS_CS)
551 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_MUMPS;
552
553 cs_param_sles_class_t ret_class =
554 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_MUMPS);
555
556 if (ret_class == CS_PARAM_SLES_N_CLASSES)
557 bft_error(__FILE__, __LINE__, 0,
558 " %s(): Eq. %s Error detected while setting \"%s\" key.\n"
559 " MUMPS is not available with your installation.\n"
560 " Please check your installation settings.\n",
561 __func__, eqname, "CS_EQKEY_ITSOL");
562 else
563 eqp->sles_param->solver_class = ret_class;
564
565 /* Sanity check */
566 assert(eqp->sles_param->solver_class != CS_PARAM_SLES_CLASS_CS &&
567 eqp->sles_param->solver_class != CS_PARAM_SLES_CLASS_HYPRE);
568
569 if (strcmp(keyval, "mumps") == 0)
570 eqp->sles_param->solver = CS_PARAM_ITSOL_MUMPS;
571 else if (strcmp(keyval, "mumps_float") == 0)
572 eqp->sles_param->solver = CS_PARAM_ITSOL_MUMPS_FLOAT;
573 else if (strcmp(keyval, "mumps_float_ldlt") == 0)
574 eqp->sles_param->solver = CS_PARAM_ITSOL_MUMPS_FLOAT_LDLT;
575 else if (strcmp(keyval, "mumps_ldlt") == 0)
576 eqp->sles_param->solver = CS_PARAM_ITSOL_MUMPS_LDLT;
577
578 }
579 else if (strcmp(keyval, "sym_gauss_seidel") == 0 ||
580 strcmp(keyval, "sgs") == 0) {
581 eqp->sles_param->solver = CS_PARAM_ITSOL_SYM_GAUSS_SEIDEL;
582 eqp->sles_param->precond = CS_PARAM_PRECOND_NONE;
583 }
584 else if (strcmp(keyval, "user") == 0) {
585 eqp->sles_param->solver = CS_PARAM_ITSOL_USER_DEFINED;
586 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_CS;
587 }
588 else if (strcmp(keyval, "none") == 0)
589 eqp->sles_param->solver = CS_PARAM_ITSOL_NONE;
590
591 else { /* keyval not found among the available keyvals */
592 const char *_val = keyval;
593 bft_error(__FILE__, __LINE__, 0,
594 emsg, __func__, eqname, _val, "CS_EQKEY_ITSOL");
595 }
596 break;
597
598 case CS_EQKEY_ITSOL_MAX_ITER:
599 eqp->sles_param->n_max_iter = atoi(keyval);
600 break;
601
602 case CS_EQKEY_ITSOL_EPS:
603 eqp->sles_param->eps = atof(keyval);
604 break;
605
606 case CS_EQKEY_ITSOL_RESNORM_TYPE:
607 if (strcmp(keyval, "none") == 0 || strcmp(keyval, "false") == 0 ||
608 strcmp(keyval, "") == 0)
609 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NONE;
610 else if (strcmp(keyval, "rhs") == 0)
611 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
612 else if (strcmp(keyval, "weighted_rhs") == 0 ||
613 strcmp(keyval, "weighted") == 0)
614 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_WEIGHTED_RHS;
615 else if (strcmp(keyval, "filtered_rhs") == 0 ||
616 strcmp(keyval, "filtered") == 0)
617 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_FILTERED_RHS;
618 else {
619 const char *_val = keyval;
620 bft_error(__FILE__, __LINE__, 0,
621 emsg, __func__, eqname, _val, "CS_EQKEY_ITSOL_RESNORM_TYPE");
622 }
623 break;
624
625 case CS_EQKEY_ITSOL_RESTART:
626 eqp->sles_param->restart = atoi(keyval);
627 break;
628
629 case CS_EQKEY_OMP_ASSEMBLY_STRATEGY:
630 if (strcmp(keyval, "critical") == 0)
631 eqp->omp_assembly_choice = CS_PARAM_ASSEMBLE_OMP_CRITICAL;
632 else if (strcmp(keyval, "atomic") == 0)
633 eqp->omp_assembly_choice = CS_PARAM_ASSEMBLE_OMP_ATOMIC;
634 else {
635 const char *_val = keyval;
636 bft_error(__FILE__, __LINE__, 0,
637 emsg, __func__, eqname, _val, "CS_EQKEY_OMP_ASSEMBLY_STRATEGY");
638 }
639 break;
640
641 case CS_EQKEY_PRECOND:
642 if (strcmp(keyval, "none") == 0) {
643 eqp->sles_param->precond = CS_PARAM_PRECOND_NONE;
644 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_NONE;
645 eqp->sles_param->amg_type = CS_PARAM_AMG_NONE;
646 }
647 else if (strcmp(keyval, "jacobi") == 0 || strcmp(keyval, "diag") == 0)
648 eqp->sles_param->precond = CS_PARAM_PRECOND_DIAG;
649
650 else if (strcmp(keyval, "block_jacobi") == 0 ||
651 strcmp(keyval, "bjacobi") == 0) {
652
653 if (eqp->sles_param->pcd_block_type == CS_PARAM_PRECOND_BLOCK_NONE)
654 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_DIAG;
655
656 /* Either with PETSc or with PETSc/HYPRE using Euclid */
657 eqp->sles_param->solver_class = _get_petsc_or_hypre(eqp->sles_param,
658 "CS_EQKEY_PRECOND");
659
660 eqp->sles_param->precond = CS_PARAM_PRECOND_BJACOB_ILU0;
661 /* Default when using PETSc */
662 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
663
664 }
665 else if (strcmp(keyval, "bjacobi_sgs") == 0 ||
666 strcmp(keyval, "bjacobi_ssor") == 0) {
667
668 if (eqp->sles_param->pcd_block_type == CS_PARAM_PRECOND_BLOCK_NONE)
669 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_DIAG;
670
671 cs_param_sles_class_t ret_class =
672 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_PETSC);
673
674 if (ret_class != CS_PARAM_SLES_CLASS_PETSC)
675 bft_error(__FILE__, __LINE__, 0,
676 " %s(): Eq. %s Error detected while setting \"%s\" key.\n"
677 " PETSc is not available with your installation.\n"
678 " Please check your installation settings.\n",
679 __func__, eqname, "CS_EQKEY_PRECOND");
680
681 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
682 eqp->sles_param->precond = CS_PARAM_PRECOND_BJACOB_SGS;
683 /* Default when using PETSc */
684 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
685
686 }
687 else if (strcmp(keyval, "lu") == 0) {
688
689 eqp->sles_param->precond = CS_PARAM_PRECOND_LU;
690
691 cs_param_sles_class_t ret_class =
692 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_PETSC);
693
694 if (ret_class != CS_PARAM_SLES_CLASS_PETSC)
695 bft_error(__FILE__, __LINE__, 0,
696 " %s(): Eq. %s Error detected while setting \"%s\" key.\n"
697 " PETSc is not available with your installation.\n"
698 " Please check your installation settings.\n",
699 __func__, eqname, "CS_EQKEY_PRECOND");
700
701 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
702
703 /* Default when using PETSc */
704 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
705
706 }
707 else if (strcmp(keyval, "ilu0") == 0) {
708
709 eqp->sles_param->precond = CS_PARAM_PRECOND_ILU0;
710
711 /* Either with PETSc or with PETSc/HYPRE using Euclid */
712 eqp->sles_param->solver_class = _get_petsc_or_hypre(eqp->sles_param,
713 "CS_EQKEY_PRECOND");
714
715 /* Default when using PETSc */
716 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
717
718 }
719 else if (strcmp(keyval, "icc0") == 0) {
720
721 eqp->sles_param->precond = CS_PARAM_PRECOND_ICC0;
722
723 /* Either with PETSc or with PETSc/HYPRE using Euclid */
724 eqp->sles_param->solver_class = _get_petsc_or_hypre(eqp->sles_param,
725 "CS_EQKEY_PRECOND");
726
727 /* Default when using PETSc */
728 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
729
730 }
731 else if (strcmp(keyval, "amg") == 0) {
732
733 eqp->sles_param->precond = CS_PARAM_PRECOND_AMG;
734
735 cs_param_sles_class_t ret_class =
736 cs_param_sles_check_class(eqp->sles_param->solver_class);
737
738 /* Set the default AMG choice according to the class of solver */
739 switch (ret_class) {
740
741 case CS_PARAM_SLES_CLASS_CS:
742 eqp->sles_param->amg_type = CS_PARAM_AMG_HOUSE_K;
743 break;
744 case CS_PARAM_SLES_CLASS_PETSC:
745 eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_GAMG_V;
746 /* Default when using PETSc */
747 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
748 break;
749 case CS_PARAM_SLES_CLASS_HYPRE:
750 eqp->sles_param->amg_type = CS_PARAM_AMG_HYPRE_BOOMER_V;
751 /* Up to now HYPRE is available only through the PETSc interface.
752 Default when using PETSc */
753 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
754 break;
755
756 default:
757 bft_error(__FILE__, __LINE__, 0,
758 "%s(): Eq. %s Invalid choice of AMG type.\n"
759 " Please modify your settings", __func__, eqname);
760
761 } /* End of switch */
762
763 }
764 else if (strcmp(keyval, "amg_block") == 0 ||
765 strcmp(keyval, "block_amg") == 0) {
766
767 eqp->sles_param->precond = CS_PARAM_PRECOND_AMG;
768 if (eqp->sles_param->pcd_block_type == CS_PARAM_PRECOND_BLOCK_NONE)
769 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_DIAG;
770
771 cs_param_sles_class_t ret_class =
772 cs_param_sles_check_class(eqp->sles_param->solver_class);
773
774 /* Set the default AMG choice according to the class of solver */
775 switch (ret_class) {
776
777 case CS_PARAM_SLES_CLASS_CS:
778 eqp->sles_param->amg_type = CS_PARAM_AMG_HOUSE_K;
779 break;
780 case CS_PARAM_SLES_CLASS_PETSC:
781 eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_GAMG_V;
782 /* Default when using PETSc */
783 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
784 break;
785 case CS_PARAM_SLES_CLASS_HYPRE:
786 eqp->sles_param->amg_type = CS_PARAM_AMG_HYPRE_BOOMER_V;
787 /* Up to now HYPRE is available only through the PETSc interface.
788 Default when using PETSc */
789 eqp->sles_param->resnorm_type = CS_PARAM_RESNORM_NORM2_RHS;
790 break;
791
792 default:
793 bft_error(__FILE__, __LINE__, 0,
794 "%s(): Eq. %s Invalid choice of block AMG type.\n"
795 " Please modify your settings", __func__, eqname);
796
797 } /* End of switch */
798
799 }
800 else if (strcmp(keyval, "poly1") == 0) {
801 eqp->sles_param->precond = CS_PARAM_PRECOND_POLY1;
802 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_CS;
803 }
804 else if (strcmp(keyval, "poly2") == 0) {
805 eqp->sles_param->precond = CS_PARAM_PRECOND_POLY2;
806 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_CS;
807 }
808 else if (strcmp(keyval, "ssor") == 0)
809 eqp->sles_param->precond = CS_PARAM_PRECOND_SSOR;
810 else {
811 const char *_val = keyval;
812 bft_error(__FILE__, __LINE__, 0,
813 emsg, __func__, eqname, _val, "CS_EQKEY_PRECOND");
814 }
815 break; /* Precond */
816
817 case CS_EQKEY_PRECOND_BLOCK_TYPE:
818 if (strcmp(keyval, "none") == 0)
819 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_NONE;
820 else if (strcmp(keyval, "diag") == 0)
821 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_DIAG;
822 else if (strcmp(keyval, "full_diag") == 0)
823 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_FULL_DIAG;
824 else if (strcmp(keyval, "full_lower") == 0)
825 eqp->sles_param->pcd_block_type =
826 CS_PARAM_PRECOND_BLOCK_FULL_LOWER_TRIANGULAR;
827 else if (strcmp(keyval, "full_symm") == 0)
828 eqp->sles_param->pcd_block_type =
829 CS_PARAM_PRECOND_BLOCK_FULL_SYM_GAUSS_SEIDEL;
830 else if (strcmp(keyval, "full_upper") == 0)
831 eqp->sles_param->pcd_block_type =
832 CS_PARAM_PRECOND_BLOCK_FULL_UPPER_TRIANGULAR;
833 else if (strcmp(keyval, "lower") == 0)
834 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_LOWER_TRIANGULAR;
835 else if (strcmp(keyval, "symm") == 0)
836 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_SYM_GAUSS_SEIDEL;
837 else if (strcmp(keyval, "upper") == 0)
838 eqp->sles_param->pcd_block_type = CS_PARAM_PRECOND_BLOCK_UPPER_TRIANGULAR;
839 else {
840 const char *_val = keyval;
841 bft_error(__FILE__, __LINE__, 0,
842 emsg, __func__, eqname, _val, "CS_EQKEY_PRECOND_BLOCK_TYPE");
843 }
844 break;
845
846 case CS_EQKEY_SLES_VERBOSITY: /* "verbosity" for SLES structures */
847 eqp->sles_param->verbosity = atoi(keyval);
848 break;
849
850 case CS_EQKEY_SOLVER_FAMILY:
851 if (strcmp(keyval, "cs") == 0)
852 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_CS;
853
854 else if (strcmp(keyval, "hypre") == 0) {
855
856 cs_param_sles_class_t ret_class =
857 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_HYPRE);
858
859 if (ret_class == CS_PARAM_SLES_N_CLASSES)
860 bft_error(__FILE__, __LINE__, 0,
861 " %s(): Eq. %s Error detected while setting \"%s\" key.\n"
862 " Neither PETSc nor HYPRE is available.\n"
863 " Please check your installation settings.\n",
864 __func__, eqname, "CS_EQKEY_SOLVER_FAMILY");
865 else if (ret_class == CS_PARAM_SLES_CLASS_PETSC)
866 bft_error(__FILE__, __LINE__, 0,
867 " %s(): Eq. %s Error detected while setting \"%s\" key.\n"
868 " PETSc with HYPRE is not available.\n"
869 " Please check your installation settings.\n",
870 __func__, eqname, "CS_EQKEY_SOLVER_FAMILY");
871
872 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_HYPRE;
873
874 }
875 else if (strcmp(keyval, "mumps") == 0) {
876
877 cs_param_sles_class_t ret_class =
878 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_MUMPS);
879
880 if (ret_class == CS_PARAM_SLES_N_CLASSES)
881 bft_error(__FILE__, __LINE__, 0,
882 " %s(): Eq. %s Error detected while setting \"%s\" key.\n"
883 " MUMPS is not available.\n"
884 " Please check your installation settings.\n",
885 __func__, eqname, "CS_EQKEY_SOLVER_FAMILY");
886
887 eqp->sles_param->solver_class = ret_class; /* PETSc or MUMPS */
888
889 }
890 else if (strcmp(keyval, "petsc") == 0) {
891
892 cs_param_sles_class_t ret_class =
893 cs_param_sles_check_class(CS_PARAM_SLES_CLASS_PETSC);
894
895 if (ret_class == CS_PARAM_SLES_N_CLASSES)
896 bft_error(__FILE__, __LINE__, 0,
897 " %s(): Eq. %s Error detected while setting \"%s\" key.\n"
898 " PETSc is not available.\n"
899 " Please check your installation settings.\n",
900 __func__, eqname, "CS_EQKEY_SOLVER_FAMILY");
901
902 eqp->sles_param->solver_class = CS_PARAM_SLES_CLASS_PETSC;
903
904 }
905 else {
906 const char *_val = keyval;
907 bft_error(__FILE__, __LINE__, 0,
908 emsg, __func__, eqname, _val, "CS_EQKEY_SOLVER_FAMILY");
909 }
910 break;
911
912 case CS_EQKEY_SPACE_SCHEME:
913 if (strcmp(keyval, "cdo_vb") == 0 ||
914 strcmp(keyval, "cdovb") == 0) {
915 eqp->space_scheme = CS_SPACE_SCHEME_CDOVB;
916 eqp->space_poly_degree = 0;
917 eqp->time_hodgep.type = CS_HODGE_TYPE_VPCD;
918 eqp->diffusion_hodgep.type = CS_HODGE_TYPE_EPFD;
919 eqp->diffusion_hodgep.algo = CS_HODGE_ALGO_BUBBLE;
920 eqp->diffusion_hodgep.coef = 2*cs_math_1ov3;
921 eqp->reaction_hodgep.type = CS_HODGE_TYPE_VPCD;
922 eqp->reaction_hodgep.algo = CS_HODGE_ALGO_WBS;
923 }
924 else if (strcmp(keyval, "cdo_vcb") == 0 ||
925 strcmp(keyval, "cdovcb") == 0) {
926 eqp->space_scheme = CS_SPACE_SCHEME_CDOVCB;
927 eqp->space_poly_degree = 0;
928 eqp->time_hodgep.type = CS_HODGE_TYPE_VPCD;
929 eqp->diffusion_hodgep.algo = CS_HODGE_ALGO_WBS;
930 eqp->diffusion_hodgep.type = CS_HODGE_TYPE_VC;
931 eqp->reaction_hodgep.type = CS_HODGE_TYPE_VPCD;
932 eqp->reaction_hodgep.algo = CS_HODGE_ALGO_WBS;
933 }
934 else if (strcmp(keyval, "cdo_fb") == 0 ||
935 strcmp(keyval, "cdofb") == 0) {
936 eqp->space_scheme = CS_SPACE_SCHEME_CDOFB;
937 eqp->space_poly_degree = 0;
938 eqp->time_hodgep.type = CS_HODGE_TYPE_CPVD;
939 eqp->time_hodgep.algo = CS_HODGE_ALGO_VORONOI;
940 eqp->reaction_hodgep.algo = CS_HODGE_ALGO_VORONOI;
941 eqp->diffusion_hodgep.type = CS_HODGE_TYPE_EDFP;
942 }
943 else if (strcmp(keyval, "cdo_eb") == 0 ||
944 strcmp(keyval, "cdoeb") == 0) {
945 eqp->space_scheme = CS_SPACE_SCHEME_CDOEB;
946 eqp->space_poly_degree = 0;
947 eqp->time_hodgep.type = CS_HODGE_TYPE_EPFD;
948 eqp->diffusion_hodgep.type = CS_HODGE_TYPE_FPED;
949 eqp->reaction_hodgep.type = CS_HODGE_TYPE_EPFD;
950 }
951
952 /* Only diffusion is implemented for HHO schemes up to now */
953 else if (strcmp(keyval, "hho_p0") == 0) {
954 eqp->space_scheme = CS_SPACE_SCHEME_HHO_P0;
955 eqp->space_poly_degree = 0;
956 eqp->time_hodgep.type = CS_HODGE_TYPE_CPVD;
957 eqp->diffusion_hodgep.type = CS_HODGE_TYPE_EDFP;
958 }
959 else if (strcmp(keyval, "hho_p1") == 0) {
960 eqp->space_scheme = CS_SPACE_SCHEME_HHO_P1;
961 eqp->space_poly_degree = 1;
962 eqp->time_hodgep.type = CS_HODGE_TYPE_CPVD;
963 eqp->diffusion_hodgep.type = CS_HODGE_TYPE_EDFP;
964 }
965 else if (strcmp(keyval, "hho_p2") == 0) {
966 eqp->space_scheme = CS_SPACE_SCHEME_HHO_P2;
967 eqp->space_poly_degree = 2;
968 eqp->time_hodgep.type = CS_HODGE_TYPE_CPVD;
969 eqp->diffusion_hodgep.type = CS_HODGE_TYPE_EDFP;
970 }
971 else {
972 const char *_val = keyval;
973 bft_error(__FILE__, __LINE__, 0,
974 emsg, __func__, eqname, _val, "CS_EQKEY_SPACE_SCHEME");
975 }
976 break;
977
978 case CS_EQKEY_TIME_SCHEME:
979
980 if (strcmp(keyval, "no") == 0 || strcmp(keyval, "steady") == 0) {
981 eqp->time_scheme = CS_TIME_SCHEME_STEADY;
982 }
983 else if (strcmp(keyval, "euler_implicit") == 0 ||
984 strcmp(keyval, "forward_euler") == 0) {
985 eqp->time_scheme = CS_TIME_SCHEME_EULER_IMPLICIT;
986 eqp->theta = 1.;
987 }
988 else if (strcmp(keyval, "euler_explicit") == 0 ||
989 strcmp(keyval, "backward_euler") == 0) {
990 eqp->time_scheme = CS_TIME_SCHEME_EULER_EXPLICIT;
991 eqp->theta = 0.;
992 }
993 else if (strcmp(keyval, "crank_nicolson") == 0) {
994 eqp->time_scheme = CS_TIME_SCHEME_CRANKNICO;
995 eqp->theta = 0.5;
996 }
997 else if (strcmp(keyval, "theta_scheme") == 0)
998 eqp->time_scheme = CS_TIME_SCHEME_THETA;
999 else if (strcmp(keyval, "bdf2") == 0) {
1000 eqp->time_scheme = CS_TIME_SCHEME_BDF2;
1001 bft_error(__FILE__, __LINE__, 0, " Eq. %s\n"
1002 " Soon available...", eqname);
1003 }
1004 else {
1005 const char *_val = keyval;
1006 bft_error(__FILE__, __LINE__, 0,
1007 emsg, __func__, eqname, _val, "CS_EQKEY_TIME_SCHEME");
1008 }
1009 break;
1010
1011 case CS_EQKEY_TIME_THETA:
1012 eqp->theta = atof(keyval);
1013 break;
1014
1015 case CS_EQKEY_VERBOSITY: /* "verbosity" */
1016 eqp->verbosity = atoi(keyval);
1017 break;
1018
1019 default:
1020 bft_error(__FILE__, __LINE__, 0,
1021 _(" %s: Invalid key for setting the equation %s."),
1022 __func__, eqname);
1023
1024 } /* Switch on keys */
1025
1026 }
1027
1028 /*============================================================================
1029 * Public function prototypes
1030 *============================================================================*/
1031
1032 /*----------------------------------------------------------------------------*/
1033 /*!
1034 * \brief Create a \ref cs_equation_param_t structure with the given
1035 * parameters. The remaining parameters are set with default values;
1036 *
1037 * \param[in] name name of the equation
1038 * \param[in] type type of equation
1039 * \param[in] dim dim of the variable associated to this equation
1040 * \param[in] default_bc type of boundary condition set by default
1041 *
1042 * \return a pointer to a new allocated \ref cs_equation_param_t structure
1043 */
1044 /*----------------------------------------------------------------------------*/
1045
1046 cs_equation_param_t *
cs_equation_param_create(const char * name,cs_equation_type_t type,int dim,cs_param_bc_type_t default_bc)1047 cs_equation_param_create(const char *name,
1048 cs_equation_type_t type,
1049 int dim,
1050 cs_param_bc_type_t default_bc)
1051 {
1052 cs_equation_param_t *eqp = NULL;
1053
1054 BFT_MALLOC(eqp, 1, cs_equation_param_t);
1055
1056 /* Store the name of the equation */
1057
1058 size_t len = strlen(name);
1059 BFT_MALLOC(eqp->name, len + 1, char);
1060 strncpy(eqp->name, name, len);
1061 eqp->name[len] = '\0';
1062
1063 /* Set additional members */
1064
1065 eqp->type = type;
1066 eqp->dim = dim;
1067
1068 /* Other default settings */
1069
1070 eqp->verbosity = 2;
1071 eqp->flag = 0;
1072 eqp->process_flag = 0;
1073
1074 /* Vertex-based schemes imply specific discrete Hodge operators for
1075 diffusion, time and reaction terms.
1076 Default initialization is made in accordance with this choice */
1077
1078 eqp->space_scheme = CS_SPACE_SCHEME_CDOVB;
1079 eqp->dof_reduction = CS_PARAM_REDUCTION_DERHAM;
1080 eqp->space_poly_degree = 0;
1081
1082 /* Default initialization for the legacy var_col_opt structure which is now
1083 * shared inside the cs_equation_param_t structure The default value used
1084 * here should be the same as the one considered in src/base/cs_parameters.c
1085 * (static variable called _equation_param_default) */
1086
1087 eqp->iconv = 1;
1088 eqp->istat = 1;
1089 eqp->idircl = 1;
1090 eqp->ndircl = 0;
1091 eqp->idiff = 1;
1092 eqp->idifft = 1;
1093 eqp->idften = CS_ISOTROPIC_DIFFUSION;
1094 eqp->iswdyn = -1;
1095 eqp->ischcv = 1;
1096 eqp->ibdtso = 1;
1097 eqp->isstpc = 1;
1098 eqp->nswrgr = 100;
1099 eqp->nswrsm = 1;
1100 eqp->imvisf = 0;
1101 eqp->imrgra = -1;
1102 eqp->imligr = -1;
1103 eqp->ircflu = 1;
1104 eqp->iwgrec = 0;
1105 eqp->icoupl = -1;
1106 eqp->thetav = 1.;
1107 eqp->blencv = 1.;
1108 eqp->blend_st = 0.;
1109 eqp->epsilo = 1.e-5;
1110 eqp->epsrsm = 1.e-4;
1111 eqp->epsrgr = 1.e-4;
1112 eqp->climgr = 1.5;
1113 eqp->relaxv = 1.;
1114
1115 /* Boundary conditions structure.
1116 One assigns a boundary condition by default */
1117
1118 eqp->default_bc = default_bc;
1119 eqp->n_bc_defs = 0;
1120 eqp->bc_defs = NULL;
1121 eqp->default_enforcement = CS_PARAM_BC_ENFORCE_ALGEBRAIC;
1122 eqp->strong_pena_bc_coeff = _strong_pena_bc_coef_by_default;
1123 eqp->weak_pena_bc_coeff = _weak_pena_bc_coef_by_default;
1124
1125 /* Initial condition (zero value by default) */
1126
1127 eqp->n_ic_defs = 0;
1128 eqp->ic_defs = NULL;
1129
1130 /* Description of the time discretization (default values) */
1131
1132 eqp->time_property = NULL;
1133 eqp->time_scheme = CS_TIME_SCHEME_EULER_IMPLICIT;
1134 eqp->theta = 1.0;
1135 eqp->do_lumping = false;
1136 eqp->time_hodgep = (cs_hodge_param_t) {
1137 .inv_pty = false,
1138 .algo = CS_HODGE_ALGO_VORONOI,
1139 .type = CS_HODGE_TYPE_VPCD,
1140 .coef = 1.,
1141 };
1142
1143 /* Description of the discetization of the diffusion term */
1144
1145 eqp->diffusion_property = NULL;
1146 eqp->diffusion_hodgep = (cs_hodge_param_t) {
1147 .inv_pty = false,
1148 .algo = CS_HODGE_ALGO_COST,
1149 .type = CS_HODGE_TYPE_EPFD,
1150 .coef = 1./3.,
1151 };
1152
1153 /* Description of the discetization of the curl-curl term */
1154
1155 eqp->curlcurl_property = NULL;
1156 eqp->curlcurl_hodgep = (cs_hodge_param_t) {
1157 .inv_pty = true,
1158 .algo = CS_HODGE_ALGO_COST,
1159 .type = CS_HODGE_TYPE_FPED,
1160 .coef = 1./3.,
1161 };
1162
1163 /* Description of the discetization of the grad-div term */
1164
1165 eqp->graddiv_property = NULL;
1166 eqp->graddiv_hodgep = (cs_hodge_param_t) {
1167 .inv_pty = false,
1168 .algo = CS_HODGE_ALGO_VORONOI,
1169 .type = CS_HODGE_TYPE_EPFD,
1170 .coef = 1./3.,
1171 };
1172
1173 /* Advection term */
1174
1175 eqp->adv_field = NULL;
1176 eqp->adv_scaling_property = NULL;
1177 eqp->adv_extrapol = CS_PARAM_ADVECTION_EXTRAPOL_NONE;
1178 eqp->adv_formulation = CS_PARAM_ADVECTION_FORM_CONSERV;
1179 eqp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_UPWIND;
1180 eqp->adv_strategy = CS_PARAM_ADVECTION_IMPLICIT_FULL;
1181 eqp->upwind_portion = 0.15;
1182
1183 /* Description of the discretization of the reaction term.
1184 No reaction term by default */
1185
1186 eqp->n_reaction_terms = 0;
1187 eqp->reaction_properties = NULL;
1188 eqp->reaction_hodgep = (cs_hodge_param_t) {
1189 .inv_pty = false,
1190 .algo = CS_HODGE_ALGO_VORONOI,
1191 .type = CS_HODGE_TYPE_VPCD,
1192 };
1193
1194 /* Source term (always in the right-hand side)
1195 No source term by default */
1196
1197 eqp->n_source_terms = 0;
1198 eqp->source_terms = NULL;
1199
1200 /* Mass injection in the volume term (always in the right-hand side)
1201 No volume mass injection term by default */
1202
1203 eqp->n_volume_mass_injections = 0;
1204 eqp->volume_mass_injections = NULL;
1205
1206 /* Members of the structure handling the enforcement of (internal) DoFs */
1207
1208 eqp->n_enforcements = 0;
1209 eqp->enforcement_params = NULL;
1210
1211 /* Settings for driving the linear algebra */
1212
1213 eqp->sles_param = cs_param_sles_create(-1, name); /* field_id, system_name */
1214
1215 /* Settings for the OpenMP strategy */
1216
1217 eqp->omp_assembly_choice = CS_PARAM_ASSEMBLE_OMP_CRITICAL;
1218
1219 return eqp;
1220 }
1221
1222 /*----------------------------------------------------------------------------*/
1223 /*!
1224 * \brief Copy the settings from one \ref cs_equation_param_t structure to
1225 * another one. The name is not copied.
1226 *
1227 * \param[in] ref pointer to the reference \ref cs_equation_param_t
1228 * \param[in, out] dst pointer to the \ref cs_equation_param_t to update
1229 * \param[in] copy_fid copy also the field id or not
1230 */
1231 /*----------------------------------------------------------------------------*/
1232
1233 void
cs_equation_param_copy_from(const cs_equation_param_t * ref,cs_equation_param_t * dst,bool copy_fid)1234 cs_equation_param_copy_from(const cs_equation_param_t *ref,
1235 cs_equation_param_t *dst,
1236 bool copy_fid)
1237 {
1238 /* Generic members */
1239
1240 dst->type = ref->type;
1241 dst->dim = ref->dim;
1242 dst->verbosity = ref->verbosity;
1243 dst->process_flag = ref->process_flag;
1244 dst->flag = ref->flag;
1245 dst->space_scheme = ref->space_scheme;
1246 dst->dof_reduction = ref->dof_reduction;
1247 dst->space_poly_degree = ref->space_poly_degree;
1248
1249 /* Members originally located in the cs_var_cal_opt_t structure */
1250
1251 dst->iconv = ref->iconv;
1252 dst->istat = ref->istat;
1253 dst->idircl = ref->idircl;
1254 dst->ndircl = ref->ndircl;
1255 dst->idiff = ref->idiff;
1256 dst->idifft = ref->idifft;
1257 dst->idften = ref->idften;
1258 dst->iswdyn = ref->iswdyn;
1259 dst->ischcv = ref->ischcv;
1260 dst->ibdtso = ref->ibdtso;
1261 dst->isstpc = ref->isstpc;
1262 dst->nswrgr = ref->nswrgr;
1263 dst->nswrsm = ref->nswrsm;
1264 dst->imvisf = ref->imvisf;
1265 dst->imrgra = ref->imrgra;
1266 dst->imligr = ref->imligr;
1267 dst->ircflu = ref->ircflu;
1268 dst->iwgrec = ref->iwgrec;
1269 dst->icoupl = ref->icoupl;
1270 dst->thetav = ref->thetav;
1271 dst->blencv = ref->blencv;
1272 dst->blend_st = ref->blend_st;
1273 dst->epsilo = ref->epsilo;
1274 dst->epsrsm = ref->epsrsm;
1275 dst->epsrgr = ref->epsrgr;
1276 dst->climgr = ref->climgr;
1277 dst->relaxv = ref->relaxv;
1278
1279 /* Boundary conditions structure */
1280
1281 dst->default_bc = ref->default_bc;
1282 dst->default_enforcement = ref->default_enforcement;
1283 dst->strong_pena_bc_coeff = ref->strong_pena_bc_coeff;
1284 dst->n_bc_defs = ref->n_bc_defs;
1285 BFT_MALLOC(dst->bc_defs, dst->n_bc_defs, cs_xdef_t *);
1286 for (int i = 0; i < ref->n_bc_defs; i++)
1287 dst->bc_defs[i] = cs_xdef_copy(ref->bc_defs[i]);
1288
1289 /* Description of the time discretization */
1290
1291 dst->time_scheme = ref->time_scheme;
1292 dst->theta = ref->theta;
1293 dst->do_lumping = ref->do_lumping;
1294 dst->time_property = ref->time_property;
1295
1296 cs_hodge_copy_parameters(&(ref->time_hodgep), &(dst->time_hodgep));
1297
1298 /* Initial condition (zero value by default) */
1299
1300 dst->n_ic_defs = ref->n_ic_defs;
1301 BFT_MALLOC(dst->ic_defs, dst->n_ic_defs, cs_xdef_t *);
1302 for (int i = 0; i < ref->n_ic_defs; i++)
1303 dst->ic_defs[i] = cs_xdef_copy(ref->ic_defs[i]);
1304
1305 /* Diffusion term */
1306
1307 dst->diffusion_property = ref->diffusion_property;
1308
1309 cs_hodge_copy_parameters(&(ref->diffusion_hodgep), &(dst->diffusion_hodgep));
1310
1311 /* Curl-curl term */
1312
1313 dst->curlcurl_property = ref->curlcurl_property;
1314
1315 cs_hodge_copy_parameters(&(ref->curlcurl_hodgep), &(dst->curlcurl_hodgep));
1316
1317 /* Grad-div term */
1318
1319 dst->graddiv_property = ref->graddiv_property;
1320
1321 cs_hodge_copy_parameters(&(ref->graddiv_hodgep), &(dst->graddiv_hodgep));
1322
1323 /* Advection term */
1324
1325 dst->adv_extrapol = ref->adv_extrapol;
1326 dst->adv_formulation = ref->adv_formulation;
1327 dst->adv_scheme = ref->adv_scheme;
1328 dst->adv_strategy = ref->adv_strategy;
1329 dst->upwind_portion = ref->upwind_portion;
1330 dst->adv_field = ref->adv_field;
1331 dst->adv_scaling_property = ref->adv_scaling_property;
1332
1333 /* Reaction term */
1334
1335 dst->n_reaction_terms = ref->n_reaction_terms;
1336 BFT_MALLOC(dst->reaction_properties, dst->n_reaction_terms, cs_property_t *);
1337 for (int i = 0; i < ref->n_reaction_terms; i++)
1338 dst->reaction_properties[i] = ref->reaction_properties[i];
1339
1340 cs_hodge_copy_parameters(&(ref->reaction_hodgep), &(dst->reaction_hodgep));
1341
1342 /* Source term */
1343
1344 dst->n_source_terms = ref->n_source_terms;
1345 BFT_MALLOC(dst->source_terms, dst->n_source_terms, cs_xdef_t *);
1346 for (int i = 0; i < dst->n_source_terms; i++)
1347 dst->source_terms[i] = cs_xdef_copy(ref->source_terms[i]);
1348
1349 /* Mass injection term */
1350
1351 dst->n_volume_mass_injections = ref->n_volume_mass_injections;
1352 BFT_MALLOC(dst->volume_mass_injections,
1353 dst->n_volume_mass_injections,
1354 cs_xdef_t *);
1355 for (int i = 0; i < dst->n_volume_mass_injections; i++)
1356 dst->volume_mass_injections[i]
1357 = cs_xdef_copy(ref->volume_mass_injections[i]);
1358
1359 /* Enforcement of internal DoFs */
1360
1361 if (ref->n_enforcements > 0) {
1362
1363 dst->n_enforcements = ref->n_enforcements;
1364 BFT_MALLOC(dst->enforcement_params, dst->n_enforcements,
1365 cs_enforcement_param_t *);
1366
1367 for (int i = 0; i < ref->n_enforcements; i++)
1368 dst->enforcement_params[i] =
1369 cs_enforcement_param_copy(ref->enforcement_params[i]);
1370
1371 }
1372
1373 /* Copy the settings driving the linear algebra algorithms */
1374
1375 if (copy_fid)
1376 cs_param_sles_copy_from(ref->sles_param, dst->sles_param);
1377
1378 else {
1379
1380 int save_field_id = dst->sles_param->field_id;
1381 cs_param_sles_copy_from(ref->sles_param, dst->sles_param);
1382 dst->sles_param->field_id = save_field_id;
1383
1384 }
1385
1386 /* Settings related to the performance */
1387
1388 dst->omp_assembly_choice = ref->omp_assembly_choice;
1389 }
1390
1391 /*----------------------------------------------------------------------------*/
1392 /*!
1393 * \brief Free the contents of a \ref cs_equation_param_t
1394 *
1395 * The cs_equation_param_t structure itself is not freed, but all its
1396 * sub-structures are freed.
1397 *
1398 * This is useful for equation parameters which are accessed through
1399 * field keywords.
1400 *
1401 * \param[in, out] eqp pointer to a \ref cs_equation_param_t
1402 */
1403 /*----------------------------------------------------------------------------*/
1404
1405 void
cs_equation_param_clear(cs_equation_param_t * eqp)1406 cs_equation_param_clear(cs_equation_param_t *eqp)
1407 {
1408 if (eqp == NULL)
1409 return;
1410
1411 /* Information related to the definition of the boundary conditions */
1412
1413 if (eqp->n_bc_defs > 0) {
1414
1415 for (int i = 0; i < eqp->n_bc_defs; i++)
1416 eqp->bc_defs[i] = cs_xdef_free(eqp->bc_defs[i]);
1417 BFT_FREE(eqp->bc_defs);
1418
1419 }
1420
1421 /* Information related to the definition of reaction terms */
1422
1423 if (eqp->n_reaction_terms > 0) {
1424
1425 BFT_FREE(eqp->reaction_properties);
1426
1427 /* Remark: properties are freed when the global cs_domain_t structure is
1428 freed thanks to a call to cs_property_free() */
1429 }
1430
1431 /* Information related to the definition of source terms */
1432
1433 if (eqp->n_source_terms > 0) {
1434
1435 for (int i = 0; i < eqp->n_source_terms; i++)
1436 eqp->source_terms[i] = cs_xdef_free(eqp->source_terms[i]);
1437 BFT_FREE(eqp->source_terms);
1438
1439 }
1440
1441 /* Information related to the definition of mass injection terms */
1442
1443 if (eqp->n_volume_mass_injections > 0) {
1444
1445 for (int i = 0; i < eqp->n_volume_mass_injections; i++)
1446 eqp->volume_mass_injections[i]
1447 = cs_xdef_free(eqp->volume_mass_injections[i]);
1448 BFT_FREE(eqp->volume_mass_injections);
1449
1450 }
1451
1452 /* Information related to the enforcement of internal DoFs */
1453
1454 if (eqp->n_enforcements > 0) {
1455
1456 for (int i = 0; i < eqp->n_enforcements; i++)
1457 cs_enforcement_param_free(&(eqp->enforcement_params[i]));
1458 BFT_FREE(eqp->enforcement_params);
1459
1460 }
1461
1462 /* Information related to the definition of initial conditions */
1463
1464 if (eqp->n_ic_defs > 0) {
1465
1466 for (int i = 0; i < eqp->n_ic_defs; i++)
1467 eqp->ic_defs[i] = cs_xdef_free(eqp->ic_defs[i]);
1468 BFT_FREE(eqp->ic_defs);
1469
1470 }
1471
1472 /* Information related to the linear algebra settings */
1473
1474 cs_param_sles_free(&(eqp->sles_param));
1475
1476 BFT_FREE(eqp->name);
1477 }
1478
1479 /*----------------------------------------------------------------------------*/
1480 /*!
1481 * \brief Free a \ref cs_equation_param_t
1482 *
1483 * \param[in, out] eqp pointer to a \ref cs_equation_param_t
1484 *
1485 * \return a NULL pointer
1486 */
1487 /*----------------------------------------------------------------------------*/
1488
1489 cs_equation_param_t *
cs_equation_param_free(cs_equation_param_t * eqp)1490 cs_equation_param_free(cs_equation_param_t *eqp)
1491 {
1492 if (eqp == NULL)
1493 return NULL;
1494
1495 cs_equation_clear_param(eqp);
1496
1497 BFT_FREE(eqp);
1498
1499 return NULL;
1500 }
1501
1502 /*----------------------------------------------------------------------------*/
1503 /*!
1504 * \brief Set a parameter attached to a keyname in a \ref cs_equation_param_t
1505 * structure
1506 *
1507 * \param[in, out] eqp pointer to a \ref cs_equation_param_t structure
1508 * \param[in] key key related to the member of eq to set
1509 * \param[in] keyval accessor to the value to set
1510 */
1511 /*----------------------------------------------------------------------------*/
1512
1513 void
cs_equation_param_set(cs_equation_param_t * eqp,cs_equation_key_t key,const char * keyval)1514 cs_equation_param_set(cs_equation_param_t *eqp,
1515 cs_equation_key_t key,
1516 const char *keyval)
1517 {
1518 /* Sanity checks */
1519 if (eqp == NULL)
1520 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
1521 if (keyval == NULL)
1522 bft_error(__FILE__, __LINE__, 0, "%s: Eq: %s: Key value is empty",
1523 __func__, eqp->name);
1524 if (eqp->flag & CS_EQUATION_LOCKED)
1525 bft_error(__FILE__, __LINE__, 0,
1526 _(" %s: Equation %s is not modifiable anymore.\n"
1527 " Please check your settings."), eqp->name, __func__);
1528
1529 /* Conversion of the string to lower case */
1530 char val[CS_BASE_STRING_LEN];
1531 for (size_t i = 0; i < strlen(keyval); i++)
1532 val[i] = tolower(keyval[i]);
1533 val[strlen(keyval)] = '\0';
1534
1535 /* Set the couple (key,keyval) */
1536 _set_key(eqp, key, val);
1537 }
1538
1539 /*----------------------------------------------------------------------------*/
1540 /*!
1541 * \brief Set parameters for initializing SLES structures used for the
1542 * resolution of the linear system.
1543 * Settings are related to this equation.
1544 *
1545 * \param[in, out] eqp pointer to a cs_equation_param_t structure
1546 */
1547 /*----------------------------------------------------------------------------*/
1548
1549 void
cs_equation_param_set_sles(cs_equation_param_t * eqp)1550 cs_equation_param_set_sles(cs_equation_param_t *eqp)
1551 {
1552 /* Define a cs_sles_t structure using the field_id related to the variable
1553 * field associated to this equation */
1554 int ierr = cs_param_sles_set(true, eqp->sles_param);
1555
1556 if (ierr == -1)
1557 bft_error(__FILE__, __LINE__, 0,
1558 "%s: The requested class of solvers is not available"
1559 " for the equation %s\n"
1560 " Please modify your settings.", __func__, eqp->name);
1561 }
1562
1563 /*----------------------------------------------------------------------------*/
1564 /*!
1565 * \brief Last modification of the cs_equation_param_t structure before
1566 * launching the computation
1567 *
1568 * \param[in, out] eqp pointer to a \ref cs_equation_param_t structure
1569 */
1570 /*----------------------------------------------------------------------------*/
1571
1572 void
cs_equation_param_last_stage(cs_equation_param_t * eqp)1573 cs_equation_param_last_stage(cs_equation_param_t *eqp)
1574 {
1575 if (eqp == NULL)
1576 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
1577
1578 if (eqp->flag & CS_EQUATION_LOCKED)
1579 bft_error(__FILE__, __LINE__, 0,
1580 _(" %s: Equation %s is not modifiable anymore.\n"
1581 " Please check your settings."), eqp->name, __func__);
1582
1583 if (eqp->do_lumping) { /* Activate a set of options */
1584
1585 eqp->reaction_hodgep.algo = CS_HODGE_ALGO_VORONOI;
1586 eqp->time_hodgep.algo = CS_HODGE_ALGO_VORONOI;
1587
1588 /* A simple barycentric quadrature rule is applied to all source terms */
1589 for (int i = 0; i < eqp->n_source_terms; i++)
1590 cs_xdef_set_quadrature(eqp->source_terms[i], CS_QUADRATURE_BARY);
1591
1592 }
1593
1594 }
1595
1596 /*----------------------------------------------------------------------------*/
1597 /*!
1598 * \brief Print the detail of a \ref cs_equation_param_t structure
1599 *
1600 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
1601 */
1602 /*----------------------------------------------------------------------------*/
1603
1604 void
cs_equation_param_log(const cs_equation_param_t * eqp)1605 cs_equation_param_log(const cs_equation_param_t *eqp)
1606 {
1607 if (eqp == NULL)
1608 return;
1609
1610 const char *eqname = eqp->name;
1611
1612 char prefix[256];
1613 assert(strlen(eqname) < 200); /* Check that prefix is large enough */
1614
1615 /* High-level settings */
1616
1617 cs_log_printf(CS_LOG_SETUP, "\n### %s | High-level settings\n", eqname);
1618 cs_log_printf(CS_LOG_SETUP, " * %s | Type: ", eqname);
1619 switch (eqp->type) {
1620 case CS_EQUATION_TYPE_GROUNDWATER:
1621 cs_log_printf(CS_LOG_SETUP, "Associated to groundwater flows\n");
1622 break;
1623 case CS_EQUATION_TYPE_MAXWELL:
1624 cs_log_printf(CS_LOG_SETUP, "Associated to the Maxwell module\n");
1625 break;
1626 case CS_EQUATION_TYPE_NAVSTO:
1627 cs_log_printf(CS_LOG_SETUP, "Associated to the Navier-Stokes system\n");
1628 break;
1629 case CS_EQUATION_TYPE_PREDEFINED:
1630 cs_log_printf(CS_LOG_SETUP, "Predefined\n");
1631 break;
1632 case CS_EQUATION_TYPE_SOLIDIFICATION:
1633 cs_log_printf(CS_LOG_SETUP, "Associated to the solidification module\n");
1634 break;
1635 case CS_EQUATION_TYPE_THERMAL:
1636 cs_log_printf(CS_LOG_SETUP, "Associated to the thermal module\n");
1637 break;
1638 case CS_EQUATION_TYPE_USER:
1639 cs_log_printf(CS_LOG_SETUP, "User-defined\n");
1640 break;
1641
1642 default:
1643 bft_error(__FILE__, __LINE__, 0,
1644 " Eq. %s has no type.\n Please check your settings.", eqname);
1645 }
1646
1647 bool unsteady = (eqp->flag & CS_EQUATION_UNSTEADY) ? true : false;
1648 bool convection = (eqp->flag & CS_EQUATION_CONVECTION) ? true : false;
1649 bool diffusion = (eqp->flag & CS_EQUATION_DIFFUSION) ? true : false;
1650 bool curlcurl = (eqp->flag & CS_EQUATION_CURLCURL) ? true : false;
1651 bool graddiv = (eqp->flag & CS_EQUATION_GRADDIV) ? true : false;
1652 bool reaction = (eqp->flag & CS_EQUATION_REACTION) ? true : false;
1653 bool source_term = (eqp->n_source_terms > 0) ? true : false;
1654 bool force_values = (eqp->flag & CS_EQUATION_FORCE_VALUES) ? true : false;
1655
1656 cs_log_printf(CS_LOG_SETUP,
1657 " * %s | Terms: unsteady:%s, convection:%s, diffusion:%s\n",
1658 eqname, cs_base_strtf(unsteady), cs_base_strtf(convection),
1659 cs_base_strtf(diffusion));
1660 cs_log_printf(CS_LOG_SETUP,
1661 " * %s | Terms: curl-curl:%s, grad-div:%s\n",
1662 eqname, cs_base_strtf(curlcurl), cs_base_strtf(graddiv));
1663 cs_log_printf(CS_LOG_SETUP,
1664 " * %s | Terms: reaction:%s, source term:%s,"
1665 " force internal values: %s\n",
1666 eqname,cs_base_strtf(reaction), cs_base_strtf(source_term),
1667 cs_base_strtf(force_values));
1668
1669 if (eqp->space_scheme < CS_SPACE_N_SCHEMES)
1670 cs_log_printf(CS_LOG_SETUP, " * %s | Space scheme: %s\n",
1671 eqname, cs_param_get_space_scheme_name(eqp->space_scheme));
1672 else
1673 bft_error(__FILE__, __LINE__, 0,
1674 " Undefined space scheme for eq. %s", eqname);
1675
1676 cs_log_printf(CS_LOG_SETUP, " * %s | Space poly degree: %d\n",
1677 eqname, eqp->space_poly_degree);
1678 cs_log_printf(CS_LOG_SETUP, " * %s | Verbosity: %d\n",
1679 eqname, eqp->verbosity);
1680
1681 if (cs_glob_n_threads > 1) {
1682 if (eqp->omp_assembly_choice == CS_PARAM_ASSEMBLE_OMP_CRITICAL)
1683 cs_log_printf(CS_LOG_SETUP, " * %s | OpenMP.Assembly.Choice: %s\n",
1684 eqname, "critical");
1685 else if (eqp->omp_assembly_choice == CS_PARAM_ASSEMBLE_OMP_ATOMIC)
1686 cs_log_printf(CS_LOG_SETUP, " * %s | OpenMP.Assembly.Choice: %s\n",
1687 eqname, "atomic");
1688 }
1689
1690 /* Boundary conditions */
1691
1692 cs_log_printf(CS_LOG_SETUP, "\n### %s | Boundary condition settings\n",
1693 eqname);
1694 cs_log_printf(CS_LOG_SETUP,
1695 " * %s | Boundary conditions | Default: %s\n",
1696 eqname, cs_param_get_bc_name(eqp->default_bc));
1697 cs_log_printf(CS_LOG_SETUP,
1698 " * %s | Boundary conditions | Enforcement: %s\n",
1699 eqname,
1700 cs_param_get_bc_enforcement_name(eqp->default_enforcement));
1701 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_PENALIZED)
1702 cs_log_printf(CS_LOG_SETUP,
1703 " * %s | Boundary conditions | Penalization coefficient:"
1704 " %5.3e\n", eqname, eqp->strong_pena_bc_coeff);
1705 else if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
1706 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM)
1707 cs_log_printf(CS_LOG_SETUP,
1708 " * %s | Boundary conditions | Penalization coefficient:"
1709 " %5.3e\n", eqname, eqp->weak_pena_bc_coeff);
1710
1711 cs_log_printf(CS_LOG_SETUP, " * %s | Boundary conditions |"
1712 " Number of definitions: %d\n", eqname, eqp->n_bc_defs);
1713
1714 if (eqp->verbosity > 0) {
1715 char desc[128];
1716 for (int id = 0; id < eqp->n_bc_defs; id++) {
1717 const cs_xdef_t *d = eqp->bc_defs[id];
1718
1719 cs_cdo_bc_get_desc(d->meta, desc);
1720 sprintf(prefix, " Definition %3d", id);
1721 cs_log_printf(CS_LOG_SETUP, "\n%s | Type: %s\n", prefix, desc);
1722 cs_xdef_log(prefix, d);
1723 }
1724 }
1725
1726 if (unsteady) {
1727
1728 cs_log_printf(CS_LOG_SETUP, "\n### %s | Time settings\n", eqname);
1729 cs_log_printf(CS_LOG_SETUP,
1730 " * %s | Initial conditions | Number of definitions: %d",
1731 eqname, eqp->n_ic_defs);
1732 if (eqp->n_ic_defs > 0)
1733 cs_log_printf(CS_LOG_SETUP, "\n\n");
1734 for (int i = 0; i < eqp->n_ic_defs; i++) {
1735 sprintf(prefix, " Definition %3d", i);
1736 cs_xdef_log(prefix, eqp->ic_defs[i]);
1737 }
1738
1739 const char *time_scheme = cs_param_get_time_scheme_name(eqp->time_scheme);
1740 if (time_scheme != NULL) {
1741 cs_log_printf(CS_LOG_SETUP, "\n * %s | Time scheme: %s",
1742 eqname, time_scheme);
1743 if (eqp->time_scheme == CS_TIME_SCHEME_THETA)
1744 cs_log_printf(CS_LOG_SETUP, " with value %f\n", eqp->theta);
1745 else
1746 cs_log_printf(CS_LOG_SETUP, "\n");
1747 }
1748 else
1749 bft_error(__FILE__, __LINE__, 0, " Invalid time scheme.");
1750
1751 cs_log_printf(CS_LOG_SETUP, " * %s | Mass.Lumping: %s\n",
1752 eqname, cs_base_strtf(eqp->do_lumping));
1753 cs_log_printf(CS_LOG_SETUP, " * %s | Time property: %s\n\n",
1754 eqname, cs_property_get_name(eqp->time_property));
1755
1756 sprintf(prefix, " Time Hodge op. ");
1757 cs_hodge_param_log(prefix, eqp->time_property, eqp->time_hodgep);
1758
1759 } /* Unsteady term */
1760
1761 if (diffusion) {
1762
1763 cs_log_printf(CS_LOG_SETUP, "\n### %s | Diffusion term settings\n", eqname);
1764 cs_log_printf(CS_LOG_SETUP, " * %s | Diffusion property: %s\n\n",
1765 eqname, cs_property_get_name(eqp->diffusion_property));
1766
1767 sprintf(prefix, " Diffusion Hodge op. ");
1768 cs_hodge_param_log(prefix, eqp->diffusion_property, eqp->diffusion_hodgep);
1769
1770 } /* Diffusion term */
1771
1772 if (curlcurl) {
1773
1774 cs_log_printf(CS_LOG_SETUP, "\n### %s | Curl-Curl term settings\n", eqname);
1775 cs_log_printf(CS_LOG_SETUP, " * %s | Curl-Curl property: %s\n\n",
1776 eqname, cs_property_get_name(eqp->curlcurl_property));
1777
1778 sprintf(prefix, " Curl-curl Hodge op. ");
1779 cs_hodge_param_log(prefix, eqp->curlcurl_property, eqp->curlcurl_hodgep);
1780
1781 } /* Curl-curl term */
1782
1783 if (graddiv) {
1784
1785 cs_log_printf(CS_LOG_SETUP, "\n### %s | Grad-Div term settings\n", eqname);
1786 cs_log_printf(CS_LOG_SETUP, " * %s | Grad-Div property: %s\n\n",
1787 eqname, cs_property_get_name(eqp->graddiv_property));
1788
1789 sprintf(prefix, " Grad-Div Hodge op. ");
1790 cs_hodge_param_log(prefix, eqp->graddiv_property, eqp->graddiv_hodgep);
1791
1792 } /* Curl-curl term */
1793
1794 if (convection) {
1795
1796 cs_log_printf(CS_LOG_SETUP, "\n### %s | Advection term settings\n", eqname);
1797 cs_log_printf(CS_LOG_SETUP, " * %s | Advection.Field: %s\n",
1798 eqname, cs_advection_field_get_name(eqp->adv_field));
1799 if (eqp->adv_scaling_property != NULL)
1800 cs_log_printf(CS_LOG_SETUP, " * %s | Scaling.Property: %s\n",
1801 eqname, cs_property_get_name(eqp->adv_scaling_property));
1802
1803 cs_log_printf(CS_LOG_SETUP, " * %s | Advection.Formulation: %s\n",
1804 eqname,
1805 cs_param_get_advection_form_name(eqp->adv_formulation));
1806
1807 cs_log_printf(CS_LOG_SETUP, " * %s | Advection.Scheme: %s\n",
1808 eqname,
1809 cs_param_get_advection_scheme_name(eqp->adv_scheme));
1810
1811 /* Piece of information specific to a scheme */
1812
1813 if (eqp->adv_scheme == CS_PARAM_ADVECTION_SCHEME_HYBRID_CENTERED_UPWIND)
1814 cs_log_printf(CS_LOG_SETUP, " * %s | Upwind.Portion: %3.2f %%\n",
1815 eqname, 100*eqp->upwind_portion);
1816 else if (eqp->adv_scheme == CS_PARAM_ADVECTION_SCHEME_CIP ||
1817 eqp->adv_scheme == CS_PARAM_ADVECTION_SCHEME_CIP_CW)
1818 cs_log_printf(CS_LOG_SETUP, " * %s | CIP.coef: %f\n",
1819 eqname, cs_cdo_advection_get_cip_coef());
1820
1821 cs_log_printf(CS_LOG_SETUP, " * %s | Advection.Strategy: %s\n",
1822 eqname,
1823 cs_param_get_advection_strategy_name(eqp->adv_strategy));
1824 cs_log_printf(CS_LOG_SETUP, " * %s | Advection.Extrapolation: %s\n",
1825 eqname,
1826 cs_param_get_advection_extrapol_name(eqp->adv_extrapol));
1827
1828 } /* Advection term */
1829
1830 if (reaction) {
1831
1832 cs_log_printf(CS_LOG_SETUP, "\n### %s | Reaction settings\n", eqname);
1833 cs_log_printf(CS_LOG_SETUP, " * %s | Reaction | Number of terms: %d\n",
1834 eqname, eqp->n_reaction_terms);
1835
1836 sprintf(prefix, " Reaction Hodge op. ");
1837 cs_hodge_param_log(prefix, NULL, eqp->reaction_hodgep);
1838
1839 } /* Reaction terms */
1840
1841 if (source_term) {
1842
1843 cs_log_printf(CS_LOG_SETUP, "\n### %s | Source term settings\n", eqname);
1844 cs_log_printf(CS_LOG_SETUP, " * %s | Source terms | Number of terms: %d\n",
1845 eqname, eqp->n_source_terms);
1846
1847 for (int s_id = 0; s_id < eqp->n_source_terms; s_id++) {
1848 sprintf(prefix, " Definition %3d", s_id);
1849 cs_xdef_log(prefix, eqp->source_terms[s_id]);
1850 }
1851
1852 } /* Source terms */
1853
1854 /* Interior enforcement(s) */
1855
1856 cs_log_printf(CS_LOG_SETUP,
1857 "\n### %s | Number of interior enforcements: %d\n",
1858 eqname, eqp->n_enforcements);
1859
1860 if (eqp->n_enforcements > 0) {
1861
1862 for (int i = 0; i < eqp->n_enforcements; i++)
1863 cs_enforcement_param_log(eqname, eqp->enforcement_params[i]);
1864
1865 }
1866
1867 /* Iterative solver information */
1868
1869 cs_param_sles_log(eqp->sles_param);
1870
1871 }
1872
1873 /*----------------------------------------------------------------------------*/
1874 /*!
1875 * \brief Ask if the parameter settings of the equation has requested the
1876 * treatment of Robin boundary conditions
1877 *
1878 * \param[in] eqp pointer to a \ref cs_equation_param_t
1879 *
1880 * \return true or false
1881 */
1882 /*----------------------------------------------------------------------------*/
1883
1884 bool
cs_equation_param_has_robin_bc(const cs_equation_param_t * eqp)1885 cs_equation_param_has_robin_bc(const cs_equation_param_t *eqp)
1886 {
1887 if (eqp == NULL)
1888 return false;
1889
1890 for (int i = 0; i < eqp->n_bc_defs; i++) {
1891 cs_xdef_t *def = eqp->bc_defs[i];
1892 if (def->meta & CS_CDO_BC_ROBIN)
1893 return true;
1894 }
1895 return false;
1896 }
1897
1898 /*----------------------------------------------------------------------------*/
1899 /*!
1900 * \brief Define the initial condition for the unknown related to this equation
1901 * This definition can be done on a specified mesh location.
1902 * By default, the unknown is set to zero everywhere.
1903 * Here a constant value is set to all the entities belonging to the
1904 * given mesh location
1905 *
1906 * \param[in, out] eqp pointer to a cs_equation_param_t structure
1907 * \param[in] z_name name of the associated zone (if NULL or
1908 * "" all cells are considered)
1909 * \param[in] val pointer to the value
1910 *
1911 * \return a pointer to the new \ref cs_xdef_t structure
1912 */
1913 /*----------------------------------------------------------------------------*/
1914
1915 cs_xdef_t *
cs_equation_add_ic_by_value(cs_equation_param_t * eqp,const char * z_name,cs_real_t * val)1916 cs_equation_add_ic_by_value(cs_equation_param_t *eqp,
1917 const char *z_name,
1918 cs_real_t *val)
1919 {
1920 if (eqp == NULL)
1921 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
1922
1923 /* Add a new cs_xdef_t structure */
1924 int z_id = cs_get_vol_zone_id(z_name);
1925
1926 cs_flag_t meta_flag = 0;
1927 if (z_id == 0)
1928 meta_flag |= CS_FLAG_FULL_LOC;
1929
1930 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
1931 eqp->dim,
1932 z_id,
1933 CS_FLAG_STATE_UNIFORM, /* state flag */
1934 meta_flag,
1935 val);
1936
1937 int new_id = eqp->n_ic_defs;
1938 eqp->n_ic_defs += 1;
1939 BFT_REALLOC(eqp->ic_defs, eqp->n_ic_defs, cs_xdef_t *);
1940 eqp->ic_defs[new_id] = d;
1941
1942 return d;
1943 }
1944
1945 /*----------------------------------------------------------------------------*/
1946 /*!
1947 * \brief Define the initial condition for the unknown related to this equation
1948 * This definition can be done on a specified mesh location.
1949 * By default, the unknown is set to zero everywhere.
1950 * Here the value related to all the entities belonging to the
1951 * given mesh location is such that the integral over these cells
1952 * returns the requested quantity
1953 *
1954 * \param[in, out] eqp pointer to a cs_equation_param_t structure
1955 * \param[in] z_name name of the associated zone (if NULL or
1956 * "" all cells are considered)
1957 * \param[in] quantity quantity to distribute over the mesh location
1958 *
1959 * \return a pointer to the new \ref cs_xdef_t structure
1960 */
1961 /*----------------------------------------------------------------------------*/
1962
1963 cs_xdef_t *
cs_equation_add_ic_by_qov(cs_equation_param_t * eqp,const char * z_name,double quantity)1964 cs_equation_add_ic_by_qov(cs_equation_param_t *eqp,
1965 const char *z_name,
1966 double quantity)
1967 {
1968 if (eqp == NULL)
1969 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
1970
1971 /* Add a new cs_xdef_t structure */
1972 int z_id = cs_get_vol_zone_id(z_name);
1973
1974 cs_flag_t meta_flag = 0;
1975 if (z_id == 0)
1976 meta_flag |= CS_FLAG_FULL_LOC;
1977
1978 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_QOV,
1979 eqp->dim,
1980 z_id,
1981 0, /* state flag */
1982 meta_flag,
1983 &quantity);
1984
1985 int new_id = eqp->n_ic_defs;
1986 eqp->n_ic_defs += 1;
1987 BFT_REALLOC(eqp->ic_defs, eqp->n_ic_defs, cs_xdef_t *);
1988 eqp->ic_defs[new_id] = d;
1989
1990 return d;
1991 }
1992
1993 /*----------------------------------------------------------------------------*/
1994 /*!
1995 * \brief Define the initial condition for the unknown related to this
1996 * equation. This definition can be done on a specified mesh location.
1997 * By default, the unknown is set to zero everywhere.
1998 * Here the initial value is set according to an analytical function
1999 *
2000 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2001 * \param[in] z_name name of the associated zone (if NULL or "" if
2002 * all cells are considered)
2003 * \param[in] analytic pointer to an analytic function
2004 * \param[in] input NULL or pointer to a structure cast on-the-fly
2005 *
2006 * \return a pointer to the new \ref cs_xdef_t structure
2007 */
2008 /*----------------------------------------------------------------------------*/
2009
2010 cs_xdef_t *
cs_equation_add_ic_by_analytic(cs_equation_param_t * eqp,const char * z_name,cs_analytic_func_t * analytic,void * input)2011 cs_equation_add_ic_by_analytic(cs_equation_param_t *eqp,
2012 const char *z_name,
2013 cs_analytic_func_t *analytic,
2014 void *input)
2015 {
2016 if (eqp == NULL)
2017 bft_error(__FILE__, __LINE__, 0, _err_empty_eqp);
2018
2019 /* Add a new cs_xdef_t structure */
2020 int z_id = cs_get_vol_zone_id(z_name);
2021
2022 cs_flag_t meta_flag = 0;
2023 if (z_id == 0)
2024 meta_flag |= CS_FLAG_FULL_LOC;
2025
2026 cs_xdef_analytic_context_t ac = { .z_id = z_id,
2027 .func = analytic,
2028 .input = input,
2029 .free_input = NULL };
2030
2031 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
2032 eqp->dim, z_id,
2033 0, /* state flag */
2034 meta_flag,
2035 &ac);
2036
2037 int new_id = eqp->n_ic_defs;
2038 eqp->n_ic_defs += 1;
2039 BFT_REALLOC(eqp->ic_defs, eqp->n_ic_defs, cs_xdef_t *);
2040 eqp->ic_defs[new_id] = d;
2041
2042 return d;
2043 }
2044
2045 /*----------------------------------------------------------------------------*/
2046 /*!
2047 * \brief Set a boundary condition from an existing \ref cs_xdef_t structure
2048 * The lifecycle of the cs_xdef_t structure is now managed by the
2049 * current \ref cs_equation_param_t structure.
2050 *
2051 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2052 * \param[in] xdef pointer to the \ref cs_xdef_t structure to transfer
2053 */
2054 /*----------------------------------------------------------------------------*/
2055
2056 void
cs_equation_add_xdef_bc(cs_equation_param_t * eqp,cs_xdef_t * xdef)2057 cs_equation_add_xdef_bc(cs_equation_param_t *eqp,
2058 cs_xdef_t *xdef)
2059 {
2060 if (eqp == NULL)
2061 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2062
2063 int new_id = eqp->n_bc_defs;
2064 eqp->n_bc_defs += 1;
2065 BFT_REALLOC(eqp->bc_defs, eqp->n_bc_defs, cs_xdef_t *);
2066 eqp->bc_defs[new_id] = xdef;
2067 }
2068
2069 /*----------------------------------------------------------------------------*/
2070 /*!
2071 * \brief Define and initialize a new structure to set a boundary condition
2072 * related to the given equation structure
2073 * z_name corresponds to the name of a pre-existing cs_zone_t
2074 *
2075 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2076 * \param[in] bc_type type of boundary condition to add
2077 * \param[in] z_name name of the related boundary zone
2078 * \param[in] values pointer to a array storing the values
2079 *
2080 * \return a pointer to the new \ref cs_xdef_t structure
2081 */
2082 /*----------------------------------------------------------------------------*/
2083
2084 cs_xdef_t *
cs_equation_add_bc_by_value(cs_equation_param_t * eqp,const cs_param_bc_type_t bc_type,const char * z_name,cs_real_t * values)2085 cs_equation_add_bc_by_value(cs_equation_param_t *eqp,
2086 const cs_param_bc_type_t bc_type,
2087 const char *z_name,
2088 cs_real_t *values)
2089 {
2090 if (eqp == NULL)
2091 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2092
2093 /* Add a new cs_xdef_t structure */
2094 int dim = eqp->dim;
2095 if (bc_type == CS_PARAM_BC_NEUMANN_FULL ||
2096 bc_type == CS_PARAM_BC_HMG_NEUMANN)
2097 dim *= 3; /* vector if scalar eq, tensor if vector eq. */
2098
2099 if (bc_type == CS_PARAM_BC_ROBIN) {
2100 /* FluxNormal + alpha * (u - u_0) = beta => Set (alpha, u_0, beta) */
2101 if (eqp->dim == 1)
2102 dim = 3;
2103 else
2104 bft_error(__FILE__, __LINE__, 0,
2105 "%s: This situation is not handled yet.\n", __func__);
2106 }
2107
2108 cs_flag_t meta_flag = (eqp-> space_scheme == CS_SPACE_SCHEME_LEGACY) ?
2109 (cs_flag_t)bc_type : cs_cdo_bc_get_flag(bc_type);
2110
2111 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2112 dim,
2113 cs_get_bdy_zone_id(z_name),
2114 CS_FLAG_STATE_UNIFORM, /* state flag */
2115 meta_flag,
2116 (void *)values);
2117
2118 int new_id = eqp->n_bc_defs;
2119 eqp->n_bc_defs += 1;
2120 BFT_REALLOC(eqp->bc_defs, eqp->n_bc_defs, cs_xdef_t *);
2121 eqp->bc_defs[new_id] = d;
2122
2123 return d;
2124 }
2125
2126 /*----------------------------------------------------------------------------*/
2127 /*!
2128 * \brief Define and initialize a new structure to set a boundary condition
2129 * related to the given equation structure
2130 * z_name corresponds to the name of a pre-existing cs_zone_t
2131 *
2132 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2133 * \param[in] bc_type type of boundary condition to add
2134 * \param[in] z_name name of the related boundary zone
2135 * \param[in] loc information to know where are located values
2136 * \param[in] array pointer to an array
2137 * \param[in] is_owner transfer the lifecycle to the cs_xdef_t structure
2138 * (true or false)
2139 * \param[in] index optional pointer to the array index
2140 *
2141 * \return a pointer to the new allocated \ref cs_xdef_t structure
2142 */
2143 /*----------------------------------------------------------------------------*/
2144
2145 cs_xdef_t *
cs_equation_add_bc_by_array(cs_equation_param_t * eqp,const cs_param_bc_type_t bc_type,const char * z_name,cs_flag_t loc,cs_real_t * array,bool is_owner,cs_lnum_t * index)2146 cs_equation_add_bc_by_array(cs_equation_param_t *eqp,
2147 const cs_param_bc_type_t bc_type,
2148 const char *z_name,
2149 cs_flag_t loc,
2150 cs_real_t *array,
2151 bool is_owner,
2152 cs_lnum_t *index)
2153 {
2154 if (eqp == NULL)
2155 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2156
2157 assert(cs_flag_test(loc, cs_flag_primal_face) ||
2158 cs_flag_test(loc, cs_flag_boundary_face) ||
2159 cs_flag_test(loc, cs_flag_primal_vtx) ||
2160 cs_flag_test(loc, cs_flag_primal_edge)); /* for circulation */
2161
2162 int z_id = cs_get_bdy_zone_id(z_name);
2163
2164 /* Add a new cs_xdef_t structure */
2165 cs_xdef_array_context_t input = {.z_id = z_id,
2166 .stride = eqp->dim,
2167 .loc = loc,
2168 .values = array,
2169 .index = index,
2170 .is_owner = is_owner};
2171
2172 cs_flag_t state_flag = 0;
2173 if (loc == cs_flag_primal_face || loc == cs_flag_boundary_face)
2174 state_flag = CS_FLAG_STATE_FACEWISE;
2175
2176 int dim = eqp->dim;
2177 if (bc_type == CS_PARAM_BC_NEUMANN_FULL ||
2178 bc_type == CS_PARAM_BC_HMG_NEUMANN)
2179 dim *= 3; /* vector if scalar eq, tensor if vector eq. */
2180
2181 if (bc_type == CS_PARAM_BC_ROBIN) {
2182 /* FluxNormal = alpha * (u_0 - u) + beta => Set (alpha, beta, u_0) */
2183 if (eqp->dim == 1)
2184 dim = 3;
2185 else
2186 bft_error(__FILE__, __LINE__, 0,
2187 "%s: This situation is not handled yet.\n", __func__);
2188 }
2189
2190 cs_flag_t meta_flag = (eqp-> space_scheme == CS_SPACE_SCHEME_LEGACY) ?
2191 (cs_flag_t)bc_type : cs_cdo_bc_get_flag(bc_type);
2192
2193 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_ARRAY,
2194 dim,
2195 z_id,
2196 state_flag,
2197 meta_flag,
2198 (void *)&input);
2199
2200 int new_id = eqp->n_bc_defs;
2201 eqp->n_bc_defs += 1;
2202 BFT_REALLOC(eqp->bc_defs, eqp->n_bc_defs, cs_xdef_t *);
2203 eqp->bc_defs[new_id] = d;
2204
2205 return d;
2206 }
2207
2208 /*----------------------------------------------------------------------------*/
2209 /*!
2210 * \brief Define and initialize a new structure to set a boundary condition
2211 * related to the given equation structure
2212 * z_name corresponds to the name of a pre-existing cs_zone_t
2213 *
2214 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2215 * \param[in] bc_type type of boundary condition to add
2216 * \param[in] z_name name of the related boundary zone
2217 * \param[in] field pointer to a cs_field_t structure
2218 *
2219 * \return a pointer to the new allocated \ref cs_xdef_t structure
2220 */
2221 /*----------------------------------------------------------------------------*/
2222
2223 cs_xdef_t *
cs_equation_add_bc_by_field(cs_equation_param_t * eqp,const cs_param_bc_type_t bc_type,const char * z_name,cs_field_t * field)2224 cs_equation_add_bc_by_field(cs_equation_param_t *eqp,
2225 const cs_param_bc_type_t bc_type,
2226 const char *z_name,
2227 cs_field_t *field)
2228 {
2229 if (eqp == NULL)
2230 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2231
2232 int z_id = cs_get_bdy_zone_id(z_name);
2233
2234 int dim = eqp->dim;
2235 if (bc_type == CS_PARAM_BC_NEUMANN_FULL ||
2236 bc_type == CS_PARAM_BC_HMG_NEUMANN)
2237 dim *= 3; /* vector if scalar eq, tensor if vector eq. */
2238
2239 if (bc_type == CS_PARAM_BC_ROBIN) {
2240 /* FluxNormal = alpha * (u_0 - u) + beta => Set (alpha, beta, u_0) */
2241 if (eqp->dim == 1)
2242 dim = 3;
2243 else
2244 bft_error(__FILE__, __LINE__, 0,
2245 "%s: This situation is not handled yet.\n", __func__);
2246 }
2247
2248 /* Sanity checks */
2249 assert(field != NULL);
2250
2251 if (dim != field->dim)
2252 bft_error(__FILE__, __LINE__, 0,
2253 "%s: Invalid dimension for field %s\n", __func__, field->name);
2254
2255 cs_flag_t state_flag = CS_FLAG_STATE_FACEWISE;
2256 cs_flag_t meta_flag = (eqp-> space_scheme == CS_SPACE_SCHEME_LEGACY) ?
2257 (cs_flag_t)bc_type : cs_cdo_bc_get_flag(bc_type);
2258
2259 /* Add a new cs_xdef_t structure */
2260 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_FIELD,
2261 dim,
2262 z_id,
2263 state_flag,
2264 meta_flag,
2265 field);
2266
2267 int new_id = eqp->n_bc_defs;
2268 eqp->n_bc_defs += 1;
2269 BFT_REALLOC(eqp->bc_defs, eqp->n_bc_defs, cs_xdef_t *);
2270 eqp->bc_defs[new_id] = d;
2271
2272 return d;
2273 }
2274
2275 /*----------------------------------------------------------------------------*/
2276 /*!
2277 * \brief Define and initialize a new structure to set a boundary condition
2278 * related to the given equation param structure
2279 * ml_name corresponds to the name of a pre-existing cs_mesh_location_t
2280 *
2281 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2282 * \param[in] bc_type type of boundary condition to add
2283 * \param[in] z_name name of the associated zone (if NULL or "" if
2284 * all cells are considered)
2285 * \param[in] analytic pointer to an analytic function defining the value
2286 * \param[in] input NULL or pointer to a structure cast on-the-fly
2287 *
2288 * \return a pointer to the new \ref cs_xdef_t structure
2289 */
2290 /*----------------------------------------------------------------------------*/
2291
2292 cs_xdef_t *
cs_equation_add_bc_by_analytic(cs_equation_param_t * eqp,const cs_param_bc_type_t bc_type,const char * z_name,cs_analytic_func_t * analytic,void * input)2293 cs_equation_add_bc_by_analytic(cs_equation_param_t *eqp,
2294 const cs_param_bc_type_t bc_type,
2295 const char *z_name,
2296 cs_analytic_func_t *analytic,
2297 void *input)
2298 {
2299 if (eqp == NULL)
2300 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2301
2302 /* Set the value for dim */
2303 int dim = eqp->dim;
2304
2305 if (bc_type == CS_PARAM_BC_NEUMANN_FULL ||
2306 bc_type == CS_PARAM_BC_HMG_NEUMANN)
2307 dim *= 3; /* vector if scalar eq, tensor if vector eq. */
2308
2309 if (bc_type == CS_PARAM_BC_CIRCULATION) {
2310 /* This is a vector-valued equation but the DoF is scalar-valued since
2311 * it is a circulation associated to each edge */
2312 if (eqp->dim == 3)
2313 dim = 1;
2314 else
2315 bft_error(__FILE__, __LINE__, 0,
2316 "%s: This situation is not handled.\n", __func__);
2317 }
2318
2319 if (bc_type == CS_PARAM_BC_ROBIN) {
2320 /* FluxNormal = alpha * (u_0 - u) + beta => Set (alpha, beta, u_0) */
2321 if (eqp->dim == 1)
2322 dim = 3;
2323 else
2324 bft_error(__FILE__, __LINE__, 0,
2325 "%s: This situation is not handled yet.\n", __func__);
2326 }
2327
2328 int z_id = cs_get_bdy_zone_id(z_name);
2329
2330 /* Add a new cs_xdef_t structure */
2331 cs_xdef_analytic_context_t ac = {.z_id = z_id,
2332 .func = analytic,
2333 .input = input,
2334 .free_input = NULL};
2335
2336 cs_flag_t meta_flag = (eqp-> space_scheme == CS_SPACE_SCHEME_LEGACY) ?
2337 (cs_flag_t)bc_type : cs_cdo_bc_get_flag(bc_type);
2338
2339 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
2340 dim,
2341 z_id,
2342 0, /* state */
2343 meta_flag,
2344 &ac);
2345
2346 int new_id = eqp->n_bc_defs;
2347 eqp->n_bc_defs += 1;
2348 BFT_REALLOC(eqp->bc_defs, eqp->n_bc_defs, cs_xdef_t *);
2349 eqp->bc_defs[new_id] = d;
2350
2351 return d;
2352 }
2353
2354 /*----------------------------------------------------------------------------*/
2355 /*!
2356 * \brief Define and initialize a new structure to set a boundary condition
2357 * related to the given cs_equation_param_t structure
2358 * ml_name corresponds to the name of a pre-existing cs_mesh_location_t
2359 *
2360 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2361 * \param[in] bc_type type of boundary condition to add
2362 * \param[in] z_name name of the associated zone (if NULL or "" if
2363 * all cells are considered)
2364 * \param[in] loc_flag location where values are computed
2365 * \param[in] func pointer to cs_dof_func_t function
2366 * \param[in] input NULL or pointer to a structure cast on-the-fly
2367 *
2368 * \return a pointer to the new \ref cs_xdef_t structure
2369 */
2370 /*----------------------------------------------------------------------------*/
2371
2372 cs_xdef_t *
cs_equation_add_bc_by_dof_func(cs_equation_param_t * eqp,const cs_param_bc_type_t bc_type,const char * z_name,cs_flag_t loc_flag,cs_dof_func_t * func,void * input)2373 cs_equation_add_bc_by_dof_func(cs_equation_param_t *eqp,
2374 const cs_param_bc_type_t bc_type,
2375 const char *z_name,
2376 cs_flag_t loc_flag,
2377 cs_dof_func_t *func,
2378 void *input)
2379 {
2380 if (eqp == NULL)
2381 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2382
2383 /* Set the value for dim */
2384 int dim = eqp->dim;
2385
2386 if (bc_type == CS_PARAM_BC_NEUMANN_FULL ||
2387 bc_type == CS_PARAM_BC_HMG_NEUMANN)
2388 dim *= 3; /* vector if scalar eq, tensor if vector eq. */
2389
2390 if (bc_type == CS_PARAM_BC_CIRCULATION) {
2391 /* This is a vector-valued equation but the DoF is scalar-valued since
2392 * it is a circulation associated to each edge */
2393 if (eqp->dim == 3)
2394 dim = 1;
2395 else
2396 bft_error(__FILE__, __LINE__, 0,
2397 "%s: This situation is not handled.\n", __func__);
2398 }
2399
2400 if (bc_type == CS_PARAM_BC_ROBIN) {
2401 /* FluxNormal = alpha * (u_0 - u) + beta => Set (alpha, beta, u_0) */
2402 if (eqp->dim == 1)
2403 dim = 3;
2404 else
2405 bft_error(__FILE__, __LINE__, 0,
2406 "%s: This situation is not handled yet.\n", __func__);
2407 }
2408
2409 int z_id = cs_get_bdy_zone_id(z_name);
2410
2411 /* Add a new cs_xdef_t structure */
2412 cs_xdef_dof_context_t cx = { .z_id = z_id,
2413 .loc = loc_flag,
2414 .func = func,
2415 .input = input,
2416 .free_input = NULL };
2417
2418 cs_flag_t meta_flag = (eqp-> space_scheme == CS_SPACE_SCHEME_LEGACY) ?
2419 (cs_flag_t)bc_type : cs_cdo_bc_get_flag(bc_type);
2420
2421 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_DOF_FUNCTION,
2422 dim,
2423 z_id,
2424 0, /* state */
2425 meta_flag,
2426 &cx);
2427
2428 int new_id = eqp->n_bc_defs;
2429 eqp->n_bc_defs += 1;
2430 BFT_REALLOC(eqp->bc_defs, eqp->n_bc_defs, cs_xdef_t *);
2431 eqp->bc_defs[new_id] = d;
2432
2433 return d;
2434 }
2435
2436 /*----------------------------------------------------------------------------*/
2437 /*!
2438 * \brief Return pointer to existing boundary condition definition structure
2439 * for the given equation param structure and zone.
2440 *
2441 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2442 * \param[in] z_name name of the associated zone (if NULL or "" if
2443 * all cells are considered)
2444 *
2445 * \return a pointer to the \ref cs_xdef_t structure if present, or NULL
2446 */
2447 /*----------------------------------------------------------------------------*/
2448
2449 cs_xdef_t *
cs_equation_find_bc(cs_equation_param_t * eqp,const char * z_name)2450 cs_equation_find_bc(cs_equation_param_t *eqp,
2451 const char *z_name)
2452 {
2453 if (eqp == NULL)
2454 return NULL;
2455
2456 int z_id = -2;
2457
2458 const cs_zone_t *z = cs_boundary_zone_by_name_try(z_name);
2459 if (z != NULL)
2460 z_id = z->id;
2461
2462 /* Search for given BC. */
2463 for (int i = 0; i < eqp->n_bc_defs; i++) {
2464 if (eqp->bc_defs[i]->z_id == z_id) {
2465 return eqp->bc_defs[i];
2466 }
2467 }
2468
2469 return NULL;
2470 }
2471
2472 /*----------------------------------------------------------------------------*/
2473 /*!
2474 * \brief Remove boundary condition from the given equation param structure
2475 * for a given zone.
2476 *
2477 * If no matching boundary condition is found, the function returns
2478 * silently.
2479 *
2480 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2481 * \param[in] z_name name of the associated zone (if NULL or "" if
2482 * all cells are considered)
2483 */
2484 /*----------------------------------------------------------------------------*/
2485
2486 void
cs_equation_remove_bc(cs_equation_param_t * eqp,const char * z_name)2487 cs_equation_remove_bc(cs_equation_param_t *eqp,
2488 const char *z_name)
2489 {
2490 if (eqp == NULL)
2491 return;
2492
2493 int z_id = -2;
2494
2495 const cs_zone_t *z = cs_boundary_zone_by_name_try(z_name);
2496 if (z != NULL)
2497 z_id = z->id;
2498
2499 /* Search for given BC. */
2500 int j = -1;
2501 for (int i = 0; i < eqp->n_bc_defs; i++) {
2502 if (eqp->bc_defs[i]->z_id == z_id) {
2503 j = i;
2504 break;
2505 }
2506 }
2507
2508 /* Remove it if found */
2509 if (j > -1) {
2510 eqp->bc_defs[j] = cs_xdef_free(eqp->bc_defs[j]);
2511 for (int i = j+1; i < eqp->n_bc_defs; i++) {
2512 eqp->bc_defs[i-1] = eqp->bc_defs[i];
2513 }
2514 eqp->n_bc_defs -= 1;
2515 BFT_REALLOC(eqp->bc_defs, eqp->n_bc_defs, cs_xdef_t *);
2516 }
2517 }
2518
2519 /*----------------------------------------------------------------------------*/
2520 /*!
2521 * \brief Define and initialize a new structure to set a sliding boundary
2522 * condition related to the given equation structure
2523 * z_name corresponds to the name of a pre-existing cs_zone_t
2524 *
2525 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2526 * \param[in] z_name name of the related boundary zone
2527 */
2528 /*----------------------------------------------------------------------------*/
2529
2530 void
cs_equation_add_sliding_condition(cs_equation_param_t * eqp,const char * z_name)2531 cs_equation_add_sliding_condition(cs_equation_param_t *eqp,
2532 const char *z_name)
2533 {
2534 if (eqp == NULL)
2535 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2536 if (eqp->dim < 3)
2537 bft_error(__FILE__, __LINE__, 0, "%s: Invalid dimension of equation\n",
2538 __func__);
2539
2540 /* Add two definitions: one for the normal component and one for the
2541 tangential component */
2542 BFT_REALLOC(eqp->bc_defs, eqp->n_bc_defs + 1, cs_xdef_t *);
2543
2544 cs_xdef_t *d = NULL;
2545 cs_real_t val = 0;
2546
2547 /* Add the homogeneous Dirichlet on the normal component */
2548 d = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2549 1,
2550 cs_get_bdy_zone_id(z_name),
2551 CS_FLAG_STATE_UNIFORM, /* state flag */
2552 CS_CDO_BC_SLIDING, /* meta */
2553 (void *)&val);
2554
2555 eqp->bc_defs[eqp->n_bc_defs] = d;
2556 eqp->n_bc_defs += 1;
2557 }
2558
2559 /*----------------------------------------------------------------------------*/
2560 /*!
2561 * \brief Associate a new term related to the Laplacian operator for the
2562 * equation associated to the given \ref cs_equation_param_t structure
2563 * Laplacian means div-grad (either for vector-valued or scalar-valued
2564 * equations)
2565 *
2566 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2567 * \param[in] property pointer to a cs_property_t structure
2568 */
2569 /*----------------------------------------------------------------------------*/
2570
2571 void
cs_equation_add_diffusion(cs_equation_param_t * eqp,cs_property_t * property)2572 cs_equation_add_diffusion(cs_equation_param_t *eqp,
2573 cs_property_t *property)
2574 {
2575 if (eqp == NULL)
2576 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2577 if (property == NULL)
2578 bft_error(__FILE__, __LINE__, 0,
2579 "%s: Eq. %s: Stop adding an empty property.",
2580 __func__, eqp->name);
2581
2582 eqp->flag |= CS_EQUATION_DIFFUSION;
2583 eqp->diffusion_property = property;
2584 }
2585
2586 /*----------------------------------------------------------------------------*/
2587 /*!
2588 * \brief Associate a new term related to the curl-curl operator for the
2589 * equation associated to the given \ref cs_equation_param_t structure
2590 *
2591 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2592 * \param[in] property pointer to a cs_property_t structure
2593 * \param[in] inversion > 0: true, false otherwise
2594 */
2595 /*----------------------------------------------------------------------------*/
2596
2597 void
cs_equation_add_curlcurl(cs_equation_param_t * eqp,cs_property_t * property,int inversion)2598 cs_equation_add_curlcurl(cs_equation_param_t *eqp,
2599 cs_property_t *property,
2600 int inversion)
2601 {
2602 if (eqp == NULL)
2603 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2604 if (property == NULL)
2605 bft_error(__FILE__, __LINE__, 0,
2606 "%s: Eq. %s: Stop adding an empty property.",
2607 __func__, eqp->name);
2608
2609 eqp->flag |= CS_EQUATION_CURLCURL;
2610 eqp->curlcurl_property = property;
2611
2612 if (inversion > 0)
2613 eqp->curlcurl_hodgep.inv_pty = true;
2614 else
2615 eqp->curlcurl_hodgep.inv_pty = false;
2616 }
2617
2618 /*----------------------------------------------------------------------------*/
2619 /*!
2620 * \brief Associate a new term related to the grad-div operator for the
2621 * equation associated to the given \ref cs_equation_param_t structure
2622 *
2623 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2624 * \param[in] property pointer to a cs_property_t structure
2625 */
2626 /*----------------------------------------------------------------------------*/
2627
2628 void
cs_equation_add_graddiv(cs_equation_param_t * eqp,cs_property_t * property)2629 cs_equation_add_graddiv(cs_equation_param_t *eqp,
2630 cs_property_t *property)
2631 {
2632 if (eqp == NULL)
2633 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2634 if (property == NULL)
2635 bft_error(__FILE__, __LINE__, 0,
2636 "%s: Eq. %s: Stop adding an empty property.",
2637 __func__, eqp->name);
2638
2639 eqp->flag |= CS_EQUATION_GRADDIV;
2640 eqp->graddiv_property = property;
2641 }
2642
2643 /*----------------------------------------------------------------------------*/
2644 /*!
2645 * \brief Associate a new term related to the time derivative operator for
2646 * the equation associated to the given \ref cs_equation_param_t
2647 * structure
2648 *
2649 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2650 * \param[in] property pointer to a cs_property_t structure
2651 */
2652 /*----------------------------------------------------------------------------*/
2653
2654 void
cs_equation_add_time(cs_equation_param_t * eqp,cs_property_t * property)2655 cs_equation_add_time(cs_equation_param_t *eqp,
2656 cs_property_t *property)
2657 {
2658 if (eqp == NULL)
2659 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2660 if (property == NULL)
2661 bft_error(__FILE__, __LINE__, 0,
2662 "%s: Eq. %s: Stop adding an empty property.",
2663 __func__, eqp->name);
2664
2665 eqp->flag |= CS_EQUATION_UNSTEADY;
2666 eqp->time_property = property;
2667 }
2668
2669 /*----------------------------------------------------------------------------*/
2670 /*!
2671 * \brief Associate a new term related to the advection operator for the
2672 * equation associated to the given \ref cs_equation_param_t structure
2673 *
2674 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2675 * \param[in] adv_field pointer to a cs_adv_field_t structure
2676 */
2677 /*----------------------------------------------------------------------------*/
2678
2679 void
cs_equation_add_advection(cs_equation_param_t * eqp,cs_adv_field_t * adv_field)2680 cs_equation_add_advection(cs_equation_param_t *eqp,
2681 cs_adv_field_t *adv_field)
2682 {
2683 if (eqp == NULL)
2684 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2685 if (adv_field == NULL)
2686 bft_error(__FILE__, __LINE__, 0,
2687 "%s: Eq: %s: Stop adding an empty advection field.",
2688 __func__, eqp->name);
2689
2690 eqp->flag |= CS_EQUATION_CONVECTION;
2691 eqp->adv_field = adv_field;
2692 }
2693
2694 /*----------------------------------------------------------------------------*/
2695 /*!
2696 * \brief Associate a scaling property to the advection
2697 *
2698 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2699 * \param[in] property pointer to a cs_property_t structure
2700 */
2701 /*----------------------------------------------------------------------------*/
2702
2703 void
cs_equation_add_advection_scaling_property(cs_equation_param_t * eqp,cs_property_t * property)2704 cs_equation_add_advection_scaling_property(cs_equation_param_t *eqp,
2705 cs_property_t *property)
2706 {
2707 if (eqp == NULL)
2708 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2709 if (property == NULL)
2710 bft_error(__FILE__, __LINE__, 0,
2711 "%s: Eq. %s: Stop adding an empty property.",
2712 __func__, eqp->name);
2713
2714 eqp->adv_scaling_property = property;
2715 }
2716
2717 /*----------------------------------------------------------------------------*/
2718 /*!
2719 * \brief Associate a new term related to the reaction operator for the
2720 * equation associated to the given \ref cs_equation_param_t structure
2721 *
2722 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2723 * \param[in] property pointer to a cs_property_t structure
2724 *
2725 * \return the id related to the reaction term
2726 */
2727 /*----------------------------------------------------------------------------*/
2728
2729 int
cs_equation_add_reaction(cs_equation_param_t * eqp,cs_property_t * property)2730 cs_equation_add_reaction(cs_equation_param_t *eqp,
2731 cs_property_t *property)
2732 {
2733 if (eqp == NULL)
2734 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2735 if (property == NULL)
2736 bft_error(__FILE__, __LINE__, 0,
2737 "%s: Eq. %s: Stop adding an empty property.",
2738 __func__, eqp->name);
2739
2740 /* Only this kind of reaction term is available up to now.
2741 Add a new reaction term */
2742 int new_id = eqp->n_reaction_terms;
2743 eqp->n_reaction_terms += 1;
2744 BFT_REALLOC(eqp->reaction_properties, eqp->n_reaction_terms, cs_property_t *);
2745 eqp->reaction_properties[new_id] = property;
2746
2747 /* Flag the equation with "reaction" */
2748 eqp->flag |= CS_EQUATION_REACTION;
2749
2750 return new_id;
2751 }
2752
2753 /*----------------------------------------------------------------------------*/
2754 /*!
2755 * \brief Add a new source term by initializing a cs_xdef_t structure.
2756 * Case of a definition by a constant value
2757 *
2758 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2759 * \param[in] z_name name of the associated zone (if NULL or
2760 * "" all cells are considered)
2761 * \param[in] val pointer to the value
2762 *
2763 * \return a pointer to the new \ref cs_xdef_t structure
2764 */
2765 /*----------------------------------------------------------------------------*/
2766
2767 cs_xdef_t *
cs_equation_add_source_term_by_val(cs_equation_param_t * eqp,const char * z_name,cs_real_t * val)2768 cs_equation_add_source_term_by_val(cs_equation_param_t *eqp,
2769 const char *z_name,
2770 cs_real_t *val)
2771 {
2772 if (eqp == NULL)
2773 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2774
2775 /* Add a new cs_xdef_t structure */
2776 int z_id = cs_get_vol_zone_id(z_name);
2777
2778 /* Define a flag according to the kind of space discretization */
2779 cs_flag_t state_flag = 0, meta_flag = 0;
2780 cs_source_term_set_default_flag(eqp->space_scheme, &state_flag, &meta_flag);
2781 state_flag |= CS_FLAG_STATE_UNIFORM;
2782
2783 if (z_id == 0)
2784 meta_flag |= CS_FLAG_FULL_LOC;
2785
2786 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
2787 eqp->dim,
2788 z_id,
2789 state_flag,
2790 meta_flag,
2791 (void *)val);
2792
2793 int new_id = eqp->n_source_terms;
2794 eqp->n_source_terms += 1;
2795 BFT_REALLOC(eqp->source_terms, eqp->n_source_terms, cs_xdef_t *);
2796 eqp->source_terms[new_id] = d;
2797
2798 return d;
2799 }
2800
2801 /*----------------------------------------------------------------------------*/
2802 /*!
2803 * \brief Add a new source term by initializing a cs_xdef_t structure.
2804 * Case of a definition by an analytical function
2805 *
2806 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2807 * \param[in] z_name name of the associated zone (if NULL or "" if
2808 * all cells are considered)
2809 * \param[in] func pointer to an analytical function
2810 * \param[in] input NULL or pointer to a structure cast on-the-fly
2811 *
2812 * \return a pointer to the new \ref cs_xdef_t structure
2813 */
2814 /*----------------------------------------------------------------------------*/
2815
2816 cs_xdef_t *
cs_equation_add_source_term_by_analytic(cs_equation_param_t * eqp,const char * z_name,cs_analytic_func_t * func,void * input)2817 cs_equation_add_source_term_by_analytic(cs_equation_param_t *eqp,
2818 const char *z_name,
2819 cs_analytic_func_t *func,
2820 void *input)
2821 {
2822 if (eqp == NULL)
2823 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2824
2825 /* Define a flag according to the kind of space discretization */
2826 cs_flag_t state_flag = 0, meta_flag = 0;
2827 cs_source_term_set_default_flag(eqp->space_scheme, &state_flag, &meta_flag);
2828
2829 int z_id = cs_get_vol_zone_id(z_name);
2830 if (z_id == 0)
2831 meta_flag |= CS_FLAG_FULL_LOC;
2832
2833 cs_xdef_analytic_context_t ac = { .z_id = z_id,
2834 .func = func,
2835 .input = input,
2836 .free_input = NULL };
2837
2838 /* Add a new cs_xdef_t structure */
2839 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
2840 eqp->dim,
2841 z_id,
2842 state_flag,
2843 meta_flag,
2844 &ac);
2845
2846 /* Default setting for quadrature is different in this case */
2847 cs_xdef_set_quadrature(d, CS_QUADRATURE_BARY_SUBDIV);
2848
2849 int new_id = eqp->n_source_terms;
2850 eqp->n_source_terms += 1;
2851 BFT_REALLOC(eqp->source_terms, eqp->n_source_terms, cs_xdef_t *);
2852 eqp->source_terms[new_id] = d;
2853
2854 return d;
2855 }
2856
2857 /*----------------------------------------------------------------------------*/
2858 /*!
2859 * \brief Add a new source term by initializing a cs_xdef_t structure.
2860 * Case of a definition by a DoF function
2861 *
2862 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2863 * \param[in] z_name name of the associated zone (if NULL or "" if
2864 * all cells are considered)
2865 * \param[in] loc_flag location of element ids given as parameter
2866 * \param[in] func pointer to a DoF function
2867 * \param[in] input NULL or pointer to a structure cast on-the-fly
2868 *
2869 * \return a pointer to the new \ref cs_xdef_t structure
2870 */
2871 /*----------------------------------------------------------------------------*/
2872
2873 cs_xdef_t *
cs_equation_add_source_term_by_dof_func(cs_equation_param_t * eqp,const char * z_name,cs_flag_t loc_flag,cs_dof_func_t * func,void * input)2874 cs_equation_add_source_term_by_dof_func(cs_equation_param_t *eqp,
2875 const char *z_name,
2876 cs_flag_t loc_flag,
2877 cs_dof_func_t *func,
2878 void *input)
2879 {
2880 if (eqp == NULL)
2881 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2882
2883 /* Add a new cs_xdef_t structure */
2884 int z_id = cs_get_vol_zone_id(z_name);
2885
2886 /* Define a flag according to the kind of space discretization */
2887 cs_flag_t state_flag = 0, meta_flag = 0;
2888 cs_source_term_set_default_flag(eqp->space_scheme, &state_flag, &meta_flag);
2889
2890 if (z_id == 0)
2891 meta_flag |= CS_FLAG_FULL_LOC;
2892
2893 cs_xdef_dof_context_t context = { .func = func,
2894 .input = input,
2895 .free_input = NULL,
2896 .loc = loc_flag };
2897
2898 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_DOF_FUNCTION,
2899 eqp->dim,
2900 z_id,
2901 state_flag,
2902 meta_flag,
2903 &context);
2904
2905 /* Default setting for quadrature is different in this case */
2906 cs_xdef_set_quadrature(d, CS_QUADRATURE_BARY_SUBDIV);
2907
2908 int new_id = eqp->n_source_terms;
2909 eqp->n_source_terms += 1;
2910 BFT_REALLOC(eqp->source_terms, eqp->n_source_terms, cs_xdef_t *);
2911 eqp->source_terms[new_id] = d;
2912
2913 return d;
2914 }
2915
2916 /*----------------------------------------------------------------------------*/
2917 /*!
2918 * \brief Add a new source term by initializing a cs_xdef_t structure.
2919 * Case of a definition by an array.
2920 *
2921 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2922 * \param[in] z_name name of the associated zone (if NULL or "" if
2923 * all cells are considered)
2924 * \param[in] loc information to know where are located values
2925 * \param[in] array pointer to an array
2926 * \param[in] is_owner transfer the lifecycle to the cs_xdef_t structure
2927 * (true or false)
2928 * \param[in] index optional pointer to the array index
2929 *
2930 * \return a pointer to the new \ref cs_xdef_t structure
2931 */
2932 /*----------------------------------------------------------------------------*/
2933
2934 cs_xdef_t *
cs_equation_add_source_term_by_array(cs_equation_param_t * eqp,const char * z_name,cs_flag_t loc,cs_real_t * array,bool is_owner,cs_lnum_t * index)2935 cs_equation_add_source_term_by_array(cs_equation_param_t *eqp,
2936 const char *z_name,
2937 cs_flag_t loc,
2938 cs_real_t *array,
2939 bool is_owner,
2940 cs_lnum_t *index)
2941 {
2942 if (eqp == NULL)
2943 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
2944
2945 /* Add a new cs_xdef_t structure */
2946 int z_id = cs_get_vol_zone_id(z_name);
2947
2948 /* Define a flag according to the kind of space discretization */
2949 cs_flag_t state_flag = 0, meta_flag = 0;
2950 cs_source_term_set_default_flag(eqp->space_scheme, &state_flag, &meta_flag);
2951
2952 if (cs_flag_test(loc, cs_flag_primal_vtx) == true)
2953 state_flag = CS_FLAG_STATE_POTENTIAL; /* erase predefined settings */
2954 else if (cs_flag_test(loc, cs_flag_primal_cell) == true)
2955 state_flag |= CS_FLAG_STATE_CELLWISE;
2956
2957 if (z_id == 0)
2958 meta_flag |= CS_FLAG_FULL_LOC;
2959
2960 cs_xdef_array_context_t ctxt = {.z_id = z_id,
2961 .stride = eqp->dim,
2962 .loc = loc,
2963 .values = array,
2964 .is_owner = is_owner,
2965 .index = index };
2966
2967 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_ARRAY,
2968 eqp->dim,
2969 z_id,
2970 state_flag,
2971 meta_flag,
2972 (void *)&ctxt);
2973
2974 int new_id = eqp->n_source_terms;
2975 eqp->n_source_terms += 1;
2976 BFT_REALLOC(eqp->source_terms, eqp->n_source_terms, cs_xdef_t *);
2977 eqp->source_terms[new_id] = d;
2978
2979 return d;
2980 }
2981
2982 /*----------------------------------------------------------------------------*/
2983 /*!
2984 * \brief Add a new volume mass injection definition source term by
2985 * initializing a cs_xdef_t structure, using a constant value.
2986 *
2987 * \param[in, out] eqp pointer to a cs_equation_param_t structure
2988 * \param[in] z_name name of the associated zone (if NULL or "" if
2989 * all cells are considered)
2990 * \param[in] val pointer to the value
2991 *
2992 * \return a pointer to the new \ref cs_xdef_t structure
2993 */
2994 /*----------------------------------------------------------------------------*/
2995
2996 cs_xdef_t *
cs_equation_add_volume_mass_injection_by_value(cs_equation_param_t * eqp,const char * z_name,double * val)2997 cs_equation_add_volume_mass_injection_by_value(cs_equation_param_t *eqp,
2998 const char *z_name,
2999 double *val)
3000 {
3001 if (eqp == NULL)
3002 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
3003
3004 /* Add a new cs_xdef_t structure */
3005 int z_id = cs_get_vol_zone_id(z_name);
3006
3007 cs_flag_t state_flag = 0, meta_flag = 0;
3008
3009 if (z_id == 0)
3010 meta_flag |= CS_FLAG_FULL_LOC;
3011
3012 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
3013 eqp->dim,
3014 z_id,
3015 state_flag,
3016 meta_flag,
3017 val);
3018
3019 int new_id = eqp->n_volume_mass_injections;
3020 eqp->n_volume_mass_injections += 1;
3021 BFT_REALLOC(eqp->volume_mass_injections,
3022 eqp->n_volume_mass_injections,
3023 cs_xdef_t *);
3024 eqp->volume_mass_injections[new_id] = d;
3025
3026 return d;
3027 }
3028
3029 /*----------------------------------------------------------------------------*/
3030 /*!
3031 * \brief Add a new volume mass injection definition source term by
3032 * initializing a cs_xdef_t structure, using a constant quantity
3033 * distributed over the associated zone's volume.
3034 *
3035 * \param[in, out] eqp pointer to a cs_equation_param_t structure
3036 * \param[in] z_name name of the associated zone (if NULL or "" if
3037 * all cells are considered)
3038 * \param[in] quantity pointer to quantity to distribute over the zone
3039 *
3040 * \return a pointer to the new \ref cs_xdef_t structure
3041 */
3042 /*----------------------------------------------------------------------------*/
3043
3044 cs_xdef_t *
cs_equation_add_volume_mass_injection_by_qov(cs_equation_param_t * eqp,const char * z_name,double * quantity)3045 cs_equation_add_volume_mass_injection_by_qov(cs_equation_param_t *eqp,
3046 const char *z_name,
3047 double *quantity)
3048 {
3049 if (eqp == NULL)
3050 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
3051
3052 /* Add a new cs_xdef_t structure */
3053 int z_id = cs_get_vol_zone_id(z_name);
3054
3055 cs_flag_t state_flag = 0, meta_flag = 0;
3056
3057 if (z_id == 0)
3058 meta_flag |= CS_FLAG_FULL_LOC;
3059
3060 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_QOV,
3061 eqp->dim,
3062 z_id,
3063 state_flag,
3064 meta_flag,
3065 quantity);
3066
3067 int new_id = eqp->n_volume_mass_injections;
3068 eqp->n_volume_mass_injections += 1;
3069 BFT_REALLOC(eqp->volume_mass_injections,
3070 eqp->n_volume_mass_injections,
3071 cs_xdef_t *);
3072 eqp->volume_mass_injections[new_id] = d;
3073
3074 return d;
3075 }
3076
3077 /*----------------------------------------------------------------------------*/
3078 /*!
3079 * \brief Add a new volume mass injection definition source term by
3080 * initializing a cs_xdef_t structure, using an analytical function.
3081 *
3082 * \param[in, out] eqp pointer to a cs_equation_param_t structure
3083 * \param[in] z_name name of the associated zone (if NULL or "" if
3084 * all cells are considered)
3085 * \param[in] func pointer to an analytical function
3086 * \param[in] input NULL or pointer to a structure cast on-the-fly
3087 *
3088 * \return a pointer to the new \ref cs_xdef_t structure
3089 */
3090 /*----------------------------------------------------------------------------*/
3091
3092 cs_xdef_t *
cs_equation_add_volume_mass_injection_by_analytic(cs_equation_param_t * eqp,const char * z_name,cs_analytic_func_t * func,void * input)3093 cs_equation_add_volume_mass_injection_by_analytic(cs_equation_param_t *eqp,
3094 const char *z_name,
3095 cs_analytic_func_t *func,
3096 void *input)
3097 {
3098 if (eqp == NULL)
3099 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
3100
3101 /* Add a new cs_xdef_t structure */
3102 int z_id = cs_get_vol_zone_id(z_name);
3103
3104 cs_flag_t state_flag = 0, meta_flag = 0;
3105 if (z_id == 0)
3106 meta_flag |= CS_FLAG_FULL_LOC;
3107
3108 cs_xdef_analytic_context_t ac = { .z_id = z_id,
3109 .func = func,
3110 .input = input,
3111 .free_input = NULL };
3112
3113 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
3114 eqp->dim,
3115 z_id,
3116 state_flag,
3117 meta_flag,
3118 &ac);
3119
3120 int new_id = eqp->n_volume_mass_injections;
3121 eqp->n_volume_mass_injections += 1;
3122 BFT_REALLOC(eqp->volume_mass_injections,
3123 eqp->n_volume_mass_injections,
3124 cs_xdef_t *);
3125 eqp->volume_mass_injections[new_id] = d;
3126
3127 return d;
3128 }
3129
3130 /*----------------------------------------------------------------------------*/
3131 /*!
3132 * \brief Add an enforcement of the value of degrees of freedom located at
3133 * the mesh vertices.
3134 * The spatial discretization scheme for the given equation has to be
3135 * CDO vertex-based or CDO vertex+cell-based schemes.
3136 *
3137 * One assumes that values are interlaced if eqp->dim > 1
3138 * ref_value or elt_values has to be defined. If both parameters are
3139 * defined, one keeps the definition in elt_values
3140 *
3141 * \param[in, out] eqp pointer to a cs_equation_param_t structure
3142 * \param[in] n_vertices number of vertices to enforce
3143 * \param[in] vertex_ids list of vertices
3144 * \param[in] ref_value default values or ignored (may be NULL)
3145 * \param[in] vtx_values list of associated values, ignored if NULL
3146 *
3147 * \return a pointer to a cs_enforcement_param_t structure
3148 */
3149 /*----------------------------------------------------------------------------*/
3150
3151 cs_enforcement_param_t *
cs_equation_add_vertex_dof_enforcement(cs_equation_param_t * eqp,cs_lnum_t n_vertices,const cs_lnum_t vertex_ids[],const cs_real_t ref_value[],const cs_real_t vtx_values[])3152 cs_equation_add_vertex_dof_enforcement(cs_equation_param_t *eqp,
3153 cs_lnum_t n_vertices,
3154 const cs_lnum_t vertex_ids[],
3155 const cs_real_t ref_value[],
3156 const cs_real_t vtx_values[])
3157 {
3158 if (eqp == NULL)
3159 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
3160 if (eqp->space_scheme != CS_SPACE_SCHEME_CDOVB &&
3161 eqp->space_scheme != CS_SPACE_SCHEME_CDOVCB)
3162 bft_error(__FILE__, __LINE__, 0,
3163 "%s: Eq: %s: Invalid space scheme.\n"
3164 "This should be a vertex-based one.", __func__, eqp->name);
3165 if (ref_value == NULL && vtx_values == NULL)
3166 bft_error(__FILE__, __LINE__, 0,
3167 "%s: Eq: %s: No enforcement value.\n", __func__, eqp->name);
3168
3169 cs_enforcement_param_t *efp = NULL;
3170 int enforcement_id = eqp->n_enforcements;
3171
3172 eqp->n_enforcements++;
3173
3174 if (vtx_values == NULL) {
3175
3176 assert(ref_value != NULL);
3177 efp = cs_enforcement_param_create(CS_ENFORCEMENT_SELECTION_VERTICES,
3178 CS_ENFORCEMENT_BY_CONSTANT,
3179 eqp->dim,
3180 n_vertices,
3181 vertex_ids,
3182 ref_value);
3183
3184 }
3185 else
3186 efp = cs_enforcement_param_create(CS_ENFORCEMENT_SELECTION_VERTICES,
3187 CS_ENFORCEMENT_BY_DOF_VALUES,
3188 eqp->dim,
3189 n_vertices,
3190 vertex_ids,
3191 vtx_values);
3192
3193 BFT_REALLOC(eqp->enforcement_params, eqp->n_enforcements,
3194 cs_enforcement_param_t);
3195 eqp->enforcement_params[enforcement_id] = efp;
3196 eqp->flag |= CS_EQUATION_FORCE_VALUES;
3197
3198 return efp;
3199 }
3200
3201 /*----------------------------------------------------------------------------*/
3202 /*!
3203 * \brief Add an enforcement of the value of degrees of freedom located at
3204 * the mesh edges.
3205 * The spatial discretization scheme for the given equation has to be
3206 * CDO edge-based schemes.
3207 *
3208 * One assumes that values are interlaced if eqp->dim > 1
3209 * ref_value or elt_values has to be defined. If both parameters are
3210 * defined, one keeps the definition in elt_values
3211 *
3212 * \param[in, out] eqp pointer to a cs_equation_param_t structure
3213 * \param[in] n_edges number of edges to enforce
3214 * \param[in] edge_ids list of edges
3215 * \param[in] ref_value default values or ignored (may be NULL)
3216 * \param[in] edge_values list of associated values, ignored if NULL
3217 *
3218 * \return a pointer to a cs_enforcement_param_t structure
3219 */
3220 /*----------------------------------------------------------------------------*/
3221
3222 cs_enforcement_param_t *
cs_equation_add_edge_dof_enforcement(cs_equation_param_t * eqp,cs_lnum_t n_edges,const cs_lnum_t edge_ids[],const cs_real_t ref_value[],const cs_real_t edge_values[])3223 cs_equation_add_edge_dof_enforcement(cs_equation_param_t *eqp,
3224 cs_lnum_t n_edges,
3225 const cs_lnum_t edge_ids[],
3226 const cs_real_t ref_value[],
3227 const cs_real_t edge_values[])
3228 {
3229 if (eqp == NULL)
3230 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
3231 if (eqp->space_scheme != CS_SPACE_SCHEME_CDOEB)
3232 bft_error(__FILE__, __LINE__, 0,
3233 "%s: Eq: %s: Invalid space scheme.\n"
3234 "This should be a edge-based one.", __func__, eqp->name);
3235 if (ref_value == NULL && edge_values == NULL)
3236 bft_error(__FILE__, __LINE__, 0,
3237 "%s: Eq: %s: No enforcement value.\n", __func__, eqp->name);
3238
3239 cs_enforcement_param_t *efp = NULL;
3240 int enforcement_id = eqp->n_enforcements;
3241
3242 eqp->n_enforcements++;
3243
3244 /* Edge-based schemes are related to a vector-valued equation but DoF are
3245 scalar-valued. They are circulation. */
3246
3247 if (edge_values == NULL) {
3248
3249 assert(ref_value != NULL);
3250 efp = cs_enforcement_param_create(CS_ENFORCEMENT_SELECTION_EDGES,
3251 CS_ENFORCEMENT_BY_CONSTANT,
3252 1,
3253 n_edges,
3254 edge_ids,
3255 ref_value);
3256
3257 }
3258 else
3259 efp = cs_enforcement_param_create(CS_ENFORCEMENT_SELECTION_EDGES,
3260 CS_ENFORCEMENT_BY_DOF_VALUES,
3261 1,
3262 n_edges,
3263 edge_ids,
3264 edge_values);
3265
3266 BFT_REALLOC(eqp->enforcement_params, eqp->n_enforcements,
3267 cs_enforcement_param_t);
3268 eqp->enforcement_params[enforcement_id] = efp;
3269 eqp->flag |= CS_EQUATION_FORCE_VALUES;
3270
3271 return efp;
3272 }
3273
3274 /*----------------------------------------------------------------------------*/
3275 /*!
3276 * \brief Add an enforcement of the value of degrees of freedom located at
3277 * the mesh faces.
3278 * The spatial discretization scheme for the given equation has to be
3279 * CDO face-based schemes.
3280 *
3281 * One assumes that values are interlaced if eqp->dim > 1
3282 * ref_value or elt_values has to be defined. If both parameters are
3283 * defined, one keeps the definition in elt_values
3284 *
3285 * \param[in, out] eqp pointer to a cs_equation_param_t structure
3286 * \param[in] n_faces number of faces to enforce
3287 * \param[in] face_ids list of faces
3288 * \param[in] ref_value default values or ignored (may be NULL)
3289 * \param[in] face_values list of associated values, ignored if NULL
3290 *
3291 * \return a pointer to a cs_enforcement_param_t structure
3292 */
3293 /*----------------------------------------------------------------------------*/
3294
3295 cs_enforcement_param_t *
cs_equation_add_face_dof_enforcement(cs_equation_param_t * eqp,cs_lnum_t n_faces,const cs_lnum_t face_ids[],const cs_real_t ref_value[],const cs_real_t face_values[])3296 cs_equation_add_face_dof_enforcement(cs_equation_param_t *eqp,
3297 cs_lnum_t n_faces,
3298 const cs_lnum_t face_ids[],
3299 const cs_real_t ref_value[],
3300 const cs_real_t face_values[])
3301 {
3302 if (eqp == NULL)
3303 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
3304 if (eqp->space_scheme != CS_SPACE_SCHEME_CDOFB &&
3305 eqp->space_scheme != CS_SPACE_SCHEME_HHO_P0)
3306 bft_error(__FILE__, __LINE__, 0,
3307 "%s: Eq: %s: Invalid space scheme.\n"
3308 "This should be a face-based one.", __func__, eqp->name);
3309 if (ref_value == NULL && face_values == NULL)
3310 bft_error(__FILE__, __LINE__, 0,
3311 "%s: Eq: %s: No enforcement value.\n", __func__, eqp->name);
3312
3313 cs_enforcement_param_t *efp = NULL;
3314 int enforcement_id = eqp->n_enforcements;
3315
3316 eqp->n_enforcements++;
3317
3318 if (face_values == NULL) {
3319
3320 assert(ref_value != NULL);
3321 efp = cs_enforcement_param_create(CS_ENFORCEMENT_SELECTION_FACES,
3322 CS_ENFORCEMENT_BY_CONSTANT,
3323 eqp->dim,
3324 n_faces,
3325 face_ids,
3326 ref_value);
3327
3328 }
3329 else
3330 efp = cs_enforcement_param_create(CS_ENFORCEMENT_SELECTION_FACES,
3331 CS_ENFORCEMENT_BY_DOF_VALUES,
3332 eqp->dim,
3333 n_faces,
3334 face_ids,
3335 face_values);
3336
3337 BFT_REALLOC(eqp->enforcement_params, eqp->n_enforcements,
3338 cs_enforcement_param_t);
3339 eqp->enforcement_params[enforcement_id] = efp;
3340 eqp->flag |= CS_EQUATION_FORCE_VALUES;
3341
3342 return efp;
3343 }
3344
3345 /*----------------------------------------------------------------------------*/
3346 /*!
3347 * \brief Add an enforcement of the value related to the degrees of freedom
3348 * associated to the list of selected cells.
3349 *
3350 * One assumes that values are interlaced if eqp->dim > 1
3351 * ref_value or elt_values has to be defined. If both parameters are
3352 * defined, one keeps the definition in elt_values
3353 *
3354 * \param[in, out] eqp pointer to a cs_equation_param_t structure
3355 * \param[in] n_cells number of selected cells
3356 * \param[in] cell_ids list of cell ids
3357 * \param[in] ref_value ignored if NULL
3358 * \param[in] cell_values list of associated values, ignored if NULL
3359 *
3360 * \return a pointer to a cs_enforcement_param_t structure
3361 */
3362 /*----------------------------------------------------------------------------*/
3363
3364 cs_enforcement_param_t *
cs_equation_add_cell_enforcement(cs_equation_param_t * eqp,cs_lnum_t n_cells,const cs_lnum_t cell_ids[],const cs_real_t ref_value[],const cs_real_t cell_values[])3365 cs_equation_add_cell_enforcement(cs_equation_param_t *eqp,
3366 cs_lnum_t n_cells,
3367 const cs_lnum_t cell_ids[],
3368 const cs_real_t ref_value[],
3369 const cs_real_t cell_values[])
3370 {
3371 if (eqp == NULL)
3372 bft_error(__FILE__, __LINE__, 0, "%s: %s\n", __func__, _err_empty_eqp);
3373 if (ref_value == NULL && cell_values == NULL)
3374 bft_error(__FILE__, __LINE__, 0,
3375 "%s: Eq: %s: No enforcement value.\n", __func__, eqp->name);
3376
3377 cs_enforcement_param_t *efp = NULL;
3378 int enforcement_id = eqp->n_enforcements;
3379
3380 eqp->n_enforcements++;
3381
3382 if (cell_values == NULL) {
3383
3384 assert(ref_value != NULL);
3385 efp = cs_enforcement_param_create(CS_ENFORCEMENT_SELECTION_CELLS,
3386 CS_ENFORCEMENT_BY_CONSTANT,
3387 eqp->dim,
3388 n_cells,
3389 cell_ids,
3390 ref_value);
3391
3392 }
3393 else
3394 efp = cs_enforcement_param_create(CS_ENFORCEMENT_SELECTION_CELLS,
3395 CS_ENFORCEMENT_BY_DOF_VALUES,
3396 eqp->dim,
3397 n_cells,
3398 cell_ids,
3399 cell_values);
3400
3401 BFT_REALLOC(eqp->enforcement_params, eqp->n_enforcements,
3402 cs_enforcement_param_t);
3403 eqp->enforcement_params[enforcement_id] = efp;
3404 eqp->flag |= CS_EQUATION_FORCE_VALUES;
3405
3406 return efp;
3407 }
3408
3409 /*----------------------------------------------------------------------------*/
3410
3411 END_C_DECLS
3412