1 /*============================================================================
2 * Functions to handle cs_navsto_param_t structure
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <stdlib.h>
35 #include <string.h>
36 #include <ctype.h>
37
38 #if defined(HAVE_MPI)
39 #include <mpi.h>
40 #endif
41
42 #if defined(HAVE_PETSC)
43 #include <petscversion.h>
44 #endif
45
46 /*----------------------------------------------------------------------------
47 * Local headers
48 *----------------------------------------------------------------------------*/
49
50 #include <bft_mem.h>
51 #include <bft_printf.h>
52
53 #include "cs_base.h"
54 #include "cs_equation.h"
55 #include "cs_log.h"
56
57 /*----------------------------------------------------------------------------
58 * Header for the current file
59 *----------------------------------------------------------------------------*/
60
61 #include "cs_navsto_param.h"
62
63 /*----------------------------------------------------------------------------*/
64
65 BEGIN_C_DECLS
66
67 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
68
69 /*=============================================================================
70 * Local Macro definitions
71 *============================================================================*/
72
73 #define CS_NAVSTO_PARAM_DBG 0
74
75 /*============================================================================
76 * Type definitions
77 *============================================================================*/
78
79 /*============================================================================
80 * Private variables
81 *============================================================================*/
82
83 static const char _err_empty_nsp[] =
84 N_(" %s: Stop setting an empty cs_navsto_param_t structure.\n"
85 " Please check your settings.\n");
86
87 static const char
88 cs_navsto_param_model_name[CS_NAVSTO_N_MODELS][CS_BASE_STRING_LEN] =
89 { N_("Stokes equations"),
90 N_("Oseen equations"),
91 N_("Incompressible Navier-Stokes equations"),
92 };
93
94 static const char
95 cs_navsto_param_coupling_name[CS_NAVSTO_N_COUPLINGS][CS_BASE_STRING_LEN] =
96 { N_("Artificial compressibility algorithm"),
97 N_("Monolithic"),
98 N_("Incremental projection algorithm"),
99 };
100
101 /* Keys to transfer settings from cs_param_navsto_t to cs_equation_param_t */
102
103 static const char
104 _space_scheme_key[CS_SPACE_N_SCHEMES][CS_BASE_STRING_LEN] =
105 { "fv",
106 "cdo_vb",
107 "cdo_vcb",
108 "cdo_eb",
109 "cdo_fb",
110 "hho_p0",
111 "hho_p1",
112 "hho_p2"
113 };
114
115 static const char
116 _time_scheme_key[CS_TIME_N_SCHEMES][CS_BASE_STRING_LEN] =
117 { "steady",
118 "euler_implicit",
119 "euler_explicit",
120 "crank_nicolson",
121 "theta_scheme",
122 "bdf2"
123 };
124
125 static const char
126 _dof_reduction_key[CS_PARAM_N_REDUCTIONS][CS_BASE_STRING_LEN] =
127 { "derham",
128 "average"
129 };
130
131 static const char
132 _quad_type_key[CS_QUADRATURE_N_TYPES][CS_BASE_STRING_LEN] =
133 { "none",
134 "bary",
135 "bary_subdiv",
136 "higher",
137 "highest"
138 };
139
140 static const char
141 _adv_extrap_key[CS_PARAM_N_ADVECTION_EXTRAPOLATIONS][CS_BASE_STRING_LEN] =
142 {
143 "none",
144 "taylor",
145 "adams_bashforth"
146 };
147
148 static const char
149 _adv_formulation_key[CS_PARAM_N_ADVECTION_FORMULATIONS][CS_BASE_STRING_LEN] =
150 {
151 "conservative",
152 "non_conservative",
153 "skew_symmetric"
154 };
155
156 static const char
157 _adv_scheme_key[CS_PARAM_N_ADVECTION_SCHEMES][CS_BASE_STRING_LEN] =
158 {
159 "centered",
160 "cip",
161 "cip_cw",
162 "hybrid_centered_upwind",
163 "samarskii",
164 "sg",
165 "upwind"
166 };
167
168 static const char
169 _adv_strategy_key[CS_PARAM_N_ADVECTION_STRATEGIES][CS_BASE_STRING_LEN] =
170 {
171 "fully_implicit",
172 "linearized",
173 "explicit"
174 };
175
176 /* scaling coefficient used in Notay's transformation devised in
177 * "Algebraic multigrid for Stokes equations" SIAM J. Sci. Comput. Vol. 39 (5),
178 * 2017 */
179 static double cs_navsto_param_notay_scaling = 1.0;
180
181 /*============================================================================
182 * Private function prototypes
183 *============================================================================*/
184
185 /*----------------------------------------------------------------------------*/
186 /*!
187 * \brief Check if the prerequisite are fullfilled when a PETSC-related type
188 * of sles strategy is requested
189 *
190 * \param[in] val keyval
191 * \param[in] sles_type type of SLES strategy
192 *
193 * \return the same sles_type if ok
194 */
195 /*----------------------------------------------------------------------------*/
196
197 static inline cs_navsto_sles_t
_check_petsc_strategy(const char * val,cs_navsto_sles_t sles_type)198 _check_petsc_strategy(const char *val,
199 cs_navsto_sles_t sles_type)
200 {
201 #if defined(HAVE_PETSC)
202 #if PETSC_VERSION_GE(3,11,0)
203 CS_UNUSED(val);
204 return sles_type;
205 #else
206 if (sles_type == CS_NAVSTO_SLES_GKB_GMRES ||
207 sles_type == CS_NAVSTO_SLES_GKB_PETSC)
208 bft_error(__FILE__, __LINE__, 0,
209 "%s: PETSc version greater or equal to 3.11 is required"
210 " when using the keyval \"%s\"\n", __func__, val);
211 return sles_type;
212 #endif
213 #else
214 bft_error(__FILE__, __LINE__, 0,
215 " %s: \"CS_NSKEY_SLES_STRATEGY\" keyval %s requires"
216 " an installation with PETSC\n", __func__, val);
217 return CS_NAVSTO_SLES_N_TYPES;
218 #endif
219 }
220
221 /*----------------------------------------------------------------------------*/
222 /*!
223 * \brief Retrieve the \ref cs_equation_param_t structure related to the
224 * momentum equation according to the type of coupling
225 *
226 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
227 *
228 * \return a pointer to the corresponding \ref cs_equation_param_t structure
229 */
230 /*----------------------------------------------------------------------------*/
231
232 static inline cs_equation_param_t *
_get_momentum_param(cs_navsto_param_t * nsp)233 _get_momentum_param(cs_navsto_param_t *nsp)
234 {
235 switch (nsp->coupling) {
236
237 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
238 case CS_NAVSTO_COUPLING_MONOLITHIC:
239 return cs_equation_param_by_name("momentum");
240
241 case CS_NAVSTO_COUPLING_PROJECTION:
242 return cs_equation_param_by_name("velocity_prediction");
243 break;
244
245 default:
246 return NULL;
247
248 } /* Switch */
249 }
250
251 /*----------------------------------------------------------------------------*/
252 /*!
253 * \brief Retrieve the \ref cs_equation_param_t structure related to the
254 * momentum equation according to the type of coupling
255 *
256 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
257 *
258 * \return a pointer to the corresponding \ref cs_equation_param_t structure
259 */
260 /*----------------------------------------------------------------------------*/
261
262 static void
_propagate_qtype(cs_navsto_param_t * nsp)263 _propagate_qtype(cs_navsto_param_t *nsp)
264 {
265 /* Loop on velocity ICs */
266
267 for (int i = 0; i < nsp->n_velocity_ic_defs; i++)
268 cs_xdef_set_quadrature(nsp->velocity_ic_defs[i], nsp->qtype);
269
270 /* Loop on pressure ICs */
271
272 for (int i = 0; i < nsp->n_pressure_ic_defs; i++)
273 cs_xdef_set_quadrature(nsp->pressure_ic_defs[i], nsp->qtype);
274
275 /* Loop on velocity BCs */
276
277 for (int i = 0; i < nsp->n_velocity_bc_defs; i++)
278 cs_xdef_set_quadrature(nsp->velocity_bc_defs[i], nsp->qtype);
279
280 /* Loop on pressure BCs */
281
282 for (int i = 0; i < nsp->n_pressure_bc_defs; i++)
283 cs_xdef_set_quadrature(nsp->pressure_bc_defs[i], nsp->qtype);
284 }
285
286
287 /*----------------------------------------------------------------------------*/
288 /*!
289 * \brief Create a new structure to store the settings related to the
290 * resolution of the Navier-Stokes (NS) system: linear or non-linear
291 * algorithms
292 *
293 * \param[in] model type of model related to the NS system
294 * \param[in] model_flag additional high-level model options
295 * \param[in] algo_coupling algorithm used for solving the NS system
296 *
297 * \return a pointer to a new allocated structure
298 */
299 /*----------------------------------------------------------------------------*/
300
301 static cs_navsto_param_sles_t *
_navsto_param_sles_create(cs_navsto_param_model_t model,cs_navsto_param_model_flag_t model_flag,cs_navsto_param_coupling_t algo_coupling)302 _navsto_param_sles_create(cs_navsto_param_model_t model,
303 cs_navsto_param_model_flag_t model_flag,
304 cs_navsto_param_coupling_t algo_coupling)
305 {
306 CS_UNUSED(model_flag);
307
308 cs_navsto_param_sles_t *nslesp = NULL;
309 BFT_MALLOC(nslesp, 1, cs_navsto_param_sles_t);
310
311 /* Set the non-linear algorithm (only useful if advection is implicit and
312 * Navier-Stokes is set) */
313
314 nslesp->nl_algo_type = CS_PARAM_NL_ALGO_PICARD;
315 nslesp->nl_algo_param.n_max_algo_iter = 25;
316 nslesp->nl_algo_param.rtol = 1e-5;
317 nslesp->nl_algo_param.atol = 1e-5;
318 nslesp->nl_algo_param.dtol = 1e3;
319 nslesp->nl_algo_param.verbosity = 1;
320
321 nslesp->anderson_param.n_max_dir = 6;
322 nslesp->anderson_param.starting_iter = 3;
323 nslesp->anderson_param.max_cond = -1; /* No test by default */
324 nslesp->anderson_param.beta = 1.0; /* No damping by default */
325 nslesp->anderson_param.dp_type = CS_PARAM_DOTPROD_EUCLIDEAN;
326
327 /* Set the default solver options for the main linear algorithm */
328
329 nslesp->il_algo_param.n_max_algo_iter = 100;
330 nslesp->il_algo_param.rtol = 1e-08;
331 nslesp->il_algo_param.atol = 1e-08;
332 nslesp->il_algo_param.dtol = 1e3;
333 nslesp->il_algo_param.verbosity = 0;
334
335 nslesp->il_algo_restart = 10;
336
337 switch (algo_coupling) {
338
339 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
340 nslesp->strategy = CS_NAVSTO_SLES_EQ_WITHOUT_BLOCK;
341 break;
342
343 case CS_NAVSTO_COUPLING_MONOLITHIC:
344 if (model == CS_NAVSTO_MODEL_STOKES)
345 nslesp->strategy = CS_NAVSTO_SLES_UZAWA_CG;
346 else
347 nslesp->strategy = CS_NAVSTO_SLES_UZAWA_AL;
348 break;
349
350 case CS_NAVSTO_COUPLING_PROJECTION:
351 nslesp->strategy = CS_NAVSTO_SLES_EQ_WITHOUT_BLOCK;
352 break;
353
354 default:
355 /* Nothing done */
356 break;
357
358 }
359
360 nslesp->schur_approximation = CS_PARAM_SCHUR_LUMPED_INVERSE;
361
362 /* Settings for driving the linear algebra related to the Schur complement
363 approximation */
364
365 cs_param_sles_t *schur_slesp = cs_param_sles_create(-1,
366 "schur_approximation");
367
368 schur_slesp->precond = CS_PARAM_PRECOND_AMG; /* preconditioner */
369 schur_slesp->solver = CS_PARAM_ITSOL_CG; /* iterative solver */
370 schur_slesp->amg_type = CS_PARAM_AMG_HOUSE_K; /* no predefined AMG type */
371 schur_slesp->eps = 1e-4; /* relative tolerance to stop
372 an iterative solver */
373
374 nslesp->schur_sles_param = schur_slesp;
375
376 return nslesp;
377 }
378
379 /*----------------------------------------------------------------------------*/
380 /*!
381 * \brief Free a cs_navsto_param_sles_t structure
382 *
383 * \param[in] p_nslesp double pointer on the cs_navsto_param_sles_t to free
384 */
385 /*----------------------------------------------------------------------------*/
386
387 static void
_navsto_param_sles_free(cs_navsto_param_sles_t ** p_nslesp)388 _navsto_param_sles_free(cs_navsto_param_sles_t **p_nslesp)
389 {
390 if (p_nslesp == NULL)
391 return;
392
393 cs_navsto_param_sles_t *nslesp = *p_nslesp;
394 if (nslesp == NULL)
395 return;
396
397 cs_param_sles_t *schur_slesp = nslesp->schur_sles_param;
398
399 cs_sles_t *schur_sles = cs_sles_find(-1, schur_slesp->name);
400 cs_sles_free(schur_sles);
401
402 cs_param_sles_free(&schur_slesp);
403
404 BFT_FREE(nslesp);
405
406 *p_nslesp = NULL;
407 }
408
409 /*----------------------------------------------------------------------------*/
410 /*!
411 * \brief Log the settings related to the resolution of the Navier-Stokes (NS)
412 * system: linear algorithms only (non-linear algo. depends on the
413 * choice of the model and is this treated in cs_navsto_param_log())
414 *
415 * \param[in] nslesp pointer to the cs_navsto_param_sles_t structure to log
416 */
417 /*----------------------------------------------------------------------------*/
418
419 static void
_navsto_param_sles_log(const cs_navsto_param_sles_t * nslesp)420 _navsto_param_sles_log(const cs_navsto_param_sles_t *nslesp)
421 {
422 if (nslesp == NULL)
423 return;
424
425 const char navsto[16] = " * NavSto |";
426
427 cs_log_printf(CS_LOG_SETUP, "%s SLES strategy: ", navsto);
428 switch (nslesp->strategy) {
429
430 case CS_NAVSTO_SLES_ADDITIVE_GMRES_BY_BLOCK:
431 cs_log_printf(CS_LOG_SETUP, "Additive block preconditioner + GMRES\n");
432 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
433 nslesp->il_algo_restart);
434 break;
435 case CS_NAVSTO_SLES_BLOCK_MULTIGRID_CG:
436 cs_log_printf(CS_LOG_SETUP, "Block AMG + CG\n");
437 break;
438 case CS_NAVSTO_SLES_BY_BLOCKS:
439 cs_log_printf(CS_LOG_SETUP,
440 "Blocks for velocity and pressure (deprecated)\n");
441 break;
442 case CS_NAVSTO_SLES_DIAG_SCHUR_GCR:
443 cs_log_printf(CS_LOG_SETUP, "Diag. block preconditioner with Schur approx."
444 " + (in-house) GCR\n");
445 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
446 nslesp->il_algo_restart);
447 cs_log_printf(CS_LOG_SETUP, "%s Schur approximation: %s\n", navsto,
448 cs_param_get_schur_approx_name(nslesp->schur_approximation));
449 break;
450 case CS_NAVSTO_SLES_DIAG_SCHUR_GMRES:
451 cs_log_printf(CS_LOG_SETUP, "Diag. block preconditioner with Schur approx."
452 " + GMRES\n");
453 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
454 nslesp->il_algo_restart);
455 break;
456 case CS_NAVSTO_SLES_DIAG_SCHUR_MINRES:
457 cs_log_printf(CS_LOG_SETUP, "Diag. block preconditioner with Schur approx."
458 " + in-house MINRES\n");
459 cs_log_printf(CS_LOG_SETUP, "%s Schur approximation: %s\n", navsto,
460 cs_param_get_schur_approx_name(nslesp->schur_approximation));
461 break;
462 case CS_NAVSTO_SLES_EQ_WITHOUT_BLOCK:
463 cs_log_printf(CS_LOG_SETUP, "Handle the full system as it is.\n");
464 break;
465 case CS_NAVSTO_SLES_GCR:
466 cs_log_printf(CS_LOG_SETUP, "in-house GCR\n");
467 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
468 nslesp->il_algo_restart);
469 break;
470 case CS_NAVSTO_SLES_GKB_PETSC:
471 cs_log_printf(CS_LOG_SETUP, "GKB algorithm (through PETSc)\n");
472 break;
473 case CS_NAVSTO_SLES_GKB_GMRES:
474 cs_log_printf(CS_LOG_SETUP, "GMRES with a GKB preconditioner\n");
475 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
476 nslesp->il_algo_restart);
477 break;
478 case CS_NAVSTO_SLES_GKB_SATURNE:
479 cs_log_printf(CS_LOG_SETUP, "in-house GKB algorithm\n");
480 break;
481 case CS_NAVSTO_SLES_LOWER_SCHUR_GCR:
482 cs_log_printf(CS_LOG_SETUP, "Lower block preconditioner with Schur approx."
483 " + (in-house) GCR\n");
484 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
485 nslesp->il_algo_restart);
486 cs_log_printf(CS_LOG_SETUP, "%s Schur approximation: %s\n", navsto,
487 cs_param_get_schur_approx_name(nslesp->schur_approximation));
488 break;
489 case CS_NAVSTO_SLES_MINRES:
490 cs_log_printf(CS_LOG_SETUP, "in-house MINRES\n");
491 break;
492 case CS_NAVSTO_SLES_MULTIPLICATIVE_GMRES_BY_BLOCK:
493 cs_log_printf(CS_LOG_SETUP, "Multiplicative block preconditioner"
494 " + GMRES\n");
495 break;
496 case CS_NAVSTO_SLES_MUMPS:
497 cs_log_printf(CS_LOG_SETUP, "LU factorization with MUMPS\n");
498 break;
499 case CS_NAVSTO_SLES_NOTAY_TRANSFORM:
500 cs_log_printf(CS_LOG_SETUP,
501 "Notay's transformation (scaling=%4.2e) -- Experimental\n",
502 cs_navsto_param_notay_scaling);
503 break;
504 case CS_NAVSTO_SLES_SGS_SCHUR_GCR:
505 cs_log_printf(CS_LOG_SETUP, "Symmetric Gauss-Seidel block preconditioner"
506 " with Schur approx. + (in-house) GCR\n");
507 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
508 nslesp->il_algo_restart);
509 cs_log_printf(CS_LOG_SETUP, "%s Schur approximation: %s\n", navsto,
510 cs_param_get_schur_approx_name(nslesp->schur_approximation));
511 break;
512 case CS_NAVSTO_SLES_UPPER_SCHUR_GMRES:
513 cs_log_printf(CS_LOG_SETUP, "Upper block preconditioner with Schur approx."
514 " + GMRES\n");
515 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
516 nslesp->il_algo_restart);
517 break;
518 case CS_NAVSTO_SLES_UPPER_SCHUR_GCR:
519 cs_log_printf(CS_LOG_SETUP, "Upper block preconditioner with Schur approx."
520 " + (in-house) GCR\n");
521 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
522 nslesp->il_algo_restart);
523 cs_log_printf(CS_LOG_SETUP, "%s Schur approximation: %s\n", navsto,
524 cs_param_get_schur_approx_name(nslesp->schur_approximation));
525 break;
526 case CS_NAVSTO_SLES_USER:
527 cs_log_printf(CS_LOG_SETUP, "User-defined\n");
528 break;
529 case CS_NAVSTO_SLES_UZAWA_AL:
530 cs_log_printf(CS_LOG_SETUP, "Augmented Lagrangian-Uzawa\n");
531 break;
532 case CS_NAVSTO_SLES_UZAWA_CG:
533 cs_log_printf(CS_LOG_SETUP, "Uzawa accelerated with Conjugate Gradient\n");
534 cs_log_printf(CS_LOG_SETUP, "%s Schur approximation: %s\n", navsto,
535 cs_param_get_schur_approx_name(nslesp->schur_approximation));
536 break;
537 case CS_NAVSTO_SLES_UZAWA_SCHUR_GCR:
538 cs_log_printf(CS_LOG_SETUP, "Uzawa block preconditioner with Schur approx."
539 " + (in-house) GCR\n");
540 cs_log_printf(CS_LOG_SETUP, "%s Restart threshold: %d\n", navsto,
541 nslesp->il_algo_restart);
542 cs_log_printf(CS_LOG_SETUP, "%s Schur approximation: %s\n", navsto,
543 cs_param_get_schur_approx_name(nslesp->schur_approximation));
544 break;
545
546 default:
547 cs_log_printf(CS_LOG_SETUP, "Not set\n");
548 break;
549 }
550
551 cs_log_printf(CS_LOG_SETUP, "%s Tolerances of inner linear algo:"
552 " rtol: %5.3e; atol: %5.3e; dtol: %5.3e; verbosity: %d\n",
553 navsto, nslesp->il_algo_param.rtol, nslesp->il_algo_param.atol,
554 nslesp->il_algo_param.dtol, nslesp->il_algo_param.verbosity);
555 cs_log_printf(CS_LOG_SETUP, "%s Max of inner-linear iterations: %d\n",
556 navsto, nslesp->il_algo_param.n_max_algo_iter);
557
558 /* Additional settings for the Schur complement solver */
559 if (nslesp->strategy == CS_NAVSTO_SLES_UZAWA_CG ||
560 nslesp->strategy == CS_NAVSTO_SLES_DIAG_SCHUR_MINRES ||
561 nslesp->strategy == CS_NAVSTO_SLES_DIAG_SCHUR_GCR ||
562 nslesp->strategy == CS_NAVSTO_SLES_LOWER_SCHUR_GCR ||
563 nslesp->strategy == CS_NAVSTO_SLES_SGS_SCHUR_GCR ||
564 nslesp->strategy == CS_NAVSTO_SLES_UPPER_SCHUR_GCR ||
565 nslesp->strategy == CS_NAVSTO_SLES_UZAWA_SCHUR_GCR)
566 cs_param_sles_log(nslesp->schur_sles_param);
567
568 }
569
570 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
571
572 /*============================================================================
573 * Public function prototypes
574 *============================================================================*/
575
576 /*----------------------------------------------------------------------------*/
577 /*!
578 * \brief Retrieve the scaling coefficient used in the Notay's transformation
579 * devised in "Algebraic multigrid for Stokes equations" SIAM
580 * J. Sci. Comput. Vol. 39 (5), 2017
581 * In this article, this scaling is denoted by alpha
582 *
583 * \return the value of the scaling coefficient
584 */
585 /*----------------------------------------------------------------------------*/
586
587 double
cs_navsto_param_get_notay_scaling(void)588 cs_navsto_param_get_notay_scaling(void)
589 {
590 return cs_navsto_param_notay_scaling;
591 }
592
593 /*----------------------------------------------------------------------------*/
594 /*!
595 * \brief Set the scaling coefficient used in the Notay's transformation
596 * devised in "Algebraic multigrid for Stokes equations" SIAM
597 * J. Sci. Comput. Vol. 39 (5), 2017
598 * In this article, this scaling is denoted by alpha
599 *
600 * \param[in] scaling_coef valued of the scaling coefficient
601 */
602 /*----------------------------------------------------------------------------*/
603
604 void
cs_navsto_param_set_notay_scaling(double scaling_coef)605 cs_navsto_param_set_notay_scaling(double scaling_coef)
606 {
607 cs_navsto_param_notay_scaling = scaling_coef;
608 }
609
610 /*----------------------------------------------------------------------------*/
611 /*!
612 * \brief Create a new structure to store all numerical parameters related
613 * to the resolution of the Navier-Stokes (NS) system
614 *
615 * \param[in] boundaries pointer to a cs_boundary_t structure
616 * \param[in] model type of model related to the NS system
617 * \param[in] model_flag additional high-level model options
618 * \param[in] algo_coupling algorithm used for solving the NS system
619 * \param[in] post_flag predefined post-processings options
620 *
621 * \return a pointer to a new allocated structure
622 */
623 /*----------------------------------------------------------------------------*/
624
625 cs_navsto_param_t *
cs_navsto_param_create(const cs_boundary_t * boundaries,cs_navsto_param_model_t model,cs_navsto_param_model_flag_t model_flag,cs_navsto_param_coupling_t algo_coupling,cs_navsto_param_post_flag_t post_flag)626 cs_navsto_param_create(const cs_boundary_t *boundaries,
627 cs_navsto_param_model_t model,
628 cs_navsto_param_model_flag_t model_flag,
629 cs_navsto_param_coupling_t algo_coupling,
630 cs_navsto_param_post_flag_t post_flag)
631 {
632 cs_navsto_param_t *nsp = NULL;
633 BFT_MALLOC(nsp, 1, cs_navsto_param_t);
634
635 /* Physical modelling */
636 /* ------------------ */
637
638 /* Which equations are solved and which terms are needed */
639
640 nsp->model = model;
641 nsp->model_flag = model_flag;
642
643 /* Turbulence modelling (pointer to global structures) */
644
645 nsp->turbulence = cs_turbulence_param_create();
646
647 /* Main set of properties */
648 /* ---------------------- */
649
650 nsp->phys_constants = cs_get_glob_physical_constants();
651
652 nsp->mass_density = cs_property_by_name(CS_PROPERTY_MASS_DENSITY);
653 if (nsp->mass_density == NULL)
654 nsp->mass_density = cs_property_add(CS_PROPERTY_MASS_DENSITY,
655 CS_PROPERTY_ISO);
656
657 nsp->lam_viscosity = cs_property_add(CS_NAVSTO_LAM_VISCOSITY,
658 CS_PROPERTY_ISO);
659
660 if (nsp->turbulence->model->iturb == CS_TURB_NONE)
661 nsp->tot_viscosity = nsp->lam_viscosity;
662 else
663 nsp->tot_viscosity = cs_property_add(CS_NAVSTO_TOTAL_VISCOSITY,
664 CS_PROPERTY_ISO);
665
666 /* Default numerical settings */
667 /* -------------------------- */
668
669 nsp->dof_reduction_mode = CS_PARAM_REDUCTION_AVERAGE;
670 nsp->coupling = algo_coupling;
671 nsp->space_scheme = CS_SPACE_SCHEME_CDOFB;
672
673 /* Advection settings */
674
675 nsp->adv_form = CS_PARAM_ADVECTION_FORM_NONCONS;
676 nsp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_UPWIND;
677 nsp->adv_strategy = CS_PARAM_ADVECTION_IMPLICIT_FULL;
678 nsp->adv_extrapol = CS_PARAM_ADVECTION_EXTRAPOL_NONE;
679
680 /* Forcing steady state in order to avoid inconsistencies */
681
682 if (model_flag & CS_NAVSTO_MODEL_STEADY)
683 nsp->time_scheme = CS_TIME_SCHEME_STEADY;
684 else
685 nsp->time_scheme = CS_TIME_SCHEME_EULER_IMPLICIT;
686 nsp->theta = 1.0;
687
688 /* Boussinesq term(s) */
689
690 nsp->n_boussinesq_terms = 0;
691 nsp->boussinesq_param = NULL;
692
693 /* Default level of quadrature */
694
695 nsp->qtype = CS_QUADRATURE_BARY;
696
697 /* Resolution parameters (inner linear system then the non-linear system )*/
698
699 nsp->sles_param = _navsto_param_sles_create(model, model_flag, algo_coupling);
700
701 /* Management of the outer resolution steps (i.e. the full system including
702 the turbulence modelling or the the thermal system) */
703
704 nsp->n_max_outer_iter = 5;
705 nsp->delta_thermal_tolerance = 1e-2;
706
707 /* Output indicators */
708
709 nsp->verbosity = 1;
710 nsp->post_flag = post_flag;
711
712 /* Set no augmentation of the linear system. This could be modified according
713 to the type of coupling or the strategy used to solved to the linear
714 system. */
715
716 nsp->gd_scale_coef = 0.0;
717
718 /* Initial conditions
719 * ------------------
720 *
721 * Remark: As velocity and pressure may not be associated to an equation
722 * directly, one stores the definition of initial conditions and boundary
723 * conditions at this level */
724
725 switch (algo_coupling) {
726
727 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
728 nsp->gd_scale_coef = 1.0;
729
730 nsp->velocity_ic_is_owner = false;
731 nsp->velocity_bc_is_owner = false;
732 nsp->pressure_ic_is_owner = true;
733 nsp->pressure_bc_is_owner = true;
734 break;
735
736 case CS_NAVSTO_COUPLING_MONOLITHIC:
737 if (model != CS_NAVSTO_MODEL_STOKES)
738 /* The default strategy when one does not solve the Stokes equations is
739 * set in _navsto_param_sles_create() which is CS_NAVSTO_SLES_UZAWA_AL.
740 * Thus, one adds a slight augmentation of the linear system */
741
742 nsp->gd_scale_coef = 1.0;
743
744 nsp->velocity_ic_is_owner = false;
745 nsp->velocity_bc_is_owner = false;
746 nsp->pressure_ic_is_owner = true;
747 nsp->pressure_bc_is_owner = true;
748 break;
749
750 case CS_NAVSTO_COUPLING_PROJECTION:
751 nsp->velocity_ic_is_owner = false;
752 nsp->velocity_bc_is_owner = false;
753 nsp->pressure_ic_is_owner = false;
754 nsp->pressure_bc_is_owner = false;
755 break;
756
757 default:
758 /* Nothing done */
759 break;
760 }
761
762 /* Initial conditions for the pressure field */
763
764 nsp->n_pressure_ic_defs = 0;
765 nsp->pressure_ic_defs = NULL;
766
767 /* Initial conditions for the velocity field */
768
769 nsp->n_velocity_ic_defs = 0;
770 nsp->velocity_ic_defs = NULL;
771
772 /* Boundary conditions */
773 /* ------------------- */
774
775 /* Physical boundaries specific to the problem at stake */
776
777 nsp->boundaries = boundaries; /* shared structure */
778
779 /* Boundary conditions for the pressure field */
780
781 nsp->n_pressure_bc_defs = 0;
782 nsp->pressure_bc_defs = NULL;
783
784 /* Boundary conditions for the velocity field */
785
786 nsp->n_velocity_bc_defs = 0;
787 nsp->velocity_bc_defs = NULL;
788
789 /* Other conditions */
790 /* ---------------- */
791
792 /* Rescaling of the pressure */
793
794 nsp->reference_pressure = 0.;
795
796 return nsp;
797 }
798
799 /*----------------------------------------------------------------------------*/
800 /*!
801 * \brief Free a \ref cs_navsto_param_t structure
802 *
803 * \param[in, out] param pointer to a \ref cs_navsto_param_t structure
804 *
805 * \return a NULL pointer
806 */
807 /*----------------------------------------------------------------------------*/
808
809 cs_navsto_param_t *
cs_navsto_param_free(cs_navsto_param_t * param)810 cs_navsto_param_free(cs_navsto_param_t *param)
811 {
812 if (param == NULL)
813 return param;
814
815 /* Turbulence modelling */
816
817 BFT_FREE(param->turbulence);
818
819 /* Boussinesq term(s) */
820
821 if (param->n_boussinesq_terms > 0)
822 BFT_FREE(param->boussinesq_param);
823
824 /* Velocity initial conditions */
825
826 if (param->n_velocity_ic_defs > 0) {
827
828 /* Otherwise this is freed inside the related equation */
829 if (param->velocity_ic_is_owner) {
830 for (int i = 0; i < param->n_velocity_ic_defs; i++)
831 param->velocity_ic_defs[i] = cs_xdef_free(param->velocity_ic_defs[i]);
832 }
833 BFT_FREE(param->velocity_ic_defs);
834 param->velocity_ic_defs = NULL;
835
836 }
837
838 /* Pressure initial conditions */
839
840 if (param->n_pressure_ic_defs > 0) {
841
842 if (param->pressure_ic_is_owner) {
843 for (int i = 0; i < param->n_pressure_ic_defs; i++)
844 param->pressure_ic_defs[i] = cs_xdef_free(param->pressure_ic_defs[i]);
845 }
846 BFT_FREE(param->pressure_ic_defs);
847 param->pressure_ic_defs = NULL;
848
849 }
850
851 /* Velocity boundary conditions */
852
853 if (param->n_velocity_bc_defs > 0) {
854 if (param->velocity_bc_is_owner) {
855 for (int i = 0; i < param->n_velocity_bc_defs; i++)
856 param->velocity_bc_defs[i] = cs_xdef_free(param->velocity_bc_defs[i]);
857 }
858 BFT_FREE(param->velocity_bc_defs);
859 param->velocity_bc_defs = NULL;
860 }
861
862 /* Pressure boundary conditions */
863
864 if (param->n_pressure_bc_defs > 0) {
865
866 if (param->pressure_bc_is_owner) {
867 for (int i = 0; i < param->n_pressure_bc_defs; i++)
868 param->pressure_bc_defs[i] = cs_xdef_free(param->pressure_bc_defs[i]);
869 }
870 BFT_FREE(param->pressure_bc_defs);
871 param->pressure_bc_defs = NULL;
872
873 }
874
875 _navsto_param_sles_free(&(param->sles_param));
876
877 /* Free the main structure */
878 BFT_FREE(param);
879
880 return NULL;
881 }
882
883 /*----------------------------------------------------------------------------*/
884 /*!
885 * \brief Set a parameter attached to a keyname in a \ref cs_navsto_param_t
886 * structure
887 *
888 * \param[in, out] nsp pointer to a \ref cs_navsto_param_t structure to set
889 * \param[in] key key related to the member of eq to set
890 * \param[in] keyval accessor to the value to set
891 */
892 /*----------------------------------------------------------------------------*/
893
894 void
cs_navsto_param_set(cs_navsto_param_t * nsp,cs_navsto_key_t key,const char * keyval)895 cs_navsto_param_set(cs_navsto_param_t *nsp,
896 cs_navsto_key_t key,
897 const char *keyval)
898 {
899 if (nsp == NULL)
900 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
901 if (keyval == NULL)
902 bft_error(__FILE__, __LINE__, 0, "%s: Empty key value.\n", __func__);
903
904 /* Conversion of the string to lower case */
905 char val[CS_BASE_STRING_LEN];
906 for (size_t i = 0; i < strlen(keyval); i++)
907 val[i] = tolower(keyval[i]);
908 val[strlen(keyval)] = '\0';
909
910 switch(key) {
911
912 case CS_NSKEY_ADVECTION_EXTRAPOL:
913 if (strcmp(val, "none") == 0)
914 nsp->adv_extrapol = CS_PARAM_ADVECTION_EXTRAPOL_NONE;
915 else if (strcmp(val, "taylor") == 0)
916 nsp->adv_extrapol = CS_PARAM_ADVECTION_EXTRAPOL_TAYLOR_2;
917 else if (strcmp(val, "adams_bashforth") == 0)
918 nsp->adv_extrapol = CS_PARAM_ADVECTION_EXTRAPOL_ADAMS_BASHFORTH_2;
919 else {
920 const char *_val = val;
921 bft_error(__FILE__, __LINE__, 0,
922 _(" %s: Invalid val %s related to key"
923 " CS_NSKEY_ADVECTION_EXTRAPOL\n"), __func__, _val);
924 }
925 break;
926
927 case CS_NSKEY_ADVECTION_FORMULATION:
928 if (strcmp(val, "conservative") == 0)
929 nsp->adv_form = CS_PARAM_ADVECTION_FORM_CONSERV;
930 else if (strcmp(val, "non_conservative") == 0)
931 nsp->adv_form = CS_PARAM_ADVECTION_FORM_NONCONS;
932 else if (strcmp(val, "skew_symmetric") == 0)
933 nsp->adv_form = CS_PARAM_ADVECTION_FORM_SKEWSYM;
934 else {
935 const char *_val = val;
936 bft_error(__FILE__, __LINE__, 0,
937 _(" %s: Invalid val %s related to key"
938 " CS_NSKEY_ADVECTION_FORMULATION\n"
939 " Choice between conservative, non_conservative"),
940 __func__, _val);
941 }
942 break;
943
944 case CS_NSKEY_ADVECTION_SCHEME:
945 if (strcmp(val, "upwind") == 0)
946 nsp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_UPWIND;
947 else if (strcmp(val, "samarskii") == 0)
948 nsp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_SAMARSKII;
949 else if (strcmp(val, "sg") == 0)
950 nsp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_SG;
951 else if (strcmp(val, "centered") == 0)
952 nsp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_CENTERED;
953 else if (strcmp(val, "mix_centered_upwind") == 0 ||
954 strcmp(val, "hybrid_centered_upwind") == 0)
955 nsp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_HYBRID_CENTERED_UPWIND;
956 else if (strcmp(val, "cip") == 0) {
957 nsp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_CIP;
958 /* Automatically switch to a non-conservative formulation */
959 nsp->adv_form = CS_PARAM_ADVECTION_FORM_NONCONS;
960 }
961 else if (strcmp(val, "cip_cw") == 0) {
962 nsp->adv_scheme = CS_PARAM_ADVECTION_SCHEME_CIP_CW;
963 /* Automatically switch to a non-conservative formulation */
964 nsp->adv_form = CS_PARAM_ADVECTION_FORM_NONCONS;
965 }
966 else {
967 const char *_val = val;
968 bft_error(__FILE__, __LINE__, 0,
969 _(" %s: Invalid val %s related to key"
970 " CS_NSKEY_ADVECTION_SCHEME\n"
971 " Choices between upwind, samarskii, sg, centered, cip,"
972 " cip_cw, hybrid_centered_upwind, mix_centered_upwind"),
973 __func__, _val);
974 }
975 break;
976
977 case CS_NSKEY_ADVECTION_STRATEGY:
978 if (strcmp(val, "fully_implicit") == 0 ||
979 strcmp(val, "implicit") == 0)
980 nsp->adv_strategy = CS_PARAM_ADVECTION_IMPLICIT_FULL;
981 else if (strcmp(val, "implicit_linear") == 0 ||
982 strcmp(val, "linearized") == 0)
983 nsp->adv_strategy = CS_PARAM_ADVECTION_IMPLICIT_LINEARIZED;
984 else if (strcmp(val, "explicit") == 0)
985 nsp->adv_strategy = CS_PARAM_ADVECTION_EXPLICIT;
986 else {
987 const char *_val = val;
988 bft_error(__FILE__, __LINE__, 0,
989 _(" %s: Invalid val %s related to key"
990 " CS_NSKEY_ADVECTION_STRATEGY\n"), __func__, _val);
991 }
992 break;
993
994 case CS_NSKEY_DOF_REDUCTION:
995 if (strcmp(val, "derham") == 0)
996 nsp->dof_reduction_mode = CS_PARAM_REDUCTION_DERHAM;
997 else if (strcmp(val, "average") == 0)
998 nsp->dof_reduction_mode = CS_PARAM_REDUCTION_AVERAGE;
999 else {
1000 const char *_val = val;
1001 bft_error(__FILE__, __LINE__, 0,
1002 _(" %s: Invalid val %s related to key CS_NSKEY_DOF_REDUCTION\n"
1003 " Choice between \"derham\" or \"average\"."),
1004 __func__, _val);
1005 }
1006 break;
1007
1008 case CS_NSKEY_GD_SCALE_COEF:
1009 switch (nsp->coupling) {
1010 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
1011 case CS_NAVSTO_COUPLING_MONOLITHIC:
1012 nsp->gd_scale_coef = atof(val);
1013 break;
1014
1015 case CS_NAVSTO_COUPLING_PROJECTION:
1016 cs_base_warn(__FILE__, __LINE__);
1017 bft_printf(" %s: Trying to set the zeta parameter with the %s\n "
1018 " although this will not have use in the algorithm.\n",
1019 __func__, cs_navsto_param_coupling_name[nsp->coupling]);
1020
1021 default:
1022 break;
1023
1024 } /* End of switch */
1025 break;
1026
1027 case CS_NSKEY_IL_ALGO_ATOL:
1028 nsp->sles_param->il_algo_param.atol = atof(val);
1029 if (nsp->sles_param->il_algo_param.rtol < 0)
1030 bft_error(__FILE__, __LINE__, 0,
1031 " %s: Invalid value for the absolute tolerance\n", __func__);
1032 break;
1033
1034 case CS_NSKEY_IL_ALGO_DTOL:
1035 nsp->sles_param->il_algo_param.dtol = atof(val);
1036 if (nsp->sles_param->il_algo_param.dtol < 0)
1037 bft_error(__FILE__, __LINE__, 0,
1038 " %s: Invalid value for the divergence tolerance\n", __func__);
1039 break;
1040
1041 case CS_NSKEY_IL_ALGO_RTOL:
1042 nsp->sles_param->il_algo_param.rtol = atof(val);
1043 if (nsp->sles_param->il_algo_param.rtol < 0)
1044 bft_error(__FILE__, __LINE__, 0,
1045 " %s: Invalid value for the residual tolerance\n", __func__);
1046 break;
1047
1048 case CS_NSKEY_IL_ALGO_RESTART:
1049 nsp->sles_param->il_algo_restart = atoi(val);
1050 break;
1051
1052 case CS_NSKEY_IL_ALGO_VERBOSITY:
1053 nsp->sles_param->il_algo_param.verbosity = atoi(val);
1054 break;
1055
1056 case CS_NSKEY_MAX_IL_ALGO_ITER:
1057 nsp->sles_param->il_algo_param.n_max_algo_iter = atoi(val);
1058 break;
1059
1060 case CS_NSKEY_MAX_NL_ALGO_ITER:
1061 nsp->sles_param->nl_algo_param.n_max_algo_iter = atoi(val);
1062 break;
1063
1064 case CS_NSKEY_MAX_OUTER_ITER:
1065 nsp->n_max_outer_iter = atoi(val);
1066 break;
1067
1068 case CS_NSKEY_NL_ALGO:
1069 {
1070 if (strcmp(val, "picard") == 0 || strcmp(val, "fixed-point") == 0)
1071 nsp->sles_param->nl_algo_type = CS_PARAM_NL_ALGO_PICARD;
1072 else if (strcmp(val, "anderson") == 0)
1073 nsp->sles_param->nl_algo_type = CS_PARAM_NL_ALGO_ANDERSON;
1074 else {
1075 const char *_val = val;
1076 bft_error(__FILE__, __LINE__, 0,
1077 " %s: Invalid value \"%s\" for key CS_NSKEY_NL_ALGO\n"
1078 " Valid choices are \"picard\", \"fixed-point\".",
1079 __func__, _val);
1080 }
1081 }
1082 break; /* Non-linear algorithm */
1083
1084 case CS_NSKEY_NL_ALGO_ATOL:
1085 nsp->sles_param->nl_algo_param.atol = atof(val);
1086 if (nsp->sles_param->nl_algo_param.atol < 0)
1087 bft_error(__FILE__, __LINE__, 0,
1088 " %s: Invalid value for the absolute tolerance of the"
1089 " non-linear algorithm\n",
1090 __func__);
1091 break;
1092
1093 case CS_NSKEY_NL_ALGO_DTOL:
1094 nsp->sles_param->nl_algo_param.dtol = atof(val);
1095 if (nsp->sles_param->nl_algo_param.dtol < 0)
1096 bft_error(__FILE__, __LINE__, 0,
1097 " %s: Invalid value for the divergence tolerance of the"
1098 " non-linear algorithm\n",
1099 __func__);
1100 break;
1101
1102 case CS_NSKEY_NL_ALGO_RTOL:
1103 nsp->sles_param->nl_algo_param.rtol = atof(val);
1104 if (nsp->sles_param->nl_algo_param.rtol < 0)
1105 bft_error(__FILE__, __LINE__, 0,
1106 " %s: Invalid value for the relative tolerance of the"
1107 " non-linear algorithm\n",
1108 __func__);
1109 break;
1110
1111 case CS_NSKEY_NL_ALGO_VERBOSITY:
1112 nsp->sles_param->nl_algo_param.verbosity = atoi(val);
1113 break;
1114
1115 case CS_NSKEY_QUADRATURE:
1116 {
1117 nsp->qtype = CS_QUADRATURE_NONE;
1118
1119 if (strcmp(val, "bary") == 0)
1120 nsp->qtype = CS_QUADRATURE_BARY;
1121 else if (strcmp(val, "bary_subdiv") == 0)
1122 nsp->qtype = CS_QUADRATURE_BARY_SUBDIV;
1123 else if (strcmp(val, "higher") == 0)
1124 nsp->qtype = CS_QUADRATURE_HIGHER;
1125 else if (strcmp(val, "highest") == 0)
1126 nsp->qtype = CS_QUADRATURE_HIGHEST;
1127 else {
1128 const char *_val = val;
1129 bft_error(__FILE__, __LINE__, 0,
1130 _(" %s: Invalid value \"%s\" for key CS_NSKEY_QUADRATURE\n"
1131 " Valid choices are \"bary\", \"bary_subdiv\", \"higher\""
1132 " and \"highest\"."), __func__, _val);
1133 }
1134
1135 _propagate_qtype(nsp);
1136 }
1137 break; /* Quadrature */
1138
1139 case CS_NSKEY_SCHUR_STRATEGY:
1140 if (strcmp(val, "diag_schur") == 0)
1141 nsp->sles_param->schur_approximation = CS_PARAM_SCHUR_DIAG_INVERSE;
1142 else if (strcmp(val, "lumped_schur") == 0)
1143 nsp->sles_param->schur_approximation = CS_PARAM_SCHUR_LUMPED_INVERSE;
1144 else if (strcmp(val, "elman") == 0)
1145 nsp->sles_param->schur_approximation = CS_PARAM_SCHUR_ELMAN;
1146 else {
1147 const char *_val = val;
1148 bft_error(__FILE__, __LINE__, 0,
1149 _(" %s: Invalid value \"%s\" not among valid choices:\n"
1150 " \"diag_schur\", \"lumped_schur\", \"elman\"."),
1151 __func__, _val);
1152 }
1153 break;
1154
1155 case CS_NSKEY_SLES_STRATEGY:
1156
1157 /* In-house strategies. Do not need a third-part library. */
1158 /* ------------------------------------------------------ */
1159
1160 if (strcmp(val, "block_amg_cg") == 0)
1161 nsp->sles_param->strategy = CS_NAVSTO_SLES_BLOCK_MULTIGRID_CG;
1162 else if (strcmp(val, "diag_schur_gcr") == 0)
1163 nsp->sles_param->strategy = CS_NAVSTO_SLES_DIAG_SCHUR_GCR;
1164 else if (strcmp(val, "diag_schur_minres") == 0)
1165 nsp->sles_param->strategy = CS_NAVSTO_SLES_DIAG_SCHUR_MINRES;
1166 else if (strcmp(val, "gcr") == 0)
1167 nsp->sles_param->strategy = CS_NAVSTO_SLES_GCR;
1168 else if (strcmp(val, "gkb_saturne") == 0 ||
1169 strcmp(val, "gkb") == 0)
1170 nsp->sles_param->strategy = CS_NAVSTO_SLES_GKB_SATURNE;
1171 else if (strcmp(val, "lower_schur_gcr") == 0)
1172 nsp->sles_param->strategy = CS_NAVSTO_SLES_LOWER_SCHUR_GCR;
1173 else if (strcmp(val, "minres") == 0)
1174 nsp->sles_param->strategy = CS_NAVSTO_SLES_MINRES;
1175 else if (strcmp(val, "no_block") == 0)
1176 nsp->sles_param->strategy = CS_NAVSTO_SLES_EQ_WITHOUT_BLOCK;
1177 else if (strcmp(val, "sgs_schur_gcr") == 0)
1178 nsp->sles_param->strategy = CS_NAVSTO_SLES_SGS_SCHUR_GCR;
1179 else if (strcmp(val, "upper_schur_gcr") == 0)
1180 nsp->sles_param->strategy = CS_NAVSTO_SLES_UPPER_SCHUR_GCR;
1181 else if (strcmp(val, "uzawa_al") == 0 || strcmp(val, "alu") == 0)
1182 nsp->sles_param->strategy = CS_NAVSTO_SLES_UZAWA_AL;
1183 else if (strcmp(val, "uzawa_cg") == 0 || strcmp(val, "uzapcg") == 0)
1184 nsp->sles_param->strategy = CS_NAVSTO_SLES_UZAWA_CG;
1185 else if (strcmp(val, "uza_schur_gcr") == 0 ||
1186 strcmp(val, "uzawa_schur_gcr") == 0)
1187 nsp->sles_param->strategy = CS_NAVSTO_SLES_UZAWA_SCHUR_GCR;
1188
1189 /* All the following options need either PETSC or MUMPS */
1190 /* ---------------------------------------------------- */
1191
1192 else if (strcmp(val, "additive_gmres") == 0)
1193 nsp->sles_param->strategy =
1194 _check_petsc_strategy(val, CS_NAVSTO_SLES_ADDITIVE_GMRES_BY_BLOCK);
1195 else if (strcmp(val, "diag_schur_gmres") == 0)
1196 nsp->sles_param->strategy =
1197 _check_petsc_strategy(val, CS_NAVSTO_SLES_DIAG_SCHUR_GMRES);
1198 else if (strcmp(val, "gkb_gmres") == 0)
1199 nsp->sles_param->strategy =
1200 _check_petsc_strategy(val, CS_NAVSTO_SLES_GKB_GMRES);
1201 else if (strcmp(val, "gkb_petsc") == 0)
1202 nsp->sles_param->strategy =
1203 _check_petsc_strategy(val, CS_NAVSTO_SLES_GKB_PETSC);
1204 else if (strcmp(val, "multiplicative_gmres") == 0)
1205 nsp->sles_param->strategy =
1206 _check_petsc_strategy(val,
1207 CS_NAVSTO_SLES_MULTIPLICATIVE_GMRES_BY_BLOCK);
1208 else if (strcmp(val, "mumps") == 0) {
1209 #if defined(HAVE_MUMPS)
1210 nsp->sles_param->strategy = CS_NAVSTO_SLES_MUMPS;
1211 #else
1212 #if defined(HAVE_PETSC)
1213 #if defined(PETSC_HAVE_MUMPS)
1214 nsp->sles_param->strategy = CS_NAVSTO_SLES_MUMPS;
1215 #else
1216 bft_error(__FILE__, __LINE__, 0,
1217 " %s: Error detected while setting \"%s\" key\n"
1218 " MUMPS is not available with your installation.\n"
1219 " Please check your installation settings.\n",
1220 __func__, "CS_NSKEY_SLES_STRATEGY");
1221 #endif /* PETSC_HAVE_MUMPS */
1222 #endif /* HAVE_PETSC */
1223 #endif /* HAVE_MUMPS */
1224 }
1225 else if (strcmp(val, "notay") == 0)
1226 nsp->sles_param->strategy = CS_NAVSTO_SLES_NOTAY_TRANSFORM;
1227 else if (strcmp(val, "upper_schur_gmres") == 0)
1228 nsp->sles_param->strategy =
1229 _check_petsc_strategy(val, CS_NAVSTO_SLES_UPPER_SCHUR_GMRES);
1230
1231 else {
1232 const char *_val = val;
1233 bft_error(__FILE__, __LINE__, 0,
1234 " %s: Invalid val %s related to key CS_NSKEY_SLES_STRATEGY\n"
1235 " Choice between: no_block, block_amg_cg, minres, gcr,\n"
1236 " {additive,multiplicative}_gmres, {diag,upper}_schur_gmres,\n"
1237 " gkb, gkb_petsc, gkb_gmres, gkb_saturne, mumps, notay\n"
1238 " uzawa_al or alu, uzawa_cg, uzawa_schur_gcr,\n"
1239 " {diag,upper,lower,sgs}_schur_gcr, diag_schur_minres",
1240 __func__, _val);
1241 }
1242 break;
1243
1244 case CS_NSKEY_SPACE_SCHEME:
1245 if (strcmp(val, "cdo_fb") == 0) {
1246 nsp->space_scheme = CS_SPACE_SCHEME_CDOFB;
1247 }
1248 else if (strcmp(val, "hho_p0") == 0) {
1249 nsp->space_scheme = CS_SPACE_SCHEME_HHO_P0;
1250 }
1251 else if (strcmp(val, "hho_p1") == 0) {
1252 nsp->space_scheme = CS_SPACE_SCHEME_HHO_P1;
1253 }
1254 else if (strcmp(val, "hho_p2") == 0) {
1255 nsp->space_scheme = CS_SPACE_SCHEME_HHO_P2;
1256 }
1257 else {
1258 const char *_val = val;
1259 bft_error(__FILE__, __LINE__, 0,
1260 _(" %s: Invalid val %s related to key CS_NSKEY_SPACE_SCHEME\n"
1261 " Choice between hho_{p0, p1, p2} or cdo_fb"),
1262 __func__, _val);
1263 }
1264 break;
1265
1266 case CS_NSKEY_THERMAL_TOLERANCE:
1267 nsp->delta_thermal_tolerance = atof(val);
1268 /* If tolerance is set to a negative value then it stops the outer
1269 iteration process after the first iteration */
1270 break;
1271
1272 case CS_NSKEY_TIME_SCHEME:
1273 if (strcmp(val, "euler_implicit") == 0) {
1274 nsp->time_scheme = CS_TIME_SCHEME_EULER_IMPLICIT;
1275 nsp->theta = 1.;
1276 }
1277 else if (strcmp(val, "euler_explicit") == 0) {
1278 nsp->time_scheme = CS_TIME_SCHEME_EULER_EXPLICIT;
1279 nsp->theta = 0.;
1280 }
1281 else if (strcmp(val, "crank_nicolson") == 0) {
1282 nsp->time_scheme = CS_TIME_SCHEME_CRANKNICO;
1283 nsp->theta = 0.5;
1284 }
1285 else if (strcmp(val, "theta_scheme") == 0)
1286 nsp->time_scheme = CS_TIME_SCHEME_THETA;
1287 else if (strcmp(val, "bdf2") == 0)
1288 nsp->time_scheme = CS_TIME_SCHEME_BDF2;
1289 else {
1290 const char *_val = val;
1291 bft_error(__FILE__, __LINE__, 0,
1292 _(" %s: Invalid value \"%s\" for CS_EQKEY_TIME_SCHEME\n"
1293 " Valid choices are \"euler_implicit\","
1294 " \"euler_explicit\"," " \"crank_nicolson\","
1295 " \"theta_scheme\" and \"bdf2\"."), __func__, _val);
1296 }
1297 break;
1298
1299 case CS_NSKEY_TIME_THETA:
1300 nsp->theta = atof(val);
1301 if (nsp->theta < 0. - cs_math_zero_threshold ||
1302 nsp->theta > 1.0 + cs_math_zero_threshold)
1303 bft_error(__FILE__, __LINE__, 0,
1304 " %s: Invalid value for theta\n", __func__);
1305 break;
1306
1307 case CS_NSKEY_VERBOSITY:
1308 nsp->verbosity = atoi(val);
1309 break;
1310
1311 default:
1312 bft_error(__FILE__, __LINE__, 0,
1313 _(" %s: Invalid key for setting the Navier-Stokes system."),
1314 __func__);
1315
1316 }
1317
1318 }
1319
1320 /*----------------------------------------------------------------------------*/
1321 /*!
1322 * \brief Apply the numerical settings defined for the Navier-Stokes system
1323 * to an equation related to this system.
1324 *
1325 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1326 * \param[in, out] eqp pointer to a \ref cs_equation_param_t structure
1327 */
1328 /*----------------------------------------------------------------------------*/
1329
1330 void
cs_navsto_param_transfer(const cs_navsto_param_t * nsp,cs_equation_param_t * eqp)1331 cs_navsto_param_transfer(const cs_navsto_param_t *nsp,
1332 cs_equation_param_t *eqp)
1333 {
1334 assert(nsp != NULL && eqp != NULL);
1335
1336 /* Set the space discretization scheme */
1337 const char *ss_key = _space_scheme_key[nsp->space_scheme];
1338
1339 cs_equation_param_set(eqp, CS_EQKEY_SPACE_SCHEME, ss_key);
1340
1341 /* Set the time discretization scheme */
1342 const char *ts_key = _time_scheme_key[nsp->time_scheme];
1343
1344 cs_equation_param_set(eqp, CS_EQKEY_TIME_SCHEME, ts_key);
1345 if (nsp->time_scheme == CS_TIME_SCHEME_THETA) {
1346 char cvalue[36]; /* include '\0' */
1347 snprintf(cvalue, 35*sizeof(char), "%g", nsp->theta);
1348 cs_equation_param_set(eqp, CS_EQKEY_TIME_THETA, cvalue);
1349 }
1350
1351 /* Set the way DoFs are defined */
1352 const char *dof_key = _dof_reduction_key[nsp->dof_reduction_mode];
1353
1354 cs_equation_param_set(eqp, CS_EQKEY_DOF_REDUCTION, dof_key);
1355
1356 /* Set quadratures type */
1357 const char *quad_key = _quad_type_key[nsp->qtype];
1358
1359 /* If requested, add advection parameters */
1360 if ((nsp->model & (CS_NAVSTO_MODEL_INCOMPRESSIBLE_NAVIER_STOKES |
1361 CS_NAVSTO_MODEL_OSEEN)) > 0) {
1362
1363 /* If different from default value */
1364 const char *extrap_key = _adv_extrap_key[nsp->adv_extrapol];
1365 cs_equation_param_set(eqp, CS_EQKEY_ADV_EXTRAPOL, extrap_key);
1366
1367 const char *stra_key = _adv_strategy_key[nsp->adv_strategy];
1368 cs_equation_param_set(eqp, CS_EQKEY_ADV_STRATEGY, stra_key);
1369
1370 const char *form_key = _adv_formulation_key[nsp->adv_form];
1371 cs_equation_param_set(eqp, CS_EQKEY_ADV_FORMULATION, form_key);
1372
1373 const char *scheme_key = _adv_scheme_key[nsp->adv_scheme];
1374 cs_equation_param_set(eqp, CS_EQKEY_ADV_SCHEME, scheme_key);
1375
1376 }
1377
1378 cs_equation_param_set(eqp, CS_EQKEY_BC_QUADRATURE, quad_key);
1379 }
1380
1381 /*----------------------------------------------------------------------------*/
1382 /*!
1383 * \brief Summary of the main cs_navsto_param_t structure
1384 *
1385 * \param[in] nsp pointer to a cs_navsto_param_t structure
1386 */
1387 /*----------------------------------------------------------------------------*/
1388
1389 void
cs_navsto_param_log(const cs_navsto_param_t * nsp)1390 cs_navsto_param_log(const cs_navsto_param_t *nsp)
1391 {
1392 if (nsp == NULL)
1393 return;
1394
1395 const char navsto[16] = " * NavSto |";
1396
1397 /* Sanity checks */
1398
1399 if (nsp->model == CS_NAVSTO_N_MODELS)
1400 bft_error(__FILE__, __LINE__, 0, "%s: Invalid model for Navier-Stokes.\n",
1401 __func__);
1402 if (nsp->coupling == CS_NAVSTO_N_COUPLINGS)
1403 bft_error(__FILE__, __LINE__, 0,
1404 "%s: Invalid way of coupling the Navier-Stokes equations.\n",
1405 __func__);
1406
1407 cs_log_printf(CS_LOG_SETUP, "%s Verbosity: %d\n", navsto, nsp->verbosity);
1408
1409 /* Describe the physical modelling */
1410
1411 cs_log_printf(CS_LOG_SETUP, "%s Model: %s\n",
1412 navsto, cs_navsto_param_get_model_name(nsp->model));
1413
1414 if (nsp->model_flag & CS_NAVSTO_MODEL_GRAVITY_EFFECTS)
1415 cs_log_printf(CS_LOG_SETUP, "%s Model: Gravity effect activated\n",
1416 navsto);
1417
1418 if (nsp->model_flag & CS_NAVSTO_MODEL_CORIOLIS_EFFECTS)
1419 cs_log_printf(CS_LOG_SETUP, "%s Model: Coriolis effect activated\n",
1420 navsto);
1421
1422 if (nsp->model_flag & CS_NAVSTO_MODEL_BOUSSINESQ) {
1423
1424 cs_log_printf(CS_LOG_SETUP, "%s Model:"
1425 " Boussinesq approximation activated (%d term(s))\n",
1426 navsto, nsp->n_boussinesq_terms);
1427
1428 for (int i = 0; i < nsp->n_boussinesq_terms; i++) {
1429
1430 cs_navsto_param_boussinesq_t bp = nsp->boussinesq_param[i];
1431
1432 cs_log_printf(CS_LOG_SETUP,
1433 "%s Dilatation coef. = %5.3e; Reference value = %5.3e\n",
1434 navsto, bp.beta, bp.var0);
1435
1436 }
1437
1438 if (nsp->n_boussinesq_terms == 0)
1439 bft_error(__FILE__, __LINE__, 0,
1440 "%s: Boussinesq model is activated but no parameter is given\n"
1441 " to define the Boussinesq term.\n", __func__);
1442
1443 } /* Boussinesq term(s) added */
1444
1445 /* Describe the space-time discretization */
1446
1447 cs_log_printf(CS_LOG_SETUP, "%s Coupling: %s\n",
1448 navsto, cs_navsto_param_coupling_name[nsp->coupling]);
1449
1450 if (cs_navsto_param_is_steady(nsp))
1451 cs_log_printf(CS_LOG_SETUP, "%s Time status: Steady\n", navsto);
1452
1453 else {
1454
1455 cs_log_printf(CS_LOG_SETUP, "%s Time status: Unsteady\n", navsto);
1456
1457 const char *time_scheme = cs_param_get_time_scheme_name(nsp->time_scheme);
1458 if (time_scheme != NULL) {
1459 cs_log_printf(CS_LOG_SETUP, "%s Time scheme: %s", navsto, time_scheme);
1460 if (nsp->time_scheme == CS_TIME_SCHEME_THETA)
1461 cs_log_printf(CS_LOG_SETUP, " with value %f\n", nsp->theta);
1462 else
1463 cs_log_printf(CS_LOG_SETUP, "\n");
1464 }
1465 else
1466 bft_error(__FILE__, __LINE__, 0, "%s: Invalid time scheme.", __func__);
1467
1468 }
1469
1470 const char *space_scheme = cs_param_get_space_scheme_name(nsp->space_scheme);
1471 if (space_scheme != NULL)
1472 cs_log_printf(CS_LOG_SETUP, "%s Space scheme: %s\n", navsto, space_scheme);
1473 else
1474 bft_error(__FILE__, __LINE__, 0, " %s: Undefined space scheme.", __func__);
1475
1476 if (nsp->model == CS_NAVSTO_MODEL_INCOMPRESSIBLE_NAVIER_STOKES) {
1477
1478 /* Advection treament */
1479
1480 cs_log_printf(CS_LOG_SETUP, "%s Advection scheme: %s\n",
1481 navsto, cs_param_get_advection_scheme_name(nsp->adv_scheme));
1482 cs_log_printf(CS_LOG_SETUP, "%s Advection formulation: %s\n",
1483 navsto, cs_param_get_advection_form_name(nsp->adv_form));
1484 cs_log_printf(CS_LOG_SETUP, "%s Advection strategy: %s\n",
1485 navsto,
1486 cs_param_get_advection_strategy_name(nsp->adv_strategy));
1487 cs_log_printf(CS_LOG_SETUP, "%s Advection extrapolation: %s\n",
1488 navsto,
1489 cs_param_get_advection_extrapol_name(nsp->adv_extrapol));
1490
1491 /* Describe if needed the SLES settings for the non-linear algorithm */
1492
1493 if (nsp->adv_strategy == CS_PARAM_ADVECTION_IMPLICIT_FULL) {
1494
1495 const cs_navsto_param_sles_t *nslesp = nsp->sles_param;
1496
1497 cs_log_printf(CS_LOG_SETUP, "%s Non-linear algo: %s\n",
1498 navsto, cs_param_get_nl_algo_name(nslesp->nl_algo_type));
1499 cs_log_printf(CS_LOG_SETUP, "%s Tolerances of non-linear algo:"
1500 " rtol: %5.3e; atol: %5.3e; dtol: %5.3e\n",
1501 navsto, nslesp->nl_algo_param.rtol,
1502 nslesp->nl_algo_param.atol, nslesp->nl_algo_param.dtol);
1503 cs_log_printf(CS_LOG_SETUP, "%s Max of non-linear iterations: %d\n",
1504 navsto, nslesp->nl_algo_param.n_max_algo_iter);
1505
1506 if (nslesp->nl_algo_type == CS_PARAM_NL_ALGO_ANDERSON) {
1507
1508 const cs_iter_algo_param_aa_t aap = nslesp->anderson_param;
1509
1510 cs_log_printf(CS_LOG_SETUP, "%s Anderson param: max. dir: %d; "
1511 " start: %d; drop. tol: %5.3e; relax: %5.3e\n",
1512 navsto, aap.n_max_dir, aap.starting_iter, aap.max_cond,
1513 aap.beta);
1514 cs_log_printf(CS_LOG_SETUP, "%s Anderson param: Dot product type: %s\n",
1515 navsto, cs_param_get_dotprod_type_name(aap.dp_type));
1516
1517 }
1518
1519 } /* A non-linear treatment is requested */
1520
1521 } /* Navier-Stokes */
1522
1523 /* Describe the strategy to inverse the (inner) linear system */
1524
1525 _navsto_param_sles_log(nsp->sles_param);
1526
1527 if (nsp->gd_scale_coef > 0)
1528 cs_log_printf(CS_LOG_SETUP, "%s Grad-div scaling %e\n",
1529 navsto, nsp->gd_scale_coef);
1530
1531 /* Default quadrature type */
1532
1533 cs_log_printf(CS_LOG_SETUP, "%s Default quadrature: %s\n",
1534 navsto, cs_quadrature_get_type_name(nsp->qtype));
1535
1536 /* Initial conditions for the velocity */
1537
1538 cs_log_printf(CS_LOG_SETUP,
1539 "%s Velocity.Init.Cond | Number of definitions %2d\n",
1540 navsto, nsp->n_velocity_ic_defs);
1541
1542 char prefix[256];
1543 for (int i = 0; i < nsp->n_velocity_ic_defs; i++) {
1544 sprintf(prefix, "%s Velocity.Init.Cond | Definition %2d", navsto, i);
1545 cs_xdef_log(prefix, nsp->velocity_ic_defs[i]);
1546 }
1547
1548 /* Initial conditions for the pressure */
1549
1550 cs_log_printf(CS_LOG_SETUP,
1551 "%s Pressure.Init.Cond | Number of definitions: %d\n",
1552 navsto, nsp->n_pressure_ic_defs);
1553 for (int i = 0; i < nsp->n_pressure_ic_defs; i++) {
1554 sprintf(prefix, "%s Pressure.Init.Cond | Definition %2d", navsto, i);
1555 cs_xdef_log(prefix, nsp->pressure_ic_defs[i]);
1556 }
1557 }
1558
1559 /*----------------------------------------------------------------------------*/
1560 /*!
1561 * \brief Add a new Boussinesq term (source term for the momemtum equation)
1562 *
1563 * \param[in, out] nsp pointer to a cs_navsto_param_t structure
1564 *
1565 * \return a pointer to the newly added structure
1566 */
1567 /*----------------------------------------------------------------------------*/
1568
1569 cs_navsto_param_boussinesq_t *
cs_navsto_param_add_boussinesq_term(cs_navsto_param_t * nsp,cs_real_t dilatation_coef,cs_real_t reference_value)1570 cs_navsto_param_add_boussinesq_term(cs_navsto_param_t *nsp,
1571 cs_real_t dilatation_coef,
1572 cs_real_t reference_value)
1573 {
1574 if (nsp == NULL)
1575 return NULL;
1576
1577 nsp->model_flag |= CS_NAVSTO_MODEL_BOUSSINESQ;
1578
1579 int b_id = nsp->n_boussinesq_terms;
1580
1581 nsp->n_boussinesq_terms += 1;
1582 BFT_REALLOC(nsp->boussinesq_param, nsp->n_boussinesq_terms + 1,
1583 cs_navsto_param_boussinesq_t);
1584
1585 cs_navsto_param_boussinesq_t *bp = nsp->boussinesq_param + b_id;
1586
1587 bp->var0 = reference_value;
1588 bp->beta = dilatation_coef;
1589
1590 /* bp->var (shared pointer) is set with another function */
1591
1592 return bp;
1593 }
1594
1595 /*----------------------------------------------------------------------------*/
1596 /*!
1597 * \brief Set the array of values to consider in the Boussinesq term
1598 *
1599 * \param[in, out] bp pointer to a cs_navsto_param_boussinesq_t structure
1600 * \param[in] var shared pointer to the array of values to consider
1601 */
1602 /*----------------------------------------------------------------------------*/
1603
1604 void
cs_navsto_param_set_boussinesq_array(cs_navsto_param_boussinesq_t * bp,const cs_real_t * var)1605 cs_navsto_param_set_boussinesq_array(cs_navsto_param_boussinesq_t *bp,
1606 const cs_real_t *var)
1607 {
1608 if (bp == NULL)
1609 bft_error(__FILE__, __LINE__, 0,
1610 "%s: Boussinesq structure is empty\n", __func__);
1611
1612 bp->var = var;
1613 }
1614
1615 /*----------------------------------------------------------------------------*/
1616 /*!
1617 * \brief Retrieve the \ref cs_equation_param_t structure related to the
1618 * velocity equation (momentum equation in most of the cases)
1619 *
1620 * \param[in] nsp pointer to a cs_navsto_param_t structure
1621 *
1622 * \return a pointer to the set of SLES parameters
1623 */
1624 /*----------------------------------------------------------------------------*/
1625
1626 cs_navsto_param_sles_t *
cs_navsto_param_get_sles_param(const cs_navsto_param_t * nsp)1627 cs_navsto_param_get_sles_param(const cs_navsto_param_t *nsp)
1628 {
1629 if (nsp == NULL)
1630 return NULL;
1631
1632 return nsp->sles_param;
1633 }
1634
1635 /*----------------------------------------------------------------------------*/
1636 /*!
1637 * \brief Retrieve the \ref cs_equation_param_t structure related to the
1638 * velocity equation (momentum equation in most of the cases)
1639 *
1640 * \param[in] nsp pointer to a cs_navsto_param_t structure
1641 *
1642 * \return a pointer to the set of parameters related to the momentum equation
1643 */
1644 /*----------------------------------------------------------------------------*/
1645
1646 cs_equation_param_t *
cs_navsto_param_get_velocity_param(const cs_navsto_param_t * nsp)1647 cs_navsto_param_get_velocity_param(const cs_navsto_param_t *nsp)
1648 {
1649 cs_equation_param_t *eqp = NULL;
1650
1651 switch (nsp->coupling) {
1652
1653 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
1654 case CS_NAVSTO_COUPLING_MONOLITHIC:
1655 eqp = cs_equation_param_by_name("momentum");
1656 break;
1657
1658 case CS_NAVSTO_COUPLING_PROJECTION:
1659 eqp = cs_equation_param_by_name("velocity_prediction");
1660 break;
1661
1662 default:
1663 bft_error(__FILE__, __LINE__, 0,
1664 "%s: Invalid coupling algorithm", __func__);
1665 break;
1666
1667 }
1668
1669 return eqp;
1670 }
1671
1672 /*----------------------------------------------------------------------------*/
1673 /*!
1674 * \brief Retrieve the name of the model system of equations
1675 *
1676 * \param[in] model a \ref cs_navsto_param_model_t
1677 *
1678 * \return the corresponding name
1679 */
1680 /*----------------------------------------------------------------------------*/
1681
1682 const char *
cs_navsto_param_get_model_name(cs_navsto_param_model_t model)1683 cs_navsto_param_get_model_name(cs_navsto_param_model_t model)
1684 {
1685 switch (model) {
1686
1687 case CS_NAVSTO_MODEL_STOKES:
1688 case CS_NAVSTO_MODEL_OSEEN:
1689 case CS_NAVSTO_MODEL_INCOMPRESSIBLE_NAVIER_STOKES:
1690 return cs_navsto_param_model_name[model];
1691
1692 default:
1693 bft_error(__FILE__, __LINE__, 0, "%s: Invalid model.", __func__);
1694 break;
1695
1696 }
1697
1698 return NULL;
1699 }
1700
1701 /*----------------------------------------------------------------------------*/
1702 /*!
1703 * \brief Retrieve the name of the coupling algorithm
1704 *
1705 * \param[in] coupling a \ref cs_navsto_param_coupling_t
1706 *
1707 * \return the name
1708 */
1709 /*----------------------------------------------------------------------------*/
1710
1711 const char *
cs_navsto_param_get_coupling_name(cs_navsto_param_coupling_t coupling)1712 cs_navsto_param_get_coupling_name(cs_navsto_param_coupling_t coupling)
1713 {
1714 switch (coupling) {
1715
1716 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
1717 case CS_NAVSTO_COUPLING_MONOLITHIC:
1718 case CS_NAVSTO_COUPLING_PROJECTION:
1719 return cs_navsto_param_coupling_name[coupling];
1720
1721 default:
1722 bft_error(__FILE__, __LINE__, 0, "%s: Invalid coupling.", __func__);
1723 break;
1724
1725 }
1726
1727 return NULL;
1728 }
1729
1730 /*----------------------------------------------------------------------------*/
1731 /*!
1732 * \brief Set the value to consider for the reference pressure
1733 *
1734 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1735 * \param[in] pref value of the reference pressure
1736 */
1737 /*----------------------------------------------------------------------------*/
1738
1739 void
cs_navsto_set_reference_pressure(cs_navsto_param_t * nsp,cs_real_t pref)1740 cs_navsto_set_reference_pressure(cs_navsto_param_t *nsp,
1741 cs_real_t pref)
1742 {
1743 if (nsp == NULL)
1744 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
1745
1746 nsp->reference_pressure = pref;
1747 }
1748
1749 /*----------------------------------------------------------------------------*/
1750 /*!
1751 * \brief Define the initial condition for the velocity unknowns.
1752 * This definition can be done on a specified mesh location.
1753 * By default, the unknown is set to zero everywhere.
1754 * Here the initial value is set to a constant value
1755 *
1756 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1757 * \param[in] z_name name of the associated zone (if NULL or "" if
1758 * all cells are considered)
1759 * \param[in] val pointer to the value
1760 *
1761 * \return a pointer to the new \ref cs_xdef_t structure
1762 */
1763 /*----------------------------------------------------------------------------*/
1764
1765 cs_xdef_t *
cs_navsto_add_velocity_ic_by_value(cs_navsto_param_t * nsp,const char * z_name,cs_real_t * val)1766 cs_navsto_add_velocity_ic_by_value(cs_navsto_param_t *nsp,
1767 const char *z_name,
1768 cs_real_t *val)
1769 {
1770 if (nsp == NULL)
1771 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
1772
1773 cs_xdef_t *d = NULL;
1774 cs_equation_param_t *eqp = _get_momentum_param(nsp);
1775
1776 if (eqp != NULL) { /* An equation related to the velocity is defined */
1777
1778 d = cs_equation_add_ic_by_value(eqp, z_name, val);
1779
1780 }
1781 else { /* No momentum equation available with the choice of velocity-pressure
1782 coupling */
1783
1784 nsp->velocity_ic_is_owner = true;
1785
1786 /* Add a new cs_xdef_t structure */
1787
1788 int z_id = 0;
1789 if (z_name != NULL && strlen(z_name) > 0)
1790 z_id = (cs_volume_zone_by_name(z_name))->id;
1791
1792 cs_flag_t meta_flag = 0;
1793 if (z_id == 0)
1794 meta_flag |= CS_FLAG_FULL_LOC;
1795
1796 d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
1797 3, /* dim */
1798 z_id,
1799 CS_FLAG_STATE_UNIFORM,
1800 meta_flag,
1801 val);
1802 }
1803
1804 int new_id = nsp->n_velocity_ic_defs;
1805 nsp->n_velocity_ic_defs += 1;
1806 BFT_REALLOC(nsp->velocity_ic_defs, nsp->n_velocity_ic_defs, cs_xdef_t *);
1807 nsp->velocity_ic_defs[new_id] = d;
1808
1809 return d;
1810 }
1811
1812 /*----------------------------------------------------------------------------*/
1813 /*!
1814 * \brief Define the initial condition for the velocity unknowns.
1815 * This definition can be done on a specified mesh location.
1816 * By default, the unknown is set to zero everywhere.
1817 * Here the initial value is set according to an analytical function
1818 *
1819 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1820 * \param[in] z_name name of the associated zone (if NULL or "" if
1821 * all cells are considered)
1822 * \param[in] analytic pointer to an analytic function
1823 * \param[in] input NULL or pointer to a structure cast on-the-fly
1824 *
1825 * \return a pointer to the new \ref cs_xdef_t structure
1826 */
1827 /*----------------------------------------------------------------------------*/
1828
1829 cs_xdef_t *
cs_navsto_add_velocity_ic_by_analytic(cs_navsto_param_t * nsp,const char * z_name,cs_analytic_func_t * analytic,void * input)1830 cs_navsto_add_velocity_ic_by_analytic(cs_navsto_param_t *nsp,
1831 const char *z_name,
1832 cs_analytic_func_t *analytic,
1833 void *input)
1834 {
1835 if (nsp == NULL)
1836 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
1837
1838 cs_xdef_t *d = NULL;
1839 cs_equation_param_t *eqp = _get_momentum_param(nsp);
1840
1841 if (eqp != NULL) { /* An equation related to the velocity is defined */
1842
1843 d = cs_equation_add_ic_by_analytic(eqp, z_name, analytic, input);
1844
1845 }
1846 else { /* No momentum equation available with the choice of velocity-pressure
1847 coupling */
1848
1849 nsp->velocity_ic_is_owner = true;
1850
1851 /* Add a new cs_xdef_t structure */
1852
1853 int z_id = cs_get_vol_zone_id(z_name);
1854
1855 cs_flag_t meta_flag = 0;
1856 if (z_id == 0)
1857 meta_flag |= CS_FLAG_FULL_LOC;
1858
1859 cs_xdef_analytic_context_t anai = { .z_id = z_id,
1860 .func = analytic,
1861 .input = input,
1862 .free_input = NULL };
1863
1864 d = cs_xdef_volume_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
1865 3, /* dim */
1866 z_id,
1867 0, /* state flag */
1868 meta_flag,
1869 &anai);
1870 }
1871
1872 /* Assign the default quadrature type of the Navier-Stokes module to this
1873 * definition (this can be modified by the user if the same call is performed
1874 * in cs_user_finalize_setup()) */
1875
1876 cs_xdef_set_quadrature(d, nsp->qtype);
1877
1878 int new_id = nsp->n_velocity_ic_defs;
1879 nsp->n_velocity_ic_defs += 1;
1880 BFT_REALLOC(nsp->velocity_ic_defs, nsp->n_velocity_ic_defs, cs_xdef_t *);
1881 nsp->velocity_ic_defs[new_id] = d;
1882
1883 return d;
1884 }
1885
1886 /*----------------------------------------------------------------------------*/
1887 /*!
1888 * \brief Define the initial condition for the pressure unknowns.
1889 * This definition can be done on a specified mesh location.
1890 * By default, the unknown is set to zero everywhere.
1891 * Here the initial value is set to a constant value
1892 *
1893 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1894 * \param[in] z_name name of the associated zone (if NULL or "" if
1895 * all cells are considered)
1896 * \param[in] val pointer to the value
1897 *
1898 * \return a pointer to the new \ref cs_xdef_t structure
1899 */
1900 /*----------------------------------------------------------------------------*/
1901
1902 cs_xdef_t *
cs_navsto_add_pressure_ic_by_value(cs_navsto_param_t * nsp,const char * z_name,cs_real_t * val)1903 cs_navsto_add_pressure_ic_by_value(cs_navsto_param_t *nsp,
1904 const char *z_name,
1905 cs_real_t *val)
1906 {
1907 if (nsp == NULL)
1908 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
1909
1910 /* Add a new cs_xdef_t structure */
1911
1912 int z_id = 0;
1913 if (z_name != NULL && strlen(z_name) > 0)
1914 z_id = (cs_volume_zone_by_name(z_name))->id;
1915
1916 cs_flag_t meta_flag = 0;
1917 if (z_id == 0)
1918 meta_flag |= CS_FLAG_FULL_LOC;
1919
1920 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
1921 1, /* dim */
1922 z_id,
1923 CS_FLAG_STATE_UNIFORM,
1924 meta_flag,
1925 val);
1926
1927 int new_id = nsp->n_pressure_ic_defs;
1928 nsp->n_pressure_ic_defs += 1;
1929 BFT_REALLOC(nsp->pressure_ic_defs, nsp->n_pressure_ic_defs, cs_xdef_t *);
1930 nsp->pressure_ic_defs[new_id] = d;
1931
1932 return d;
1933 }
1934
1935 /*----------------------------------------------------------------------------*/
1936 /*!
1937 * \brief Define the initial condition for the pressure unknowns.
1938 * This definition can be done on a specified mesh location.
1939 * By default, the unknown is set to zero everywhere.
1940 * Here the initial value is set according to an analytical function
1941 *
1942 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1943 * \param[in] z_name name of the associated zone (if NULL or "" if
1944 * all cells are considered)
1945 * \param[in] analytic pointer to an analytic function
1946 * \param[in] input NULL or pointer to a structure cast on-the-fly
1947 *
1948 * \return a pointer to the new \ref cs_xdef_t structure
1949 */
1950 /*----------------------------------------------------------------------------*/
1951
1952 cs_xdef_t *
cs_navsto_add_pressure_ic_by_analytic(cs_navsto_param_t * nsp,const char * z_name,cs_analytic_func_t * analytic,void * input)1953 cs_navsto_add_pressure_ic_by_analytic(cs_navsto_param_t *nsp,
1954 const char *z_name,
1955 cs_analytic_func_t *analytic,
1956 void *input)
1957 {
1958 if (nsp == NULL)
1959 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
1960
1961 /* Add a new cs_xdef_t structure */
1962
1963 int z_id = cs_get_vol_zone_id(z_name);
1964
1965 cs_flag_t meta_flag = 0;
1966 if (z_id == 0)
1967 meta_flag |= CS_FLAG_FULL_LOC;
1968
1969 cs_xdef_analytic_context_t ac = { .z_id = z_id,
1970 .func = analytic,
1971 .input = input,
1972 .free_input = NULL };
1973
1974 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
1975 1, /* dim */
1976 z_id,
1977 0, /* state flag */
1978 meta_flag,
1979 &ac);
1980
1981 /* Assign the default quadrature type of the Navier-Stokes module to this
1982 * definition (this can be modified by the user if the same call is
1983 * performed in cs_user_finalize_setup()) */
1984
1985 cs_xdef_set_quadrature(d, nsp->qtype);
1986
1987 int new_id = nsp->n_pressure_ic_defs;
1988 nsp->n_pressure_ic_defs += 1;
1989 BFT_REALLOC(nsp->pressure_ic_defs, nsp->n_pressure_ic_defs, cs_xdef_t *);
1990 nsp->pressure_ic_defs[new_id] = d;
1991
1992 return d;
1993 }
1994
1995 /*----------------------------------------------------------------------------*/
1996 /*!
1997 * \brief Add the definition of boundary conditions related to a fixed wall
1998 * into the set of parameters for the management of the Navier-Stokes
1999 * system of equations
2000 *
2001 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2002 */
2003 /*----------------------------------------------------------------------------*/
2004
2005 void
cs_navsto_set_fixed_walls(cs_navsto_param_t * nsp)2006 cs_navsto_set_fixed_walls(cs_navsto_param_t *nsp)
2007 {
2008 if (nsp == NULL)
2009 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2010 assert(nsp->boundaries != NULL);
2011
2012 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2013 cs_real_3_t zero = {0, 0, 0};
2014
2015 const cs_boundary_t *bdy = nsp->boundaries;
2016
2017 for (int i = 0; i < bdy->n_boundaries; i++) {
2018 if ( bdy->types[i] & CS_BOUNDARY_WALL
2019 && !(bdy->types[i] & CS_BOUNDARY_SLIDING_WALL)) {
2020
2021 /* Homogeneous Dirichlet on the velocity field. Nothing to enforce on the
2022 pressure field */
2023
2024 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2025 3, /* dim */
2026 bdy->zone_ids[i],
2027 CS_FLAG_STATE_UNIFORM, /* state */
2028 CS_CDO_BC_HMG_DIRICHLET,
2029 (void *)zero);
2030 int new_id = nsp->n_velocity_bc_defs;
2031
2032 nsp->n_velocity_bc_defs += 1;
2033 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2034 nsp->velocity_bc_defs[new_id] = d;
2035
2036 cs_equation_add_xdef_bc(eqp, d);
2037
2038 /* Homogeneous Neumann on the pressure field --> default BC */
2039
2040 } /* Fixed wall */
2041 } /* Loop on domain boundaries */
2042 }
2043
2044 /*----------------------------------------------------------------------------*/
2045 /*!
2046 * \brief Add the definition of boundary conditions related to a symmetry
2047 * into the set of parameters for the management of the Navier-Stokes
2048 * system of equations
2049 *
2050 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2051 */
2052 /*----------------------------------------------------------------------------*/
2053
2054 void
cs_navsto_set_symmetries(cs_navsto_param_t * nsp)2055 cs_navsto_set_symmetries(cs_navsto_param_t *nsp)
2056 {
2057 if (nsp == NULL)
2058 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2059 assert(nsp->boundaries != NULL);
2060
2061 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2062 cs_real_t zero = 0;
2063
2064 const cs_boundary_t *bdy = nsp->boundaries;
2065
2066 for (int i = 0; i < bdy->n_boundaries; i++) {
2067 if (bdy->types[i] & CS_BOUNDARY_SYMMETRY) {
2068
2069 /* Homogeneous Dirichlet on the normal component of the velocity field
2070 and homogeneous Neumann on the normal stress (balance between the
2071 pressure gradient and the viscous term) */
2072
2073 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2074 1, /* dim */
2075 bdy->zone_ids[i],
2076 CS_FLAG_STATE_UNIFORM, /* state */
2077 CS_CDO_BC_SLIDING,
2078 (void *)&zero);
2079
2080 cs_equation_add_xdef_bc(eqp, d);
2081
2082 int new_id = nsp->n_velocity_bc_defs;
2083
2084 nsp->n_velocity_bc_defs += 1;
2085 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2086 nsp->velocity_bc_defs[new_id] = d;
2087
2088 /* Homogeneous Neumann on the pressure field --> default BC (Nothing to
2089 do) */
2090
2091 } /* Symmetry */
2092 } /* Loop on domain boundaries */
2093 }
2094
2095 /*----------------------------------------------------------------------------*/
2096 /*!
2097 * \brief Add the definition of boundary conditions related to outlets
2098 * into the set of parameters for the management of the Navier-Stokes
2099 * system of equations
2100 *
2101 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2102 */
2103 /*----------------------------------------------------------------------------*/
2104
2105 void
cs_navsto_set_outlets(cs_navsto_param_t * nsp)2106 cs_navsto_set_outlets(cs_navsto_param_t *nsp)
2107 {
2108 if (nsp == NULL)
2109 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2110 assert(nsp->boundaries != NULL);
2111
2112 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2113 cs_real_33_t zero = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
2114
2115 const cs_boundary_t *bdy = nsp->boundaries;
2116
2117 int exclude_filter = CS_BOUNDARY_IMPOSED_P | CS_BOUNDARY_IMPOSED_VEL;
2118
2119 for (int i = 0; i < bdy->n_boundaries; i++) {
2120 if ( bdy->types[i] & CS_BOUNDARY_OUTLET
2121 && ! (bdy->types[i] & exclude_filter)) {
2122
2123 /* Add the homogeneous Neumann on the normal component */
2124
2125 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2126 9, /* dim */
2127 bdy->zone_ids[i],
2128 CS_FLAG_STATE_UNIFORM, /* state */
2129 CS_CDO_BC_HMG_NEUMANN,
2130 (void *)&zero);
2131 cs_equation_add_xdef_bc(eqp, d);
2132
2133 int new_id = nsp->n_velocity_bc_defs;
2134
2135 nsp->n_velocity_bc_defs += 1;
2136 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2137 nsp->velocity_bc_defs[new_id] = d;
2138
2139 /* Homogeneous Neumann on the pressure field --> default BC.
2140 At the end of the day, we end up with a balance between the pressure
2141 gradient and the viscous term (and advection term in Navier-Stokes) */
2142
2143 } /* Symmetry */
2144 } /* Loop on domain boundaries */
2145 }
2146
2147 /*----------------------------------------------------------------------------*/
2148 /*!
2149 * \brief Set the pressure field on a boundary using a uniform value.
2150 *
2151 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2152 * \param[in] z_name name of the associated zone (if NULL or "" all
2153 * boundary faces are considered)
2154 * \param[in] value value to set
2155 *
2156 * \return a pointer to the new \ref cs_xdef_t structure
2157 */
2158 /*----------------------------------------------------------------------------*/
2159
2160 cs_xdef_t *
cs_navsto_set_pressure_bc_by_value(cs_navsto_param_t * nsp,const char * z_name,cs_real_t * values)2161 cs_navsto_set_pressure_bc_by_value(cs_navsto_param_t *nsp,
2162 const char *z_name,
2163 cs_real_t *values)
2164 {
2165 if (nsp == NULL)
2166 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2167
2168 int z_id = cs_get_bdy_zone_id(z_name);
2169 if (z_id < 0)
2170 bft_error(__FILE__, __LINE__, 0,
2171 " %s: Zone \"%s\" does not exist.\n"
2172 " Please check your settings.", __func__, z_name);
2173
2174 int bdy_id = cs_boundary_id_by_zone_id(nsp->boundaries, z_id);
2175 if (bdy_id < 0)
2176 bft_error(__FILE__, __LINE__, 0,
2177 " %s: Zone \"%s\" does not belong to an existing boundary.\n"
2178 " Please check your settings.", __func__, z_name);
2179
2180 if (!(nsp->boundaries->types[bdy_id] & CS_BOUNDARY_IMPOSED_P))
2181 bft_error(__FILE__, __LINE__, 0,
2182 " %s: Zone \"%s\" is not related to a pressure boundary.\n"
2183 " Please check your settings.", __func__, z_name);
2184
2185 /* Set the boundary condition for the pressure field */
2186
2187 cs_xdef_t *dp = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2188 1, /* dim */
2189 z_id,
2190 CS_FLAG_STATE_UNIFORM, /* state */
2191 CS_CDO_BC_DIRICHLET,
2192 (void *)values);
2193
2194 int pnew_id = nsp->n_pressure_bc_defs;
2195
2196 nsp->n_pressure_bc_defs += 1;
2197 BFT_REALLOC(nsp->pressure_bc_defs, nsp->n_pressure_bc_defs, cs_xdef_t *);
2198 nsp->pressure_bc_defs[pnew_id] = dp;
2199
2200 if (!nsp->pressure_bc_is_owner) {
2201 bft_error(__FILE__, __LINE__, 0, "%s: Not implemented yet", __func__);
2202
2203 #if 0 /* TODO */
2204 /* Retrieve the equation related to the pressure */
2205 cs_equation_param_t *p_eqp = NULL;
2206 assert(p_eqp != NULL);
2207 #endif
2208 }
2209
2210 /* Add a new cs_xdef_t structure. For the momentum equation, this is a
2211 * homogeneous Neumann BC for the velocity */
2212
2213 cs_real_33_t zero = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
2214
2215 cs_xdef_t *du = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2216 9, /* dim */
2217 z_id,
2218 CS_FLAG_STATE_UNIFORM, /* state */
2219 CS_CDO_BC_HMG_NEUMANN,
2220 (void *)zero);
2221
2222 int unew_id = nsp->n_velocity_bc_defs;
2223
2224 nsp->n_velocity_bc_defs += 1;
2225 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2226 nsp->velocity_bc_defs[unew_id] = du;
2227
2228 cs_equation_param_t *u_eqp = _get_momentum_param(nsp);
2229 assert(u_eqp != NULL);
2230 cs_equation_add_xdef_bc(u_eqp, du);
2231
2232 return dp;
2233 }
2234
2235 /*----------------------------------------------------------------------------*/
2236 /*!
2237 * \brief Define the velocity field for a sliding wall boundary using a
2238 * uniform value
2239 *
2240 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2241 * \param[in] z_name name of the associated zone (if NULL or "" all
2242 * boundary faces are considered)
2243 * \param[in] values array of three real values
2244 *
2245 * \return a pointer to the new \ref cs_xdef_t structure
2246 */
2247 /*----------------------------------------------------------------------------*/
2248
2249 cs_xdef_t *
cs_navsto_set_velocity_wall_by_value(cs_navsto_param_t * nsp,const char * z_name,cs_real_t * values)2250 cs_navsto_set_velocity_wall_by_value(cs_navsto_param_t *nsp,
2251 const char *z_name,
2252 cs_real_t *values)
2253 {
2254 if (nsp == NULL)
2255 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2256
2257 int z_id = cs_get_bdy_zone_id(z_name);
2258 if (z_id < 0)
2259 bft_error(__FILE__, __LINE__, 0,
2260 " %s: Zone \"%s\" does not exist.\n"
2261 " Please check your settings.", __func__, z_name);
2262
2263 int bdy_id = cs_boundary_id_by_zone_id(nsp->boundaries, z_id);
2264 if (bdy_id < 0)
2265 bft_error(__FILE__, __LINE__, 0,
2266 " %s: Zone \"%s\" does not belong to an existing boundary.\n"
2267 " Please check your settings.", __func__, z_name);
2268
2269 if (! (nsp->boundaries->types[bdy_id] & CS_BOUNDARY_SLIDING_WALL))
2270 bft_error(__FILE__, __LINE__, 0,
2271 " %s: Zone \"%s\" is not related to a sliding wall boundary.\n"
2272 " Please check your settings.", __func__, z_name);
2273
2274 /* Add a new cs_xdef_t structure */
2275 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2276 3, /* dim */
2277 z_id,
2278 CS_FLAG_STATE_UNIFORM, /* state */
2279 CS_CDO_BC_DIRICHLET,
2280 (void *)values);
2281
2282 int new_id = nsp->n_velocity_bc_defs;
2283
2284 nsp->n_velocity_bc_defs += 1;
2285 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2286 nsp->velocity_bc_defs[new_id] = d;
2287
2288 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2289 assert(eqp != NULL);
2290 cs_equation_add_xdef_bc(eqp, d);
2291
2292 return d;
2293 }
2294
2295 /*----------------------------------------------------------------------------*/
2296 /*!
2297 * \brief Define the velocity field for an inlet boundary using a uniform
2298 * value
2299 *
2300 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2301 * \param[in] z_name name of the associated zone (if NULL or "" all
2302 * boundary faces are considered)
2303 * \param[in] values array of three real values
2304 *
2305 * \return a pointer to the new \ref cs_xdef_t structure
2306 */
2307 /*----------------------------------------------------------------------------*/
2308
2309 cs_xdef_t *
cs_navsto_set_velocity_inlet_by_value(cs_navsto_param_t * nsp,const char * z_name,cs_real_t * values)2310 cs_navsto_set_velocity_inlet_by_value(cs_navsto_param_t *nsp,
2311 const char *z_name,
2312 cs_real_t *values)
2313 {
2314 if (nsp == NULL)
2315 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2316
2317 int z_id = cs_get_bdy_zone_id(z_name);
2318 if (z_id < 0)
2319 bft_error(__FILE__, __LINE__, 0,
2320 " %s: Zone \"%s\" does not exist.\n"
2321 " Please check your settings.", __func__, z_name);
2322
2323 int bdy_id = cs_boundary_id_by_zone_id(nsp->boundaries, z_id);
2324 if (bdy_id < 0)
2325 bft_error(__FILE__, __LINE__, 0,
2326 " %s: Zone \"%s\" does not belong to an existing boundary.\n"
2327 " Please check your settings.", __func__, z_name);
2328
2329 if (!(nsp->boundaries->types[bdy_id] & CS_BOUNDARY_IMPOSED_VEL))
2330 bft_error
2331 (__FILE__, __LINE__, 0,
2332 " %s: Zone \"%s\" is not related to an imposed velocity boundary.\n"
2333 " Please check your settings.", __func__, z_name);
2334
2335 /* Add a new cs_xdef_t structure */
2336
2337 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
2338 3, /* dim */
2339 z_id,
2340 CS_FLAG_STATE_UNIFORM, /* state */
2341 CS_CDO_BC_DIRICHLET,
2342 (void *)values);
2343
2344 int new_id = nsp->n_velocity_bc_defs;
2345
2346 nsp->n_velocity_bc_defs += 1;
2347 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2348 nsp->velocity_bc_defs[new_id] = d;
2349
2350 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2351 assert(eqp != NULL);
2352 cs_equation_add_xdef_bc(eqp, d);
2353
2354 return d;
2355 }
2356
2357 /*----------------------------------------------------------------------------*/
2358 /*!
2359 * \brief Define the velocity field for an inlet boundary using an analytical
2360 * function
2361 *
2362 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2363 * \param[in] z_name name of the associated zone (if NULL or "" all
2364 * boundary faces are considered)
2365 * \param[in] ana pointer to an analytical 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_navsto_set_velocity_inlet_by_analytic(cs_navsto_param_t * nsp,const char * z_name,cs_analytic_func_t * ana,void * input)2373 cs_navsto_set_velocity_inlet_by_analytic(cs_navsto_param_t *nsp,
2374 const char *z_name,
2375 cs_analytic_func_t *ana,
2376 void *input)
2377 {
2378 if (nsp == NULL)
2379 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2380
2381 int z_id = cs_get_bdy_zone_id(z_name);
2382 if (z_id < 0)
2383 bft_error(__FILE__, __LINE__, 0,
2384 " %s: Zone \"%s\" does not exist.\n"
2385 " Please check your settings.", __func__, z_name);
2386
2387 int bdy_id = cs_boundary_id_by_zone_id(nsp->boundaries, z_id);
2388 if (bdy_id < 0)
2389 bft_error(__FILE__, __LINE__, 0,
2390 " %s: Zone \"%s\" does not belong to an existing boundary.\n"
2391 " Please check your settings.", __func__, z_name);
2392
2393 if (!(nsp->boundaries->types[bdy_id] & CS_BOUNDARY_IMPOSED_VEL))
2394 bft_error
2395 (__FILE__, __LINE__, 0,
2396 " %s: Zone \"%s\" is not related to an imposed velocity boundary.\n"
2397 " Please check your settings.", __func__, z_name);
2398
2399 /* Add a new cs_xdef_t structure */
2400
2401 cs_xdef_analytic_context_t ac = { .z_id = z_id,
2402 .func = ana,
2403 .input = input,
2404 .free_input = NULL };
2405
2406 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
2407 3, /* dim */
2408 z_id,
2409 0, /* state */
2410 CS_CDO_BC_DIRICHLET,
2411 &ac);
2412
2413 /* Assign the default quadrature type of the Navier-Stokes module to this
2414 * definition (this can be modified by the user if the same call is
2415 * performed in cs_user_finalize_setup()) */
2416
2417 cs_xdef_set_quadrature(d, nsp->qtype);
2418
2419 int new_id = nsp->n_velocity_bc_defs;
2420 nsp->n_velocity_bc_defs += 1;
2421 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2422 nsp->velocity_bc_defs[new_id] = d;
2423
2424 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2425 cs_equation_add_xdef_bc(eqp, d);
2426
2427 return d;
2428 }
2429
2430 /*----------------------------------------------------------------------------*/
2431 /*!
2432 * \brief Define the velocity field for an inlet boundary using an array
2433 * of values
2434 *
2435 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2436 * \param[in] z_name name of the associated zone (if NULL or "" all
2437 * boundary faces are considered)
2438 * \param[in] loc information to know where are located values
2439 * \param[in] array pointer to an array
2440 * \param[in] is_owner transfer the lifecycle to the cs_xdef_t structure
2441 * (true or false)
2442 * \param[in] index optional pointer to the array index
2443 *
2444 * \return a pointer to the new \ref cs_xdef_t structure
2445 */
2446 /*----------------------------------------------------------------------------*/
2447
2448 cs_xdef_t *
cs_navsto_set_velocity_inlet_by_array(cs_navsto_param_t * nsp,const char * z_name,cs_flag_t loc,cs_real_t * array,bool is_owner,cs_lnum_t * index)2449 cs_navsto_set_velocity_inlet_by_array(cs_navsto_param_t *nsp,
2450 const char *z_name,
2451 cs_flag_t loc,
2452 cs_real_t *array,
2453 bool is_owner,
2454 cs_lnum_t *index)
2455 {
2456 if (nsp == NULL)
2457 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2458
2459 int z_id = cs_get_bdy_zone_id(z_name);
2460 if (z_id < 0)
2461 bft_error(__FILE__, __LINE__, 0,
2462 " %s: Zone \"%s\" does not exist.\n"
2463 " Please check your settings.", __func__, z_name);
2464
2465 int bdy_id = cs_boundary_id_by_zone_id(nsp->boundaries, z_id);
2466 if (bdy_id < 0)
2467 bft_error(__FILE__, __LINE__, 0,
2468 " %s: Zone \"%s\" does not belong to an existing boundary.\n"
2469 " Please check your settings.", __func__, z_name);
2470
2471 if (!(nsp->boundaries->types[bdy_id] & CS_BOUNDARY_IMPOSED_VEL))
2472 bft_error
2473 (__FILE__, __LINE__, 0,
2474 " %s: Zone \"%s\" is not related to an imposed velocity boundary.\n"
2475 " Please check your settings.", __func__, z_name);
2476
2477 cs_xdef_array_context_t context = {.z_id = z_id,
2478 .stride = 3,
2479 .loc = loc,
2480 .values = array,
2481 .is_owner = is_owner,
2482 .index = index };
2483
2484 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_ARRAY,
2485 3,
2486 z_id,
2487 CS_FLAG_STATE_FACEWISE,
2488 CS_CDO_BC_DIRICHLET,
2489 (void *)&context);
2490
2491 int new_id = nsp->n_velocity_bc_defs;
2492
2493 nsp->n_velocity_bc_defs += 1;
2494 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2495 nsp->velocity_bc_defs[new_id] = d;
2496
2497 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2498 assert(eqp != NULL);
2499 cs_equation_add_xdef_bc(eqp, d);
2500
2501 return d;
2502 }
2503
2504 /*----------------------------------------------------------------------------*/
2505 /*!
2506 * \brief Define the velocity field for an inlet boundary using a DoF function
2507 *
2508 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2509 * \param[in] z_name name of the associated zone (if NULL or "" all
2510 * boundary faces are considered)
2511 * \param[in] func pointer to a \ref cs_dof_function_t
2512 * \param[in] func_input NULL or pointer to a structure cast on-the-fly
2513 *
2514 * \return a pointer to the new \ref cs_xdef_t structure
2515 */
2516 /*----------------------------------------------------------------------------*/
2517
2518 cs_xdef_t *
cs_navsto_set_velocity_inlet_by_dof_func(cs_navsto_param_t * nsp,const char * z_name,cs_dof_func_t * func,void * func_input)2519 cs_navsto_set_velocity_inlet_by_dof_func(cs_navsto_param_t *nsp,
2520 const char *z_name,
2521 cs_dof_func_t *func,
2522 void *func_input)
2523 {
2524 if (nsp == NULL)
2525 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2526
2527 int z_id = cs_get_bdy_zone_id(z_name);
2528 if (z_id < 0)
2529 bft_error(__FILE__, __LINE__, 0,
2530 " %s: Zone \"%s\" does not exist.\n"
2531 " Please check your settings.", __func__, z_name);
2532
2533 int bdy_id = cs_boundary_id_by_zone_id(nsp->boundaries, z_id);
2534 if (bdy_id < 0)
2535 bft_error(__FILE__, __LINE__, 0,
2536 " %s: Zone \"%s\" does not belong to an existing boundary.\n"
2537 " Please check your settings.", __func__, z_name);
2538
2539 if (!(nsp->boundaries->types[bdy_id] & CS_BOUNDARY_IMPOSED_VEL))
2540 bft_error
2541 (__FILE__, __LINE__, 0,
2542 " %s: Zone \"%s\" is not related to an imposed velocity boundary.\n"
2543 " Please check your settings.", __func__, z_name);
2544
2545 /* Add a new cs_xdef_t structure */
2546 cs_xdef_dof_context_t dc = { .z_id = z_id,
2547 .func = func,
2548 .input = func_input,
2549 .free_input = NULL };
2550
2551 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_DOF_FUNCTION,
2552 3, /* dim */
2553 z_id,
2554 0, /* state */
2555 CS_CDO_BC_DIRICHLET,
2556 &dc);
2557
2558 /* Assign the default quadrature type of the Navier-Stokes module to this
2559 * definition (this can be modified by the user if the same call is
2560 * performed in cs_user_finalize_setup()) */
2561
2562 cs_xdef_set_quadrature(d, nsp->qtype);
2563
2564 int new_id = nsp->n_velocity_bc_defs;
2565 nsp->n_velocity_bc_defs += 1;
2566 BFT_REALLOC(nsp->velocity_bc_defs, nsp->n_velocity_bc_defs, cs_xdef_t *);
2567 nsp->velocity_bc_defs[new_id] = d;
2568
2569 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2570 cs_equation_add_xdef_bc(eqp, d);
2571
2572 return d;
2573 }
2574
2575 /*----------------------------------------------------------------------------*/
2576 /*!
2577 * \brief Define a new source term structure defined by an analytical function
2578 *
2579 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2580 * \param[in] z_name name of the associated zone (if NULL or "" all
2581 * cells are considered)
2582 * \param[in] ana pointer to an analytical function
2583 * \param[in] input NULL or pointer to a structure cast on-the-fly
2584 *
2585 * \return a pointer to the new \ref cs_xdef_t structure
2586 */
2587 /*----------------------------------------------------------------------------*/
2588
2589 cs_xdef_t *
cs_navsto_add_source_term_by_analytic(cs_navsto_param_t * nsp,const char * z_name,cs_analytic_func_t * ana,void * input)2590 cs_navsto_add_source_term_by_analytic(cs_navsto_param_t *nsp,
2591 const char *z_name,
2592 cs_analytic_func_t *ana,
2593 void *input)
2594 {
2595 if (nsp == NULL)
2596 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2597
2598 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2599 cs_xdef_t *d = cs_equation_add_source_term_by_analytic(eqp,
2600 z_name, ana, input);
2601
2602 /* Assign the default quadrature type of the Navier-Stokes module to this
2603 * definition (this can be modified by the user if the same call is
2604 * performed in cs_user_finalize_setup()) */
2605
2606 cs_xdef_set_quadrature(d, nsp->qtype);
2607
2608 return d;
2609 }
2610
2611 /*----------------------------------------------------------------------------*/
2612 /*!
2613 * \brief Define a new source term structure defined by a constant value
2614 *
2615 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2616 * \param[in] z_name name of the associated zone (if NULL or "" all
2617 * cells are considered)
2618 * \param[in] val pointer to the value to set
2619 *
2620 * \return a pointer to the new \ref cs_xdef_t structure
2621 */
2622 /*----------------------------------------------------------------------------*/
2623
2624 cs_xdef_t *
cs_navsto_add_source_term_by_val(cs_navsto_param_t * nsp,const char * z_name,cs_real_t * val)2625 cs_navsto_add_source_term_by_val(cs_navsto_param_t *nsp,
2626 const char *z_name,
2627 cs_real_t *val)
2628 {
2629 if (nsp == NULL)
2630 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2631
2632 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2633
2634 return cs_equation_add_source_term_by_val(eqp, z_name, val);
2635 }
2636
2637 /*----------------------------------------------------------------------------*/
2638 /*!
2639 * \brief Define a new source term structure defined by an array
2640 *
2641 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
2642 * \param[in] z_name name of the associated zone (if NULL or "" all
2643 * cells are considered)
2644 * \param[in] loc information to know where are located values
2645 * \param[in] array pointer to an array
2646 * \param[in] is_owner transfer the lifecycle to the cs_xdef_t structure
2647 * (true or false)
2648 * \param[in] index optional pointer to the array index
2649 *
2650 * \return a pointer to the new \ref cs_xdef_t structure
2651 */
2652 /*----------------------------------------------------------------------------*/
2653
2654 cs_xdef_t *
cs_navsto_add_source_term_by_array(cs_navsto_param_t * nsp,const char * z_name,cs_flag_t loc,cs_real_t * array,bool is_owner,cs_lnum_t * index)2655 cs_navsto_add_source_term_by_array(cs_navsto_param_t *nsp,
2656 const char *z_name,
2657 cs_flag_t loc,
2658 cs_real_t *array,
2659 bool is_owner,
2660 cs_lnum_t *index)
2661 {
2662 if (nsp == NULL)
2663 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2664
2665 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2666
2667 return cs_equation_add_source_term_by_array(eqp, z_name, loc,
2668 array, is_owner, index);
2669 }
2670
2671 /*----------------------------------------------------------------------------*/
2672 /*!
2673 * \brief Add a advection field for the Oseen problem
2674 *
2675 * \param[in, out] nsp pointer to a \ref cs_navsto_param_t
2676 * \param[in, out] adv_fld pointer to a \ref cs_adv_field_t
2677 */
2678 /*----------------------------------------------------------------------------*/
2679
2680 void
cs_navsto_add_oseen_field(cs_navsto_param_t * nsp,cs_adv_field_t * adv_fld)2681 cs_navsto_add_oseen_field(cs_navsto_param_t *nsp,
2682 cs_adv_field_t *adv_fld)
2683 {
2684 if (nsp == NULL)
2685 bft_error(__FILE__, __LINE__, 0, _err_empty_nsp, __func__);
2686
2687 if (nsp->model != CS_NAVSTO_MODEL_OSEEN)
2688 bft_error(__FILE__, __LINE__, 0, " %s: Trying to set an external advection"
2689 " where there should not be one. Stopping",
2690 __func__);
2691
2692 cs_equation_param_t *eqp = _get_momentum_param(nsp);
2693
2694 cs_equation_add_advection(eqp, adv_fld);
2695 }
2696
2697 /*----------------------------------------------------------------------------*/
2698
2699 END_C_DECLS
2700