1 /*============================================================================
2 * User functions for input of calculation parameters.
3 *============================================================================*/
4
5 /* VERS */
6
7 /*
8 This file is part of Code_Saturne, a general-purpose CFD tool.
9
10 Copyright (C) 1998-2021 EDF S.A.
11
12 This program is free software; you can redistribute it and/or modify it under
13 the terms of the GNU General Public License as published by the Free Software
14 Foundation; either version 2 of the License, or (at your option) any later
15 version.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
20 details.
21
22 You should have received a copy of the GNU General Public License along with
23 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
24 Street, Fifth Floor, Boston, MA 02110-1301, USA.
25 */
26
27 /*----------------------------------------------------------------------------*/
28
29 #include "cs_defs.h"
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <math.h>
37 #include <string.h>
38
39 #if defined(HAVE_MPI)
40 #include <mpi.h>
41 #endif
42
43 /*----------------------------------------------------------------------------
44 * Local headers
45 *----------------------------------------------------------------------------*/
46
47 #include "cs_headers.h"
48
49 /*----------------------------------------------------------------------------*/
50
51 BEGIN_C_DECLS
52
53 /*----------------------------------------------------------------------------*/
54 /*!
55 * \file cs_user_parameters-cdo-condif.c
56 *
57 * \brief User functions for input of calculation parameters.
58 *
59 * See \ref parameters for examples.
60 */
61 /*----------------------------------------------------------------------------*/
62
63 /*============================================================================
64 * Private function prototypes
65 *============================================================================*/
66
67 /*----------------------------------------------------------------------------*/
68 /*!
69 * \brief Give the explicit definition of the advection field.
70 * Generic function pointer for an evaluation relying on an analytic
71 * function
72 * pt_ids is optional. If not NULL, it enables to access to the coords
73 * array with an indirection. The same indirection can be applied to
74 * fill retval if dense_output is set to false.
75 *
76 * \param[in] time when ?
77 * \param[in] n_pts number of elements to consider
78 * \param[in] pt_ids list of elements ids (in coords and retval)
79 * \param[in] xyz where ? Coordinates array
80 * \param[in] dense_output perform an indirection in res or not
81 * \param[in] input NULL or pointer to a structure cast on-the-fly
82 * \param[in, out] res resulting value(s). Must be allocated.
83 */
84 /*----------------------------------------------------------------------------*/
85
86 static void
_define_adv_field(cs_real_t time,cs_lnum_t n_pts,const cs_lnum_t * pt_ids,const cs_real_t * xyz,bool dense_output,void * input,cs_real_t * res)87 _define_adv_field(cs_real_t time,
88 cs_lnum_t n_pts,
89 const cs_lnum_t *pt_ids,
90 const cs_real_t *xyz,
91 bool dense_output,
92 void *input,
93 cs_real_t *res)
94 {
95 CS_UNUSED(time);
96 CS_UNUSED(input);
97
98 for (cs_lnum_t p = 0; p < n_pts; p++) {
99
100 const cs_lnum_t id = (pt_ids == NULL) ? p : pt_ids[p];
101 const cs_lnum_t ii = dense_output ? p : id;
102 const cs_real_t *pxyz = xyz + 3*id;
103 cs_real_t *pres = res + 3*ii;
104
105 pres[0] = pxyz[1] - 0.5;
106 pres[1] = 0.5 - pxyz[0];
107 pres[2] = pxyz[2] + 1;
108
109 }
110 }
111
112 /*----------------------------------------------------------------------------*/
113 /*!
114 * \brief Give the explicit definition of the dirichlet boundary conditions
115 * Generic function pointer for an evaluation relying on an analytic
116 * function
117 * pt_ids is optional. If not NULL, it enables to access to the coords
118 * array with an indirection. The same indirection can be applied to
119 * fill retval if dense_output is set to false.
120 *
121 * \param[in] time when ?
122 * \param[in] n_pts number of elements to consider
123 * \param[in] pt_ids list of elements ids (in coords and retval)
124 * \param[in] xyz where ? Coordinates array
125 * \param[in] dense_output perform an indirection in res or not
126 * \param[in] input NULL or pointer to a structure cast on-the-fly
127 * \param[in, out] res resulting value(s). Must be allocated.
128 */
129 /*----------------------------------------------------------------------------*/
130
131 static void
_define_bcs(cs_real_t time,cs_lnum_t n_pts,const cs_lnum_t * pt_ids,const cs_real_t * xyz,bool dense_output,void * input,cs_real_t * res)132 _define_bcs(cs_real_t time,
133 cs_lnum_t n_pts,
134 const cs_lnum_t *pt_ids,
135 const cs_real_t *xyz,
136 bool dense_output,
137 void *input,
138 cs_real_t *res)
139 {
140 CS_UNUSED(time);
141 CS_UNUSED(input);
142
143 const double pi = cs_math_pi;
144 for (cs_lnum_t p = 0; p < n_pts; p++) {
145
146 const cs_lnum_t id = (pt_ids == NULL) ? p : pt_ids[p];
147 const cs_lnum_t ii = dense_output ? p : id;
148 const cs_real_t *_xyz = xyz + 3*id;
149 const double x = _xyz[0], y = _xyz[1], z = _xyz[2];
150
151 res[ii] = 1 + sin(pi*x)*sin(pi*(y+0.5))*sin(pi*(z+0.25));
152
153 }
154 }
155
156 /*============================================================================
157 * User function definitions
158 *============================================================================*/
159
160 /*----------------------------------------------------------------------------*/
161 /*!
162 * \brief Select physical model options, including user fields.
163 *
164 * This function is called at the earliest stages of the data setup,
165 * so field ids are not available yet.
166 */
167 /*----------------------------------------------------------------------------*/
168
169 void
cs_user_model(void)170 cs_user_model(void)
171 {
172 /* Activate CDO/HHO module so that main additional structure are built */
173
174 /*! [param_cdo_activation] */
175 {
176 cs_domain_t *domain = cs_glob_domain;
177
178 cs_domain_set_cdo_mode(domain, CS_DOMAIN_CDO_MODE_ONLY);
179 }
180 /*! [param_cdo_activation] */
181
182 /* ======================
183 Boundary of the domain
184 ====================== */
185
186 /*! [param_cdo_domain_boundary] */
187 {
188 cs_domain_t *domain = cs_glob_domain;
189 cs_boundary_t *bdy = domain->boundaries;
190
191 /* Choose a boundary by default */
192 cs_boundary_set_default(bdy, CS_BOUNDARY_WALL);
193
194 /* Add a new boundary
195 * >> cs_domain_add_boundary(boundary, type_of_boundary, zone_name);
196 *
197 * zone_name is either a predefined one or user-defined one
198 * type_of_boundary is one of the following keywords:
199 * CS_BOUNDARY_WALL,
200 * CS_BOUNDARY_INLET,
201 * CS_BOUNDARY_OUTLET,
202 * CS_BOUNDARY_SYMMETRY
203 */
204
205 cs_boundary_add(bdy, CS_BOUNDARY_INLET, "in");
206 cs_boundary_add(bdy, CS_BOUNDARY_OUTLET, "out");
207 }
208 /*! [param_cdo_domain_boundary] */
209
210 /* =========================
211 Generic output management
212 ========================= */
213
214 /*! [param_cdo_domain_output] */
215 {
216 cs_domain_t *domain = cs_glob_domain;
217
218 cs_domain_set_output_param(domain,
219 -1, // restart frequency
220 10, // output log frequency
221 2); // verbosity (-1: no, 0, ...)
222
223 }
224 /*! [param_cdo_domain_output] */
225
226 /* ====================
227 Time step management
228 ==================== */
229
230 /*! [param_cdo_time_step] */
231 {
232 cs_domain_t *domain = cs_glob_domain;
233
234 /*
235 If there is an inconsistency between the max. number of iteration in
236 time and the final physical time, the first condition encountered stops
237 the calculation.
238 */
239
240 cs_domain_set_time_param(domain,
241 100, // nt_max or -1 (automatic)
242 -1.); // t_max or < 0. (automatic)
243
244 /* Define the value of the time step
245 >> cs_domain_def_time_step_by_value(domain, dt_val);
246 >> cs_domain_def_time_step_by_func(domain, dt_func);
247
248 The second way to define the time step enable complex definitions.
249 */
250
251 cs_domain_def_time_step_by_value(domain, 1.0);
252
253 }
254 /*! [param_cdo_time_step] */
255
256 /*! [param_cdo_wall_distance] */
257 {
258 /* Activate predefined module as the computation of the wall distance */
259 cs_walldistance_activate();
260 }
261 /*! [param_cdo_wall_distance] */
262
263 /*! [param_cdo_add_user_equation] */
264 {
265 /* Add a new user equation.
266 * The default boundary condition has to be chosen among:
267 * CS_PARAM_BC_HMG_DIRICHLET
268 * CS_PARAM_BC_HMG_NEUMANN
269 */
270
271 cs_equation_add_user("AdvDiff.Upw", // equation name
272 "Pot.Upw", // associated variable field name
273 1, // dimension of the unknown
274 CS_PARAM_BC_HMG_DIRICHLET); // default boundary
275
276 cs_equation_add_user("AdvDiff.SG", // equation name
277 "Pot.SG", // associated variable field name
278 1, // dimension of the unknown
279 CS_PARAM_BC_HMG_DIRICHLET); // default boundary
280 }
281 /*! [param_cdo_add_user_equation] */
282
283 /* ========================================
284 Add material properties
285 ======================================== */
286
287 /*! [param_cdo_add_user_properties] */
288 {
289 cs_property_add("conductivity", /* property name */
290 CS_PROPERTY_ANISO); /* type of material property */
291 cs_property_add("rho.cp", /* property name */
292 CS_PROPERTY_ISO); /* type of material property */
293
294 }
295 /*! [param_cdo_add_user_properties] */
296
297 /*! [param_cdo_add_user_properties_opt] */
298 {
299 /* Retrieve a property named "conductivity" */
300 cs_property_t *pty = cs_property_by_name("conductivity");
301
302 /* Activate the computation of the Fourier number for this property */
303 cs_property_set_option(pty, CS_PTYKEY_POST_FOURIER);
304 }
305 /*! [param_cdo_add_user_properties_opt] */
306
307 /* ========================================
308 Add advection fields
309 ======================================== */
310
311 /*! [param_cdo_add_user_adv_field] */
312 {
313 /* Add a user-defined advection field named "adv_field" */
314 cs_adv_field_t *adv = cs_advection_field_add_user("adv_field");
315
316 CS_UNUSED(adv); /* adv can be used to set options */
317 }
318 /*! [param_cdo_add_user_adv_field] */
319
320 /*! [param_cdo_add_adv_field] */
321 {
322 /* Add a user-defined advection field named "adv_field" */
323 cs_advection_field_status_t adv_status =
324 CS_ADVECTION_FIELD_USER | /* = user-defined */
325 CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR | /* = define by a vector field */
326 CS_ADVECTION_FIELD_DEFINE_AT_VERTICES | /* = add a field at vertices */
327 CS_ADVECTION_FIELD_DEFINE_AT_BOUNDARY_FACES; /* = add boundary fluxes */
328
329 cs_adv_field_t *adv = cs_advection_field_add("adv_field", adv_status);
330
331 CS_UNUSED(adv); /* adv can be used to set options */
332 }
333 /*! [param_cdo_add_adv_field] */
334
335 /*! [param_cdo_add_user_adv_field_post] */
336 {
337 /* Retrieve an advection field named "adv_field" */
338 cs_adv_field_t *adv = cs_advection_field_by_name("adv_field");
339
340 /* Compute the Courant number (if unsteady simulation) */
341 cs_advection_field_set_postprocess(adv, CS_ADVECTION_FIELD_POST_COURANT);
342
343 }
344 /*! [param_cdo_add_user_adv_field_post] */
345
346 }
347
348 /*----------------------------------------------------------------------------*/
349 /*!
350 * \brief Define or modify general numerical and physical user parameters.
351 *
352 * At the calling point of this function, most model-related most variables
353 * and other fields have been defined, so specific settings related to those
354 * fields may be set here.
355 *
356 * \param[in, out] domain pointer to a cs_domain_t structure
357 */
358 /*----------------------------------------------------------------------------*/
359
360 void
cs_user_parameters(cs_domain_t * domain)361 cs_user_parameters(cs_domain_t *domain)
362 {
363 CS_UNUSED(domain);
364
365 /* ─────────────────────────────────────────────────────────────────────────
366 * Modify the setting of an equation using a generic process
367 * cf. the DOXYGEN documentation (website or in your installation path)
368 *
369 * cs_equation_set_param(eqp, CS_EQKEY_*, "key_val")
370 *
371 * ─────────────────────────────────────────────────────────────────────────*/
372
373 /*! [param_cdo_numerics] */
374 {
375 cs_equation_param_t *eqp = cs_equation_param_by_name("AdvDiff.Upw");
376
377 /* The modification of the space discretization should be apply first */
378 cs_equation_param_set(eqp, CS_EQKEY_SPACE_SCHEME, "cdo_vb");
379 cs_equation_param_set(eqp, CS_EQKEY_ADV_SCHEME, "upwind");
380
381 /* Modify other parameters than the space discretization */
382 cs_equation_param_set(eqp, CS_EQKEY_VERBOSITY, "2");
383 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_ALGO, "cost");
384 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_COEF, "dga");
385
386 /* Linear algebra settings */
387 #if defined(HAVE_PETSC)
388 cs_equation_param_set(eqp, CS_EQKEY_SOLVER_FAMILY, "petsc");
389 cs_equation_param_set(eqp, CS_EQKEY_ITSOL, "cg");
390 cs_equation_param_set(eqp, CS_EQKEY_PRECOND, "amg");
391 #else
392 cs_equation_param_set(eqp, CS_EQKEY_SOLVER_FAMILY, "cs");
393 cs_equation_param_set(eqp, CS_EQKEY_PRECOND, "jacobi");
394 cs_equation_param_set(eqp, CS_EQKEY_ITSOL, "cg");
395 #endif
396 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_MAX_ITER, "2500");
397 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_EPS, "1e-12");
398 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_RESNORM_TYPE, "false");
399
400 }
401 /*! [param_cdo_numerics] */
402 }
403
404 /*----------------------------------------------------------------------------*/
405 /*!
406 * \brief Advanced user-defined settings for the linear algebra related
407 * to CDO equations
408 * This is closed to cs_user_linear_solvers() but called once the fields
409 * and equations have been created (this happens at a different stage
410 * in the CDO framework)
411 */
412 /*----------------------------------------------------------------------------*/
413
414 void
cs_user_linear_solvers(void)415 cs_user_linear_solvers(void)
416 {
417 /*! [param_cdo_mg_aggreg] */
418 {
419 cs_equation_t *eq = cs_equation_by_name("AdvDiff.Upw");
420 cs_equation_param_t *eqp = cs_equation_get_param(eq);
421 cs_field_t *fld = cs_equation_get_field(eq);
422
423 /* In case of a in-house K-cylcle multigrid as a preconditioner of a
424 linear iterative solver */
425 if (eqp->sles_param->precond == CS_PARAM_PRECOND_AMG) {
426 /* If multigrid is the chosen preconditioner */
427 if (eqp->sles_param->amg_type == CS_PARAM_AMG_HOUSE_K) {
428 /* If this is a K-cycle multigrid */
429
430 /* Retrieve the different context structures to apply additional
431 settings */
432 cs_sles_t *sles = cs_sles_find_or_add(fld->id, NULL);
433 cs_sles_it_t *itsol = cs_sles_get_context(sles);
434 cs_sles_pc_t *pc = cs_sles_it_get_pc(itsol);
435 cs_multigrid_t *mg = cs_sles_pc_get_context(pc);
436
437 /* Available settings:
438 * - max. number of elements in an aggregation
439 * - type of algorithm to perform the aggregation
440 * - max. number of levels (i.e. grids)
441 * - max globalnumber of rows at the coarsest level
442 * - type of relaxation (weighting between a P_0 and P_1). For K-cycle,
443 * this should be equal to 0.
444 * - Activation of the postprocessing for the aggregation if > 0.
445 * Aggregation set is numbered by its coarse row number modulo this
446 * value
447 */
448 cs_multigrid_set_coarsening_options(mg,
449 8, /* aggregation_limit*/
450 CS_GRID_COARSENING_SPD_PW,
451 10, /* n_max_levels */
452 30, /* min_g_cells (default 30) */
453 0., /* P0P1 relaxation */
454 12); /* postprocess (default 0) */
455
456 } /* K-cycle */
457 } /* Multigrid as preconditioner */
458 }
459 /*! [param_cdo_mg_aggreg] */
460 }
461
462 /*----------------------------------------------------------------------------*/
463 /*!
464 * \brief Specify the elements such as properties, advection fields,
465 * user-defined equations and modules which have been previously
466 * added.
467 *
468 * \param[in, out] domain pointer to a cs_domain_t structure
469 */
470 /*----------------------------------------------------------------------------*/
471
472 void
cs_user_finalize_setup(cs_domain_t * domain)473 cs_user_finalize_setup(cs_domain_t *domain)
474 {
475 CS_UNUSED(domain);
476
477 /* =======================
478 User-defined properties
479 ======================= */
480
481 /*! [param_cdo_setup_property] */
482 {
483 cs_property_t *conductivity = cs_property_by_name("conductivity");
484 cs_real_33_t tensor = {{1.0, 0.5, 0.0}, {0.5, 1.0, 0.5}, {0.0, 0.5, 1.0}};
485
486 cs_property_def_aniso_by_value(conductivity, // property structure
487 "cells", // name of the volume zone
488 tensor); // values of the property
489
490 cs_property_t *rhocp = cs_property_by_name("rho.cp");
491 cs_real_t iso_val = 2.0;
492
493 cs_property_def_iso_by_value(rhocp, // property structure
494 "cells", // name of the volume zone
495 iso_val); // value of the property
496
497 }
498 /*! [param_cdo_setup_property] */
499
500 /* =============================
501 User-defined advection fields
502 ============================= */
503
504 /*! [param_cdo_setup_advfield] */
505 {
506 cs_adv_field_t *adv = cs_advection_field_by_name("adv_field");
507
508 cs_advection_field_def_by_analytic(adv, _define_adv_field, NULL);
509
510 /* Activate the post-processing of the related Courant number */
511 cs_advection_field_set_postprocess(adv, CS_ADVECTION_FIELD_POST_COURANT);
512 }
513 /*! [param_cdo_setup_advfield] */
514
515 /* ======================
516 User-defined equations
517 ====================== */
518
519 /* Define the boundary conditions
520 >> cs_equation_add_bc_by_analytic(eqp,
521 bc_type,
522 "zone_name",
523 analytic_function);
524
525 -> eq is the structure related to the equation to set
526 -> type of boundary condition:
527 CS_PARAM_BC_DIRICHLET, CS_PARAM_BC_HMG_DIRICHLET,
528 CS_PARAM_BC_NEUMANN_FULL, CS_PARAM_BC_HMG_NEUMANN, CS_PARAM_BC_ROBIN
529
530 >> cs_equation_add_bc_by_value(eqp,
531 bc_type,
532 "mesh_location_name",
533 values); // pointer
534 */
535
536 /*! [param_cdo_setup_bcs] */
537 {
538 cs_equation_param_t *eqp = cs_equation_param_by_name("AdvDiff.Upw");
539
540 cs_equation_add_bc_by_analytic(eqp,
541 CS_PARAM_BC_DIRICHLET,
542 "boundary_faces", // zone name
543 _define_bcs, // pointer to the function
544 NULL); // input structure
545
546 }
547 /*! [param_cdo_setup_bcs] */
548
549 /*! [param_cdo_add_terms] */
550 {
551 cs_equation_param_t *eqp = cs_equation_param_by_name("AdvDiff.Upw");
552 cs_property_t *rhocp = cs_property_by_name("rho.cp");
553 cs_property_t *conductivity = cs_property_by_name("conductivity");
554 cs_adv_field_t *adv = cs_advection_field_by_name("adv_field");
555
556 /* Activate the unsteady term */
557 cs_equation_add_time(eqp, rhocp);
558
559 /* Activate the diffusion term */
560 cs_equation_add_diffusion(eqp, conductivity);
561
562 /* Activate advection effect */
563 cs_equation_add_advection(eqp, adv);
564
565 /* Simple definition with cs_equation_add_source_term_by_val
566 where the value of the source term is given by m^3
567 */
568 cs_real_t st_val = -0.1;
569 cs_xdef_t *st = cs_equation_add_source_term_by_val(eqp,
570 "cells",
571 &st_val);
572
573 CS_UNUSED(st); /* st can be used for advanced settings like quadrature
574 rules */
575 }
576 /*! [param_cdo_add_terms] */
577
578 /*! [param_cdo_copy_settings] */
579 {
580 /* Copy the settings from AdvDiff.Upw to AdvDiff.SG */
581 cs_equation_param_t *eqp_ref = cs_equation_param_by_name("AdvDiff.Upw");
582
583 /* Another example to retrieve the cs_equation_param structure */
584 cs_equation_t *eq = cs_equation_by_name("AdvDiff.SG");
585 cs_equation_param_t *eqp = cs_equation_get_param(eq);
586
587 /* Copy the settings to start from a common state as in the reference
588 * equation (in this case, do not copy the associated field id since these
589 * are two different equations with their own variable field) */
590 bool copy_field_id = false;
591 cs_equation_copy_param_from(eqp_ref, eqp, copy_field_id);
592
593 /* Keep all the settings from "AdvDiff.Upw and then only change the
594 advection scheme for the second equation */
595 cs_equation_param_set(eqp, CS_EQKEY_ADV_SCHEME, "sg");
596
597 /* Call this function to be sure that the linear solver is set to what
598 one wants (if there is no modification of the SLES parameters, this
599 step can be skipped) */
600 cs_equation_param_set_sles(eqp);
601
602 }
603 /*! [param_cdo_copy_settings] */
604 }
605
606 /*----------------------------------------------------------------------------*/
607
608 END_C_DECLS
609