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