1 /*============================================================================
2 * Functions to handle cs_navsto_system_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
37 #if defined(HAVE_MPI)
38 #include <mpi.h>
39 #endif
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include <bft_mem.h>
46
47 #include "cs_cdofb_ac.h"
48 #include "cs_cdofb_monolithic.h"
49 #include "cs_cdofb_monolithic_sles.h"
50 #include "cs_cdofb_navsto.h"
51 #include "cs_cdofb_predco.h"
52 #include "cs_hho_stokes.h"
53 #include "cs_evaluate.h"
54 #include "cs_flag.h"
55 #include "cs_log.h"
56 #include "cs_navsto_coupling.h"
57 #include "cs_post.h"
58 #include "cs_volume_zone.h"
59
60 /*----------------------------------------------------------------------------
61 * Header for the current file
62 *----------------------------------------------------------------------------*/
63
64 #include "cs_navsto_system.h"
65
66 /*----------------------------------------------------------------------------*/
67
68 BEGIN_C_DECLS
69
70 /*=============================================================================
71 * Additional doxygen documentation
72 *============================================================================*/
73
74 /*!
75 * \file cs_navsto_system.c
76 *
77 * \brief Functions to handle the cs_navsto_system_t structure which is
78 * the high-level structure to manage the Navier-Stokes system of
79 * equations.
80 */
81
82 /*=============================================================================
83 * Local Macro definitions
84 *============================================================================*/
85
86 #define CS_NAVSTO_SYSTEM_DBG 0
87
88 /*============================================================================
89 * Type definitions
90 *============================================================================*/
91
92 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
93
94 /*============================================================================
95 * Private variables
96 *============================================================================*/
97
98 static const char _err_empty_ns[] =
99 " Stop execution. The structure related to the Navier-Stokes system is"
100 " empty.\n Please check your settings.\n";
101
102 static const char _err_invalid_coupling[] =
103 " %s: Invalid case for the coupling algorithm.\n";
104
105 static cs_navsto_system_t *cs_navsto_system = NULL;
106
107 /*============================================================================
108 * Private function prototypes
109 *============================================================================*/
110
111 /*----------------------------------------------------------------------------*/
112 /*!
113 * \brief Retrieve the \ref cs_equation_param_t structure related to the
114 * momentum equation according to the type of coupling
115 *
116 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
117 *
118 * \return a pointer to the corresponding \ref cs_equation_param_t structure
119 */
120 /*----------------------------------------------------------------------------*/
121
122 static inline bool
_handle_non_linearities(cs_navsto_param_t * nsp)123 _handle_non_linearities(cs_navsto_param_t *nsp)
124 {
125 if (nsp == NULL)
126 return false;
127
128 switch (nsp->model) {
129
130 case CS_NAVSTO_MODEL_OSEEN:
131 case CS_NAVSTO_MODEL_STOKES:
132 return false;
133
134 case CS_NAVSTO_MODEL_INCOMPRESSIBLE_NAVIER_STOKES:
135 if (nsp->adv_strategy == CS_PARAM_ADVECTION_IMPLICIT_FULL)
136 return true;
137 else {
138 assert(nsp->adv_strategy == CS_PARAM_ADVECTION_IMPLICIT_LINEARIZED ||
139 nsp->adv_strategy == CS_PARAM_ADVECTION_EXPLICIT);
140 return false;
141 }
142 break;
143
144 default:
145 return true;
146
147 }
148 }
149
150 /*----------------------------------------------------------------------------*/
151 /*!
152 * \brief Allocate an empty Navier-Stokes system
153 *
154 * \return a pointer to a new allocated groundwater flow structure
155 */
156 /*----------------------------------------------------------------------------*/
157
158 static cs_navsto_system_t *
_allocate_navsto_system(void)159 _allocate_navsto_system(void)
160 {
161 cs_navsto_system_t *navsto = NULL;
162
163 BFT_MALLOC(navsto, 1, cs_navsto_system_t);
164
165 navsto->param = NULL;
166
167 /* Array of boundary type */
168
169 navsto->bf_type = NULL;
170
171 /* Main set of variables */
172
173 navsto->velocity = NULL;
174 navsto->pressure = NULL;
175
176 /* Advection field */
177
178 navsto->adv_field = NULL;
179 navsto->mass_flux_array = NULL;
180 navsto->mass_flux_array_pre = NULL;
181
182 /* Related modules */
183
184 navsto->turbulence = NULL;
185
186 /* Post-processing fields */
187
188 navsto->plot_writer = NULL;
189 navsto->velocity_divergence = NULL;
190 navsto->pressure_gradient = NULL;
191 navsto->mass_density = NULL;
192 navsto->mass_flux_balance = NULL;
193 navsto->kinetic_energy = NULL;
194 navsto->velocity_gradient = NULL;
195 navsto->vorticity = NULL;
196 navsto->helicity = NULL;
197 navsto->enstrophy = NULL;
198
199 /* Stream function is associated to the variable field of an equation
200 So the treatment is different */
201
202 navsto->stream_function_eq = NULL;
203
204 /* Additional data fitting the choice of the coupling model */
205
206 navsto->coupling_context = NULL;
207 navsto->scheme_context = NULL;
208
209 /* Function pointers */
210
211 navsto->init_scheme_context = NULL;
212 navsto->free_scheme_context = NULL;
213 navsto->init_velocity = NULL;
214 navsto->init_pressure = NULL;
215 navsto->compute_steady = NULL;
216 navsto->compute= NULL;
217
218 return navsto;
219 }
220
221 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
222
223 /*============================================================================
224 * Public function prototypes
225 *============================================================================*/
226
227 /*----------------------------------------------------------------------------*/
228 /*!
229 * \brief Check if the resolution of the Navier-Stokes system has been
230 * activated
231 *
232 * \return true or false
233 */
234 /*----------------------------------------------------------------------------*/
235
236 bool
cs_navsto_system_is_activated(void)237 cs_navsto_system_is_activated(void)
238 {
239 if (cs_navsto_system == NULL)
240 return false;
241 else
242 return true;
243 }
244
245 /*----------------------------------------------------------------------------*/
246 /*!
247 * \brief Update the flag associated to the modelling options
248 *
249 * \param[in] with_thermal true or false
250 */
251 /*----------------------------------------------------------------------------*/
252
253 void
cs_navsto_system_update_model(bool with_thermal)254 cs_navsto_system_update_model(bool with_thermal)
255 {
256 if (cs_navsto_system == NULL)
257 bft_error(__FILE__, __LINE__, 0,
258 " %s: The main structure for the Navier-Stokes equations is not"
259 " allocated.\n"
260 " Please check your settings", __func__);
261
262 cs_navsto_param_t *nsp = cs_navsto_system->param;
263
264 if (with_thermal) { /* Thermal system is switch on and relies on the mass flux
265 for the advection */
266
267 if ((nsp->model_flag & (CS_NAVSTO_MODEL_PASSIVE_THERMAL_TRACER |
268 CS_NAVSTO_MODEL_BOUSSINESQ)) == 0) {
269
270 /* Thermal system is linked to the Navier-Stokes one but nothing has been
271 * set. Add the "minimal" flag. */
272
273 nsp->model_flag |= CS_NAVSTO_MODEL_PASSIVE_THERMAL_TRACER;
274
275 }
276
277 }
278 }
279
280 /*----------------------------------------------------------------------------*/
281 /*!
282 * \brief Allocate and initialize the Navier-Stokes (NS) system
283 *
284 * \param[in] boundaries pointer to the domain boundaries
285 * \param[in] model type of model related to the NS system
286 * \param[in] model_flag additional high-level model options
287 * \param[in] algo_coupling algorithm used for solving the NS system
288 * \param[in] post_flag predefined post-processings options
289 *
290 * \return a pointer to a new allocated cs_navsto_system_t structure
291 */
292 /*----------------------------------------------------------------------------*/
293
294 cs_navsto_system_t *
cs_navsto_system_activate(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)295 cs_navsto_system_activate(const cs_boundary_t *boundaries,
296 cs_navsto_param_model_t model,
297 cs_navsto_param_model_flag_t model_flag,
298 cs_navsto_param_coupling_t algo_coupling,
299 cs_navsto_param_post_flag_t post_flag)
300 {
301 if (model == CS_NAVSTO_N_MODELS)
302 bft_error(__FILE__, __LINE__, 0, "%s: Invalid model for Navier-Stokes.\n",
303 __func__);
304
305 /* Allocate an empty structure */
306
307 cs_navsto_system_t *navsto = _allocate_navsto_system();
308
309 /* Initialize the set of parameters */
310
311 navsto->param = cs_navsto_param_create(boundaries,
312 model,
313 model_flag,
314 algo_coupling,
315 post_flag);
316
317 /* Set the default boundary condition for the equations of the Navier-Stokes
318 system according to the default domain boundary */
319
320 cs_param_bc_type_t default_bc = CS_PARAM_N_BC_TYPES;
321 switch (boundaries->default_type) {
322
323 case CS_BOUNDARY_WALL:
324 default_bc = CS_PARAM_BC_HMG_DIRICHLET;
325 break;
326 case CS_BOUNDARY_SYMMETRY:
327 default_bc = CS_PARAM_BC_SLIDING;
328 break;
329
330 default:
331 bft_error(__FILE__, __LINE__, 0, " %s: Invalid boundary default type\n",
332 __func__);
333 break;
334
335 } /* End of switch */
336
337 /* Advection field related to the resolved velocity */
338
339 cs_advection_field_status_t adv_status =
340 CS_ADVECTION_FIELD_NAVSTO | CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX;
341
342 if (cs_navsto_param_is_steady(navsto->param))
343 adv_status |= CS_ADVECTION_FIELD_STEADY;
344
345 navsto->adv_field = cs_advection_field_add("mass_flux", adv_status);
346
347 /* Additional initialization fitting the choice of model */
348
349 switch (navsto->param->coupling) {
350
351 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
352 navsto->coupling_context =
353 cs_navsto_ac_create_context(default_bc, navsto->param);
354 break;
355 case CS_NAVSTO_COUPLING_MONOLITHIC:
356 navsto->coupling_context =
357 cs_navsto_monolithic_create_context(default_bc, navsto->param);
358 break;
359 case CS_NAVSTO_COUPLING_PROJECTION:
360 navsto->coupling_context =
361 cs_navsto_projection_create_context(default_bc, navsto->param);
362 break;
363
364 default:
365 bft_error(__FILE__, __LINE__, 0, _err_invalid_coupling, __func__);
366 break;
367
368 }
369
370 if (post_flag & CS_NAVSTO_POST_STREAM_FUNCTION) {
371
372 navsto->stream_function_eq = cs_equation_add(CS_NAVSTO_STREAM_EQNAME,
373 "stream_function",
374 CS_EQUATION_TYPE_NAVSTO,
375 1,
376 CS_PARAM_BC_HMG_NEUMANN);
377
378 cs_equation_param_t *eqp =
379 cs_equation_get_param(navsto->stream_function_eq);
380 assert(eqp != NULL);
381
382 /* Default settings for this equation */
383
384 cs_equation_param_set(eqp, CS_EQKEY_SPACE_SCHEME, "cdo_vb");
385 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_COEF, "dga");
386 cs_equation_param_set(eqp, CS_EQKEY_PRECOND, "amg");
387 cs_equation_param_set(eqp, CS_EQKEY_AMG_TYPE, "k_cycle");
388 cs_equation_param_set(eqp, CS_EQKEY_ITSOL, "cg");
389
390 /* This is for post-processing purpose, so, there is no need to have
391 * a restrictive convergence tolerance on the resolution of the linear
392 * system */
393
394 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_EPS, "1e-6");
395
396 }
397
398 /* Create the main structure to handle the turbulence modelling */
399
400 navsto->turbulence = cs_turbulence_create(navsto->param->turbulence);
401
402 /* Set the static variable */
403
404 cs_navsto_system = navsto;
405
406 return navsto;
407 }
408
409 /*----------------------------------------------------------------------------*/
410 /*!
411 * \brief Free the main structure related to the Navier-Stokes system
412 */
413 /*----------------------------------------------------------------------------*/
414
415 void
cs_navsto_system_destroy(void)416 cs_navsto_system_destroy(void)
417 {
418 cs_navsto_system_t *navsto = cs_navsto_system;
419
420 if (navsto == NULL)
421 return;
422
423 BFT_FREE(navsto->bf_type);
424
425 /*
426 * Properties, advection fields, equations and fields are all destroyed
427 * respectively inside cs_property_destroy_all(),
428 * cs_advection_field_destroy_all(), cs_equation_destroy_all() and
429 * cs_field_destroy_all()
430 */
431
432 BFT_FREE(navsto->mass_flux_array);
433 BFT_FREE(navsto->mass_flux_array_pre);
434
435 /* Free the plot writer */
436
437 if (navsto->plot_writer != NULL)
438 cs_time_plot_finalize(&navsto->plot_writer);
439
440 cs_navsto_param_t *nsp = navsto->param;
441
442 /* Free the context according to the model choice */
443
444 switch (nsp->coupling) {
445
446 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
447 navsto->coupling_context =
448 cs_navsto_ac_free_context(navsto->coupling_context);
449 break;
450 case CS_NAVSTO_COUPLING_MONOLITHIC:
451 navsto->coupling_context =
452 cs_navsto_monolithic_free_context(navsto->coupling_context);
453 if (nsp->space_scheme == CS_SPACE_SCHEME_CDOFB)
454 cs_cdofb_monolithic_finalize_common(nsp);
455 break;
456 case CS_NAVSTO_COUPLING_PROJECTION:
457 navsto->coupling_context =
458 cs_navsto_projection_free_context(navsto->coupling_context);
459 break;
460
461 default:
462 bft_error(__FILE__, __LINE__, 0, _err_invalid_coupling, __func__);
463 break;
464 }
465
466 /* Destroy the context related to the discretization scheme */
467
468 navsto->free_scheme_context(navsto->scheme_context);
469
470 /* Free the context and the structure related to the turbulence modelling */
471
472 cs_turbulence_free(&(navsto->turbulence));
473
474 /* Set of numerical parameters */
475
476 navsto->param = cs_navsto_param_free(nsp);
477
478 BFT_FREE(navsto);
479 cs_navsto_system = NULL;
480 }
481
482 /*----------------------------------------------------------------------------*/
483 /*!
484 * \brief Retrieve the structure storing the parameters for the Navier--Stokes
485 * system
486 *
487 * \return NULL or the pointer to a \ref cs_navsto_param_t structure
488 */
489 /*----------------------------------------------------------------------------*/
490
491 cs_navsto_param_t *
cs_navsto_system_get_param(void)492 cs_navsto_system_get_param(void)
493 {
494 cs_navsto_system_t *navsto = cs_navsto_system;
495
496 if (navsto == NULL)
497 return NULL;
498
499 return navsto->param;
500 }
501
502 /*----------------------------------------------------------------------------*/
503 /*!
504 * \brief Retrieve a pointer to the equation related to the momentum equation
505 *
506 * \return NULL or the pointer
507 */
508 /*----------------------------------------------------------------------------*/
509
510 cs_equation_t *
cs_navsto_system_get_momentum_eq(void)511 cs_navsto_system_get_momentum_eq(void)
512 {
513 cs_navsto_system_t *navsto = cs_navsto_system;
514
515 if (navsto == NULL)
516 return NULL;
517
518 cs_navsto_param_t *nsp = navsto->param;
519 cs_equation_t *eq = NULL;
520
521 switch (nsp->coupling) {
522
523 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
524 eq = cs_navsto_ac_get_momentum_eq(navsto->coupling_context);
525 break;
526 case CS_NAVSTO_COUPLING_MONOLITHIC:
527 eq = cs_navsto_monolithic_get_momentum_eq(navsto->coupling_context);
528 break;
529 case CS_NAVSTO_COUPLING_PROJECTION:
530 eq = cs_navsto_projection_get_momentum_eq(navsto->coupling_context);
531 break;
532
533 default:
534 bft_error(__FILE__, __LINE__, 0, _err_invalid_coupling, __func__);
535 break;
536
537 }
538
539 return eq;
540 }
541
542 /*----------------------------------------------------------------------------*/
543 /*!
544 * \brief Retrieve the advection field structure (the mass flux) related to
545 * the Navier-Stokes system.
546 *
547 * \return a pointer to the advection field structure
548 */
549 /*----------------------------------------------------------------------------*/
550
551 cs_adv_field_t *
cs_navsto_get_adv_field(void)552 cs_navsto_get_adv_field(void)
553 {
554 cs_navsto_system_t *ns = cs_navsto_system;
555
556 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
557
558 return ns->adv_field;
559 }
560
561 /*----------------------------------------------------------------------------*/
562 /*!
563 * \brief Retrieve the mass flux array related to the Navier-Stokes system.
564 *
565 * \param[in] previous if true return the previous state otherwise the
566 * current state.
567 *
568 * \return a pointer to an array of cs_real_t
569 */
570 /*----------------------------------------------------------------------------*/
571
572 cs_real_t *
cs_navsto_get_mass_flux(bool previous)573 cs_navsto_get_mass_flux(bool previous)
574 {
575 cs_navsto_system_t *ns = cs_navsto_system;
576
577 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
578
579 cs_real_t *mass_flux = ns->mass_flux_array;
580 if (previous)
581 mass_flux = ns->mass_flux_array_pre;
582
583 return mass_flux;
584 }
585
586 /*----------------------------------------------------------------------------*/
587 /*!
588 * \brief Start setting-up the Navier-Stokes system
589 * At this stage, numerical settings should be completely determined
590 * but connectivity and geometrical information is not yet available.
591 */
592 /*----------------------------------------------------------------------------*/
593
594 void
cs_navsto_system_init_setup(void)595 cs_navsto_system_init_setup(void)
596 {
597 cs_navsto_system_t *ns = cs_navsto_system;
598
599 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
600
601 cs_navsto_param_t *nsp = ns->param;
602
603 /* Plotter output:
604 * By default: there is at least the norm of the velocity field divergence */
605
606 const char **labels;
607 int n_plotter_outputs = 1;
608
609 /* Set field metadata */
610
611 const int log_key = cs_field_key_id("log");
612 const int post_key = cs_field_key_id("post_vis");
613 const bool has_previous = cs_navsto_param_is_steady(nsp) ? false : true;
614 int field_mask = CS_FIELD_INTENSIVE | CS_FIELD_VARIABLE | CS_FIELD_CDO;
615
616 /* Set the location id to define a mesh location support */
617
618 int location_id = -1;
619 switch (nsp->space_scheme) {
620
621 case CS_SPACE_SCHEME_CDOFB:
622 case CS_SPACE_SCHEME_HHO_P0:
623 case CS_SPACE_SCHEME_HHO_P1:
624 case CS_SPACE_SCHEME_HHO_P2:
625 location_id = cs_mesh_location_get_id_by_name("cells");
626 break; /* Face-based scheme family */
627
628 default:
629 bft_error(__FILE__, __LINE__, 0,
630 "%s: Invalid space discretization scheme.", __func__);
631 }
632
633 /* Create if needed velocity and pressure fields */
634
635 const int field_post_flag = CS_POST_ON_LOCATION | CS_POST_MONITOR;
636
637 /* Handle the velocity field */
638
639 ns->velocity = cs_field_find_or_create("velocity",
640 field_mask,
641 location_id,
642 3, /* dimension */
643 has_previous);
644
645 /* Set the default value for keys related to log and post-processing */
646
647 cs_field_set_key_int(ns->velocity, log_key, 1);
648 cs_field_set_key_int(ns->velocity, post_key, field_post_flag);
649
650 /* Handle the pressure field */
651
652 ns->pressure = cs_field_find_or_create("pressure",
653 field_mask,
654 location_id,
655 1, /* dimension */
656 has_previous);
657
658 /* Set the default value for keys related to log and post-processing */
659
660 cs_field_set_key_int(ns->pressure, log_key, 1);
661 cs_field_set_key_int(ns->pressure, post_key, field_post_flag);
662
663 /* Handle the divergence of the velocity field.
664 * Up to now, always defined the divergence of the velocity field. This
665 * should be changed in the future */
666
667 int p_mask = CS_FIELD_INTENSIVE | CS_FIELD_PROPERTY | CS_FIELD_CDO;
668
669 nsp->post_flag |= CS_NAVSTO_POST_VELOCITY_DIVERGENCE;
670
671 /* Always post-process the velocity divergence */
672
673 ns->velocity_divergence = cs_field_find_or_create("velocity_divergence",
674 p_mask,
675 location_id,
676 1, /* dimension */
677 has_previous);
678
679 /* Set default value for keys related to log and post-processing */
680
681 cs_field_set_key_int(ns->velocity_divergence, log_key, 1);
682 cs_field_set_key_int(ns->velocity_divergence, post_key, field_post_flag);
683
684 if (nsp->post_flag & CS_NAVSTO_POST_KINETIC_ENERGY) {
685
686 n_plotter_outputs += 1; /* Integral of the kinetic energy is monitored */
687
688 ns->kinetic_energy = cs_field_find_or_create("kinetic_energy",
689 p_mask,
690 location_id,
691 1, /* dimension */
692 has_previous);
693
694 /* Set the default value for keys related to log and post-processing */
695
696 cs_field_set_key_int(ns->kinetic_energy, log_key, 1);
697 cs_field_set_key_int(ns->kinetic_energy, post_key, field_post_flag);
698
699 }
700
701 if (nsp->post_flag & CS_NAVSTO_POST_MASS_DENSITY) {
702
703 n_plotter_outputs += 1; /* Integral of the mass density is monitored */
704
705 ns->mass_density = cs_field_find_or_create("mass_density",
706 p_mask,
707 location_id,
708 1, /* dimension */
709 has_previous);
710
711 /* Set the default value for keys related to log and post-processing */
712
713 cs_field_set_key_int(ns->mass_density, log_key, 1);
714 cs_field_set_key_int(ns->mass_density, post_key, field_post_flag);
715
716 }
717
718 if (nsp->post_flag & CS_NAVSTO_POST_CELL_MASS_FLUX_BALANCE) {
719
720 ns->mass_flux_balance = cs_field_find_or_create("mass_flux_balance",
721 p_mask,
722 location_id,
723 1, /* dimension */
724 has_previous);
725
726 /* Set the default value for keys related to log and post-processing */
727
728 cs_field_set_key_int(ns->mass_flux_balance, log_key, 1);
729 cs_field_set_key_int(ns->mass_flux_balance, post_key, field_post_flag);
730
731 }
732
733 if (nsp->post_flag & CS_NAVSTO_POST_PRESSURE_GRADIENT) {
734
735 ns->pressure_gradient = cs_field_find_or_create("pressure_gradient",
736 p_mask,
737 location_id,
738 3, /* dimension */
739 has_previous);
740
741 /* Set the default value for keys related to post-processing */
742
743 cs_field_set_key_int(ns->pressure_gradient, post_key, field_post_flag);
744
745 }
746
747 if (nsp->post_flag & CS_NAVSTO_POST_STREAM_FUNCTION)
748 nsp->post_flag |= CS_NAVSTO_POST_VORTICITY; /* automatic */
749
750 if (nsp->post_flag & CS_NAVSTO_POST_HELICITY) {
751
752 n_plotter_outputs += 1; /* Integral of the helicity is monitored */
753
754 nsp->post_flag |= CS_NAVSTO_POST_VORTICITY; /* automatic */
755 ns->helicity = cs_field_find_or_create("helicity",
756 p_mask,
757 location_id,
758 1, /* dimension */
759 has_previous);
760
761 /* Set default value for keys related to log and post-processing */
762
763 cs_field_set_key_int(ns->helicity, log_key, 1);
764 cs_field_set_key_int(ns->helicity, post_key, field_post_flag);
765
766 }
767
768 if (nsp->post_flag & CS_NAVSTO_POST_ENSTROPHY) {
769
770 n_plotter_outputs += 1; /* Integral of the enstrophy is monitored */
771
772 nsp->post_flag |= CS_NAVSTO_POST_VORTICITY; /* automatic */
773 ns->enstrophy = cs_field_find_or_create("enstrophy",
774 p_mask,
775 location_id,
776 1, /* dimension */
777 has_previous);
778
779 /* Set default value for keys related to log and post-processing */
780
781 cs_field_set_key_int(ns->enstrophy, log_key, 1);
782 cs_field_set_key_int(ns->enstrophy, post_key, field_post_flag);
783
784 }
785
786 if (nsp->post_flag & CS_NAVSTO_POST_VORTICITY) {
787
788 ns->vorticity = cs_field_find_or_create("vorticity",
789 p_mask,
790 location_id,
791 3, /* dimension */
792 has_previous);
793
794 /* Set default value for keys related to log and post-processing */
795
796 cs_field_set_key_int(ns->vorticity, log_key, 1);
797 cs_field_set_key_int(ns->vorticity, post_key, field_post_flag);
798
799 }
800
801 if (nsp->post_flag & CS_NAVSTO_POST_VELOCITY_GRADIENT) {
802
803 ns->velocity_gradient = cs_field_find_or_create("velocity_gradient",
804 p_mask,
805 location_id,
806 9, /* dimension */
807 has_previous);
808
809 /* Set default value for keys related to log and post-processing */
810
811 cs_field_set_key_int(ns->velocity_gradient, log_key, 1);
812 cs_field_set_key_int(ns->velocity_gradient, post_key, field_post_flag);
813
814 }
815
816 /* Time plot monitor of global quantities predefined by the automatic
817 post-processing */
818
819 if (cs_glob_rank_id < 1) {
820
821 assert(n_plotter_outputs > 0);
822 BFT_MALLOC(labels, n_plotter_outputs, const char *);
823
824 int n_cols = 0;
825
826 if (nsp->post_flag & CS_NAVSTO_POST_VELOCITY_DIVERGENCE)
827 labels[n_cols++] = "vel_div_norm";
828 if (nsp->post_flag & CS_NAVSTO_POST_MASS_DENSITY)
829 labels[n_cols++] = "mass";
830 if (nsp->post_flag & CS_NAVSTO_POST_KINETIC_ENERGY)
831 labels[n_cols++] = "kinetic_energy";
832 if (nsp->post_flag & CS_NAVSTO_POST_ENSTROPHY)
833 labels[n_cols++] = "enstrophy";
834 if (nsp->post_flag & CS_NAVSTO_POST_HELICITY)
835 labels[n_cols++] = "helicity";
836
837 ns->plot_writer = cs_time_plot_init_probe("navsto_monitor",
838 "",
839 CS_TIME_PLOT_DAT,
840 false, /* use iteration */
841 300, /* flush time */
842 -1,
843 n_cols,
844 NULL,
845 NULL,
846 labels);
847
848 BFT_FREE(labels);
849
850 } /* monitoring */
851
852 /* Setup data according to the type of coupling (add the advection field in
853 case of Navier-Stokes equations) */
854
855 switch (nsp->coupling) {
856
857 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
858 cs_navsto_ac_init_setup(nsp, ns->adv_field, ns->coupling_context);
859 break;
860 case CS_NAVSTO_COUPLING_MONOLITHIC:
861 cs_navsto_monolithic_init_setup(nsp, ns->adv_field, ns->coupling_context);
862 break;
863 case CS_NAVSTO_COUPLING_PROJECTION:
864 cs_navsto_projection_init_setup(nsp,
865 ns->adv_field,
866 location_id,
867 has_previous,
868 ns->coupling_context);
869 break;
870
871 default:
872 bft_error(__FILE__, __LINE__, 0, _err_invalid_coupling, __func__);
873 break;
874
875 }
876
877 /* Initialize the turbulence modelling */
878
879 cs_equation_t *mom_eq = cs_navsto_system_get_momentum_eq();
880 cs_turbulence_init_setup(ns->turbulence, mom_eq);
881 }
882
883 /*----------------------------------------------------------------------------*/
884 /*!
885 * \brief Define the settings for SLES related to the Navier-Stokes system
886 */
887 /*----------------------------------------------------------------------------*/
888
889 void
cs_navsto_system_set_sles(void)890 cs_navsto_system_set_sles(void)
891 {
892 cs_navsto_system_t *ns = cs_navsto_system;
893
894 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
895
896 void *nscc = ns->coupling_context;
897 cs_navsto_param_t *nsp = ns->param;
898
899 switch (nsp->space_scheme) {
900
901 case CS_SPACE_SCHEME_CDOFB:
902 case CS_SPACE_SCHEME_HHO_P0:
903 switch (nsp->coupling) {
904
905 case CS_NAVSTO_COUPLING_MONOLITHIC:
906 cs_cdofb_monolithic_set_sles(nsp, nscc);
907 break;
908
909 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
910 cs_cdofb_ac_set_sles(nsp, nscc);
911 break;
912
913 case CS_NAVSTO_COUPLING_PROJECTION:
914 cs_cdofb_predco_set_sles(nsp, nscc);
915 break;
916
917 default:
918 bft_error(__FILE__, __LINE__, 0, _err_invalid_coupling, __func__);
919 break;
920
921 } /* Switch algo. for coupling velocity/pressure */
922 break; /* Face-based scheme family */
923
924 default:
925 bft_error(__FILE__, __LINE__, 0,
926 "%s: Invalid space discretization scheme.", __func__);
927
928 } /* Switch space scheme */
929
930 if (nsp->post_flag & CS_NAVSTO_POST_STREAM_FUNCTION) {
931
932 cs_equation_param_t *eqp = cs_equation_get_param(ns->stream_function_eq);
933 assert(eqp != NULL);
934
935 /* Equation related to Navier-Stokes do not follow the classical setup
936 stage */
937 cs_equation_param_set_sles(eqp);
938
939 }
940 }
941
942 /*----------------------------------------------------------------------------*/
943 /*!
944 * \brief Last step of the setup of the Navier-Stokes system
945 *
946 * \param[in] mesh pointer to a cs_mesh_t structure
947 * \param[in] connect pointer to a cs_cdo_connect_t structure
948 * \param[in] quant pointer to a cs_cdo_quantities_t structure
949 * \param[in] time_step pointer to a cs_time_step_t structure
950 */
951 /*----------------------------------------------------------------------------*/
952
953 void
cs_navsto_system_finalize_setup(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)954 cs_navsto_system_finalize_setup(const cs_mesh_t *mesh,
955 const cs_cdo_connect_t *connect,
956 const cs_cdo_quantities_t *quant,
957 const cs_time_step_t *time_step)
958 {
959 cs_navsto_system_t *ns = cs_navsto_system;
960
961 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
962
963 cs_navsto_param_t *nsp = ns->param;
964 assert(connect != NULL && quant != NULL && nsp != NULL);
965
966 /* Setup checkings */
967
968 if ((nsp->model_flag & (CS_NAVSTO_MODEL_BOUSSINESQ |
969 CS_NAVSTO_MODEL_PASSIVE_THERMAL_TRACER)) > 0) {
970
971 if (cs_thermal_system_is_activated() == false)
972 bft_error(__FILE__, __LINE__, 0,
973 " %s: The Navier-Stokes module is activated with options"
974 " that imply a thermal module.\n"
975 " Please check that cs_thermal_system_activate() has been"
976 " called.\n", __func__);
977
978 /* Check that an advection term has been set */
979
980 cs_equation_param_t *eqp = cs_equation_param_by_name(CS_THERMAL_EQNAME);
981 if (cs_equation_param_has_convection(eqp) == false)
982 bft_error(__FILE__, __LINE__, 0,
983 " %s: No advection field is associated with the thermal"
984 " equation\n whereas the Navier-Stokes is associated with"
985 " the thermal equation.\n"
986 " Please check your settings.", __func__);
987
988 }
989
990 /* Remaining boundary conditions:
991 * 1. Walls
992 * 2. Symmetries
993 * 3. Outlets
994 */
995
996 cs_navsto_set_fixed_walls(nsp);
997 cs_navsto_set_symmetries(nsp);
998 cs_navsto_set_outlets(nsp);
999
1000 /* Define the advection field which relies on the mass flux. */
1001
1002 switch (nsp->space_scheme) {
1003
1004 case CS_SPACE_SCHEME_CDOFB:
1005 case CS_SPACE_SCHEME_HHO_P0:
1006 {
1007 BFT_MALLOC(ns->mass_flux_array, quant->n_faces, cs_real_t);
1008 memset(ns->mass_flux_array, 0, sizeof(cs_real_t)*quant->n_faces);
1009
1010 BFT_MALLOC(ns->mass_flux_array_pre, quant->n_faces, cs_real_t);
1011 memset(ns->mass_flux_array_pre, 0, sizeof(cs_real_t)*quant->n_faces);
1012
1013 cs_flag_t loc_flag =
1014 CS_FLAG_FULL_LOC | CS_FLAG_SCALAR | cs_flag_primal_face;
1015
1016 cs_advection_field_def_by_array(ns->adv_field,
1017 loc_flag,
1018 ns->mass_flux_array,
1019 false, /* advection field is not owner */
1020 NULL); /* index (not useful here) */
1021
1022 }
1023 break;
1024
1025 default:
1026 bft_error(__FILE__, __LINE__, 0,
1027 "%s: Invalid space discretization scheme.", __func__);
1028
1029 } /* End of switch on space_scheme */
1030
1031 /* Last setup stage according to the type of coupling (not related to
1032 space discretization scheme */
1033
1034 switch (nsp->coupling) {
1035
1036 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
1037 cs_navsto_ac_last_setup(nsp, ns->coupling_context);
1038 break;
1039 case CS_NAVSTO_COUPLING_MONOLITHIC:
1040 cs_navsto_monolithic_last_setup(nsp, ns->coupling_context);
1041 break;
1042 case CS_NAVSTO_COUPLING_PROJECTION:
1043 cs_navsto_projection_last_setup(quant, nsp, ns->coupling_context);
1044 break;
1045
1046 default:
1047 bft_error(__FILE__, __LINE__, 0, _err_invalid_coupling, __func__);
1048 break;
1049
1050 }
1051
1052 /* Settings with respect to the discretization scheme */
1053
1054 switch (nsp->space_scheme) {
1055
1056 case CS_SPACE_SCHEME_CDOFB:
1057 case CS_SPACE_SCHEME_HHO_P0:
1058
1059 /* Setup functions and data according to the type of coupling */
1060
1061 switch (nsp->coupling) {
1062
1063 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
1064 /* ============================================= */
1065
1066 ns->init_scheme_context = cs_cdofb_ac_init_scheme_context;
1067 ns->free_scheme_context = cs_cdofb_ac_free_scheme_context;
1068 ns->init_velocity = NULL;
1069 ns->init_pressure = cs_cdofb_navsto_init_pressure;
1070 ns->compute_steady = NULL;
1071
1072 if (_handle_non_linearities(nsp)) {
1073
1074 /* Same function. The difference will be taken care of by mean of a
1075 * build sub-function */
1076
1077 if (nsp->time_scheme == CS_TIME_SCHEME_EULER_IMPLICIT)
1078 ns->compute = cs_cdofb_ac_compute_implicit_nl;
1079 else
1080 bft_error(__FILE__, __LINE__, 0,
1081 "%s: Invalid time scheme for the AC and Navier-Stokes",
1082 __func__);
1083 }
1084 else {
1085
1086 switch (nsp->time_scheme) {
1087
1088 case CS_TIME_SCHEME_EULER_IMPLICIT:
1089 ns->compute = cs_cdofb_ac_compute_implicit;
1090 break;
1091
1092 case CS_TIME_SCHEME_THETA:
1093 case CS_TIME_SCHEME_CRANKNICO:
1094 case CS_TIME_SCHEME_BDF2:
1095 bft_error(__FILE__, __LINE__, 0,
1096 "%s: Time scheme not implemented for the AC coupling",
1097 __func__);
1098 break;
1099
1100 case CS_TIME_SCHEME_STEADY:
1101 default:
1102 bft_error(__FILE__, __LINE__, 0,
1103 "%s: Invalid time scheme for the AC coupling", __func__);
1104 break;
1105
1106 } /* Switch */
1107
1108 } /* If nonlinear */
1109
1110 cs_cdofb_ac_init_common(quant, connect, time_step);
1111 break;
1112
1113 case CS_NAVSTO_COUPLING_MONOLITHIC:
1114 /* ============================= */
1115
1116 ns->init_scheme_context = cs_cdofb_monolithic_init_scheme_context;
1117 ns->free_scheme_context = cs_cdofb_monolithic_free_scheme_context;
1118 ns->init_velocity = NULL;
1119 ns->init_pressure = cs_cdofb_navsto_init_pressure;
1120 if (_handle_non_linearities(nsp))
1121 ns->compute_steady = cs_cdofb_monolithic_steady_nl;
1122 else
1123 ns->compute_steady = cs_cdofb_monolithic_steady;
1124
1125 switch (nsp->time_scheme) {
1126
1127 case CS_TIME_SCHEME_STEADY:
1128 if (_handle_non_linearities(nsp))
1129 ns->compute = cs_cdofb_monolithic_steady_nl;
1130 else
1131 ns->compute = cs_cdofb_monolithic_steady;
1132 break; /* Nothing to set */
1133
1134 case CS_TIME_SCHEME_EULER_IMPLICIT:
1135 case CS_TIME_SCHEME_THETA:
1136 case CS_TIME_SCHEME_CRANKNICO:
1137 if (_handle_non_linearities(nsp))
1138 ns->compute = cs_cdofb_monolithic_nl;
1139 else
1140 ns->compute = cs_cdofb_monolithic;
1141 break;
1142
1143 case CS_TIME_SCHEME_BDF2:
1144 default:
1145 bft_error(__FILE__, __LINE__, 0,
1146 "%s: Invalid time scheme for the monolithic coupling",
1147 __func__);
1148 break;
1149
1150 } /* Switch */
1151
1152 cs_cdofb_monolithic_init_common(nsp, mesh, quant, connect, time_step);
1153 break;
1154
1155 case CS_NAVSTO_COUPLING_PROJECTION:
1156 /* ============================= */
1157
1158 ns->init_scheme_context = cs_cdofb_predco_init_scheme_context;
1159 ns->free_scheme_context = cs_cdofb_predco_free_scheme_context;
1160 ns->init_velocity = NULL;
1161 ns->init_pressure = cs_cdofb_navsto_init_pressure;
1162 ns->compute_steady = NULL;
1163
1164 switch (nsp->time_scheme) {
1165
1166 case CS_TIME_SCHEME_STEADY:
1167 bft_error(__FILE__, __LINE__, 0,
1168 "%s: The projection coupling algorithm "
1169 "can be used only in unsteady problems", __func__);
1170 break;
1171
1172 case CS_TIME_SCHEME_EULER_IMPLICIT:
1173 ns->compute = cs_cdofb_predco_compute_implicit;
1174 break;
1175
1176 case CS_TIME_SCHEME_THETA:
1177 case CS_TIME_SCHEME_CRANKNICO:
1178 case CS_TIME_SCHEME_BDF2:
1179 default:
1180 bft_error(__FILE__, __LINE__, 0,
1181 "%s: Invalid time scheme for the projection coupling"
1182 " algorithm", __func__);
1183 break;
1184
1185 } /* Switch */
1186
1187 cs_cdofb_predco_init_common(quant, connect, time_step);
1188 break;
1189
1190 default:
1191 bft_error(__FILE__, __LINE__, 0, _err_invalid_coupling, __func__);
1192 break;
1193
1194 }
1195 break; /* Lowest-order face-based schemes */
1196
1197 case CS_SPACE_SCHEME_HHO_P1:
1198 case CS_SPACE_SCHEME_HHO_P2:
1199 /* TODO: set function pointers */
1200 break; /* HHO schemes */
1201
1202 default:
1203 bft_error(__FILE__, __LINE__, 0,
1204 "%s: Invalid space discretization scheme.", __func__);
1205 }
1206
1207 if (fabs(nsp->reference_pressure) > 0 && nsp->n_pressure_ic_defs == 0) {
1208
1209 /* Initialize the initial pressure to the reference pressure */
1210
1211 cs_navsto_add_pressure_ic_by_value(nsp, NULL, &(nsp->reference_pressure));
1212
1213 }
1214
1215 /* Add default post-processing related to the Navier-Stokes system */
1216
1217 cs_post_add_time_mesh_dep_output(cs_navsto_system_extra_post, ns);
1218
1219 if (nsp->post_flag & CS_NAVSTO_POST_STREAM_FUNCTION) {
1220
1221 cs_equation_param_t *eqp = cs_equation_get_param(ns->stream_function_eq);
1222 assert(eqp != NULL);
1223 cs_field_t *w = cs_field_by_name("vorticity");
1224
1225 /* Add a laplacian term: -div.grad */
1226
1227 cs_equation_add_diffusion(eqp, cs_property_by_name("unity"));
1228
1229 /* Add source term as the vorticity w.r.t. the z-axis */
1230
1231 cs_equation_add_source_term_by_dof_func(eqp,
1232 NULL,
1233 cs_flag_primal_cell,
1234 cs_cdofb_navsto_stream_source_term,
1235 (void *)w->val);
1236
1237 } /* Post-processing of the stream function is requested */
1238
1239 /* Finalize the setting of the turbulence modelling */
1240
1241 cs_turbulence_finalize_setup(mesh, connect, quant, time_step,
1242 ns->turbulence);
1243 }
1244
1245 /*----------------------------------------------------------------------------*/
1246 /*!
1247 * \brief Initialize the scheme context structure used to build the algebraic
1248 * system. This is done after the setup step.
1249 * Set an initial value for the velocity and pressure field if needed
1250 *
1251 * \param[in] mesh pointer to a cs_mesh_t structure
1252 * \param[in] connect pointer to a cs_cdo_connect_t structure
1253 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1254 * \param[in] time_step pointer to a cs_time_step_t structure
1255 */
1256 /*----------------------------------------------------------------------------*/
1257
1258 void
cs_navsto_system_initialize(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)1259 cs_navsto_system_initialize(const cs_mesh_t *mesh,
1260 const cs_cdo_connect_t *connect,
1261 const cs_cdo_quantities_t *quant,
1262 const cs_time_step_t *time_step)
1263 {
1264 cs_navsto_system_t *ns = cs_navsto_system;
1265
1266 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
1267
1268 const cs_navsto_param_t *nsp = ns->param;
1269 assert(nsp != NULL);
1270 if (nsp->space_scheme != CS_SPACE_SCHEME_CDOFB)
1271 bft_error(__FILE__, __LINE__, 0,
1272 "%s: Invalid space discretization scheme.", __func__);
1273
1274 /* Allocate then define an array of boundary types for each boundary face */
1275
1276 BFT_MALLOC(ns->bf_type, mesh->n_b_faces, cs_boundary_type_t);
1277 cs_boundary_build_type_array(nsp->boundaries, mesh->n_b_faces, ns->bf_type);
1278
1279 /* Allocate and initialize the scheme context structure */
1280
1281 ns->scheme_context = ns->init_scheme_context(nsp,
1282 ns->adv_field,
1283 ns->mass_flux_array,
1284 ns->mass_flux_array_pre,
1285 ns->bf_type,
1286 ns->coupling_context);
1287
1288 if (time_step->nt_cur < 1) {
1289
1290 /* Initial conditions for the velocity */
1291
1292 if (ns->init_velocity != NULL)
1293 ns->init_velocity(nsp, quant, time_step, ns->scheme_context);
1294
1295 /* Initial conditions for the pressure */
1296
1297 if (ns->init_pressure != NULL)
1298 ns->init_pressure(nsp, quant, time_step, ns->pressure);
1299 }
1300
1301 if (nsp->space_scheme == CS_SPACE_SCHEME_CDOFB) {
1302
1303 if (nsp->coupling == CS_NAVSTO_COUPLING_PROJECTION) {
1304
1305 /* The call to the initialization of the cell pressure should be done
1306 before */
1307
1308 cs_real_t *pr_f = cs_cdofb_predco_get_face_pressure(ns->scheme_context);
1309
1310 cs_cdofb_navsto_init_face_pressure(nsp, connect, time_step, pr_f);
1311
1312 }
1313
1314 /* Initialize the mass flux */
1315
1316 const cs_equation_t *mom_eq = cs_navsto_system_get_momentum_eq();
1317 const cs_real_t *face_vel = cs_equation_get_face_values(mom_eq, false);
1318 cs_cdofb_navsto_mass_flux(nsp, quant, face_vel, ns->mass_flux_array);
1319
1320 } /* Face-based schemes */
1321
1322 cs_turbulence_initialize(mesh, connect, quant, time_step, ns->turbulence);
1323 }
1324
1325 /*----------------------------------------------------------------------------*/
1326 /*!
1327 * \brief Set a solid zone related to the Navier-Stokes equations
1328 *
1329 * \param[in] n_solid_cells number of solid cells
1330 * \param[in] solid_cell_ids list of cell ids (local numbering)
1331 */
1332 /*----------------------------------------------------------------------------*/
1333
1334 void
cs_navsto_system_set_solid_cells(cs_lnum_t n_solid_cells,cs_lnum_t solid_cell_ids[])1335 cs_navsto_system_set_solid_cells(cs_lnum_t n_solid_cells,
1336 cs_lnum_t solid_cell_ids[])
1337 {
1338 cs_navsto_system_t *ns = cs_navsto_system;
1339
1340 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
1341
1342 cs_navsto_param_t *nsp = ns->param;
1343 assert(nsp != NULL);
1344
1345 /* Do not exit the function if n_elts = 0 since there is a parallel
1346 synchronization to count the global number of cells to enforce in the next
1347 call */
1348
1349 /* The momentum equation has to enforce a zero-velocity */
1350
1351 cs_equation_t *mom_eq = cs_navsto_system_get_momentum_eq();
1352 cs_equation_param_t *mom_eqp = cs_equation_get_param(mom_eq);
1353 cs_real_t zero_velocity[3] = {0, 0, 0};
1354
1355 cs_equation_add_cell_enforcement(mom_eqp,
1356 n_solid_cells,
1357 solid_cell_ids,
1358 zero_velocity,
1359 NULL);
1360 }
1361
1362 /*----------------------------------------------------------------------------*/
1363 /*!
1364 * \brief Update variables and related quantities when a new state of the
1365 * Navier-Stokes system has been computed
1366 *
1367 * \param[in] mesh pointer to a cs_mesh_t structure
1368 * \param[in] connect pointer to a cs_cdo_connect_t structure
1369 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1370 * \param[in] time_step structure managing the time stepping
1371 */
1372 /*----------------------------------------------------------------------------*/
1373
1374 void
cs_navsto_system_update(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)1375 cs_navsto_system_update(const cs_mesh_t *mesh,
1376 const cs_cdo_connect_t *connect,
1377 const cs_cdo_quantities_t *quant,
1378 const cs_time_step_t *time_step)
1379 {
1380 CS_UNUSED(mesh);
1381
1382 cs_navsto_system_t *ns = cs_navsto_system;
1383
1384 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
1385
1386 /* Synchronize the cell pressure before doing something */
1387
1388 cs_halo_sync_var(mesh->halo, CS_HALO_STANDARD, ns->pressure->val);
1389
1390 /* Update quantities relative to turbulence (e.g. turbulence
1391 * viscosity...) */
1392
1393 cs_turbulence_t *tbs = ns->turbulence;
1394 if (tbs->update != NULL)
1395 tbs->update(mesh,
1396 connect,
1397 quant,
1398 time_step,
1399 tbs);
1400 }
1401
1402 /*----------------------------------------------------------------------------*/
1403 /*!
1404 * \brief Build, solve and update the Navier-Stokes system in case of a
1405 * steady-state approach
1406 *
1407 * \param[in] mesh pointer to a cs_mesh_t structure
1408 * \param[in] connect pointer to a cs_cdo_connect_t structure
1409 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1410 * \param[in] time_step structure managing the time stepping
1411 */
1412 /*----------------------------------------------------------------------------*/
1413
1414 void
cs_navsto_system_compute_steady_state(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)1415 cs_navsto_system_compute_steady_state(const cs_mesh_t *mesh,
1416 const cs_cdo_connect_t *connect,
1417 const cs_cdo_quantities_t *quant,
1418 const cs_time_step_t *time_step)
1419 {
1420 cs_navsto_system_t *ns = cs_navsto_system;
1421
1422 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
1423
1424 cs_navsto_param_t *nsp = ns->param;
1425 if (!cs_navsto_param_is_steady(nsp))
1426 return;
1427
1428 cs_turbulence_t *tbs = ns->turbulence;
1429
1430 /* First resolve the thermal system if needed */
1431
1432 cs_equation_t *th_eq = cs_thermal_system_get_equation();
1433
1434 if (nsp->model_flag & CS_NAVSTO_MODEL_PASSIVE_THERMAL_TRACER) {
1435
1436 assert(th_eq != NULL);
1437
1438 /* Update variable, properties according to the computed variables */
1439
1440 cs_navsto_system_update(mesh, connect, quant, time_step);
1441
1442 /* Build and solve the Navier-Stokes system */
1443
1444 ns->compute_steady(mesh, nsp, ns->scheme_context);
1445
1446 /* Build and solve the turbulence variable system */
1447
1448 if (tbs->compute_steady != NULL)
1449 tbs->compute_steady(mesh,
1450 connect,
1451 quant,
1452 time_step,
1453 tbs);
1454
1455 /* Solve the thermal equation */
1456
1457 if (cs_equation_param_has_time(cs_equation_get_param(th_eq)) == false)
1458 cs_thermal_system_compute_steady_state(mesh, connect, quant, time_step);
1459
1460 }
1461 else if (cs_flag_test(nsp->model_flag, CS_NAVSTO_MODEL_BOUSSINESQ) &&
1462 cs_flag_test(nsp->model_flag,
1463 CS_NAVSTO_MODEL_WITH_SOLIDIFICATION) == false) {
1464
1465 /* Remark: The "solidification" case is handled in a dedicated module */
1466
1467 assert(th_eq != NULL);
1468 if (cs_equation_param_has_time(cs_equation_get_param(th_eq)))
1469 bft_error(__FILE__, __LINE__, 0,
1470 " %s: Steady-state computation for Navier-Stokes with a"
1471 " Boussinesq approximation\n"
1472 " whereas an unsteady thermal equation is set.\n"
1473 " Please check your settings.", __func__);
1474
1475 cs_real_t *th_var = cs_equation_get_cell_values(th_eq, false);
1476
1477 cs_real_t *th_var_iter_prev = NULL;
1478
1479 BFT_MALLOC(th_var_iter_prev, quant->n_cells, cs_real_t);
1480 memcpy(th_var_iter_prev, th_var, quant->n_cells*sizeof(cs_real_t));
1481
1482 cs_real_t inv_tref = cs_thermal_system_get_reference_temperature();
1483
1484 if (fabs(inv_tref) > cs_math_zero_threshold)
1485 inv_tref = 1./inv_tref;
1486 else
1487 inv_tref = 1.;
1488
1489 cs_real_t delta_th_tolerance = 1 + nsp->delta_thermal_tolerance;
1490 int iter = 0;
1491
1492 do {
1493
1494 /* Update variable, properties according to the computed variables */
1495
1496 cs_navsto_system_update(mesh, connect, quant, time_step);
1497
1498 /* Build and solve the thermal system */
1499
1500 cs_thermal_system_compute_steady_state(mesh, connect, quant, time_step);
1501
1502 /* Build and solve the Navier-Stokes system */
1503
1504 ns->compute_steady(mesh, nsp, ns->scheme_context);
1505
1506 /* Build and solve the turbulence variable system */
1507
1508 if (tbs->compute_steady != NULL)
1509 tbs->compute_steady(mesh,
1510 connect,
1511 quant,
1512 time_step,
1513 tbs);
1514
1515 /* Check convergence */
1516
1517 cs_real_t delta = -1;
1518 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1519 cs_real_t _delta = fabs(th_var[c_id] - th_var_iter_prev[c_id]);
1520 th_var_iter_prev[c_id] = th_var[c_id];
1521 if (_delta > delta)
1522 delta = _delta;
1523 }
1524
1525 /* Update the stopping criteria */
1526
1527 delta_th_tolerance = delta * inv_tref;
1528 iter++;
1529
1530 if (nsp->verbosity > 0)
1531 cs_log_printf(CS_LOG_DEFAULT,
1532 "### Boussinesq.Iteration: %2d | delta_th_var= %5.3e\n",
1533 iter, delta_th_tolerance);
1534
1535 }
1536 while ( (delta_th_tolerance > nsp->delta_thermal_tolerance) &&
1537 (iter < nsp->n_max_outer_iter) );
1538
1539 cs_log_printf(CS_LOG_DEFAULT, " Steady algorithm exits.\n"
1540 " Number of iterations: %d\n"
1541 " Tolerance on the thermal variable: %5.3e\n",
1542 iter, delta_th_tolerance);
1543
1544 BFT_FREE(th_var_iter_prev);
1545
1546 }
1547 else {
1548
1549 /* Build and solve the Navier-Stokes system */
1550
1551 ns->compute_steady(mesh, nsp, ns->scheme_context);
1552
1553 /* Build and solve the turbulence variable system
1554 * TODO, we should loop over NS+turbulence and add a
1555 * stopping criteria on the NS matrix for instance
1556 * */
1557
1558 if (tbs->compute_steady != NULL)
1559 tbs->compute_steady(mesh,
1560 connect,
1561 quant,
1562 time_step,
1563 tbs);
1564 }
1565 }
1566
1567 /*----------------------------------------------------------------------------*/
1568 /*!
1569 * \brief Build, solve and update the Navier-Stokes system
1570 *
1571 * \param[in] mesh pointer to a cs_mesh_t structure
1572 * \param[in] connect pointer to a cs_cdo_connect_t structure
1573 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1574 * \param[in] time_step structure managing the time stepping
1575 */
1576 /*----------------------------------------------------------------------------*/
1577
1578 void
cs_navsto_system_compute(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)1579 cs_navsto_system_compute(const cs_mesh_t *mesh,
1580 const cs_cdo_connect_t *connect,
1581 const cs_cdo_quantities_t *quant,
1582 const cs_time_step_t *time_step)
1583 {
1584 cs_navsto_system_t *ns = cs_navsto_system;
1585
1586 if (ns == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
1587
1588 const cs_navsto_param_t *nsp = ns->param;
1589
1590 cs_turbulence_t *tbs = ns->turbulence;
1591
1592 /* Update variable, properties according to the new computed variables */
1593
1594 cs_navsto_system_update(mesh, connect, quant, time_step);
1595
1596 if (nsp->model_flag & CS_NAVSTO_MODEL_PASSIVE_THERMAL_TRACER) {
1597
1598 /* First resolve the thermal system if needed */
1599
1600 cs_equation_t *th_eq = cs_thermal_system_get_equation();
1601
1602 assert(th_eq != NULL);
1603
1604 /* Build and solve the Navier-Stokes system */
1605
1606 if (cs_navsto_param_is_steady(nsp) == false) {
1607 ns->compute(mesh, nsp, ns->scheme_context);
1608
1609 /* Build and solve the turbulence variable system */
1610
1611 if (tbs->compute != NULL)
1612 tbs->compute(mesh,
1613 connect,
1614 quant,
1615 time_step,
1616 tbs);
1617
1618 }
1619
1620 /* Solve the thermal equation */
1621
1622 if (cs_equation_param_has_time(cs_equation_get_param(th_eq)))
1623 cs_thermal_system_compute(true, mesh, connect, quant, time_step);
1624
1625 }
1626 else if (cs_flag_test(nsp->model_flag, CS_NAVSTO_MODEL_BOUSSINESQ) &&
1627 cs_flag_test(nsp->model_flag,
1628 CS_NAVSTO_MODEL_WITH_SOLIDIFICATION) == false) {
1629
1630 /* Remark: The "solidification" case is handled in a dedicated module */
1631
1632 if (cs_navsto_param_is_steady(nsp))
1633 return;
1634
1635 /* First resolve the thermal system if needed */
1636
1637 cs_equation_t *th_eq = cs_thermal_system_get_equation();
1638
1639 assert(th_eq != NULL);
1640
1641 if (cs_equation_param_has_time(cs_equation_get_param(th_eq))) {
1642
1643 /* Build and solve the thermal system */
1644
1645 cs_thermal_system_compute(true, mesh, connect, quant, time_step);
1646
1647 /* Build and solve the Navier-Stokes system */
1648
1649 ns->compute(mesh, nsp, ns->scheme_context);
1650
1651 /* Build and solve the turbulence variable system */
1652
1653 if (tbs->compute != NULL)
1654 tbs->compute(mesh,
1655 connect,
1656 quant,
1657 time_step,
1658 tbs);
1659
1660
1661 }
1662 else { /* Thermal system is declared as steady. So, this is equivalent to a
1663 standard Navier-Stokes iteration */
1664
1665 /* Build and solve the Navier-Stokes system */
1666
1667 ns->compute(mesh, nsp, ns->scheme_context);
1668
1669 }
1670
1671 }
1672 else {
1673
1674 if (cs_navsto_param_is_steady(nsp))
1675 return;
1676
1677 /* Build and solve the Navier-Stokes system */
1678
1679 ns->compute(mesh, nsp, ns->scheme_context);
1680
1681 /* Build and solve the turbulence variable system */
1682
1683 if (tbs->compute != NULL)
1684 tbs->compute(mesh,
1685 connect,
1686 quant,
1687 time_step,
1688 tbs);
1689
1690 }
1691 }
1692
1693 /*----------------------------------------------------------------------------*/
1694 /*!
1695 * \brief Predefined extra-operations for the Navier-Stokes system
1696 *
1697 * \param[in] mesh pointer to a cs_mesh_t structure
1698 * \param[in] connect pointer to a cs_cdo_connect_t structure
1699 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1700 * \param[in] time_step pointer to a cs_time_step_t structure
1701 */
1702 /*----------------------------------------------------------------------------*/
1703
1704 void
cs_navsto_system_extra_op(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)1705 cs_navsto_system_extra_op(const cs_mesh_t *mesh,
1706 const cs_cdo_connect_t *connect,
1707 const cs_cdo_quantities_t *quant,
1708 const cs_time_step_t *time_step)
1709 {
1710 cs_navsto_system_t *navsto = cs_navsto_system;
1711
1712 if (navsto == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_ns));
1713
1714 const cs_navsto_param_t *nsp = navsto->param;
1715
1716 switch (nsp->space_scheme) {
1717
1718 case CS_SPACE_SCHEME_CDOFB:
1719 {
1720 /* Get the current values not the previous one */
1721
1722 bool need_prev = false;
1723
1724 /* Mass flux is a scalar associated to each face (first interior faces
1725 then border faces */
1726
1727 const cs_real_t *mass_flux = cs_navsto_get_mass_flux(need_prev);
1728
1729 const cs_adv_field_t *adv = cs_navsto_get_adv_field();
1730 const cs_equation_t *eq = cs_navsto_system_get_momentum_eq();
1731 const cs_real_t *u_face = cs_equation_get_face_values(eq, need_prev);
1732 const cs_real_t *u_cell = navsto->velocity->val;
1733 const cs_real_t *p_cell = navsto->pressure->val;
1734
1735 /* Synchronize the cell pressure before doing something */
1736
1737 if (nsp->post_flag & CS_NAVSTO_POST_PRESSURE_GRADIENT)
1738 cs_halo_sync_var(mesh->halo, CS_HALO_STANDARD, p_cell);
1739
1740 cs_cdofb_navsto_extra_op(nsp, mesh, quant, connect, time_step,
1741 navsto->plot_writer,
1742 adv, mass_flux, p_cell, u_cell, u_face);
1743 }
1744 break;
1745
1746 default:
1747 bft_error(__FILE__, __LINE__, 0,
1748 "%s: Invalid space discretization scheme.", __func__);
1749 break;
1750
1751 } /* End of switch */
1752 }
1753
1754 /*----------------------------------------------------------------------------*/
1755 /*!
1756 * \brief Predefined post-processing output for the Navier-Stokes system.
1757 * The prototype of this function is fixed since it is a function
1758 * pointer defined in cs_post.h (\ref cs_post_time_mesh_dep_output_t)
1759 *
1760 * \param[in, out] input pointer to a optional structure (here a
1761 * cs_navsto_system_t structure)
1762 * \param[in] mesh_id id of the output mesh for the current call
1763 * \param[in] cat_id category id of the output mesh for this call
1764 * \param[in] ent_flag indicate global presence of cells (ent_flag[0]),
1765 * interior faces (ent_flag[1]), boundary faces
1766 * (ent_flag[2]), particles (ent_flag[3]) or probes
1767 * (ent_flag[4])
1768 * \param[in] n_cells local number of cells of post_mesh
1769 * \param[in] n_i_faces local number of interior faces of post_mesh
1770 * \param[in] n_b_faces local number of boundary faces of post_mesh
1771 * \param[in] cell_ids list of cells (0 to n-1)
1772 * \param[in] i_face_ids list of interior faces (0 to n-1)
1773 * \param[in] b_face_ids list of boundary faces (0 to n-1)
1774 * \param[in] time_step pointer to a cs_time_step_t struct.
1775 */
1776 /*----------------------------------------------------------------------------*/
1777
1778 void
cs_navsto_system_extra_post(void * input,int mesh_id,int cat_id,int ent_flag[5],cs_lnum_t n_cells,cs_lnum_t n_i_faces,cs_lnum_t n_b_faces,const cs_lnum_t cell_ids[],const cs_lnum_t i_face_ids[],const cs_lnum_t b_face_ids[],const cs_time_step_t * time_step)1779 cs_navsto_system_extra_post(void *input,
1780 int mesh_id,
1781 int cat_id,
1782 int ent_flag[5],
1783 cs_lnum_t n_cells,
1784 cs_lnum_t n_i_faces,
1785 cs_lnum_t n_b_faces,
1786 const cs_lnum_t cell_ids[],
1787 const cs_lnum_t i_face_ids[],
1788 const cs_lnum_t b_face_ids[],
1789 const cs_time_step_t *time_step)
1790 {
1791 CS_UNUSED(n_cells);
1792 CS_UNUSED(n_i_faces);
1793 CS_UNUSED(n_b_faces);
1794 CS_UNUSED(cell_ids);
1795 CS_UNUSED(i_face_ids);
1796 CS_UNUSED(b_face_ids);
1797
1798 cs_navsto_system_t *ns = (cs_navsto_system_t *)input;
1799 if (ns == NULL)
1800 return;
1801
1802 const cs_navsto_param_t *nsp = ns->param;
1803
1804 /* Post-processing of the boundary mass flux */
1805 if (cat_id == CS_POST_MESH_BOUNDARY && ent_flag[2] > 0) {
1806
1807 switch (nsp->space_scheme) {
1808
1809 case CS_SPACE_SCHEME_CDOFB:
1810 case CS_SPACE_SCHEME_HHO_P0:
1811 {
1812 /* Get the current values not the previous one */
1813
1814 bool need_prev = false;
1815
1816 /* In case of postprocessing of the border faces, one has to check if
1817 there is a mesh modification. In particular, a removal of 2D
1818 extruded border faces*/
1819
1820 bool use_parent =
1821 (cs_glob_mesh->n_g_b_faces_all > cs_glob_mesh->n_g_b_faces) ?
1822 false : true;
1823
1824 /* Mass flux is a scalar associated to each face (first interior faces
1825 then border faces */
1826 const cs_real_t *mass_flux = cs_navsto_get_mass_flux(need_prev);
1827 const cs_real_t *mass_bflux = mass_flux + cs_glob_mesh->n_i_faces;
1828
1829 cs_post_write_var(mesh_id,
1830 CS_POST_WRITER_DEFAULT,
1831 "boundary_mass_flux",
1832 1,
1833 false, /* interlace */
1834 use_parent, /* true = original mesh */
1835 CS_POST_TYPE_cs_real_t,
1836 NULL, /* values on cells */
1837 NULL, /* values at internal faces */
1838 mass_bflux, /* values at border faces */
1839 time_step); /* time step structure */
1840
1841 }
1842 break;
1843
1844 default:
1845 bft_error(__FILE__, __LINE__, 0,
1846 "%s: Invalid space scheme\n", __func__);
1847 break;
1848
1849 }
1850
1851 } /* Boundary mesh type */
1852
1853 if (cat_id == CS_POST_MESH_VOLUME && ent_flag[0] > 0) {
1854 /* Velocity-Pressure coupling algorithm */
1855 switch (nsp->coupling) {
1856
1857 case CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY:
1858 case CS_NAVSTO_COUPLING_MONOLITHIC:
1859 /* Nothing to do up to now */
1860 break;
1861
1862 case CS_NAVSTO_COUPLING_PROJECTION:
1863 {
1864 cs_navsto_projection_t *cc
1865 = (cs_navsto_projection_t *)ns->coupling_context;
1866
1867 const cs_field_t *velp = cc->predicted_velocity;
1868
1869 /* Post-process the predicted velocity */
1870 cs_post_write_var(mesh_id,
1871 CS_POST_WRITER_DEFAULT,
1872 velp->name,
1873 3,
1874 true, // interlace
1875 true, // true = original mesh
1876 CS_POST_TYPE_cs_real_t,
1877 velp->val, // values on cells
1878 NULL, // values at internal faces
1879 NULL, // values at border faces
1880 time_step); // time step management struct.
1881
1882 /* Post-process the source term of the correction equation on the
1883 pressure increment (-div(velp_f) */
1884 cs_post_write_var(mesh_id,
1885 CS_POST_WRITER_DEFAULT,
1886 "-DivVelPred",
1887 1,
1888 true, // interlace
1889 true, // true = original mesh
1890 CS_POST_TYPE_cs_real_t,
1891 cc->div_st, // values on cells
1892 NULL, // values at internal faces
1893 NULL, // values at border faces
1894 time_step); // time step management struct.
1895 }
1896 break;
1897
1898 default:
1899 bft_error(__FILE__, __LINE__, 0, _err_invalid_coupling, __func__);
1900 break;
1901 }
1902
1903 } /* cat_id == CS_POST_MESH_VOLUME */
1904
1905 }
1906
1907 /*----------------------------------------------------------------------------*/
1908 /*!
1909 * \brief Summary of the main cs_navsto_system_t structure
1910 */
1911 /*----------------------------------------------------------------------------*/
1912
1913 void
cs_navsto_system_log_setup(void)1914 cs_navsto_system_log_setup(void)
1915 {
1916 cs_navsto_system_t *ns = cs_navsto_system;
1917
1918 if (ns == NULL)
1919 return;
1920
1921 cs_log_printf(CS_LOG_SETUP, "\n");
1922 cs_log_printf(CS_LOG_SETUP, "%s", cs_sep_h1);
1923 cs_log_printf(CS_LOG_SETUP, "\tSummary of the Navier-Stokes system\n");
1924 cs_log_printf(CS_LOG_SETUP, "%s", cs_sep_h1);
1925
1926 /* Main set of numerical parameters */
1927 cs_navsto_param_log(ns->param);
1928
1929 }
1930
1931 /*----------------------------------------------------------------------------*/
1932
1933 END_C_DECLS
1934