1 /*============================================================================
2 * Build an algebraic CDO face-based system for the Navier-Stokes equations
3 * and solved it with a prediction/correction 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_cdofb_priv.h"
54 #include "cs_cdofb_scaleq.h"
55 #include "cs_cdofb_vecteq.h"
56 #include "cs_cdofb_navsto.h"
57 #if defined(DEBUG) && !defined(NDEBUG)
58 #include "cs_dbg.h"
59 #endif
60 #include "cs_equation_bc.h"
61 #include "cs_equation_common.h"
62 #include "cs_equation_priv.h"
63 #include "cs_evaluate.h"
64 #include "cs_log.h"
65 #include "cs_math.h"
66 #include "cs_navsto_sles.h"
67 #include "cs_post.h"
68 #include "cs_source_term.h"
69 #include "cs_static_condensation.h"
70 #include "cs_timer.h"
71
72 #if defined(HAVE_PETSC)
73 #include "cs_sles_petsc.h"
74 #endif
75
76 /*----------------------------------------------------------------------------
77 * Header for the current file
78 *----------------------------------------------------------------------------*/
79
80 #include "cs_cdofb_predco.h"
81
82 /*----------------------------------------------------------------------------*/
83
84 BEGIN_C_DECLS
85
86 /*=============================================================================
87 * Additional doxygen documentation
88 *============================================================================*/
89
90 /*!
91 * \file cs_cdofb_predco.c
92 *
93 * \brief Build an algebraic CDO face-based system for the Navier-Stokes
94 * equations and solved it with a prediction/correction algorithm. A first
95 * equation related to the velocity prediction is solved and then a second
96 * equation related to the pressure correction is solved to project the velocity
97 * field into the space of divergence-free field.
98 */
99
100 /*=============================================================================
101 * Local structure definitions
102 *============================================================================*/
103
104 /*! \struct cs_cdofb_predco_t
105 * \brief Context related to CDO face-based discretization when dealing with
106 * Navier-Stokes equations with a prediction/correction algorithm
107 */
108
109 typedef struct {
110
111 /*! \var coupling_context
112
113 * Pointer to a \ref cs_navsto_projection_t_t (owned by
114 * \ref cs_navsto_system_t) containing the settings related to a prjection
115 * or prediction/correction algorithm.
116 */
117
118 cs_navsto_projection_t *coupling_context;
119
120 /*!
121 * @name Main field variables
122 * Fields for every main variable of the equation. Got from cs_navsto_system_t
123 */
124
125 /*! \var velocity
126 * Pointer to \ref cs_field_t (owned by \ref cs_navsto_system_t) containing
127 * the cell DoFs of the velocity
128 */
129
130 cs_field_t *velocity;
131
132 /*! \var pressure
133 * Pointer to \ref cs_field_t (owned by \ref cs_navsto_system_t) containing
134 * the cell DoFs of the pressure
135 */
136
137 cs_field_t *pressure;
138
139 /*! \var divergence
140 * Pointer to \ref cs_real_t containing the values of the divergence on the
141 * cells
142 */
143
144 cs_field_t *divergence;
145
146 /*! \var predicted_velocity_f
147 * Values of the predicted velocity at faces.
148 * This values may not be divergence-free
149 */
150 cs_real_t *predicted_velocity_f;
151
152 /*! \var pressure_f
153 * Values of the pressure at faces.
154 */
155 cs_real_t *pressure_f;
156
157 /*!
158 * @}
159 * @name Advection quantities
160 * Members related to the advection
161 * @{
162 *
163 * \var adv_field
164 * Pointer to the cs_adv_field_t related to the Navier-Stokes eqs (Shared)
165 */
166 cs_adv_field_t *adv_field;
167
168 /*! \var mass_flux_array
169 * Current values of the mass flux at primal faces (Shared)
170 */
171 cs_real_t *mass_flux_array;
172
173 /*! \var mass_flux_array_pre
174 * Previous values of the mass flux at primal faces (Shared)
175 */
176 cs_real_t *mass_flux_array_pre;
177
178 /*!
179 * @}
180 * @name Boundary conditions (BC) management
181 * Functions and elements used for enforcing the BCs
182 * @{
183 *
184 * \var bf_type
185 * Array of boundary type for each boundary face. (Shared)
186 */
187
188 const cs_boundary_type_t *bf_type;
189
190 /*!
191 * \var pressure_bc
192 * Structure storing the metadata after processing the user-defined boundary
193 * conditions related to the pressure field
194 */
195
196 cs_cdo_bc_face_t *pressure_bc;
197 int pressure_rescaling;
198
199 /*! \var apply_fixed_wall
200 * \ref cs_cdo_apply_boundary_t function pointer defining how to apply a
201 * wall boundary (no slip boundary)
202 *
203 * \var apply_sliding_wall
204 * \ref cs_cdo_apply_boundary_t function pointer defining how to apply a
205 * wall boundary (a tangential velocity is specified at the wall)
206 *
207 * \var apply_velocity_inlet
208 * \ref cs_cdo_apply_boundary_t function pointer defining how to apply a
209 * boundary with a fixed velocity at the inlet
210 *
211 * \var apply_symmetry
212 * \ref cs_cdo_apply_boundary_t function pointer defining how to apply a
213 * symmetry boundary
214 */
215
216 cs_cdo_apply_boundary_t *apply_fixed_wall;
217 cs_cdo_apply_boundary_t *apply_sliding_wall;
218 cs_cdo_apply_boundary_t *apply_velocity_inlet;
219 cs_cdo_apply_boundary_t *apply_symmetry;
220
221 /*!
222 * @}
223 * @name Performance monitoring
224 * Monitoring the efficiency of the algorithm used to solve the Navier-Stokes
225 * system
226 * @{
227 */
228
229 /*! \var timer
230 * Cumulated elapsed time for building and solving the Navier--Stokes system
231 */
232 cs_timer_counter_t timer;
233
234 /*! @} */
235
236 } cs_cdofb_predco_t;
237
238 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
239
240 /*=============================================================================
241 * Local Macro definitions and structure definitions
242 *============================================================================*/
243
244 #define CS_CDOFB_PREDCO_DBG 0
245
246 /*============================================================================
247 * Private variables
248 *============================================================================*/
249
250 /* Pointer to shared structures */
251
252 static const cs_cdo_quantities_t *cs_shared_quant;
253 static const cs_cdo_connect_t *cs_shared_connect;
254 static const cs_time_step_t *cs_shared_time_step;
255 static const cs_matrix_structure_t *cs_shared_mom_ms;
256 static const cs_matrix_structure_t *cs_shared_pre_ms;
257
258 /*============================================================================
259 * Private function prototypes
260 *============================================================================*/
261
262 /*----------------------------------------------------------------------------*/
263 /*!
264 * \brief Apply the boundary conditions to the local system when this should
265 * be done before the static condensation
266 *
267 * \param[in] sc pointer to a cs_cdofb_predco_t structure
268 * \param[in] eqp pointer to a cs_equation_param_t structure
269 * \param[in] cm pointer to a cellwise view of the mesh
270 * \param[in] bf_type type of boundary for the boundary face
271 * \param[in] diff_pty pointer to \ref cs_property_data_t for diffusion
272 * \param[in, out] csys pointer to a cellwise view of the system
273 * \param[in, out] cb pointer to a cellwise builder
274 */
275 /*----------------------------------------------------------------------------*/
276
277 static void
_predco_apply_bc_partly(const cs_cdofb_predco_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)278 _predco_apply_bc_partly(const cs_cdofb_predco_t *sc,
279 const cs_equation_param_t *eqp,
280 const cs_cell_mesh_t *cm,
281 const cs_boundary_type_t *bf_type,
282 const cs_property_data_t *diff_pty,
283 cs_cell_sys_t *csys,
284 cs_cell_builder_t *cb)
285 {
286 assert(cs_equation_param_has_diffusion(eqp) == true);
287
288 if (cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
289
290 /* Update the velocity-block and the right-hand side (part related to the
291 * momentum equation). */
292
293 /* Neumann boundary conditions:
294 * The common practice is to define Phi_neu = - lambda * grad(u) . n_fc
295 * An outward flux is a positive flux whereas an inward flux is negative
296 * The minus just above implies the minus just below */
297
298 if (csys->has_nhmg_neumann)
299 for (short int i = 0; i < 3*cm->n_fc; i++)
300 csys->rhs[i] -= csys->neu_values[i];
301
302 const cs_real_t pc = sc->pressure->val[cm->c_id];
303
304 for (short int i = 0; i < csys->n_bc_faces; i++) {
305
306 /* Get the boundary face in the cell numbering */
307
308 const short int f = csys->_f_ids[i];
309 const cs_quant_t pfq = cm->face[f];
310 const cs_real_t f_prs = pfq.meas * pc;
311 cs_real_t *f_rhs = csys->rhs + 3*f;
312
313 if (bf_type[i] & CS_BOUNDARY_IMPOSED_VEL) {
314 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
315 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM) {
316 sc->apply_velocity_inlet(f, eqp, cm, diff_pty, cb, csys);
317 f_rhs[0] -= f_prs * pfq.unitv[0];
318 f_rhs[1] -= f_prs * pfq.unitv[1];
319 f_rhs[2] -= f_prs * pfq.unitv[2];
320 }
321 }
322
323 else if (bf_type[i] & CS_BOUNDARY_WALL) {
324 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
325 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM) {
326 if (bf_type[i] & CS_BOUNDARY_SLIDING_WALL)
327 sc->apply_sliding_wall(f, eqp, cm, diff_pty, cb, csys);
328 else
329 sc->apply_fixed_wall(f, eqp, cm, diff_pty, cb, csys);
330 f_rhs[0] -= f_prs * pfq.unitv[0];
331 f_rhs[1] -= f_prs * pfq.unitv[1];
332 f_rhs[2] -= f_prs * pfq.unitv[2];
333 }
334 }
335
336 else if (bf_type[i] & CS_BOUNDARY_SYMMETRY) {
337
338 /* Always weakly enforce the symmetric constraint on the
339 velocity-block */
340
341 sc->apply_symmetry(f, eqp, cm, diff_pty, cb, csys);
342 f_rhs[0] -= f_prs * pfq.unitv[0];
343 f_rhs[1] -= f_prs * pfq.unitv[1];
344 f_rhs[2] -= f_prs * pfq.unitv[2];
345 }
346
347 /* default: nothing to do (case of a "natural" outlet) */
348
349 } /* Loop on boundary faces */
350
351 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_PREDCO_DBG > 1
352 if (cs_dbg_cw_test(eqp, cm, csys))
353 cs_cell_sys_dump(">> Cell system matrix after partial BC enforcement",
354 csys);
355 #endif
356 } /* Boundary cell */
357 }
358
359 /*----------------------------------------------------------------------------*/
360 /*!
361 * \brief Apply the boundary conditions to the local system when this should
362 * be done after the static condensation
363 * Case of CDO-Fb schemes with Prediction-Correction coupling
364 *
365 * \param[in] sc pointer to a cs_cdofb_predco_t structure
366 * \param[in] eqp pointer to a cs_equation_param_t structure
367 * \param[in] eqb pointer to a cs_equation_builder_t structure
368 * \param[in] cm pointer to a cellwise view of the mesh
369 * \param[in] bf_type type of boundary for the boundary face
370 * \param[in] diff_pty pointer to \ref cs_property_data_t for diffusion
371 * \param[in, out] csys pointer to a cellwise view of the system
372 * \param[in, out] cb pointer to a cellwise builder
373 */
374 /*----------------------------------------------------------------------------*/
375
376 static void
_predco_apply_remaining_bc(const cs_cdofb_predco_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)377 _predco_apply_remaining_bc(const cs_cdofb_predco_t *sc,
378 const cs_equation_param_t *eqp,
379 const cs_equation_builder_t *eqb,
380 const cs_cell_mesh_t *cm,
381 const cs_boundary_type_t *bf_type,
382 const cs_property_data_t *diff_pty,
383 cs_cell_sys_t *csys,
384 cs_cell_builder_t *cb)
385 {
386 /* Internal enforcement of DoFs: Update csys (matrix and rhs) */
387
388 if (cs_equation_param_has_internal_enforcement(eqp)) {
389
390 cs_equation_enforced_internal_dofs(eqb, cb, csys);
391
392 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_PREDCO_DBG > 2
393 if (cs_dbg_cw_test(eqp, cm, csys))
394 cs_cell_sys_dump("\n>> Cell system after the internal enforcement",
395 csys);
396 #endif
397 }
398
399 if (cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
400
401 /* Update the divergence operator and the right-hand side related to the
402 * mass equation.
403 * Enforcement of Dirichlet BC in a stronger way if this is the choice
404 */
405
406 for (short int i = 0; i < csys->n_bc_faces; i++) {
407
408 /* Get the boundary face in the cell numbering */
409
410 const short int f = csys->_f_ids[i];
411
412 if (bf_type[i] & CS_BOUNDARY_IMPOSED_VEL) {
413
414 /* Enforcement of the velocity for the velocity-block
415 * Dirichlet on the three components of the velocity field */
416
417 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_PENALIZED ||
418 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_ALGEBRAIC) {
419 sc->apply_velocity_inlet(f, eqp, cm, diff_pty, cb, csys);
420 }
421 }
422
423 else if (bf_type[i] & CS_BOUNDARY_WALL) {
424
425 /* Enforcement of the velocity for the velocity-block
426 * Dirichlet on the three components of the velocity field */
427
428 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_PENALIZED ||
429 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_ALGEBRAIC) {
430 if (bf_type[i] & CS_BOUNDARY_SLIDING_WALL)
431 sc->apply_sliding_wall(f, eqp, cm, diff_pty, cb, csys);
432 else
433 sc->apply_fixed_wall(f, eqp, cm, diff_pty, cb, csys);
434 }
435 }
436
437 #if 0
438 else if (bf_type[i] & CS_BOUNDARY_SYMMETRY) {
439 /* Weak-enforcement for the velocity-block
440 (cf. _predco_apply_bc_partly) */
441 }
442 #endif
443
444 /* default: nothing to do (case of a "natural" outlet) */
445
446 } /* Loop boundary faces */
447
448 } /* Boundary cell */
449 }
450
451 /*----------------------------------------------------------------------------*/
452 /*!
453 * \brief Prepare and solve the pressure correction step
454 *
455 * \param[in] mesh pointer to a \ref cs_mesh_t structure
456 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
457 * \param[in, out] sc pointer to a scheme context structure
458 */
459 /*----------------------------------------------------------------------------*/
460
461 static void
_solve_pressure_correction(const cs_mesh_t * mesh,const cs_navsto_param_t * nsp,cs_cdofb_predco_t * sc)462 _solve_pressure_correction(const cs_mesh_t *mesh,
463 const cs_navsto_param_t *nsp,
464 cs_cdofb_predco_t *sc)
465 {
466 cs_navsto_projection_t *cc = sc->coupling_context;
467 cs_equation_t *pre_eq = cc->correction;
468 cs_equation_param_t *pre_eqp = cs_equation_get_param(pre_eq);
469 void *pre_eqc = cs_equation_get_scheme_context(pre_eq);
470 cs_equation_builder_t *pre_eqb = pre_eq->builder;
471
472 const cs_time_step_t *ts = cs_shared_time_step;
473 const cs_cdo_connect_t *connect = cs_shared_connect;
474 const cs_cdo_quantities_t *quant = cs_shared_quant;
475 const cs_real_t *velp_f = sc->predicted_velocity_f;
476 const cs_real_t *const pressure_f = sc->pressure_f;
477
478 /* Boundary conditions are always evaluated at t + dt */
479
480 const cs_real_t time_eval = ts->t_cur + ts->dt[0];
481
482 cs_timer_t t_bld = cs_timer_time();
483
484 /* Compute the source term */
485
486 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
487
488 /* Compute the divergence of the predicted velocity */
489
490 const cs_real_t div_c
491 = cs_cdofb_navsto_cell_divergence(c_id, quant, connect->c2f, velp_f);
492
493 cc->div_st[c_id] = -div_c * quant->cell_vol[c_id];
494
495 }
496
497 /* Set the boundary conditions on the pressure increment if needed */
498
499 for (int id = 0; id < nsp->n_pressure_bc_defs; id++) {
500
501 const cs_xdef_t *pdef = nsp->pressure_bc_defs[id];
502 const cs_zone_t *z = cs_boundary_zone_by_id(pdef->z_id);
503
504 if (pdef->meta & CS_CDO_BC_DIRICHLET) {
505
506 switch (pdef->type) {
507
508 case CS_XDEF_BY_ANALYTIC_FUNCTION:
509 /* Evaluate the boundary condition at each boundary face */
510 switch(pre_eqp->dof_reduction) {
511
512 case CS_PARAM_REDUCTION_DERHAM:
513 cs_xdef_eval_at_b_faces_by_analytic(z->n_elts,
514 z->elt_ids,
515 false, /* dense output */
516 mesh,
517 connect,
518 quant,
519 time_eval,
520 pdef->context,
521 cc->bdy_pressure_incr);
522 break;
523
524 case CS_PARAM_REDUCTION_AVERAGE:
525 cs_xdef_eval_avg_at_b_faces_by_analytic(z->n_elts,
526 z->elt_ids,
527 false, /* dense output */
528 mesh,
529 connect,
530 quant,
531 time_eval,
532 pdef->context,
533 pdef->qtype,
534 pdef->dim,
535 cc->bdy_pressure_incr);
536 break;
537
538 default:
539 bft_error(__FILE__, __LINE__, 0,
540 _(" %s: Invalid type of reduction.\n"
541 " Stop computing the Dirichlet value.\n"), __func__);
542
543 } /* switch on reduction */
544
545 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
546 const cs_lnum_t f_id = z->elt_ids[i];
547 cc->bdy_pressure_incr[f_id] -= pressure_f[f_id + quant->n_i_faces];
548 }
549 break;
550
551 case CS_XDEF_BY_VALUE:
552 /* This corresponds to a homogeneous Dirichlet BC */
553 break;
554
555 default:
556 bft_error(__FILE__, __LINE__, 0, "%s: Not implemented yet.", __func__);
557 break;
558 }
559
560 }
561 else
562 bft_error(__FILE__, __LINE__, 0, "%s: Not implemented yet.", __func__);
563
564 } /* Loop on pressure definitions */
565
566 cs_timer_t t_tmp = cs_timer_time();
567 cs_timer_counter_add_diff(&(pre_eqb->tcb), &t_bld, &t_tmp);
568
569 /* Solve the equation related to the pressure increment */
570
571 cs_cdofb_scaleq_solve_steady_state(true, /* cur2prev */
572 mesh,
573 pre_eq->field_id,
574 pre_eqp,
575 pre_eqb,
576 pre_eqc);
577 }
578
579 /*----------------------------------------------------------------------------*/
580 /*!
581 * \brief Update velocity and pressure-related variables after the correction
582 * step
583 *
584 * \param[in, out] sc pointer to a scheme context structure
585 */
586 /*----------------------------------------------------------------------------*/
587
588 static void
_update_variables(cs_cdofb_predco_t * sc)589 _update_variables(cs_cdofb_predco_t *sc)
590 {
591 cs_navsto_projection_t *cc = sc->coupling_context;
592 cs_equation_t *pre_eq = cc->correction;
593 void *pre_eqc = cs_equation_get_scheme_context(pre_eq);
594 cs_equation_t *mom_eq = cc->prediction;
595 cs_equation_builder_t *mom_eqb = mom_eq->builder;
596 cs_cdofb_vecteq_t *mom_eqc= (cs_cdofb_vecteq_t *)mom_eq->scheme_context;
597
598 const cs_cdo_connect_t *connect = cs_shared_connect;
599 const cs_cdo_quantities_t *quant = cs_shared_quant;
600 const cs_lnum_t n_faces = quant->n_faces;
601 const cs_field_t *velp_fld = cc->predicted_velocity;
602 const cs_real_t *const velp_c = velp_fld->val;
603 const cs_real_t *const velp_f = sc->predicted_velocity_f;
604 const cs_real_t *const dp_f = cs_cdofb_scaleq_get_face_values(pre_eqc,
605 false);
606 const cs_real_t *const dp_c = cs_cdofb_scaleq_get_cell_values(pre_eqc,
607 false);
608 const cs_real_t dt_cur = cs_shared_time_step->dt[0];
609
610 /* Variables to update */
611
612 cs_real_t *pr_f = sc->pressure_f;
613 cs_real_t *pr_c = sc->pressure->val; /* cell DoFs for the pressure */
614 cs_real_t *vel_c = sc->velocity->val;
615 cs_real_t *vel_f = mom_eqc->face_values;
616 cs_real_t *div = sc->divergence->val;
617
618 cs_timer_t t_upd = cs_timer_time();
619
620 cs_field_current_to_previous(sc->velocity);
621 cs_field_current_to_previous(sc->pressure);
622 cs_field_current_to_previous(sc->divergence);
623
624 # pragma omp parallel if (quant->n_cells > CS_THR_MIN)
625 {
626 #if defined(HAVE_OPENMP) /* Retrieve the cell mesh structure w.r.t. the mode */
627 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(omp_get_thread_num());
628 #else
629 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(0);
630 #endif
631 cs_cell_sys_t *csys = NULL;
632 cs_cell_builder_t *cb = NULL;
633 cs_cdofb_vecteq_get(&csys, &cb);
634
635 /* Reset the velocity at faces */
636
637 # pragma omp for
638 for (cs_lnum_t i = 0; i < 3*n_faces; i++) vel_f[i] = 0;
639
640 # pragma omp for CS_CDO_OMP_SCHEDULE
641 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
642
643 /* Set the local mesh structure for the current cell */
644
645 cs_cell_mesh_build(c_id,
646 CS_FLAG_COMP_PF | CS_FLAG_COMP_PFQ,
647 connect, quant, cm);
648
649 /* Update the cell pressure */
650
651 pr_c[c_id] += dt_cur*dp_c[c_id];
652
653 /* Evaluate the cell gradient for the pressure increment */
654
655 cs_real_t grd_dp[3] = {0, 0, 0};
656 for (short int f = 0; f < cm->n_fc; f++) {
657 const cs_real_t f_coef = cm->f_sgn[f]*cm->face[f].meas;
658 for (int k = 0; k < 3; k++)
659 grd_dp[k] += f_coef*dp_f[cm->f_ids[f]]*cm->face[f].unitv[k];
660 }
661
662 /* Update the cell velocity */
663
664 for (int k = 0; k < 3; k++)
665 vel_c[3*c_id + k] = velp_c[3*c_id + k] + dt_cur*grd_dp[k];
666
667 /* Partial update of the face velocity:
668 * v_f^(n+1) = vp_f^(n+1) + dt*grd(incr_p)
669 * Now: 0.5*dt_cur*grd_cell(incr_p) for an interior face
670 * or dt_cur*grd_cell(incr_p) for a border face */
671
672 for (short int f = 0; f < cm->n_fc; f++) {
673
674 const cs_lnum_t f_id = cm->f_ids[f];
675
676 cs_real_t f_coef = dt_cur;
677 if (f_id < quant->n_i_faces)
678 f_coef *= 0.5;
679
680 cs_real_t *_vel = vel_f + 3*f_id;
681 for (int k = 0; k < 3; k++)
682 _vel[k] += f_coef*grd_dp[k];
683
684 }
685
686 } /* Loop on cells */
687
688 } /* OpenMP block */
689
690 /* Parallel or periodic sum */
691
692 if (connect->interfaces[CS_CDO_CONNECT_FACE_SP0] != NULL)
693 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_FACE_SP0],
694 n_faces,
695 3,
696 true,
697 CS_REAL_TYPE,
698 vel_f);
699
700 /* Update face-related unknowns */
701
702 # pragma omp parallel for if (n_faces > CS_THR_MIN)
703 for (cs_lnum_t f = 0; f < n_faces; f++) {
704
705 /* p^(n+1) = p^n + delta_p */
706
707 pr_f[f] += dp_f[f]*dt_cur;
708
709 /* v_f^(n+1) = vp_f^(n+1) + dt*grd(incr_p) */
710
711 for (int k = 0; k < 3; k++)
712 vel_f[3*f+k] += velp_f[3*f+k];
713
714 } /* Loop on faces */
715
716 const cs_adjacency_t *c2f = connect->c2f;
717
718 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
719 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
720
721 /* Reset the divergence of the velocity before computing the updated
722 value (vel_f has to be updated first) */
723
724 div[c_id] = 0;
725 for (cs_lnum_t jf = c2f->idx[c_id]; jf < c2f->idx[c_id+1]; jf++) {
726
727 const cs_lnum_t f_id = c2f->ids[jf];
728 const cs_real_t *_vel = vel_f + 3*f_id;
729 const cs_real_t *nf = cs_quant_get_face_vector_area(f_id, quant);
730
731 div[c_id] += c2f->sgn[jf]*cs_math_3_dot_product(_vel, nf);
732
733 } /* Loop on cell faces */
734
735 const cs_real_t ovc = 1./quant->cell_vol[c_id];
736 div[c_id] *= ovc;
737
738 } /* Loop on cells */
739
740 cs_timer_t t_tmp = cs_timer_time();
741 cs_timer_counter_add_diff(&(mom_eqb->tce), &t_upd, &t_tmp);
742
743 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_PREDCO_DBG > 2
744 cs_dbg_darray_to_listing("PRED_VELOCITY_FACE", n_faces, velp_f, 9);
745 cs_dbg_darray_to_listing("VELOCITY_FACE", n_faces, vel_f, 9);
746 cs_dbg_darray_to_listing("PRESSURE", quant->n_cells, pr_c, 9);
747 cs_dbg_darray_to_listing("VELOCITY_DIV", quant->n_cells, div, 9);
748 #endif
749 }
750
751 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
752
753 /*============================================================================
754 * Public function prototypes
755 *============================================================================*/
756
757 /*----------------------------------------------------------------------------*/
758 /*!
759 * \brief Retrieve the values of the pressure at faces
760 *
761 * \param[in] context pointer to a scheme context structure
762 *
763 * \return a pointer to the pressure values
764 */
765 /*----------------------------------------------------------------------------*/
766
767 cs_real_t *
cs_cdofb_predco_get_face_pressure(void * context)768 cs_cdofb_predco_get_face_pressure(void *context)
769 {
770 cs_cdofb_predco_t *sc = (cs_cdofb_predco_t *)context;
771
772 if (sc == NULL)
773 return NULL;
774
775 return sc->pressure_f;
776 }
777
778 /*----------------------------------------------------------------------------*/
779 /*!
780 * \brief Set shared pointers from the main domain members
781 *
782 * \param[in] quant additional mesh quantities struct.
783 * \param[in] connect pointer to a \ref cs_cdo_connect_t struct.
784 * \param[in] time_step pointer to a \ref cs_time_step_t structure
785 */
786 /*----------------------------------------------------------------------------*/
787
788 void
cs_cdofb_predco_init_common(const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect,const cs_time_step_t * time_step)789 cs_cdofb_predco_init_common(const cs_cdo_quantities_t *quant,
790 const cs_cdo_connect_t *connect,
791 const cs_time_step_t *time_step)
792 {
793 /* Assign static const pointers */
794
795 cs_shared_quant = quant;
796 cs_shared_connect = connect;
797 cs_shared_time_step = time_step;
798
799 /*
800 * Matrix structure related to the algebraic system for
801 * the momentum equation --> vector-valued equation
802 * the pressure equation --> scalar-valued equation
803 */
804
805 cs_shared_mom_ms = cs_cdofb_vecteq_matrix_structure();
806 cs_shared_pre_ms = cs_cdofb_scaleq_matrix_structure();
807 }
808
809 /*----------------------------------------------------------------------------*/
810 /*!
811 * \brief Initialize a \ref cs_cdofb_predco_t structure
812 *
813 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
814 * \param[in] adv_field pointer to \ref cs_adv_field_t structure
815 * \param[in] mflux current values of the mass flux across primal faces
816 * \param[in] mflux_pre previous values of the mass flux across primal faces
817 * \param[in] fb_type type of boundary for each boundary face
818 * \param[in] nsc_input pointer to a \ref cs_navsto_predco_t structure
819 *
820 * \return a pointer to a new allocated \ref cs_cdofb_predco_t structure
821 */
822 /*----------------------------------------------------------------------------*/
823
824 void *
cs_cdofb_predco_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)825 cs_cdofb_predco_init_scheme_context(const cs_navsto_param_t *nsp,
826 cs_adv_field_t *adv_field,
827 cs_real_t *mflux,
828 cs_real_t *mflux_pre,
829 cs_boundary_type_t *fb_type,
830 void *nsc_input)
831 {
832 assert(nsp != NULL && nsc_input != NULL);
833
834 const cs_cdo_quantities_t *quant = cs_shared_quant;
835
836 /* Cast the coupling context (CC) */
837
838 cs_navsto_projection_t *cc = (cs_navsto_projection_t *)nsc_input;
839
840 /* Navier-Stokes scheme context (SC) */
841
842 cs_cdofb_predco_t *sc = NULL;
843
844 if (nsp->space_scheme != CS_SPACE_SCHEME_CDOFB)
845 bft_error(__FILE__, __LINE__, 0, " %s: Invalid space scheme.\n", __func__);
846
847 BFT_MALLOC(sc, 1, cs_cdofb_predco_t);
848
849 /* Quantities shared with the cs_navsto_system_t structure */
850
851 sc->coupling_context = cc;
852 sc->adv_field = adv_field;
853 sc->mass_flux_array = mflux;
854 sc->mass_flux_array_pre = mflux_pre;
855
856 /* Quick access to the main fields */
857
858 sc->velocity = cs_field_by_name("velocity");
859 sc->pressure = cs_field_by_name("pressure");
860
861 if (nsp->post_flag & CS_NAVSTO_POST_VELOCITY_DIVERGENCE)
862 sc->divergence = cs_field_by_name("velocity_divergence");
863 else
864 sc->divergence = NULL;
865
866 /* Values of the predicted velocity at faces */
867
868 BFT_MALLOC(sc->predicted_velocity_f, 3*quant->n_faces, cs_real_t);
869 memset(sc->predicted_velocity_f, 0, 3*quant->n_faces*sizeof(cs_real_t));
870
871 /* Values of the pressure at faces */
872
873 BFT_MALLOC(sc->pressure_f, quant->n_faces, cs_real_t);
874 memset(sc->pressure_f, 0, quant->n_faces*sizeof(cs_real_t));
875
876 /* Boundary treatment */
877
878 sc->bf_type = fb_type;
879
880 /* Processing of the pressure boundary condition */
881
882 sc->pressure_bc = cs_cdo_bc_face_define(CS_CDO_BC_HMG_NEUMANN, /* Default */
883 true, /* Steady BC up to now */
884 1, /* Dimension */
885 nsp->n_pressure_bc_defs,
886 nsp->pressure_bc_defs,
887 cs_shared_quant->n_b_faces);
888
889 sc->pressure_rescaling =
890 cs_boundary_need_pressure_rescaling(quant->n_b_faces, fb_type);
891
892 cs_equation_param_t *mom_eqp = cc->prediction->param;
893 cs_equation_builder_t *mom_eqb = cc->prediction->builder;
894
895 mom_eqb->bd_msh_flag |= CS_FLAG_COMP_PFC;
896
897 /* Set the way to enforce the Dirichlet BC on the velocity
898 * "fixed_wall" means a no-slip BC */
899
900 sc->apply_symmetry = cs_cdofb_symmetry;
901
902 switch (mom_eqp->default_enforcement) {
903
904 case CS_PARAM_BC_ENFORCE_ALGEBRAIC:
905 sc->apply_velocity_inlet = cs_cdofb_block_dirichlet_alge;
906 sc->apply_sliding_wall = cs_cdofb_block_dirichlet_alge;
907 sc->apply_fixed_wall = cs_cdofb_block_dirichlet_alge;
908 break;
909
910 case CS_PARAM_BC_ENFORCE_PENALIZED:
911 sc->apply_velocity_inlet = cs_cdofb_block_dirichlet_pena;
912 sc->apply_sliding_wall = cs_cdofb_block_dirichlet_pena;
913 sc->apply_fixed_wall = cs_cdofb_block_dirichlet_pena;
914 break;
915
916 case CS_PARAM_BC_ENFORCE_WEAK_NITSCHE:
917 sc->apply_velocity_inlet = cs_cdofb_block_dirichlet_weak;
918 sc->apply_sliding_wall = cs_cdofb_block_dirichlet_weak;
919 sc->apply_fixed_wall = cs_cdofb_block_dirichlet_weak;
920 break;
921
922 case CS_PARAM_BC_ENFORCE_WEAK_SYM:
923 sc->apply_velocity_inlet = cs_cdofb_block_dirichlet_wsym;
924 sc->apply_sliding_wall = cs_cdofb_block_dirichlet_wsym;
925 sc->apply_fixed_wall = cs_cdofb_block_dirichlet_wsym;
926 break;
927
928 default:
929 bft_error(__FILE__, __LINE__, 0,
930 " %s: Invalid type of algorithm to enforce Dirichlet BC.",
931 __func__);
932
933 }
934
935 /* Monitoring */
936
937 CS_TIMER_COUNTER_INIT(sc->timer);
938
939 return sc;
940 }
941
942 /*----------------------------------------------------------------------------*/
943 /*!
944 * \brief Destroy a \ref cs_cdofb_predco_t structure
945 *
946 * \param[in] scheme_context pointer to a scheme context structure to free
947 *
948 * \return a NULL pointer
949 */
950 /*----------------------------------------------------------------------------*/
951
952 void *
cs_cdofb_predco_free_scheme_context(void * scheme_context)953 cs_cdofb_predco_free_scheme_context(void *scheme_context)
954 {
955 cs_cdofb_predco_t *sc = (cs_cdofb_predco_t *)scheme_context;
956
957 if (sc == NULL)
958 return sc;
959
960 /* Free BC structure */
961
962 sc->pressure_bc = cs_cdo_bc_free(sc->pressure_bc);
963
964 BFT_FREE(sc->predicted_velocity_f);
965 BFT_FREE(sc->pressure_f);
966
967 /* Other pointers are only shared (i.e. not owner) */
968
969 BFT_FREE(sc);
970
971 return NULL;
972 }
973
974 /*----------------------------------------------------------------------------*/
975 /*!
976 * \brief Start setting-up the Navier-Stokes equations when a projection
977 * algorithm is used to couple the system.
978 * No mesh information is available at this stage
979 *
980 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
981 * \param[in, out] context pointer to a context structure cast on-the-fly
982 */
983 /*----------------------------------------------------------------------------*/
984
985 void
cs_cdofb_predco_set_sles(const cs_navsto_param_t * nsp,void * context)986 cs_cdofb_predco_set_sles(const cs_navsto_param_t *nsp,
987 void *context)
988 {
989 cs_navsto_projection_t *nsc = (cs_navsto_projection_t *)context;
990
991 assert(nsp != NULL && nsc != NULL);
992
993 cs_navsto_param_sles_t *nslesp = nsp->sles_param;
994 cs_equation_param_t *mom_eqp = cs_equation_get_param(nsc->prediction);
995 int field_id = cs_equation_get_field_id(nsc->prediction);
996
997 mom_eqp->sles_param->field_id = field_id;
998
999 switch (nslesp->strategy) {
1000
1001 case CS_NAVSTO_SLES_EQ_WITHOUT_BLOCK: /* "Classical" way to set SLES */
1002 cs_equation_param_set_sles(mom_eqp);
1003 break;
1004
1005 case CS_NAVSTO_SLES_BLOCK_MULTIGRID_CG:
1006 #if defined(HAVE_PETSC)
1007 if (mom_eqp->sles_param->amg_type == CS_PARAM_AMG_NONE) {
1008 #if defined(PETSC_HAVE_HYPRE)
1009 mom_eqp->sles_param->amg_type = CS_PARAM_AMG_HYPRE_BOOMER_V;
1010 #else
1011 mom_eqp->sles_param->amg_type = CS_PARAM_AMG_PETSC_GAMG_V;
1012 #endif
1013 }
1014
1015 cs_sles_petsc_init();
1016 cs_sles_petsc_define(field_id,
1017 NULL,
1018 MATMPIAIJ,
1019 cs_navsto_sles_amg_block_hook,
1020 (void *)nsp);
1021 #else
1022 bft_error(__FILE__, __LINE__, 0,
1023 "%s: Invalid strategy for solving the linear system %s\n"
1024 " PETSc is required with this option.\n"
1025 " Please build a version of Code_Saturne with the PETSc support.",
1026 __func__, mom_eqp->name);
1027 #endif /* HAVE_PETSC */
1028 break;
1029
1030 default:
1031 bft_error(__FILE__, __LINE__, 0,
1032 "%s: Invalid strategy for solving the linear system %s\n",
1033 __func__, mom_eqp->name);
1034 }
1035
1036 /* For the correction step, use the generic way to setup the SLES */
1037
1038 cs_equation_param_t *corr_eqp = cs_equation_get_param(nsc->correction);
1039
1040 corr_eqp->sles_param->field_id = cs_equation_get_field_id(nsc->correction);
1041 cs_equation_param_set_sles(corr_eqp);
1042
1043 }
1044
1045 /*----------------------------------------------------------------------------*/
1046 /*!
1047 * \brief Solve the unsteady Navier-Stokes system with a CDO face-based scheme
1048 * using a prediction/correction approach and an implicit Euler time
1049 * scheme
1050 *
1051 * \param[in] mesh pointer to a \ref cs_mesh_t structure
1052 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
1053 * \param[in, out] scheme_context pointer to a structure cast on-the-fly
1054 */
1055 /*----------------------------------------------------------------------------*/
1056
1057 void
cs_cdofb_predco_compute_implicit(const cs_mesh_t * mesh,const cs_navsto_param_t * nsp,void * scheme_context)1058 cs_cdofb_predco_compute_implicit(const cs_mesh_t *mesh,
1059 const cs_navsto_param_t *nsp,
1060 void *scheme_context)
1061 {
1062 CS_UNUSED(nsp);
1063
1064 cs_timer_t t_cmpt = cs_timer_time();
1065
1066 const cs_time_step_t *ts = cs_shared_time_step;
1067 const cs_cdo_connect_t *connect = cs_shared_connect;
1068 const cs_range_set_t *mom_rs = connect->range_sets[CS_CDO_CONNECT_FACE_VP0];
1069 const cs_cdo_quantities_t *quant = cs_shared_quant;
1070 const cs_lnum_t n_faces = quant->n_faces;
1071
1072 /* Retrieve high-level structures */
1073
1074 cs_cdofb_predco_t *sc = (cs_cdofb_predco_t *)scheme_context;
1075 cs_navsto_projection_t *cc = sc->coupling_context;
1076 cs_equation_t *mom_eq = cc->prediction;
1077 cs_cdofb_vecteq_t *mom_eqc= (cs_cdofb_vecteq_t *)mom_eq->scheme_context;
1078 cs_equation_param_t *mom_eqp = mom_eq->param;
1079 cs_equation_builder_t *mom_eqb = mom_eq->builder;
1080
1081 assert(cs_equation_param_has_time(mom_eqp) == true);
1082 assert(mom_eqp->time_scheme == CS_TIME_SCHEME_EULER_IMPLICIT);
1083
1084 /* Retrieve fields */
1085
1086 cs_real_t *pr_c = sc->pressure->val; /* cell DoFs for the pressure */
1087 cs_field_t *vel_fld = sc->velocity;
1088 cs_real_t *vel_c = vel_fld->val;
1089
1090 /*--------------------------------------------------------------------------
1091 * BUILD: START
1092 *--------------------------------------------------------------------------*/
1093
1094 cs_timer_t t_bld = cs_timer_time();
1095
1096 /* Build an array storing the Dirichlet values at faces and if needed one
1097 * computes the enforced values. Evaluation should be performed at t_cur +
1098 * dt_cur
1099 */
1100
1101 cs_cdofb_vecteq_setup(ts->t_cur + ts->dt[0], mesh, mom_eqp, mom_eqb);
1102
1103 /* Initialize the local system: matrix and rhs */
1104
1105 cs_matrix_t *matrix = cs_matrix_create(cs_shared_mom_ms);
1106 cs_real_t *rhs = NULL;
1107
1108 BFT_MALLOC(rhs, 3*n_faces, cs_real_t);
1109 # pragma omp parallel for if (3*n_faces > CS_THR_MIN)
1110 for (cs_lnum_t i = 0; i < 3*n_faces; i++) rhs[i] = 0.0;
1111
1112 /* Initialize the structure to assemble values */
1113
1114 cs_matrix_assembler_values_t *mav =
1115 cs_matrix_assembler_values_init(matrix, NULL, NULL);
1116
1117 # pragma omp parallel if (quant->n_cells > CS_THR_MIN) \
1118 shared(quant, connect, ts, mom_eqp, mom_eqb, mom_eqc, sc, rhs, matrix, \
1119 mav, mom_rs, vel_c, pr_c)
1120 {
1121 #if defined(HAVE_OPENMP) /* Determine the default number of OpenMP threads */
1122 int t_id = omp_get_thread_num();
1123 #else
1124 int t_id = 0;
1125 #endif
1126
1127 /* Each thread get back its related structures:
1128 Get the cell-wise view of the mesh and the algebraic system */
1129
1130 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
1131 cs_cdofb_navsto_builder_t nsb = cs_cdofb_navsto_create_builder(nsp,
1132 connect);
1133 cs_equation_assemble_t *eqa = cs_equation_assemble_get(t_id);
1134 cs_hodge_t *diff_hodge =
1135 (mom_eqc->diffusion_hodge == NULL) ? NULL:mom_eqc->diffusion_hodge[t_id];
1136 cs_hodge_t *mass_hodge =
1137 (mom_eqc->mass_hodge == NULL) ? NULL : mom_eqc->mass_hodge[t_id];
1138
1139 cs_cell_sys_t *csys = NULL;
1140 cs_cell_builder_t *cb = NULL;
1141
1142 cs_cdofb_vecteq_get(&csys, &cb);
1143
1144 const cs_real_t t_cur = ts->t_cur;
1145 const cs_real_t dt_cur = ts->dt[0];
1146 const cs_real_t t_eval = t_cur + dt_cur;
1147 const cs_real_t inv_dtcur = 1./dt_cur;
1148
1149 /* Set times at which one evaluates quantities when needed */
1150
1151 cb->t_pty_eval = t_eval;
1152 cb->t_bc_eval = t_eval;
1153 cb->t_st_eval = t_eval;
1154
1155 /* Initialization of the values of properties */
1156
1157 cs_equation_init_properties(mom_eqp, mom_eqb, diff_hodge, cb);
1158
1159 /* ---------------------------------------------
1160 * Main loop on cells to build the linear system
1161 * --------------------------------------------- */
1162
1163 # pragma omp for CS_CDO_OMP_SCHEDULE
1164 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1165
1166 /* Set the current cell flag */
1167
1168 cb->cell_flag = connect->cell_flag[c_id];
1169
1170 /* Set the local mesh structure for the current cell */
1171
1172 cs_cell_mesh_build(c_id,
1173 cs_equation_cell_mesh_flag(cb->cell_flag, mom_eqb),
1174 connect, quant, cm);
1175
1176 /* Starts from the stationary Stokes problem where
1177 *
1178 * | | |
1179 * | A | Bt | B is the divergence (Bt the gradient)
1180 * | | | A is csys->mat in what follows
1181 * |--------|---------| The viscous part arising from the CDO-Fb
1182 * | | | schemes for vector-valued variables
1183 * | B | 0 |
1184 * | | |
1185 */
1186
1187 /* Set the local (i.e. cellwise) structures for the current cell */
1188
1189 cs_cdofb_vecteq_init_cell_system(cm, mom_eqp, mom_eqb,
1190 mom_eqc->face_values, vel_c,
1191 NULL, NULL, /* no n-1 state is given */
1192 csys, cb);
1193
1194 /* 1- SETUP THE NAVSTO LOCAL BUILDER *
1195 * ================================= *
1196 * - Set the type of boundary
1197 * - Set the pressure boundary conditions (if required)
1198 * - Define the divergence operator used in the linear system (div_op is
1199 * equal to minus the divergence)
1200 */
1201
1202 cs_cdofb_navsto_define_builder(cb->t_bc_eval, nsp, cm, csys,
1203 sc->pressure_bc, sc->bf_type, &nsb);
1204
1205 /* 2- VELOCITY (VECTORIAL) EQUATION */
1206 /* ================================ */
1207
1208 cs_cdofb_vecteq_conv_diff_reac(mom_eqp, mom_eqb, mom_eqc, cm,
1209 mass_hodge, diff_hodge, csys, cb);
1210
1211 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_PREDCO_DBG > 1
1212 if (cs_dbg_cw_test(mom_eqp, cm, csys))
1213 cs_cell_sys_dump(">> Cell system after conv/diff/reac", csys);
1214 #endif
1215
1216 /* 3- SOURCE TERM COMPUTATION (for the momentum equation) */
1217 /* ====================================================== */
1218
1219 const bool has_sourceterm = cs_equation_param_has_sourceterm(mom_eqp);
1220 if (has_sourceterm)
1221 cs_cdofb_vecteq_sourceterm(cm, mom_eqp, cb->t_st_eval,
1222 1., /* scaling */
1223 mass_hodge,
1224 cb, mom_eqb, csys);
1225
1226 /* OTHER RHS CONTRIBUTIONS
1227 * =======================
1228 * Apply the operator gradient to the pressure field and add it to the
1229 * rhs */
1230
1231 const short int n_fc = cm->n_fc, nf_dofs = 3*n_fc;
1232
1233 cs_sdm_add_scalvect(nf_dofs, -pr_c[c_id], nsb.div_op, csys->rhs);
1234
1235 /* First part of the BOUNDARY CONDITIONS
1236 * ===================
1237 * Apply a part of BC before the time scheme */
1238
1239 _predco_apply_bc_partly(sc, mom_eqp, cm, nsb.bf_type,
1240 diff_hodge->pty_data, csys, cb);
1241
1242 /* 4- UNSTEADY TERM + TIME SCHEME
1243 * ============================== */
1244
1245 if (mom_eqb->sys_flag & CS_FLAG_SYS_TIME_DIAG) {
1246
1247 /* Mass lumping or Hodge-Voronoi */
1248
1249 const double ptyc = cb->tpty_val * cm->vol_c * inv_dtcur;
1250
1251 /* Get the cell-cell block */
1252
1253 cs_sdm_t *acc = cs_sdm_get_block(csys->mat, n_fc, n_fc);
1254
1255 for (short int k = 0; k < 3; k++) {
1256 csys->rhs[3*n_fc + k] += ptyc * csys->val_n[3*n_fc+k];
1257 /* Simply add an entry in mat[cell, cell] */
1258 acc->val[4*k] += ptyc;
1259 }
1260
1261 }
1262 else
1263 bft_error(__FILE__, __LINE__, 0,
1264 "Only diagonal time treatment available so far.");
1265
1266 /* 5- STATIC CONDENSATION
1267 * ======================
1268 * Static condensation of the local system matrix of size n_fc + 1 into
1269 * a matrix of size n_fc.
1270 * Store data in rc_tilda and acf_tilda to compute the values at cell
1271 * centers after solving the system */
1272
1273 cs_static_condensation_vector_eq(connect->c2f,
1274 mom_eqc->rc_tilda, mom_eqc->acf_tilda,
1275 cb, csys);
1276
1277 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_PREDCO_DBG > 1
1278 if (cs_dbg_cw_test(mom_eqp, cm, csys))
1279 cs_cell_sys_dump(">> Cell system matrix after static condensation",
1280 csys);
1281 #endif
1282
1283 /* 6- Remaining part of BOUNDARY CONDITIONS
1284 * ======================================== */
1285
1286 _predco_apply_remaining_bc(sc, mom_eqp, mom_eqb, cm, nsb.bf_type,
1287 diff_hodge->pty_data, csys, cb);
1288
1289 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_PREDCO_DBG > 0
1290 if (cs_dbg_cw_test(mom_eqp, cm, csys))
1291 cs_cell_sys_dump(">> (FINAL) Cell system matrix", csys);
1292 #endif
1293
1294 /* ASSEMBLY PROCESS */
1295 /* ================ */
1296
1297 cs_cdofb_vecteq_assembly(csys, mom_rs, cm, has_sourceterm,
1298 mom_eqc, eqa, mav, rhs);
1299
1300 } /* Main loop on cells */
1301
1302 /* Free temporary buffers */
1303
1304 cs_cdofb_navsto_free_builder(&nsb);
1305
1306 } /* OPENMP Block */
1307
1308 cs_matrix_assembler_values_done(mav); /* optional */
1309
1310 /* Free temporary buffers and structures */
1311
1312 cs_equation_builder_reset(mom_eqb);
1313 cs_matrix_assembler_values_finalize(&mav);
1314
1315 /* End of the system building */
1316
1317 cs_timer_t t_tmp = cs_timer_time();
1318 cs_timer_counter_add_diff(&(mom_eqb->tcb), &t_bld, &t_tmp);
1319
1320 /*--------------------------------------------------------------------------
1321 * BUILD: END
1322 *--------------------------------------------------------------------------*/
1323
1324 /* Copy current field values to previous values */
1325
1326 cs_field_t *velp_fld = cc->predicted_velocity;
1327
1328 cs_field_current_to_previous(velp_fld);
1329
1330 cs_real_t *velp_c = velp_fld->val;
1331 cs_real_t *velp_f = sc->predicted_velocity_f;
1332
1333 /* Solve the linear system (treated as a scalar-valued system
1334 * with 3 times more DoFs) */
1335
1336 cs_real_t normalization = 1.0; /* TODO */
1337 cs_sles_t *sles = cs_sles_find_or_add(mom_eqp->sles_param->field_id, NULL);
1338
1339 cs_equation_solve_scalar_system(3*n_faces,
1340 mom_eqp->sles_param,
1341 matrix,
1342 mom_rs,
1343 normalization,
1344 true, /* rhs_redux */
1345 sles,
1346 velp_f,
1347 rhs);
1348
1349 /* Update pressure, velocity and divergence fields */
1350
1351 cs_timer_t t_upd = cs_timer_time();
1352
1353 /* Compute values at cells pc from values at faces pf
1354 pc = acc^-1*(RHS - Acf*pf) */
1355
1356 cs_static_condensation_recover_vector(connect->c2f,
1357 mom_eqc->rc_tilda, mom_eqc->acf_tilda,
1358 velp_f, velp_c);
1359
1360 BFT_FREE(rhs);
1361 cs_sles_free(sles);
1362 cs_matrix_destroy(&matrix);
1363
1364 t_tmp = cs_timer_time();
1365 cs_timer_counter_add_diff(&(mom_eqb->tce), &t_upd, &t_tmp);
1366
1367 t_tmp = cs_timer_time();
1368 cs_timer_counter_add_diff(&(sc->timer), &t_cmpt, &t_tmp);
1369
1370 /* SECOND MAJOR STEP: Solve the equation related to the correction step
1371 * Correction step = Solve a diffusion equation for the pressure increment */
1372
1373 _solve_pressure_correction(mesh, nsp, sc);
1374
1375 /* LAST MAJOR STEP: Update the pressure and the velocity */
1376
1377 _update_variables(sc);
1378
1379 if (sc->pressure_rescaling == CS_BOUNDARY_PRESSURE_RESCALING)
1380 cs_cdofb_navsto_rescale_pressure_to_ref(nsp, quant, pr_c);
1381 }
1382
1383 /*----------------------------------------------------------------------------*/
1384
1385 END_C_DECLS
1386