1 /*============================================================================
2 * Build an algebraic CDO face-based system for the Navier-Stokes equations
3 * and solved it with an artificial compressibility algorithm
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <math.h>
37 #include <float.h>
38 #include <assert.h>
39 #include <string.h>
40
41 #if defined(HAVE_OPENMP)
42 #include <omp.h>
43 #endif
44
45 /*----------------------------------------------------------------------------
46 * Local headers
47 *----------------------------------------------------------------------------*/
48
49 #include <bft_mem.h>
50
51 #include "cs_blas.h"
52 #include "cs_cdo_bc.h"
53 #include "cs_cdo_blas.h"
54 #include "cs_cdofb_priv.h"
55 #include "cs_cdofb_scaleq.h"
56 #include "cs_cdofb_vecteq.h"
57 #include "cs_cdofb_navsto.h"
58 #if defined(DEBUG) && !defined(NDEBUG)
59 #include "cs_dbg.h"
60 #endif
61 #include "cs_equation_bc.h"
62 #include "cs_equation_common.h"
63 #include "cs_equation_priv.h"
64 #include "cs_evaluate.h"
65 #include "cs_iter_algo.h"
66 #include "cs_log.h"
67 #include "cs_math.h"
68 #include "cs_navsto_sles.h"
69 #include "cs_parall.h"
70 #include "cs_post.h"
71 #include "cs_sles.h"
72 #include "cs_source_term.h"
73 #include "cs_static_condensation.h"
74 #include "cs_timer.h"
75
76 #if defined(HAVE_PETSC)
77 #include "cs_sles_petsc.h"
78 #endif
79
80 /*----------------------------------------------------------------------------
81 * Header for the current file
82 *----------------------------------------------------------------------------*/
83
84 #include "cs_cdofb_ac.h"
85
86 /*----------------------------------------------------------------------------*/
87
88 BEGIN_C_DECLS
89
90 /*=============================================================================
91 * Additional doxygen documentation
92 *============================================================================*/
93
94 /*!
95 * \file cs_cdofb_ac.c
96 *
97 * \brief Build an algebraic CDO face-based system for the Navier-Stokes
98 * equations and solved it with an artificial compressibility algorithm
99 */
100
101 /*=============================================================================
102 * Local structure definitions
103 *============================================================================*/
104
105 /*! \struct cs_cdofb_ac_t
106 * \brief Context related to CDO face-based discretization when dealing with
107 * the Navier-Stokes equations with an Artificial compressibility
108 * algorithm
109 */
110
111 typedef struct {
112
113 /*! \var coupling_context
114 * Pointer to a \ref cs_navsto_ac_t (owned by \ref cs_navsto_system_t)
115 * containing the settings related to an artificial compressibility (AC)
116 * algorithm
117 */
118
119 cs_navsto_ac_t *coupling_context;
120
121 /*!
122 * @name Main field variables
123 * Fields for every main variable of the equation. Got from cs_navsto_system_t
124 */
125
126 /*! \var velocity
127 * Pointer to \ref cs_field_t (owned by \ref cs_navsto_system_t) containing
128 * the cell DoFs of the velocity
129 */
130
131 cs_field_t *velocity;
132
133 /*! \var pressure
134 * Pointer to \ref cs_field_t (owned by \ref cs_navsto_system_t) containing
135 * the cell DoFs of the pressure
136 */
137
138 cs_field_t *pressure;
139
140 /*! \var divergence
141 * Pointer to \ref cs_real_t containing the values of the divergence on the
142 * cells
143 */
144
145 cs_field_t *divergence;
146
147 /*!
148 * @}
149 * @name Parameters of the algorithm
150 * Easy access to useful features and parameters of the algorithm
151 * @{
152 */
153
154 /*! \var is_zeta_uniform
155 * Bool telling if the auxiliary parameter zeta is uniform. Not always
156 * necessary: zeta is typically used in Artificial Compressibility algos
157 */
158
159 bool is_zeta_uniform;
160
161 /*!
162 * @}
163 * @name Advection quantities
164 * Members related to the advection
165 * @{
166 *
167 * \var adv_field
168 * Pointer to the cs_adv_field_t related to the Navier-Stokes eqs (Shared)
169 */
170
171 cs_adv_field_t *adv_field;
172
173 /*! \var mass_flux_array
174 * Current values of the mass flux at primal faces (Shared)
175 */
176
177 cs_real_t *mass_flux_array;
178
179 /*! \var mass_flux_array_pre
180 * Previous values of the mass flux at primal faces (Shared)
181 */
182
183 cs_real_t *mass_flux_array_pre;
184
185 /*!
186 * @}
187 * @name Boundary conditions (BC) management
188 * Functions and elements used for enforcing the BCs
189 * @{
190 *
191 * \var bf_type
192 * Array of boundary type for each boundary face. (Shared)
193 */
194
195 const cs_boundary_type_t *bf_type;
196
197 /*!
198 * \var pressure_bc
199 * Structure storing the metadata after processing the user-defined boundary
200 * conditions related to the pressure field
201 */
202
203 cs_cdo_bc_face_t *pressure_bc;
204
205 /*! \var apply_fixed_wall
206 * \ref cs_cdo_apply_boundary_t function pointer defining how to apply a
207 * wall boundary (no slip boundary)
208 *
209 * \var apply_sliding_wall
210 * \ref cs_cdo_apply_boundary_t function pointer defining how to apply a
211 * wall boundary (a tangential velocity is specified at the wall)
212 *
213 * \var apply_velocity_inlet
214 * \ref cs_cdo_apply_boundary_t function pointer defining how to apply a
215 * boundary with a fixed velocity at the inlet
216 *
217 * \var apply_symmetry
218 * \ref cs_cdo_apply_boundary_t function pointer defining how to apply a
219 * symmetry boundary
220 */
221
222 cs_cdo_apply_boundary_t *apply_fixed_wall;
223 cs_cdo_apply_boundary_t *apply_sliding_wall;
224 cs_cdo_apply_boundary_t *apply_velocity_inlet;
225 cs_cdo_apply_boundary_t *apply_symmetry;
226
227 /*!
228 * \var add_gravity_term
229 * Compute and add the source term related to the gravity vector
230 * This can be the Boussinesq term or the hydrostatic term (rho*g)
231 */
232
233 cs_cdofb_navsto_source_t *add_gravity_term;
234
235 /*!
236 * @}
237 * @name Convergence monitoring
238 * Structure used to drive the convergence of high-level iterative algorithms
239 * @{
240 */
241
242 /*!
243 * \var nl_algo
244 * Structure driving the convergence of the non-linear algorithm
245 */
246
247 cs_iter_algo_t *nl_algo;
248
249 /*!
250 * @}
251 * @name Performance monitoring
252 * Monitoring the efficiency of the algorithm used to solve the Navier-Stokes
253 * system
254 * @{
255 */
256
257 /*! \var timer
258 * Cumulated elapsed time for building and solving the Navier--Stokes system
259 */
260
261 cs_timer_counter_t timer;
262
263 /*! @} */
264
265 } cs_cdofb_ac_t;
266
267 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
268
269 /*=============================================================================
270 * Local Macro definitions and structure definitions
271 *============================================================================*/
272
273 #define CS_CDOFB_AC_DBG 0
274
275 /*============================================================================
276 * Private variables
277 *============================================================================*/
278
279 /* Pointer to shared structures */
280 static const cs_cdo_quantities_t *cs_shared_quant;
281 static const cs_cdo_connect_t *cs_shared_connect;
282 static const cs_time_step_t *cs_shared_time_step;
283 static const cs_matrix_structure_t *cs_shared_ms;
284
285 /*============================================================================
286 * Private function prototypes
287 *============================================================================*/
288
289 /*----------------------------------------------------------------------------*/
290 /*!
291 * \brief Copy current content of the Navier-Stokes related variables-fields
292 * to previous values
293 * Case of the artificial compressibility coupling algorithm.
294 *
295 * \param[in, out] sc scheme context (\ref cs_cdofb_ac_t struct.)
296 * \param[in, out] cc coupling context (\ref cs_navsto_ac_t struct.)
297 */
298 /*----------------------------------------------------------------------------*/
299
300 static inline void
_ac_fields_to_previous(cs_cdofb_ac_t * sc,cs_navsto_ac_t * cc)301 _ac_fields_to_previous(cs_cdofb_ac_t *sc,
302 cs_navsto_ac_t *cc)
303 {
304 const cs_cdo_quantities_t *quant = cs_shared_quant;
305
306 cs_field_current_to_previous(sc->velocity);
307 cs_field_current_to_previous(sc->pressure);
308 cs_field_current_to_previous(sc->divergence);
309
310 /* Face velocity arrays */
311 cs_cdofb_vecteq_t *eqc = (cs_cdofb_vecteq_t *)cc->momentum->scheme_context;
312
313 if (eqc->face_values_pre != NULL)
314 memcpy(eqc->face_values_pre, eqc->face_values,
315 3 * quant->n_faces * sizeof(cs_real_t));
316
317 /* Mass flux arrays */
318 memcpy(sc->mass_flux_array_pre, sc->mass_flux_array,
319 quant->n_faces * sizeof(cs_real_t));
320 }
321
322 /*----------------------------------------------------------------------------*/
323 /*!
324 * \brief Compute the divergence of the velocity in each cell
325 *
326 * \param[in] vel_f velocity DoFs on faces
327 * \param[in, out] div divergence of the given velocity in each cell
328 */
329 /*----------------------------------------------------------------------------*/
330
331 static void
_ac_compute_div(const cs_real_t vel_f[],cs_real_t div[])332 _ac_compute_div(const cs_real_t vel_f[],
333 cs_real_t div[])
334 {
335 const cs_cdo_quantities_t *quant = cs_shared_quant;
336
337 /* Update the divergence of the velocity field */
338
339 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
340 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++)
341 div[c_id] = cs_cdofb_navsto_cell_divergence(c_id,
342 quant,
343 cs_shared_connect->c2f,
344 vel_f);
345
346 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_AC_DBG > 2
347 cs_dbg_darray_to_listing("VELOCITY_DIV", quant->n_cells, div, 9);
348 #endif
349 }
350
351 /*----------------------------------------------------------------------------*/
352 /*!
353 * \brief Performs the update of the pressure after the resolution of the
354 * momentum equation
355 *
356 * \param[in] t_eval time at which the property/functions are evaluated
357 * \param[in] dt_cur current value of the time step
358 * \param[in] zeta scaling of the div operator
359 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
360 * \param[in] eqb pointer to a \ref cs_equation_builder_t structure
361 * \param[in] div_c mean velocity divergence in each cell
362 * \param[in, out] pr_c pressure DoFs (on cells)
363 */
364 /*----------------------------------------------------------------------------*/
365
366 static inline void
_ac_update_pr(const cs_real_t t_eval,const cs_real_t dt_cur,const cs_property_t * zeta,const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_real_t div_c[],cs_real_t pr_c[])367 _ac_update_pr(const cs_real_t t_eval,
368 const cs_real_t dt_cur,
369 const cs_property_t *zeta,
370 const cs_equation_param_t *eqp,
371 const cs_equation_builder_t *eqb,
372 const cs_real_t div_c[],
373 cs_real_t pr_c[])
374 {
375 const cs_cdo_quantities_t *quant = cs_shared_quant;
376
377 if (cs_property_is_uniform(zeta)) {
378
379 const cs_real_t o_zeta_c = 1./cs_property_get_cell_value(0, t_eval, zeta);
380
381 if (eqb->time_pty_uniform) {
382
383 /* Best case scenario: everything is uniform */
384 const cs_real_t t_pty = cs_property_get_cell_value(0, t_eval,
385 eqp->time_property);
386 const cs_real_t coef_mult = t_pty * dt_cur * o_zeta_c;
387
388 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
389 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++)
390 pr_c[c_id] -= coef_mult * div_c[c_id];
391
392 }
393 else { /* time_pty non uniform */
394
395 const cs_real_t coef_mult = dt_cur * o_zeta_c;
396
397 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
398 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
399 cs_real_t t_pty = cs_property_get_cell_value(c_id, t_eval,
400 eqp->time_property);
401 pr_c[c_id] -= coef_mult * t_pty * div_c[c_id];
402 }
403
404 }
405
406 }
407 else { /* zeta not uniform */
408
409 if (eqb->time_pty_uniform) {
410
411 const cs_real_t coef_mult =
412 dt_cur * cs_property_get_cell_value(0, t_eval, eqp->time_property);
413
414 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
415 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
416 cs_real_t o_zeta_c = 1./cs_property_get_cell_value(c_id, t_eval, zeta);
417 pr_c[c_id] -= coef_mult * o_zeta_c * div_c[c_id];
418 }
419
420 }
421 else { /* Neither zeta nor time_pty is uniform */
422
423 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
424 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
425 cs_real_t o_zeta_c = 1./cs_property_get_cell_value(c_id, t_eval, zeta);
426 cs_real_t t_pty = cs_property_get_cell_value(c_id, t_eval,
427 eqp->time_property);
428 pr_c[c_id] -= t_pty * dt_cur * o_zeta_c * div_c[c_id];
429 }
430
431 } /* time_pty non uniform */
432
433 } /* zeta non uniform */
434 }
435
436 /*----------------------------------------------------------------------------*/
437 /*!
438 * \brief Apply the boundary conditions to the local system when this should
439 * be done before the static condensation
440 * Case of CDO Fb schemes with AC coupling
441 *
442 * \param[in] sc pointer to a cs_cdofb_ac_t structure
443 * \param[in] eqp pointer to a cs_equation_param_t structure
444 * \param[in] cm pointer to a cellwise view of the mesh
445 * \param[in] bf_type type of boundary for the boundary faces
446 * \param[in] diff_pty pointer to \ref cs_property_data_t for diffusion
447 * \param[in, out] csys pointer to a cellwise view of the system
448 * \param[in, out] cb pointer to a cellwise builder
449 */
450 /*----------------------------------------------------------------------------*/
451
452 static void
_ac_apply_bc_partly(const cs_cdofb_ac_t * sc,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_boundary_type_t bf_type[],const cs_property_data_t * diff_pty,cs_cell_sys_t * csys,cs_cell_builder_t * cb)453 _ac_apply_bc_partly(const cs_cdofb_ac_t *sc,
454 const cs_equation_param_t *eqp,
455 const cs_cell_mesh_t *cm,
456 const cs_boundary_type_t bf_type[],
457 const cs_property_data_t *diff_pty,
458 cs_cell_sys_t *csys,
459 cs_cell_builder_t *cb)
460 {
461 assert(cs_equation_param_has_diffusion(eqp) == true);
462
463 if (cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
464
465 /* Update the velocity-block and the right-hand side (part related to the
466 * momentum equation). */
467
468 /* Neumann boundary conditions:
469 * The common practice is to define Phi_neu = - lambda * grad(u) . n_fc
470 * An outward flux is a positive flux whereas an inward flux is negative
471 * The minus just above implies the minus just below */
472
473 if (csys->has_nhmg_neumann)
474 for (short int i = 0; i < 3*cm->n_fc; i++)
475 csys->rhs[i] -= csys->neu_values[i];
476
477 const cs_real_t pc = sc->pressure->val[cm->c_id];
478
479 for (short int i = 0; i < csys->n_bc_faces; i++) {
480
481 /* Get the boundary face in the cell numbering */
482 const short int f = csys->_f_ids[i];
483 const cs_quant_t pfq = cm->face[f];
484 const cs_real_t f_prs = pfq.meas * pc;
485 cs_real_t *f_rhs = csys->rhs + 3*f;
486
487 if (bf_type[i] & CS_BOUNDARY_IMPOSED_VEL) {
488 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
489 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM) {
490 sc->apply_velocity_inlet(f, eqp, cm, diff_pty, cb, csys);
491 f_rhs[0] -= f_prs * pfq.unitv[0];
492 f_rhs[1] -= f_prs * pfq.unitv[1];
493 f_rhs[2] -= f_prs * pfq.unitv[2];
494 }
495 }
496
497 else if (bf_type[i] & CS_BOUNDARY_WALL) {
498 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
499 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM) {
500 if (bf_type[i] & CS_BOUNDARY_SLIDING_WALL)
501 sc->apply_sliding_wall(f, eqp, cm, diff_pty, cb, csys);
502 else
503 sc->apply_fixed_wall(f, eqp, cm, diff_pty, cb, csys);
504 f_rhs[0] -= f_prs * pfq.unitv[0];
505 f_rhs[1] -= f_prs * pfq.unitv[1];
506 f_rhs[2] -= f_prs * pfq.unitv[2];
507 }
508 }
509
510 else if (bf_type[i] & CS_BOUNDARY_SYMMETRY) {
511 /* Always weakly enforce the symmetric constraint on the
512 velocity-block */
513 sc->apply_symmetry(f, eqp, cm, diff_pty, cb, csys);
514 f_rhs[0] -= f_prs * pfq.unitv[0];
515 f_rhs[1] -= f_prs * pfq.unitv[1];
516 f_rhs[2] -= f_prs * pfq.unitv[2];
517 }
518
519 /* default: nothing to do (case of a "natural" outlet) */
520
521 } /* Loop on boundary faces */
522
523 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_AC_DBG > 1
524 if (cs_dbg_cw_test(eqp, cm, csys))
525 cs_cell_sys_dump(">> Cell system matrix after partial BC enforcement",
526 csys);
527 #endif
528 } /* Boundary cell */
529 }
530
531 /*----------------------------------------------------------------------------*/
532 /*!
533 * \brief Apply the boundary conditions to the local system when this should be
534 * done after the static condensation. Apply the internal enforcement
535 * if needed.
536 * Case of CDO-Fb schemes with AC coupling
537 *
538 * \param[in] sc pointer to a cs_cdofb_ac_t structure
539 * \param[in] eqp pointer to a cs_equation_param_t structure
540 * \param[in] eqb pointer to a cs_equation_builder_t structure
541 * \param[in] cm pointer to a cellwise view of the mesh
542 * \param[in] bf_type type of boundary for the boundary faces
543 * \param[in] diff_pty pointer to \ref cs_property_data_t for diffusion
544 * \param[in, out] csys pointer to a cellwise view of the system
545 * \param[in, out] cb pointer to a cellwise builder
546 */
547 /*----------------------------------------------------------------------------*/
548
549 static void
_ac_apply_remaining_bc(const cs_cdofb_ac_t * sc,const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_cell_mesh_t * cm,const cs_boundary_type_t * bf_type,const cs_property_data_t * diff_pty,cs_cell_sys_t * csys,cs_cell_builder_t * cb)550 _ac_apply_remaining_bc(const cs_cdofb_ac_t *sc,
551 const cs_equation_param_t *eqp,
552 const cs_equation_builder_t *eqb,
553 const cs_cell_mesh_t *cm,
554 const cs_boundary_type_t *bf_type,
555 const cs_property_data_t *diff_pty,
556 cs_cell_sys_t *csys,
557 cs_cell_builder_t *cb)
558 {
559 if (cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
560
561 /* Update the divergence operator and the right-hand side related to the
562 * mass equation.
563 * Enforcement of Dirichlet BC in a stronger way if this is the choice
564 */
565
566 for (short int i = 0; i < csys->n_bc_faces; i++) {
567
568 /* Get the boundary face in the cell numbering */
569
570 const short int f = csys->_f_ids[i];
571
572 if (bf_type[i] & CS_BOUNDARY_IMPOSED_VEL) {
573
574 /* Enforcement of the velocity for the velocity-block
575 * Dirichlet on the three components of the velocity field */
576
577 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_PENALIZED ||
578 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_ALGEBRAIC) {
579 sc->apply_velocity_inlet(f, eqp, cm, diff_pty, cb, csys);
580 }
581 }
582
583 else if (bf_type[i] & CS_BOUNDARY_WALL) {
584
585 /* Enforcement of the velocity for the velocity-block
586 * Dirichlet on the three components of the velocity field */
587
588 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_PENALIZED ||
589 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_ALGEBRAIC) {
590 if (bf_type[i] & CS_BOUNDARY_SLIDING_WALL)
591 sc->apply_sliding_wall(f, eqp, cm, diff_pty, cb, csys);
592 else
593 sc->apply_fixed_wall(f, eqp, cm, diff_pty, cb, csys);
594 }
595 }
596
597 #if 0
598 else if (bf_type[i] & CS_BOUNDARY_SYMMETRY) {
599 /* Weak-enforcement for the velocity-block (cf. _ac_apply_bc_partly) */
600 }
601 #endif
602
603 /* default: nothing to do (case of a "natural" outlet) */
604
605 } /* Loop on boundary faces */
606
607 } /* Boundary cell */
608
609 if (cs_equation_param_has_internal_enforcement(eqp)) {
610
611 /* Internal enforcement of DoFs: Update csys (matrix and rhs) */
612
613 cs_equation_enforced_internal_block_dofs(eqb, cb, csys);
614
615 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_MONOLITHIC_DBG > 2
616 if (cs_dbg_cw_test(eqp, cm, csys))
617 cs_cell_sys_dump("\n>> Cell system after the internal enforcement",
618 csys);
619 #endif
620 }
621 }
622
623 /*----------------------------------------------------------------------------*/
624 /*!
625 * \brief Build a linear system for Stokes, Oseen or Navier-Stokes in the
626 * case of an implicit Euler time scheme
627 *
628 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
629 * \param[in] vel_f_pre velocity face DoFs of the previous time step
630 * \param[in] vel_c_pre velocity cell DoFs of the previous time step
631 * \param[in] prs_c_pre pressure cell DoFs of the previous time step
632 * \param[in, out] sc void cast into to a \ref cs_cdofb_ac_t pointer
633 * \param[in, out] matrix pointer to a \ref cs_matrix_t structure
634 * \param[in, out] rhs rhs array related to the momentum eq.
635 */
636 /*----------------------------------------------------------------------------*/
637
638 static void
_implicit_euler_build(const cs_navsto_param_t * nsp,const cs_real_t vel_f_pre[],const cs_real_t vel_c_pre[],const cs_real_t prs_c_pre[],cs_cdofb_ac_t * sc,cs_matrix_t * matrix,cs_real_t * rhs)639 _implicit_euler_build(const cs_navsto_param_t *nsp,
640 const cs_real_t vel_f_pre[],
641 const cs_real_t vel_c_pre[],
642 const cs_real_t prs_c_pre[],
643 cs_cdofb_ac_t *sc,
644 cs_matrix_t *matrix,
645 cs_real_t *rhs)
646 {
647 /* Retrieve high-level structures */
648
649 cs_navsto_ac_t *cc = sc->coupling_context;
650 cs_equation_t *mom_eq = cc->momentum;
651 cs_cdofb_vecteq_t *mom_eqc= (cs_cdofb_vecteq_t *)mom_eq->scheme_context;
652 cs_equation_param_t *mom_eqp = mom_eq->param;
653 cs_equation_builder_t *mom_eqb = mom_eq->builder;
654
655 /* Retrieve shared structures */
656
657 const cs_time_step_t *ts = cs_shared_time_step;
658 const cs_cdo_quantities_t *quant = cs_shared_quant;
659 const cs_cdo_connect_t *connect = cs_shared_connect;
660 const cs_range_set_t *rs = connect->range_sets[CS_CDO_CONNECT_FACE_VP0];
661
662 assert(cs_equation_param_has_time(mom_eqp) == true);
663 assert(mom_eqp->time_scheme == CS_TIME_SCHEME_EULER_IMPLICIT);
664 assert(matrix != NULL && rhs != NULL);
665
666 #if defined(DEBUG) && !defined(NDEBUG)
667 if (quant->n_b_faces > 0)
668 assert(mom_eqb->dir_values != NULL);
669 #endif
670
671 /* Initialize the structure to assemble values */
672
673 cs_matrix_assembler_values_t *mav =
674 cs_matrix_assembler_values_init(matrix, NULL, NULL);
675
676 # pragma omp parallel if (quant->n_cells > CS_THR_MIN)
677 {
678 #if defined(HAVE_OPENMP) /* Determine the default number of OpenMP threads */
679 int t_id = omp_get_thread_num();
680 #else
681 int t_id = 0;
682 #endif
683
684 /* Each thread get back its related structures:
685 Get the cell-wise view of the mesh and the algebraic system */
686
687 cs_cdofb_navsto_builder_t nsb = cs_cdofb_navsto_create_builder(nsp,
688 connect);
689 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
690 cs_equation_assemble_t *eqa = cs_equation_assemble_get(t_id);
691 cs_hodge_t *diff_hodge =
692 (mom_eqc->diffusion_hodge == NULL) ? NULL:mom_eqc->diffusion_hodge[t_id];
693 cs_hodge_t *mass_hodge =
694 (mom_eqc->mass_hodge == NULL) ? NULL:mom_eqc->mass_hodge[t_id];
695
696 cs_cell_sys_t *csys = NULL;
697 cs_cell_builder_t *cb = NULL;
698
699 cs_cdofb_vecteq_get(&csys, &cb);
700
701 /* Set times at which one evaluates quantities when needed */
702
703 const cs_real_t t_cur = ts->t_cur;
704 const cs_real_t dt_cur = ts->dt[0];
705 const cs_real_t inv_dtcur = 1./dt_cur;
706
707 cb->t_pty_eval = t_cur + dt_cur;
708 cb->t_bc_eval = t_cur + dt_cur;
709 cb->t_st_eval = t_cur + dt_cur;
710
711 /* Initialization of the values of properties */
712
713 cs_equation_init_properties(mom_eqp, mom_eqb, diff_hodge, cb);
714
715 cs_real_t o_zeta_c = 1./cs_property_get_cell_value(0, cb->t_pty_eval,
716 cc->zeta);
717
718 /* ---------------------------------------------
719 * Main loop on cells to build the linear system
720 * --------------------------------------------- */
721
722 # pragma omp for CS_CDO_OMP_SCHEDULE
723 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
724
725 /* Set the current cell flag */
726
727 cb->cell_flag = connect->cell_flag[c_id];
728
729 /* Set the local mesh structure for the current cell */
730
731 cs_cell_mesh_build(c_id,
732 cs_equation_cell_mesh_flag(cb->cell_flag, mom_eqb),
733 connect, quant, cm);
734
735 /* Starts from the stationary Stokes problem where
736 *
737 * | | |
738 * | A | Bt | B is the divergence (Bt the gradient)
739 * | | | A is csys->mat in what follows
740 * |--------|---------| The viscous part arising from the CDO-Fb
741 * | | | schemes for vector-valued variables
742 * | B | 0 |
743 * | | |
744 */
745
746 /* Set the local (i.e. cellwise) structures for the current cell */
747
748 cs_cdofb_vecteq_init_cell_system(cm, mom_eqp, mom_eqb,
749 vel_f_pre, vel_c_pre,
750 NULL, NULL, /* no n-1 state is given */
751 csys, cb);
752
753 /* 1- SETUP THE NAVSTO LOCAL BUILDER *
754 * ================================= *
755 * - Set the type of boundary
756 * - Set the pressure boundary conditions (if required)
757 * - Define the divergence operator used in the linear system (div_op is
758 * equal to minus the divergence)
759 */
760
761 cs_cdofb_navsto_define_builder(cb->t_bc_eval, nsp, cm, csys,
762 sc->pressure_bc, sc->bf_type, &nsb);
763
764 /* 2- VELOCITY (VECTORIAL) EQUATION */
765 /* ================================ */
766
767 cs_cdofb_vecteq_conv_diff_reac(mom_eqp, mom_eqb, mom_eqc, cm,
768 mass_hodge, diff_hodge, csys, cb);
769
770 /* Update the property related to grad-div */
771
772 if ( !(sc->is_zeta_uniform) )
773 o_zeta_c = 1./cs_property_value_in_cell(cm, cc->zeta, cb->t_pty_eval);
774
775 /* Update the property related to the unsteady term */
776
777 if (!(mom_eqb->time_pty_uniform))
778 cb->tpty_val = cs_property_value_in_cell(cm, mom_eqp->time_property,
779 cb->t_pty_eval);
780
781 const short int n_fc = cm->n_fc;
782
783 cs_cdofb_navsto_add_grad_div(n_fc,
784 cb->tpty_val * dt_cur * o_zeta_c / cm->vol_c,
785 nsb.div_op, csys->mat);
786
787 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_AC_DBG > 1
788 if (cs_dbg_cw_test(mom_eqp, cm, csys))
789 cs_cell_sys_dump(">> Local system after conv/diff/reac and grad-div"
790 " terms", csys);
791 #endif
792
793 /* 3- SOURCE TERM COMPUTATION (for the momentum equation) */
794 /* ====================================================== */
795
796 const bool has_sourceterm = cs_equation_param_has_sourceterm(mom_eqp);
797 if (has_sourceterm)
798 cs_cdofb_vecteq_sourceterm(cm, mom_eqp, cb->t_st_eval,
799 1., /* time scaling */
800 mass_hodge,
801 cb, mom_eqb, csys);
802
803 /* 3b- OTHER RHS CONTRIBUTIONS
804 * ===========================
805 * Apply the operator gradient to the pressure field and add it to the
806 * rhs */
807
808 cs_sdm_add_scalvect(3*n_fc, -prs_c_pre[c_id], nsb.div_op, csys->rhs);
809
810 /* Gravity effects and/or Boussinesq approximation rely on another
811 strategy than classical source term. The treatment is more compatible
812 with the pressure gradient by doing so. */
813
814 if (sc->add_gravity_term != NULL)
815 sc->add_gravity_term(nsp, cm, &nsb, csys);
816
817 /* First part of the BOUNDARY CONDITIONS
818 * ===================
819 * Apply a part of BC before the time scheme */
820
821 _ac_apply_bc_partly(sc, mom_eqp, cm, nsb.bf_type,
822 diff_hodge->pty_data, csys, cb);
823
824 /* 4- TIME CONTRIBUTION */
825 /* ==================== */
826
827 if (mom_eqb->sys_flag & CS_FLAG_SYS_TIME_DIAG) {
828
829 /* Mass lumping or Hodge-Voronoi */
830
831 const double ptyc = cb->tpty_val * cm->vol_c * inv_dtcur;
832
833 /* Get cell-cell block */
834
835 cs_sdm_t *acc = cs_sdm_get_block(csys->mat, n_fc, n_fc);
836
837 for (short int k = 0; k < 3; k++) {
838 csys->rhs[3*n_fc + k] += ptyc * csys->val_n[3*n_fc+k];
839 /* Simply add an entry in mat[cell, cell] */
840 acc->val[4*k] += ptyc;
841 }
842
843 }
844 else
845 bft_error(__FILE__, __LINE__, 0,
846 "Only diagonal time treatment available so far.");
847
848 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_AC_DBG > 1
849 if (cs_dbg_cw_test(mom_eqp, cm, csys))
850 cs_cell_sys_dump(">> Local system matrix before condensation", csys);
851 #endif
852
853 /* 5- STATIC CONDENSATION
854 * ======================
855 * Static condensation of the local system matrix of size n_fc + 1 into
856 * a matrix of size n_fc.
857 * Store data in rc_tilda and acf_tilda to compute the values at cell
858 * centers after solving the system */
859
860 cs_static_condensation_vector_eq(connect->c2f,
861 mom_eqc->rc_tilda, mom_eqc->acf_tilda,
862 cb, csys);
863
864 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_AC_DBG > 1
865 if (cs_dbg_cw_test(mom_eqp, cm, csys))
866 cs_cell_sys_dump(">> Local system matrix after static condensation",
867 csys);
868 #endif
869
870 /* 6- Remaining part of BOUNDARY CONDITIONS
871 * ======================================== */
872
873 _ac_apply_remaining_bc(sc, mom_eqp, mom_eqb, cm, nsb.bf_type,
874 diff_hodge->pty_data, csys, cb);
875
876 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_AC_DBG > 0
877 if (cs_dbg_cw_test(mom_eqp, cm, csys))
878 cs_cell_sys_dump(">> (FINAL) Local system matrix", csys);
879 #endif
880
881 /* ASSEMBLY PROCESS */
882 /* ================ */
883
884 cs_cdofb_vecteq_assembly(csys, rs, cm, has_sourceterm,
885 mom_eqc, eqa, mav, rhs);
886
887 } /* Main loop on cells */
888
889 /* Free temporary buffer */
890
891 cs_cdofb_navsto_free_builder(&nsb);
892
893 } /* End of the OpenMP Block */
894
895 cs_matrix_assembler_values_done(mav); /* optional */
896 cs_matrix_assembler_values_finalize(&mav);
897 }
898
899 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
900
901 /*============================================================================
902 * Public function prototypes
903 *============================================================================*/
904
905 /*----------------------------------------------------------------------------*/
906 /*!
907 * \brief Set shared pointers from the main domain members
908 *
909 * \param[in] quant additional mesh quantities struct.
910 * \param[in] connect pointer to a \ref cs_cdo_connect_t struct.
911 * \param[in] time_step pointer to a \ref cs_time_step_t structure
912 */
913 /*----------------------------------------------------------------------------*/
914
915 void
cs_cdofb_ac_init_common(const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect,const cs_time_step_t * time_step)916 cs_cdofb_ac_init_common(const cs_cdo_quantities_t *quant,
917 const cs_cdo_connect_t *connect,
918 const cs_time_step_t *time_step)
919 {
920 /* Assign static const pointers */
921
922 cs_shared_quant = quant;
923 cs_shared_connect = connect;
924 cs_shared_time_step = time_step;
925
926 /* Matrix structure related to the algebraic system for vector-valued eq. */
927
928 cs_shared_ms = cs_cdofb_vecteq_matrix_structure();
929 }
930
931 /*----------------------------------------------------------------------------*/
932 /*!
933 * \brief Initialize a \ref cs_cdofb_ac_t structure
934 *
935 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
936 * \param[in] adv_field pointer to \ref cs_adv_field_t structure
937 * \param[in] mflux current values of the mass flux across primal faces
938 * \param[in] mflux_pre current values of the mass flux across primal faces
939 * \param[in] fb_type type of boundary for each boundary face
940 * \param[in] nsc_input pointer to a \ref cs_navsto_ac_t structure
941 *
942 * \return a pointer to a new allocated \ref cs_cdofb_ac_t structure
943 */
944 /*----------------------------------------------------------------------------*/
945
946 void *
cs_cdofb_ac_init_scheme_context(const cs_navsto_param_t * nsp,cs_adv_field_t * adv_field,cs_real_t * mflux,cs_real_t * mflux_pre,cs_boundary_type_t * fb_type,void * nsc_input)947 cs_cdofb_ac_init_scheme_context(const cs_navsto_param_t *nsp,
948 cs_adv_field_t *adv_field,
949 cs_real_t *mflux,
950 cs_real_t *mflux_pre,
951 cs_boundary_type_t *fb_type,
952 void *nsc_input)
953 {
954 assert(nsp != NULL && nsc_input != NULL); /* Sanity checks */
955
956 /* Cast the coupling context (CC) */
957
958 cs_navsto_ac_t *cc = (cs_navsto_ac_t *)nsc_input;
959
960 /* Navier-Stokes scheme context (SC) */
961
962 cs_cdofb_ac_t *sc = NULL;
963
964 if (nsp->space_scheme != CS_SPACE_SCHEME_CDOFB)
965 bft_error(__FILE__, __LINE__, 0, " %s: Invalid space scheme.\n",
966 __func__);
967
968 BFT_MALLOC(sc, 1, cs_cdofb_ac_t);
969
970 /* Quantities shared with the cs_navsto_system_t structure */
971
972 sc->coupling_context = cc;
973 sc->adv_field = adv_field;
974 sc->mass_flux_array = mflux;
975 sc->mass_flux_array_pre = mflux_pre;
976
977 /* Quick access to the main fields */
978
979 sc->velocity = cs_field_by_name("velocity");
980 sc->pressure = cs_field_by_name("pressure");
981
982 if (nsp->post_flag & CS_NAVSTO_POST_VELOCITY_DIVERGENCE)
983 sc->divergence = cs_field_by_name("velocity_divergence");
984 else
985 sc->divergence = NULL;
986
987 /* Parameters related to the AC algorithm */
988
989 sc->is_zeta_uniform = true; /* s_property_is_uniform(cc->zeta); */
990
991 /* Boundary treatment */
992
993 sc->bf_type = fb_type;
994
995 /* Processing of the pressure boundary condition */
996
997 sc->pressure_bc = cs_cdo_bc_face_define(CS_CDO_BC_HMG_NEUMANN, /* Default */
998 true, /* Steady BC up to now */
999 1, /* Dimension */
1000 nsp->n_pressure_bc_defs,
1001 nsp->pressure_bc_defs,
1002 cs_shared_quant->n_b_faces);
1003
1004 cs_equation_param_t *mom_eqp = cc->momentum->param;
1005 cs_equation_builder_t *mom_eqb = cc->momentum->builder;
1006
1007 mom_eqb->bd_msh_flag |= CS_FLAG_COMP_PFC;
1008
1009 /* Set the way to enforce the Dirichlet BC on the velocity
1010 * "fixed_wall" means a no-slip BC */
1011
1012 sc->apply_symmetry = cs_cdofb_symmetry;
1013
1014 switch (mom_eqp->default_enforcement) {
1015
1016 case CS_PARAM_BC_ENFORCE_ALGEBRAIC:
1017 sc->apply_velocity_inlet = cs_cdofb_block_dirichlet_alge;
1018 sc->apply_sliding_wall = cs_cdofb_block_dirichlet_alge;
1019 sc->apply_fixed_wall = cs_cdofb_block_dirichlet_alge;
1020 break;
1021
1022 case CS_PARAM_BC_ENFORCE_PENALIZED:
1023 sc->apply_velocity_inlet = cs_cdofb_block_dirichlet_pena;
1024 sc->apply_sliding_wall = cs_cdofb_block_dirichlet_pena;
1025 sc->apply_fixed_wall = cs_cdofb_block_dirichlet_pena;
1026 break;
1027
1028 case CS_PARAM_BC_ENFORCE_WEAK_NITSCHE:
1029 sc->apply_velocity_inlet = cs_cdofb_block_dirichlet_weak;
1030 sc->apply_sliding_wall = cs_cdofb_block_dirichlet_weak;
1031 sc->apply_fixed_wall = cs_cdofb_block_dirichlet_weak;
1032 break;
1033
1034 case CS_PARAM_BC_ENFORCE_WEAK_SYM:
1035 sc->apply_velocity_inlet = cs_cdofb_block_dirichlet_wsym;
1036 sc->apply_sliding_wall = cs_cdofb_block_dirichlet_wsym;
1037 sc->apply_fixed_wall = cs_cdofb_block_dirichlet_wsym;
1038 break;
1039
1040 default:
1041 bft_error(__FILE__, __LINE__, 0,
1042 " %s: Invalid type of algorithm to enforce Dirichlet BC.",
1043 __func__);
1044
1045 }
1046
1047 /* Take into account the gravity effect if needed */
1048
1049 cs_cdofb_navsto_set_gravity_func(nsp, &(sc->add_gravity_term));
1050
1051 /* Iterative algorithm to handle the non-linearity (Picard by default) */
1052
1053 const cs_navsto_param_sles_t *nslesp = nsp->sles_param;
1054
1055 sc->nl_algo = cs_iter_algo_create(nslesp->nl_algo_param);
1056
1057 if (nslesp->nl_algo_type == CS_PARAM_NL_ALGO_ANDERSON)
1058 sc->nl_algo->context = cs_iter_algo_aa_create(nslesp->anderson_param,
1059 cs_shared_quant->n_faces);
1060
1061 /* Monitoring */
1062
1063 CS_TIMER_COUNTER_INIT(sc->timer);
1064
1065 return sc;
1066 }
1067
1068 /*----------------------------------------------------------------------------*/
1069 /*!
1070 * \brief Destroy a \ref cs_cdofb_ac_t structure
1071 *
1072 * \param[in] scheme_context pointer to a scheme context structure to free
1073 *
1074 * \return a NULL pointer
1075 */
1076 /*----------------------------------------------------------------------------*/
1077
1078 void *
cs_cdofb_ac_free_scheme_context(void * scheme_context)1079 cs_cdofb_ac_free_scheme_context(void *scheme_context)
1080 {
1081 cs_cdofb_ac_t *sc = (cs_cdofb_ac_t *)scheme_context;
1082
1083 if (sc == NULL)
1084 return sc;
1085
1086 /* Free BC structure */
1087
1088 sc->pressure_bc = cs_cdo_bc_free(sc->pressure_bc);
1089
1090 /* If the context is not NULL, this means that an Anderson algorithm has been
1091 activated otherwise nothing to do */
1092
1093 cs_iter_algo_aa_free(sc->nl_algo);
1094
1095 BFT_FREE(sc->nl_algo);
1096
1097 /* Other pointers are only shared (i.e. not owner) */
1098
1099 BFT_FREE(sc);
1100
1101 return NULL;
1102 }
1103
1104 /*----------------------------------------------------------------------------*/
1105 /*!
1106 * \brief Start setting-up the Navier-Stokes equations when an AC algorithm
1107 * is used to couple the system.
1108 * No mesh information is available at this stage
1109 *
1110 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1111 * \param[in, out] context pointer to a context structure cast on-the-fly
1112 */
1113 /*----------------------------------------------------------------------------*/
1114
1115 void
cs_cdofb_ac_set_sles(const cs_navsto_param_t * nsp,void * context)1116 cs_cdofb_ac_set_sles(const cs_navsto_param_t *nsp,
1117 void *context)
1118 {
1119 cs_navsto_ac_t *nsc = (cs_navsto_ac_t *)context;
1120
1121 assert(nsp != NULL && nsc != NULL);
1122
1123 const cs_navsto_param_sles_t *nslesp = nsp->sles_param;
1124 cs_equation_param_t *mom_eqp = cs_equation_get_param(nsc->momentum);
1125 int field_id = cs_equation_get_field_id(nsc->momentum);
1126
1127 mom_eqp->sles_param->field_id = field_id;
1128
1129 switch (nslesp->strategy) {
1130
1131 case CS_NAVSTO_SLES_EQ_WITHOUT_BLOCK: /* "Classical" way to set SLES */
1132 cs_equation_param_set_sles(mom_eqp);
1133 break;
1134
1135 case CS_NAVSTO_SLES_BLOCK_MULTIGRID_CG:
1136 #if defined(HAVE_PETSC)
1137 if (mom_eqp->sles_param->amg_type == CS_PARAM_AMG_NONE) {
1138 #if defined(PETSC_HAVE_HYPRE)
1139 mom_eqp->sles_param->amg_type = CS_PARAM_AMG_HYPRE_BOOMER_V;
1140 #else
1141 mom_eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_GAMG_V;
1142 #endif
1143 }
1144
1145 cs_sles_petsc_init();
1146 cs_sles_petsc_define(field_id,
1147 NULL,
1148 MATMPIAIJ,
1149 cs_navsto_sles_amg_block_hook,
1150 (void *)nsp);
1151 #else
1152 bft_error(__FILE__, __LINE__, 0,
1153 "%s: Invalid strategy for solving the linear system %s\n"
1154 " PETSc is required with this option.\n"
1155 " Please build a version of Code_Saturne with the PETSc support.",
1156 __func__, mom_eqp->name);
1157 #endif /* HAVE_PETSC */
1158 break;
1159
1160 default:
1161 bft_error(__FILE__, __LINE__, 0,
1162 "%s: Invalid strategy for solving the linear system %s\n",
1163 __func__, mom_eqp->name);
1164 }
1165
1166 }
1167
1168 /*----------------------------------------------------------------------------*/
1169 /*!
1170 * \brief Solve the unsteady Navier-Stokes system with a CDO face-based scheme
1171 * using a Artificial Compressibility approach and an implicit Euler
1172 * time scheme
1173 *
1174 * \param[in] mesh pointer to a \ref cs_mesh_t structure
1175 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1176 * \param[in, out] scheme_context pointer to a structure cast on-the-fly
1177 */
1178 /*----------------------------------------------------------------------------*/
1179
1180 void
cs_cdofb_ac_compute_implicit(const cs_mesh_t * mesh,const cs_navsto_param_t * nsp,void * scheme_context)1181 cs_cdofb_ac_compute_implicit(const cs_mesh_t *mesh,
1182 const cs_navsto_param_t *nsp,
1183 void *scheme_context)
1184 {
1185 cs_timer_t t_cmpt = cs_timer_time();
1186
1187 /* Retrieve high-level structures */
1188
1189 cs_cdofb_ac_t *sc = (cs_cdofb_ac_t *)scheme_context;
1190 cs_navsto_ac_t *cc = sc->coupling_context;
1191 cs_equation_t *mom_eq = cc->momentum;
1192 cs_cdofb_vecteq_t *mom_eqc= (cs_cdofb_vecteq_t *)mom_eq->scheme_context;
1193 cs_equation_param_t *mom_eqp = mom_eq->param;
1194 cs_equation_builder_t *mom_eqb = mom_eq->builder;
1195
1196 const cs_time_step_t *ts = cs_shared_time_step;
1197 const cs_cdo_connect_t *connect = cs_shared_connect;
1198 const cs_cdo_quantities_t *quant = cs_shared_quant;
1199 const cs_lnum_t n_faces = quant->n_faces;
1200 const cs_range_set_t *rs = connect->range_sets[CS_CDO_CONNECT_FACE_VP0];
1201
1202 /* Retrieve fields */
1203
1204 cs_real_t *vel_c = sc->velocity->val;
1205 cs_real_t *vel_f = mom_eqc->face_values;
1206 cs_real_t *div = sc->divergence->val;
1207 cs_real_t *pr = sc->pressure->val;
1208
1209 /*--------------------------------------------------------------------------
1210 * BUILD: START
1211 *--------------------------------------------------------------------------*/
1212
1213 cs_timer_t t_bld = cs_timer_time();
1214
1215 /* Build an array storing the Dirichlet values at faces.
1216 Evaluation should be performed at t_cur + dt_cur */
1217
1218 cs_cdofb_vecteq_setup(ts->t_cur + ts->dt[0], mesh, mom_eqp, mom_eqb);
1219
1220 /* Initialize the linear system: matrix and rhs */
1221
1222 cs_matrix_t *matrix = cs_matrix_create(cs_shared_ms);
1223 cs_real_t *rhs = NULL;
1224
1225 BFT_MALLOC(rhs, 3*n_faces, cs_real_t);
1226 #if defined(HAVE_OPENMP)
1227 # pragma omp parallel for if (3*n_faces > CS_THR_MIN)
1228 for (cs_lnum_t i = 0; i < 3*n_faces; i++) rhs[i] = 0.0;
1229 #else
1230 memset(rhs, 0, 3*n_faces*sizeof(cs_real_t));
1231 #endif
1232
1233 /* Main function for the building stage */
1234
1235 _implicit_euler_build(nsp,
1236 vel_f, /* previous values for the velocity at faces */
1237 vel_c, /* previous values for the velocity at cells */
1238 pr, /* previous values for the pressure */
1239 sc,
1240 matrix, rhs);
1241
1242 /* Free temporary buffers and structures */
1243
1244 cs_equation_builder_reset(mom_eqb);
1245
1246 /* End of the system building */
1247
1248 cs_timer_t t_tmp = cs_timer_time();
1249 cs_timer_counter_add_diff(&(mom_eqb->tcb), &t_bld, &t_tmp);
1250
1251 /*--------------------------------------------------------------------------
1252 * BUILD: END
1253 *--------------------------------------------------------------------------*/
1254
1255 /* Copy current field values to previous values */
1256
1257 cs_timer_t t_upd = cs_timer_time();
1258
1259 /* Current to previous for main variable fields */
1260
1261 _ac_fields_to_previous(sc, cc);
1262
1263 t_tmp = cs_timer_time();
1264 cs_timer_counter_add_diff(&(mom_eqb->tce), &t_upd, &t_tmp);
1265
1266 /* Solve the linear system (treated as a scalar-valued system
1267 * with 3 times more DoFs) */
1268
1269 cs_real_t normalization = 1.0; /* TODO */
1270 cs_sles_t *sles = cs_sles_find_or_add(mom_eqp->sles_param->field_id, NULL);
1271
1272 int n_solver_iter = cs_equation_solve_scalar_system(3*n_faces,
1273 mom_eqp->sles_param,
1274 matrix,
1275 rs,
1276 normalization,
1277 true, /* rhs_redux */
1278 sles,
1279 vel_f, /* updated here */
1280 rhs);
1281
1282 /* Update pressure, velocity at cells and divergence velocity fields */
1283
1284 t_upd = cs_timer_time();
1285
1286 /* Updates after the resolution:
1287 * 1. the divergence: div = B.u_f
1288 * 2. the cell velocity field
1289 * 2. the mass flux
1290 * 2. the pressure field: pr -= dt / zeta * div(u_f)
1291 */
1292
1293 _ac_compute_div(vel_f, div);
1294
1295 /* Compute values at cells vel_c from values at faces vel_f
1296 vel_c = acc^-1*(RHS - Acf*vel_f) */
1297
1298 cs_static_condensation_recover_vector(connect->c2f,
1299 mom_eqc->rc_tilda, mom_eqc->acf_tilda,
1300 vel_f, vel_c);
1301
1302 /* Compute the new mass flux used as the advection field */
1303
1304 cs_cdofb_navsto_mass_flux(nsp, quant, vel_f, sc->mass_flux_array);
1305
1306 /* Update the pressure knowing the new divergence of the velocity */
1307
1308 _ac_update_pr(ts->t_cur, ts->dt[0], cc->zeta, mom_eqp, mom_eqb, div, pr);
1309
1310 if (nsp->verbosity > 1) {
1311
1312 cs_log_printf(CS_LOG_DEFAULT, " -cvg- NavSto: cumulated_inner_iters: %d\n",
1313 n_solver_iter);
1314 cs_log_printf_flush(CS_LOG_DEFAULT);
1315
1316 }
1317
1318 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_AC_DBG > 2
1319 cs_dbg_darray_to_listing("PRESSURE", n_cells, pr, 9);
1320 cs_dbg_darray_to_listing("VELOCITY_DIV", n_cells, div, 9);
1321 #endif
1322
1323 /* Frees */
1324
1325 BFT_FREE(rhs);
1326 cs_sles_free(sles);
1327 cs_matrix_destroy(&matrix);
1328
1329 t_tmp = cs_timer_time();
1330 cs_timer_counter_add_diff(&(sc->timer), &t_cmpt, &t_tmp);
1331 }
1332
1333 /*----------------------------------------------------------------------------*/
1334 /*!
1335 * \brief Solve the unsteady Navier-Stokes system with a CDO face-based scheme
1336 * using a Artificial Compressibility approach and an implicit Euler
1337 * time scheme
1338 * For Picard - Navier--Stokes problems
1339 *
1340 * \param[in] mesh pointer to a \ref cs_mesh_t structure
1341 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1342 * \param[in, out] scheme_context pointer to a structure cast on-the-fly
1343 */
1344 /*----------------------------------------------------------------------------*/
1345
1346 void
cs_cdofb_ac_compute_implicit_nl(const cs_mesh_t * mesh,const cs_navsto_param_t * nsp,void * scheme_context)1347 cs_cdofb_ac_compute_implicit_nl(const cs_mesh_t *mesh,
1348 const cs_navsto_param_t *nsp,
1349 void *scheme_context)
1350 {
1351 cs_timer_t t_cmpt = cs_timer_time();
1352
1353 /* Retrieve high-level structures */
1354
1355 cs_cdofb_ac_t *sc = (cs_cdofb_ac_t *)scheme_context;
1356 cs_navsto_ac_t *cc = sc->coupling_context;
1357 cs_equation_t *mom_eq = cc->momentum;
1358 cs_cdofb_vecteq_t *mom_eqc= (cs_cdofb_vecteq_t *)mom_eq->scheme_context;
1359 cs_equation_param_t *mom_eqp = mom_eq->param;
1360 cs_equation_builder_t *mom_eqb = mom_eq->builder;
1361 cs_iter_algo_t *nl_algo = sc->nl_algo;
1362
1363 /*--------------------------------------------------------------------------
1364 * INITIAL BUILD: START
1365 *--------------------------------------------------------------------------*/
1366
1367 const cs_time_step_t *ts = cs_shared_time_step;
1368 const cs_cdo_connect_t *connect = cs_shared_connect;
1369 const cs_range_set_t *rs = connect->range_sets[CS_CDO_CONNECT_FACE_VP0];
1370 const cs_cdo_quantities_t *quant = cs_shared_quant;
1371 const cs_lnum_t n_faces = quant->n_faces;
1372
1373 /* Retrieve fields */
1374
1375 cs_real_t *vel_f = mom_eqc->face_values;
1376 cs_real_t *vel_c = sc->velocity->val;
1377 cs_real_t *div = sc->divergence->val;
1378 cs_real_t *pr = sc->pressure->val;
1379
1380 /*--------------------------------------------------------------------------
1381 * BUILD: START
1382 *--------------------------------------------------------------------------*/
1383
1384 cs_timer_t t_bld = cs_timer_time();
1385
1386 /* Build an array storing the Dirichlet values at faces.
1387 Evaluation should be performed at t_cur + dt_cur */
1388
1389 cs_cdofb_vecteq_setup(ts->t_cur + ts->dt[0], mesh, mom_eqp, mom_eqb);
1390
1391 /* Initialize the linear system: matrix and rhs */
1392
1393 cs_matrix_t *matrix = cs_matrix_create(cs_shared_ms);
1394 cs_real_t *rhs = NULL;
1395
1396 BFT_MALLOC(rhs, 3*n_faces, cs_real_t);
1397 #if defined(HAVE_OPENMP)
1398 # pragma omp parallel for if (3*n_faces > CS_THR_MIN)
1399 for (cs_lnum_t i = 0; i < 3*n_faces; i++) rhs[i] = 0.0;
1400 #else
1401 memset(rhs, 0, 3*n_faces*sizeof(cs_real_t));
1402 #endif
1403
1404 /* Main function for the building stage */
1405
1406 _implicit_euler_build(nsp,
1407 vel_f, /* previous values for the velocity at faces */
1408 vel_c, /* previous values for the velocity at cells */
1409 pr, /* previous values for the pressure */
1410 sc,
1411 matrix, rhs);
1412
1413 /* End of the system building */
1414
1415 cs_timer_t t_tmp = cs_timer_time();
1416 cs_timer_counter_add_diff(&(mom_eqb->tcb), &t_bld, &t_tmp);
1417
1418 /*--------------------------------------------------------------------------
1419 * BUILD: END
1420 *--------------------------------------------------------------------------*/
1421
1422 cs_timer_t t_upd = cs_timer_time();
1423
1424 /* Copy current field values to previous values */
1425
1426 _ac_fields_to_previous(sc, cc);
1427
1428 t_tmp = cs_timer_time();
1429 cs_timer_counter_add_diff(&(mom_eqb->tce), &t_upd, &t_tmp);
1430
1431 cs_timer_t t_solve_start = cs_timer_time();
1432
1433 cs_iter_algo_reset(nl_algo);
1434
1435 /* Solve the linear system (treated as a scalar-valued system
1436 * with 3 times more DoFs) */
1437
1438 cs_real_t normalization = 1.0; /* TODO */
1439 cs_sles_t *sles = cs_sles_find_or_add(mom_eqp->sles_param->field_id, NULL);
1440
1441 nl_algo->n_inner_iter = (nl_algo->last_inner_iter =
1442 cs_equation_solve_scalar_system(3*n_faces,
1443 mom_eqp->sles_param,
1444 matrix,
1445 rs,
1446 normalization,
1447 true, /* rhs_redux */
1448 sles,
1449 vel_f,
1450 rhs));
1451
1452 cs_timer_t t_solve_end = cs_timer_time();
1453 cs_timer_counter_add_diff(&(mom_eqb->tcs), &t_solve_start, &t_solve_end);
1454
1455 t_upd = cs_timer_time();
1456
1457 /* Updates after the resolution:
1458 * 1. the divergence: div = B.u_f
1459 * 2. the mass flux
1460 */
1461
1462 _ac_compute_div(vel_f, div);
1463
1464 /* Compute the new mass flux used as the advection field */
1465
1466 cs_cdofb_navsto_mass_flux(nsp, quant, vel_f, sc->mass_flux_array);
1467
1468 t_tmp = cs_timer_time();
1469 cs_timer_counter_add_diff(&(mom_eqb->tce), &t_upd, &t_tmp);
1470
1471 const cs_real_t *pr_c_pre = sc->pressure->val_pre;
1472 const cs_real_t *vel_c_pre = sc->velocity->val_pre;
1473 const cs_real_t *vel_f_pre = mom_eqc->face_values_pre;
1474 assert(vel_c_pre != NULL && vel_f_pre != NULL && pr_c_pre != NULL);
1475
1476 /*--------------------------------------------------------------------------
1477 * PICARD ITERATIONS: START
1478 *--------------------------------------------------------------------------*/
1479
1480 /* Set the normalization of the non-linear algo to the value of the first
1481 mass flux norm */
1482
1483 nl_algo->normalization =
1484 sqrt(cs_cdo_blas_square_norm_pfsf(sc->mass_flux_array));
1485
1486 /* Check the convergence status and update the nl_algo structure related
1487 * to the convergence monitoring */
1488
1489 while (cs_cdofb_navsto_nl_algo_cvg(nsp->sles_param->nl_algo_type,
1490 sc->mass_flux_array_pre,
1491 sc->mass_flux_array,
1492 nl_algo) == CS_SLES_ITERATING) {
1493
1494 /* Main loop on cells to define the linear system to solve */
1495
1496 cs_timer_t t_build_start = cs_timer_time();
1497
1498 /* rhs set to zero and cs_sles_t structure is freed in order to do
1499 * the setup once again since the matrix should be modified */
1500
1501 #if defined(HAVE_OPENMP)
1502 # pragma omp parallel for if (3*n_faces > CS_THR_MIN)
1503 for (cs_lnum_t i = 0; i < 3*n_faces; i++) rhs[i] = 0.0;
1504 #else
1505 memset(rhs, 0, 3*n_faces*sizeof(cs_real_t));
1506 #endif
1507 cs_sles_free(sles); sles = NULL;
1508
1509 /* Main loop on cells to define the linear system to solve */
1510
1511 _implicit_euler_build(nsp,
1512 vel_f_pre, /* velocity at faces: previous values */
1513 vel_c_pre, /* velocity at cells: previous values */
1514 pr_c_pre, /* pressure at cells: previous values */
1515 sc,
1516 matrix, rhs);
1517
1518 /* End of the system building */
1519
1520 cs_timer_t t_build_end = cs_timer_time();
1521 cs_timer_counter_add_diff(&(mom_eqb->tcb), &t_build_start, &t_build_end);
1522
1523 /* Solve the new system */
1524
1525 t_solve_start = cs_timer_time();
1526
1527 /* Redo the setup since the matrix is modified */
1528
1529 sles = cs_sles_find_or_add(mom_eqp->sles_param->field_id, NULL);
1530 cs_sles_setup(sles, matrix);
1531
1532 nl_algo->n_inner_iter += (nl_algo->last_inner_iter =
1533 cs_equation_solve_scalar_system(3*n_faces,
1534 mom_eqp->sles_param,
1535 matrix,
1536 rs,
1537 normalization,
1538 true, /* rhs_redux */
1539 sles,
1540 vel_f,
1541 rhs));
1542
1543 t_solve_end = cs_timer_time();
1544 cs_timer_counter_add_diff(&(mom_eqb->tcs), &t_solve_start, &t_solve_end);
1545
1546 /* Compute the velocity divergence and retrieve its L2-norm */
1547
1548 _ac_compute_div(vel_f, div);
1549
1550 /* Compute the new mass flux used as the advection field */
1551
1552 memcpy(sc->mass_flux_array_pre, sc->mass_flux_array,
1553 n_faces*sizeof(cs_real_t));
1554
1555 cs_cdofb_navsto_mass_flux(nsp, quant, vel_f, sc->mass_flux_array);
1556
1557 } /* Loop on non linear iterations */
1558
1559 /*--------------------------------------------------------------------------
1560 * PICARD ITERATIONS: END
1561 *--------------------------------------------------------------------------*/
1562
1563 if (nsp->verbosity > 1) {
1564
1565 cs_log_printf(CS_LOG_DEFAULT, " -cvg- NavSto: cumulated_inner_iters: %d\n",
1566 nl_algo->n_inner_iter);
1567 cs_log_printf_flush(CS_LOG_DEFAULT);
1568
1569 }
1570
1571 if (nsp->sles_param->nl_algo_type == CS_PARAM_NL_ALGO_PICARD)
1572 cs_iter_algo_post_check(__func__, mom_eqp->name, "Picard", nl_algo);
1573
1574 else {
1575
1576 assert(nsp->sles_param->nl_algo_type == CS_PARAM_NL_ALGO_ANDERSON);
1577
1578 cs_iter_algo_post_check(__func__, mom_eqp->name, "Anderson", nl_algo);
1579
1580 cs_iter_algo_aa_free_arrays(nl_algo->context);
1581
1582 }
1583
1584 /* Update pressure and the cell velocity */
1585
1586 t_upd = cs_timer_time();
1587
1588 /* Updates after the resolution the pressure field:
1589 * pr -= dt / zeta * div(u_f) */
1590
1591 _ac_update_pr(ts->t_cur, ts->dt[0], cc->zeta, mom_eqp, mom_eqb, div, pr);
1592
1593 /* Compute values at cells vel_c from values at faces vel_f
1594 vel_c = acc^-1*(RHS - Acf*vel_f) */
1595
1596 cs_static_condensation_recover_vector(connect->c2f,
1597 mom_eqc->rc_tilda, mom_eqc->acf_tilda,
1598 vel_f, vel_c);
1599
1600 t_tmp = cs_timer_time();
1601 cs_timer_counter_add_diff(&(mom_eqb->tce), &t_upd, &t_tmp);
1602
1603 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_AC_DBG > 2
1604 cs_dbg_darray_to_listing("PRESSURE", n_cells, pr, 9);
1605 cs_dbg_darray_to_listing("VELOCITY_DIV", n_cells, div, 9);
1606 #endif
1607
1608 BFT_FREE(rhs);
1609 cs_equation_builder_reset(mom_eqb);
1610 cs_sles_free(sles);
1611 cs_matrix_destroy(&matrix);
1612
1613 t_tmp = cs_timer_time();
1614 cs_timer_counter_add_diff(&(sc->timer), &t_cmpt, &t_tmp);
1615 }
1616
1617 /*----------------------------------------------------------------------------*/
1618
1619 END_C_DECLS
1620