1 /*============================================================================
2 * Functions shared among all face-based schemes for the discretization of the
3 * Navier-Stokes system
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 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include <bft_mem.h>
46
47 #include "cs_blas.h"
48 #include "cs_cdo_bc.h"
49 #include "cs_cdo_blas.h"
50 #if defined(DEBUG) && !defined(NDEBUG)
51 #include "cs_dbg.h"
52 #endif
53 #include "cs_equation_bc.h"
54 #include "cs_equation_priv.h"
55 #include "cs_evaluate.h"
56 #include "cs_log.h"
57 #include "cs_math.h"
58 #include "cs_navsto_coupling.h"
59 #include "cs_navsto_param.h"
60 #include "cs_parall.h"
61 #include "cs_post.h"
62 #include "cs_reco.h"
63 #include "cs_sdm.h"
64 #include "cs_timer.h"
65
66 /*----------------------------------------------------------------------------
67 * Header for the current file
68 *----------------------------------------------------------------------------*/
69
70 #include "cs_cdofb_navsto.h"
71
72 /*----------------------------------------------------------------------------*/
73
74 BEGIN_C_DECLS
75
76 /*=============================================================================
77 * Additional doxygen documentation
78 *============================================================================*/
79
80 /*!
81 * \file cs_cdofb_navsto.c
82 *
83 * \brief Shared functions among all face-based schemes for building and
84 * solving Stokes and Navier-Stokes problem
85 *
86 */
87
88 /*=============================================================================
89 * Local structure definitions
90 *============================================================================*/
91
92 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
93
94 /*=============================================================================
95 * Local Macro definitions and structure definitions
96 *============================================================================*/
97
98 #define CS_CDOFB_NAVSTO_DBG 0
99
100 /* Redefined the name of functions from cs_math to get shorter names */
101
102 #define _dp3 cs_math_3_dot_product
103
104 /*============================================================================
105 * Private variables
106 *============================================================================*/
107
108 static cs_cdofb_navsto_boussinesq_type_t cs_cdofb_navsto_boussinesq_type =
109 CS_CDOFB_NAVSTO_BOUSSINESQ_FACE_DOF;
110
111 /*----------------------------------------------------------------------------*/
112 /*!
113 * \brief Print header before dumping information
114 *
115 * \param[in] algo_name name of the algorithm
116 */
117 /*----------------------------------------------------------------------------*/
118
119 static inline void
_nl_algo_print_header(const char * algo_name)120 _nl_algo_print_header(const char *algo_name)
121 {
122 assert(algo_name != NULL);
123 cs_log_printf(CS_LOG_DEFAULT,
124 "%12s.It Algo.Res Inner Cumul Tolerance\n",
125 algo_name);
126 }
127
128 /*----------------------------------------------------------------------------*/
129 /*!
130 * \brief Print a new line in the iterative process
131 *
132 * \param[in] algo_name name of the algorithm
133 * \param[in] algo pointer to cs_iter_algo_t structure
134 */
135 /*----------------------------------------------------------------------------*/
136
137 static inline void
_nl_algo_print_entry(const char * algo_name,const cs_iter_algo_t * algo)138 _nl_algo_print_entry(const char *algo_name,
139 const cs_iter_algo_t *algo)
140 {
141 assert(algo_name != NULL);
142 cs_log_printf(CS_LOG_DEFAULT,
143 "%12s.It%02d %5.3e %5d %5d %6.4e\n",
144 algo_name, algo->n_algo_iter, algo->res, algo->last_inner_iter,
145 algo->n_inner_iter, algo->tol);
146 cs_log_printf_flush(CS_LOG_DEFAULT);
147 }
148
149 /*============================================================================
150 * Private function prototypes
151 *============================================================================*/
152
153 /*----------------------------------------------------------------------------*/
154 /*!
155 * \brief Compute \f$ \int_{fb} \nabla (u) \cdot \nu_{fb} v \f$ where \p fb
156 * is a boundary face (Co+St algorithm)
157 *
158 * \param[in] fb index of the boundary face on the local numbering
159 * \param[in] beta value of coefficient in front of the STAB part
160 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
161 * \param[in] kappa_f diffusion property against face vector for all
162 * faces
163 * \param[in, out] ntrgrd pointer to a local matrix structure
164 */
165 /*----------------------------------------------------------------------------*/
166
167 static void
_normal_flux_reco(short int fb,const double beta,const cs_cell_mesh_t * cm,const cs_real_3_t * kappa_f,cs_sdm_t * ntrgrd)168 _normal_flux_reco(short int fb,
169 const double beta,
170 const cs_cell_mesh_t *cm,
171 const cs_real_3_t *kappa_f,
172 cs_sdm_t *ntrgrd)
173 {
174 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
175 CS_FLAG_COMP_PFC));
176 assert(cm->f_sgn[fb] == 1); /* +1 because it's a boundary face */
177
178 const short int nfc = cm->n_fc;
179 const double inv_volc = 1./cm->vol_c;
180 const cs_quant_t pfbq = cm->face[fb];
181 const cs_nvec3_t debq = cm->dedge[fb];
182
183 /* |fb|^2 * nu_{fb}^T.kappa.nu_{fb} */
184
185 const double stab_scaling =
186 beta * pfbq.meas * _dp3(kappa_f[fb], pfbq.unitv) / cm->pvol_f[fb];
187
188 /* Get the 'fb' row */
189
190 cs_real_t *ntrgrd_fb = ntrgrd->val + fb * (nfc + 1);
191 double row_sum = 0.0;
192
193 for (short int f = 0; f < nfc; f++) {
194
195 const cs_quant_t pfq = cm->face[f];
196
197 /* consistent part */
198
199 const double consist_scaling = cm->f_sgn[f] * pfq.meas * inv_volc;
200 const double consist_part = consist_scaling * _dp3(kappa_f[fb], pfq.unitv);
201
202 /* stabilization part */
203
204 double stab_part = -consist_scaling*debq.meas*_dp3(debq.unitv, pfq.unitv);
205 if (f == fb) stab_part += 1;
206 stab_part *= stab_scaling;
207
208 const double fb_f_part = consist_part + stab_part;
209
210 ntrgrd_fb[f] -= fb_f_part; /* Minus because -du/dn */
211 row_sum += fb_f_part;
212
213 } /* Loop on f */
214
215 /* Cell column */
216
217 ntrgrd_fb[nfc] += row_sum;
218 }
219
220 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
221
222 /*============================================================================
223 * Public function prototypes
224 *============================================================================*/
225
226 /*----------------------------------------------------------------------------*/
227 /*!
228 * \brief Set the way to compute the Boussinesq approximation
229 *
230 * \param[in] type type of algorithm to use
231 */
232 /*----------------------------------------------------------------------------*/
233
234 void
cs_cdofb_navsto_set_boussinesq_algo(cs_cdofb_navsto_boussinesq_type_t type)235 cs_cdofb_navsto_set_boussinesq_algo(cs_cdofb_navsto_boussinesq_type_t type)
236 {
237 cs_cdofb_navsto_boussinesq_type = type;
238 }
239
240 /*----------------------------------------------------------------------------*/
241 /*!
242 * \brief Create and allocate a local NavSto builder when Fb schemes are used
243 *
244 * \param[in] nsp set of parameters to define the NavSto system
245 * \param[in] connect pointer to a cs_cdo_connect_t structure
246 *
247 * \return a cs_cdofb_navsto_builder_t structure
248 */
249 /*----------------------------------------------------------------------------*/
250
251 cs_cdofb_navsto_builder_t
cs_cdofb_navsto_create_builder(const cs_navsto_param_t * nsp,const cs_cdo_connect_t * connect)252 cs_cdofb_navsto_create_builder(const cs_navsto_param_t *nsp,
253 const cs_cdo_connect_t *connect)
254 {
255 cs_cdofb_navsto_builder_t nsb = {.rho_c = 1.,
256 .div_op = NULL,
257 .bf_type = NULL,
258 .pressure_bc_val = NULL};
259
260 assert(nsp != NULL);
261 nsb.rho_c = nsp->mass_density->ref_value;
262
263 if (connect == NULL)
264 return nsb;
265
266 BFT_MALLOC(nsb.div_op, 3*connect->n_max_fbyc, cs_real_t);
267 BFT_MALLOC(nsb.bf_type, connect->n_max_fbyc, cs_boundary_type_t);
268 BFT_MALLOC(nsb.pressure_bc_val, connect->n_max_fbyc, cs_real_t);
269
270 return nsb;
271 }
272
273 /*----------------------------------------------------------------------------*/
274 /*!
275 * \brief Destroy the given cs_cdofb_navsto_builder_t structure
276 *
277 * \param[in, out] nsb pointer to the cs_cdofb_navsto_builder_t to free
278 */
279 /*----------------------------------------------------------------------------*/
280
281 void
cs_cdofb_navsto_free_builder(cs_cdofb_navsto_builder_t * nsb)282 cs_cdofb_navsto_free_builder(cs_cdofb_navsto_builder_t *nsb)
283 {
284 if (nsb != NULL) {
285 BFT_FREE(nsb->div_op);
286 BFT_FREE(nsb->bf_type);
287 BFT_FREE(nsb->pressure_bc_val);
288 }
289 }
290
291 /*----------------------------------------------------------------------------*/
292 /*!
293 * \brief Set the members of the cs_cdofb_navsto_builder_t structure
294 *
295 * \param[in] t_eval time at which one evaluates the pressure BC
296 * \param[in] nsp set of parameters to define the NavSto system
297 * \param[in] cm cellwise view of the mesh
298 * \param[in] csys cellwise view of the algebraic system
299 * \param[in] pr_bc set of definitions for the presuure BCs
300 * \param[in] bf_type type of boundaries for all boundary faces
301 * \param[in, out] nsb builder to update
302 */
303 /*----------------------------------------------------------------------------*/
304
305 void
cs_cdofb_navsto_define_builder(cs_real_t t_eval,const cs_navsto_param_t * nsp,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,const cs_cdo_bc_face_t * pr_bc,const cs_boundary_type_t * bf_type,cs_cdofb_navsto_builder_t * nsb)306 cs_cdofb_navsto_define_builder(cs_real_t t_eval,
307 const cs_navsto_param_t *nsp,
308 const cs_cell_mesh_t *cm,
309 const cs_cell_sys_t *csys,
310 const cs_cdo_bc_face_t *pr_bc,
311 const cs_boundary_type_t *bf_type,
312 cs_cdofb_navsto_builder_t *nsb)
313 {
314 assert(cm != NULL && csys != NULL && nsp != NULL); /* sanity checks */
315
316 const short int n_fc = cm->n_fc;
317
318 /* Update the value of the mass density for the current cell if needed */
319 /* TODO: Case of a uniform but not constant in time */
320
321 if (!cs_property_is_uniform(nsp->mass_density))
322 nsb->rho_c = cs_property_value_in_cell(cm, nsp->mass_density, t_eval);
323
324 /* Build the divergence operator:
325 * D(\hat{u}) = \frac{1}{|c|} \sum_{f_c} \iota_{fc} u_f.f
326 * But in the linear system what appears is after integration
327 * [[ -div(u), q ]]_{P_c} = -|c| div(u)_c q_c
328 * Thus, the volume in the divergence drops
329 */
330
331 for (short int f = 0; f < n_fc; f++) {
332
333 const cs_quant_t pfq = cm->face[f];
334 const cs_real_t sgn_f = -cm->f_sgn[f] * pfq.meas;
335
336 cs_real_t *_div_f = nsb->div_op + 3*f;
337 _div_f[0] = sgn_f * pfq.unitv[0];
338 _div_f[1] = sgn_f * pfq.unitv[1];
339 _div_f[2] = sgn_f * pfq.unitv[2];
340
341 } /* Loop on cell faces */
342
343 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOFB_NAVSTO_DBG > 2
344 if (cs_dbg_cw_test(NULL, cm, csys)) {
345 const cs_real_t sovc = -1. / cm->vol_c;
346 # pragma omp critical
347 {
348 cs_log_printf(CS_LOG_DEFAULT, ">> Divergence:\n");
349 for (short int f = 0; f < cm->n_fc; f++)
350 cs_log_printf(CS_LOG_DEFAULT, " f%2d: %- .4e, %- .4e, %- .4e\n",
351 f, nsb->div_op[3*f]*sovc, nsb->div_op[3*f+1]*sovc,
352 nsb->div_op[3*f+2]*sovc);
353 } /* Critical section */
354 }
355 #endif
356
357 /* Build local arrays related to the boundary conditions */
358
359 for (short int i = 0; i < csys->n_bc_faces; i++) {
360
361 /* Get the boundary face in the cell numbering and the boundary face id in
362 the mesh numbering */
363
364 const short int f = csys->_f_ids[i];
365 const cs_lnum_t bf_id = cm->f_ids[f] - cm->bface_shift;
366
367 /* Set the type of boundary */
368
369 nsb->bf_type[i] = bf_type[bf_id];
370
371 /* Set the pressure BC if required */
372
373 if (nsb->bf_type[i] & CS_BOUNDARY_IMPOSED_P) {
374
375 assert(nsb->bf_type[i] & (CS_BOUNDARY_INLET | CS_BOUNDARY_OUTLET));
376
377 /* Add a Dirichlet for the pressure field */
378
379 const short int def_id = pr_bc->def_ids[bf_id];
380 const cs_xdef_t *def = nsp->pressure_bc_defs[def_id];
381 assert(pr_bc != NULL);
382
383 switch(def->type) {
384 case CS_XDEF_BY_VALUE:
385 {
386 const cs_real_t *constant_val = (cs_real_t *)def->context;
387 nsb->pressure_bc_val[i] = constant_val[0];
388 }
389 break;
390
391 case CS_XDEF_BY_ARRAY:
392 {
393 cs_xdef_array_context_t *c = (cs_xdef_array_context_t *)def->context;
394 assert(c->stride == 1);
395 assert(cs_flag_test(c->loc, cs_flag_primal_face));
396 nsb->pressure_bc_val[i] = c->values[bf_id];
397 }
398 break;
399
400 case CS_XDEF_BY_ANALYTIC_FUNCTION:
401 switch(nsp->dof_reduction_mode) {
402
403 case CS_PARAM_REDUCTION_DERHAM:
404 cs_xdef_cw_eval_at_xyz_by_analytic(cm, 1, cm->face[f].center,
405 t_eval,
406 def->context,
407 nsb->pressure_bc_val + i);
408 break;
409
410 case CS_PARAM_REDUCTION_AVERAGE:
411 cs_xdef_cw_eval_scalar_face_avg_by_analytic(cm, f, t_eval,
412 def->context,
413 def->qtype,
414 nsb->pressure_bc_val + i);
415 break;
416
417 default:
418 bft_error(__FILE__, __LINE__, 0,
419 _(" %s: Invalid type of reduction.\n"
420 " Stop computing the Dirichlet value.\n"), __func__);
421
422 } /* switch on reduction */
423 break;
424
425 default:
426 bft_error(__FILE__, __LINE__, 0,
427 _(" %s: Invalid type of definition.\n"
428 " Stop computing the Dirichlet value.\n"), __func__);
429 break;
430
431 } /* def->type */
432
433 }
434 else
435 nsb->pressure_bc_val[i] = 0.;
436
437 } /* Loop on boundary faces */
438
439 }
440
441 /*----------------------------------------------------------------------------*/
442 /*!
443 * \brief Compute the mass flux playing the role of the advection field in
444 * the Navier-Stokes equations
445 * One considers the mass flux across primal faces which relies on the
446 * velocity vector defined on each face.
447 *
448 * \param[in] nsp set of parameters to define the NavSto system
449 * \param[in] quant set of additional geometrical quantities
450 * \param[in] face_vel velocity vectors for each face
451 * \param[in, out] mass_flux array of mass flux values to update (allocated)
452 */
453 /*----------------------------------------------------------------------------*/
454
455 void
cs_cdofb_navsto_mass_flux(const cs_navsto_param_t * nsp,const cs_cdo_quantities_t * quant,const cs_real_t * face_vel,cs_real_t * mass_flux)456 cs_cdofb_navsto_mass_flux(const cs_navsto_param_t *nsp,
457 const cs_cdo_quantities_t *quant,
458 const cs_real_t *face_vel,
459 cs_real_t *mass_flux)
460 {
461 if (mass_flux == NULL)
462 return;
463
464 assert(face_vel != NULL);
465 assert(nsp->space_scheme == CS_SPACE_SCHEME_CDOFB);
466 assert(cs_property_is_uniform(nsp->mass_density));
467 assert(nsp->mass_density->n_definitions == 1);
468
469 const cs_real_t rho_val = nsp->mass_density->ref_value;
470
471 /* Define the mass flux */
472
473 # pragma omp parallel for if (quant->n_faces > CS_THR_MIN)
474 for (cs_lnum_t f_id = 0; f_id < quant->n_faces; f_id++) {
475
476 const cs_real_t *fq = cs_quant_get_face_vector_area(f_id, quant);
477 mass_flux[f_id] = rho_val*cs_math_3_dot_product(face_vel + 3*f_id, fq);
478
479 }
480 }
481
482 /*----------------------------------------------------------------------------*/
483 /*!
484 * \brief Compute the divergence in a cell of a vector-valued array defined at
485 * faces (values are defined both at interior and border faces).
486 * Variant based on the usage of \ref cs_cdo_quantities_t structure.
487 *
488 * \param[in] c_id cell id
489 * \param[in] quant pointer to a \ref cs_cdo_quantities_t
490 * \param[in] c2f pointer to cell-to-face \ref cs_adjacency_t
491 * \param[in] f_vals values of the face DoFs
492 *
493 * \return the divergence for the corresponding cell
494 */
495 /*----------------------------------------------------------------------------*/
496
497 double
cs_cdofb_navsto_cell_divergence(const cs_lnum_t c_id,const cs_cdo_quantities_t * quant,const cs_adjacency_t * c2f,const cs_real_t * f_vals)498 cs_cdofb_navsto_cell_divergence(const cs_lnum_t c_id,
499 const cs_cdo_quantities_t *quant,
500 const cs_adjacency_t *c2f,
501 const cs_real_t *f_vals)
502 {
503 const cs_lnum_t thd = 3 * quant->n_i_faces;
504
505 double div = 0.0;
506 for (cs_lnum_t f = c2f->idx[c_id]; f < c2f->idx[c_id+1]; f++) {
507
508 const cs_lnum_t shift = 3*c2f->ids[f];
509 const cs_real_t *_val = f_vals + shift;
510 const cs_real_t *fnorm = (shift < thd) ?
511 quant->i_face_normal + shift : quant->b_face_normal + (shift - thd);
512
513 div += c2f->sgn[f]*cs_math_3_dot_product(_val, fnorm);
514
515 } /* Loop on cell faces */
516
517 div /= quant->cell_vol[c_id];
518
519 return div;
520 }
521
522 /*----------------------------------------------------------------------------*/
523 /*!
524 * \brief Compute an estimation of the pressure at faces
525 *
526 * \param[in] mesh pointer to a cs_mesh_t structure
527 * \param[in] connect pointer to a cs_cdo_connect_t structure
528 * \param[in] quant pointer to a cs_cdo_quantities_t structure
529 * \param[in] time_step pointer to a cs_time_step_t structure
530 * \param[in] nsp pointer to a \ref cs_navsto_param_t struct.
531 * \param[in] p_cell value of the pressure inside each cell
532 * \param[in, out] p_face value of the pressure at each face
533 */
534 /*----------------------------------------------------------------------------*/
535
536 void
cs_cdofb_navsto_compute_face_pressure(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts,const cs_navsto_param_t * nsp,const cs_real_t * p_cell,cs_real_t * p_face)537 cs_cdofb_navsto_compute_face_pressure(const cs_mesh_t *mesh,
538 const cs_cdo_connect_t *connect,
539 const cs_cdo_quantities_t *quant,
540 const cs_time_step_t *ts,
541 const cs_navsto_param_t *nsp,
542 const cs_real_t *p_cell,
543 cs_real_t *p_face)
544 {
545 /* Interior faces */
546
547 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++) {
548
549 cs_lnum_t iid = mesh->i_face_cells[i][0];
550 cs_lnum_t jid = mesh->i_face_cells[i][1];
551
552 p_face[i] = 0.5*(p_cell[iid] + p_cell[jid]);
553
554 }
555
556 /* Border faces
557 * Use the knowledge from the BCs if possible, otherwise assume a homogeneous
558 * Neumann (i.e. p_face = p_cell). */
559
560 cs_real_t *p_bface = p_face + mesh->n_i_faces;
561
562 for (cs_lnum_t i = 0; i < mesh->n_b_faces; i++)
563 p_bface[i] = p_cell[mesh->b_face_cells[i]];
564
565 for (int def_id = 0; def_id < nsp->n_pressure_bc_defs; def_id++) {
566
567 cs_xdef_t *pbc_def = nsp->pressure_bc_defs[def_id];
568 const cs_zone_t *z = cs_boundary_zone_by_id(pbc_def->z_id);
569
570 assert(pbc_def->meta & CS_CDO_BC_DIRICHLET);
571
572 switch(pbc_def->type) {
573
574 case CS_XDEF_BY_VALUE:
575 cs_xdef_eval_scalar_by_val(z->n_elts, z->elt_ids,
576 false, /* dense output */
577 mesh,
578 connect,
579 quant,
580 ts->t_cur,
581 pbc_def->context,
582 p_bface);
583 break;
584
585 default:
586 bft_error(__FILE__, __LINE__, 0,
587 _(" %s: Type of definition not handle.\n"
588 " Stop computing the pressure BC value.\n"), __func__);
589 break;
590
591 } /* def->type */
592
593 } /* Loop on pressure BCs */
594 }
595
596 /*----------------------------------------------------------------------------*/
597 /*!
598 * \brief Add the grad-div part to the local matrix (i.e. for the current
599 * cell)
600 *
601 * \param[in] n_fc local number of faces for the current cell
602 * \param[in] zeta scalar coefficient for the grad-div operator
603 * \param[in] div divergence
604 * \param[in, out] mat local system matrix to update
605 */
606 /*----------------------------------------------------------------------------*/
607
608 void
cs_cdofb_navsto_add_grad_div(short int n_fc,const cs_real_t zeta,const cs_real_t div[],cs_sdm_t * mat)609 cs_cdofb_navsto_add_grad_div(short int n_fc,
610 const cs_real_t zeta,
611 const cs_real_t div[],
612 cs_sdm_t *mat)
613 {
614 cs_sdm_t *b = NULL;
615
616 /* Avoid dealing with cell DoFs which are not impacted */
617 for (short int bi = 0; bi < n_fc; bi++) {
618
619 const cs_real_t *divi = div + 3*bi;
620 const cs_real_t zt_di[3] = {zeta*divi[0], zeta*divi[1], zeta*divi[2]};
621
622 /* Begin with the diagonal block */
623 b = cs_sdm_get_block(mat, bi, bi);
624 assert(b->n_rows == b->n_cols && b->n_rows == 3);
625 for (short int l = 0; l < 3; l++) {
626 cs_real_t *m_l = b->val + 3*l;
627 for (short int m = 0; m < 3; m++)
628 m_l[m] += zt_di[l] * divi[m];
629 }
630
631 /* Continue with the extra-diag. blocks */
632 for (short int bj = bi+1; bj < n_fc; bj++) {
633
634 b = cs_sdm_get_block(mat, bi, bj);
635 assert(b->n_rows == b->n_cols && b->n_rows == 3);
636 cs_real_t *mij = b->val;
637 b = cs_sdm_get_block(mat, bj, bi);
638 assert(b->n_rows == b->n_cols && b->n_rows == 3);
639 cs_real_t *mji = b->val;
640
641 const cs_real_t *divj = div + 3*bj;
642
643 for (short int l = 0; l < 3; l++) {
644
645 /* Diagonal: 3*l+l = 4*l */
646 const cs_real_t gd_coef_ll = zt_di[l]*divj[l];
647 mij[4*l] += gd_coef_ll;
648 mji[4*l] += gd_coef_ll;
649
650 /* Extra-diagonal: Use the symmetry of the grad-div */
651 for (short int m = l+1; m < 3; m++) {
652 const short int lm = 3*l+m, ml = 3*m+l;
653 const cs_real_t gd_coef_lm = zt_di[l]*divj[m];
654 mij[lm] += gd_coef_lm;
655 mji[ml] += gd_coef_lm;
656 const cs_real_t gd_coef_ml = zt_di[m]*divj[l];
657 mij[ml] += gd_coef_ml;
658 mji[lm] += gd_coef_ml;
659 }
660 }
661
662 } /* Loop on column blocks: bj */
663 } /* Loop on row blocks: bi */
664
665 }
666
667 /*----------------------------------------------------------------------------*/
668 /*!
669 * \brief Initialize the pressure values
670 *
671 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
672 * \param[in] quant pointer to a \ref cs_cdo_quantities_t structure
673 * \param[in] ts pointer to a \ref cs_time_step_t structure
674 * \param[in, out] pr pointer to the pressure \ref cs_field_t structure
675 */
676 /*----------------------------------------------------------------------------*/
677
678 void
cs_cdofb_navsto_init_pressure(const cs_navsto_param_t * nsp,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts,cs_field_t * pr)679 cs_cdofb_navsto_init_pressure(const cs_navsto_param_t *nsp,
680 const cs_cdo_quantities_t *quant,
681 const cs_time_step_t *ts,
682 cs_field_t *pr)
683 {
684 /* Sanity checks */
685 assert(nsp != NULL);
686
687 /* Initial conditions for the pressure */
688 if (nsp->n_pressure_ic_defs == 0)
689 return; /* Nothing to do */
690
691 assert(nsp->pressure_ic_defs != NULL);
692
693 const cs_real_t t_cur = ts->t_cur;
694 const cs_flag_t dof_flag = CS_FLAG_SCALAR | cs_flag_primal_cell;
695
696 cs_real_t *values = pr->val;
697
698 for (int def_id = 0; def_id < nsp->n_pressure_ic_defs; def_id++) {
699
700 /* Get and then set the definition of the initial condition */
701 cs_xdef_t *def = nsp->pressure_ic_defs[def_id];
702
703 /* Initialize face-based array */
704 switch (def->type) {
705
706 /* Evaluating the integrals: the averages will be taken care of at the
707 * end when ensuring zero-mean valuedness */
708 case CS_XDEF_BY_VALUE:
709 /* When constant mean-value or the value at cell center is the same */
710 cs_evaluate_potential_at_cells_by_value(def, values);
711 break;
712
713 case CS_XDEF_BY_ANALYTIC_FUNCTION:
714 {
715 const cs_param_dof_reduction_t red = nsp->dof_reduction_mode;
716
717 switch (red) {
718 case CS_PARAM_REDUCTION_DERHAM:
719 cs_xdef_set_quadrature(def, nsp->qtype);
720 cs_evaluate_density_by_analytic(dof_flag, def, t_cur, values);
721 break;
722 case CS_PARAM_REDUCTION_AVERAGE:
723 cs_xdef_set_quadrature(def, nsp->qtype);
724 cs_evaluate_average_on_cells_by_analytic(def, t_cur, values);
725 break;
726
727 default:
728 bft_error(__FILE__, __LINE__, 0,
729 _(" %s: Incompatible reduction for the field %s.\n"),
730 __func__, pr->name);
731
732 } /* Switch on possible reduction types */
733
734 }
735 break;
736
737 default:
738 bft_error(__FILE__, __LINE__, 0,
739 _(" %s: Incompatible way to initialize the field %s.\n"),
740 __func__, pr->name);
741 break;
742
743 } /* Switch on possible type of definition */
744
745 } /* Loop on definitions */
746
747 /* We should ensure that the mean value of the pressure is zero. Thus we
748 * compute it and subtract it from every value.
749 *
750 * NOTES:
751 * - It could be useful to stored this average somewhere
752 * - The procedure is not optimized (we can avoid setting the average if
753 * it's a value), but it is the only way to allow multiple definitions
754 * and definitions that do not cover all the domain. Moreover, we need
755 * information (e.g. cs_cdo_quantities_t) which we do not know here
756 */
757 cs_cdofb_navsto_rescale_pressure_to_ref(nsp, quant, values);
758 }
759
760 /*----------------------------------------------------------------------------*/
761 /*!
762 * \brief Initialize the pressure values when the pressure is defined at
763 * faces
764 *
765 * \param[in] nsp pointer to a \ref cs_navsto_param_t structure
766 * \param[in] connect pointer to a \ref cs_cdo_connect_t structure
767 * \param[in] ts pointer to a \ref cs_time_step_t structure
768 * \param[in, out] pr_f pointer to the pressure values at faces
769 */
770 /*----------------------------------------------------------------------------*/
771
772 void
cs_cdofb_navsto_init_face_pressure(const cs_navsto_param_t * nsp,const cs_cdo_connect_t * connect,const cs_time_step_t * ts,cs_real_t * pr_f)773 cs_cdofb_navsto_init_face_pressure(const cs_navsto_param_t *nsp,
774 const cs_cdo_connect_t *connect,
775 const cs_time_step_t *ts,
776 cs_real_t *pr_f)
777 {
778 /* Sanity checks */
779 assert(nsp != NULL && pr_f != NULL);
780
781 /* Initial conditions for the pressure */
782 if (nsp->n_pressure_ic_defs == 0)
783 return; /* Nothing to do */
784
785 assert(nsp->pressure_ic_defs != NULL);
786
787 cs_lnum_t *def2f_ids = (cs_lnum_t *)cs_equation_get_tmpbuf();
788 cs_lnum_t *def2f_idx = NULL;
789 BFT_MALLOC(def2f_idx, nsp->n_pressure_ic_defs + 1, cs_lnum_t);
790
791 cs_equation_sync_vol_def_at_faces(connect,
792 nsp->n_pressure_ic_defs,
793 nsp->pressure_ic_defs,
794 def2f_idx,
795 def2f_ids);
796
797 const cs_real_t t_cur = ts->t_cur;
798
799 for (int def_id = 0; def_id < nsp->n_pressure_ic_defs; def_id++) {
800
801 /* Get and then set the definition of the initial condition */
802 cs_xdef_t *def = nsp->pressure_ic_defs[def_id];
803 const cs_lnum_t n_f_selected = def2f_idx[def_id+1] - def2f_idx[def_id];
804 const cs_lnum_t *selected_lst = def2f_ids + def2f_idx[def_id];
805
806 /* Initialize face-based array */
807 switch (def->type) {
808
809 /* Evaluating the integrals: the averages will be taken care of at the
810 * end when ensuring zero-mean valuedness */
811 case CS_XDEF_BY_VALUE:
812 cs_evaluate_potential_at_faces_by_value(def,
813 n_f_selected,
814 selected_lst,
815 pr_f);
816 break;
817
818 case CS_XDEF_BY_ANALYTIC_FUNCTION:
819 {
820 const cs_param_dof_reduction_t red = nsp->dof_reduction_mode;
821
822 switch (red) {
823 case CS_PARAM_REDUCTION_DERHAM:
824 cs_evaluate_potential_at_faces_by_analytic(def,
825 t_cur,
826 n_f_selected,
827 selected_lst,
828 pr_f);
829 break;
830
831 case CS_PARAM_REDUCTION_AVERAGE:
832 cs_xdef_set_quadrature(def, nsp->qtype);
833 cs_evaluate_average_on_faces_by_analytic(def,
834 t_cur,
835 n_f_selected,
836 selected_lst,
837 pr_f);
838 break;
839
840 default:
841 bft_error(__FILE__, __LINE__, 0,
842 _(" %s: Incompatible reduction for the pressure field\n"),
843 __func__);
844
845 } /* Switch on possible reduction types */
846
847 }
848 break;
849
850 default:
851 bft_error(__FILE__, __LINE__, 0,
852 _(" %s: Incompatible way to initialize the pressure field.\n"),
853 __func__);
854 break;
855
856 } /* Switch on possible type of definition */
857
858 } /* Loop on definitions */
859
860 BFT_FREE(def2f_idx);
861 }
862
863 /*----------------------------------------------------------------------------*/
864 /*!
865 * \brief Update the pressure field in order to get a field with a mean-value
866 * equal to the reference value
867 *
868 * \param[in] nsp pointer to a cs_navsto_param_t structure
869 * \param[in] quant pointer to a cs_cdo_quantities_t structure
870 * \param[in, out] values pressure field values
871 */
872 /*----------------------------------------------------------------------------*/
873
874 void
cs_cdofb_navsto_rescale_pressure_to_ref(const cs_navsto_param_t * nsp,const cs_cdo_quantities_t * quant,cs_real_t values[])875 cs_cdofb_navsto_rescale_pressure_to_ref(const cs_navsto_param_t *nsp,
876 const cs_cdo_quantities_t *quant,
877 cs_real_t values[])
878 {
879 const cs_lnum_t n_cells = quant->n_cells;
880
881 /*
882 * The algorithm used for summing is l3superblock60, based on the article:
883 * "Reducing Floating Point Error in Dot Product Using the Superblock Family
884 * of Algorithms" by Anthony M. Castaldo, R. Clint Whaley, and Anthony
885 * T. Chronopoulos, SIAM J. SCI. COMPUT., Vol. 31, No. 2, pp. 1156--1174
886 * 2008 Society for Industrial and Applied Mathematics
887 */
888
889 cs_real_t intgr = cs_weighted_sum(n_cells, quant->cell_vol, values);
890
891 cs_parall_sum(1, CS_REAL_TYPE, &intgr);
892
893 assert(quant->vol_tot > 0.);
894 const cs_real_t g_avg = intgr / quant->vol_tot;
895 const cs_real_t p_shift = nsp->reference_pressure - g_avg;
896
897 # pragma omp parallel for if (n_cells > CS_THR_MIN)
898 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++)
899 values[c_id] = values[c_id] + p_shift;
900 }
901
902 /*----------------------------------------------------------------------------*/
903 /*!
904 * \brief Update the pressure field in order to get a field with a zero-mean
905 * average
906 *
907 * \param[in] quant pointer to a cs_cdo_quantities_t structure
908 * \param[in, out] values pressure field values
909 */
910 /*----------------------------------------------------------------------------*/
911
912 void
cs_cdofb_navsto_set_zero_mean_pressure(const cs_cdo_quantities_t * quant,cs_real_t values[])913 cs_cdofb_navsto_set_zero_mean_pressure(const cs_cdo_quantities_t *quant,
914 cs_real_t values[])
915 {
916 /* We should ensure that the mean of the pressure is zero. Thus we compute
917 * it and subtract it from every value. */
918 /* NOTES:
919 * - It could be useful to stored this average somewhere
920 * - The procedure is not optimized (we can avoid setting the average if
921 * it's a value), but it is the only way to allow multiple definitions
922 * and definitions that do not cover all the domain. */
923
924 const cs_lnum_t n_cells = quant->n_cells;
925
926 /*
927 * The algorithm used for summing is l3superblock60, based on the article:
928 * "Reducing Floating Point Error in Dot Product Using the Superblock Family
929 * of Algorithms" by Anthony M. Castaldo, R. Clint Whaley, and Anthony
930 * T. Chronopoulos, SIAM J. SCI. COMPUT., Vol. 31, No. 2, pp. 1156--1174
931 * 2008 Society for Industrial and Applied Mathematics
932 */
933
934 cs_real_t intgr = cs_weighted_sum(n_cells, quant->cell_vol, values);
935
936 cs_parall_sum(1, CS_REAL_TYPE, &intgr);
937
938 assert(quant->vol_tot > 0.);
939 const cs_real_t g_avg = intgr / quant->vol_tot;
940
941 # pragma omp parallel for if (n_cells > CS_THR_MIN)
942 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++)
943 values[c_id] = values[c_id] - g_avg;
944 }
945
946 /*----------------------------------------------------------------------------*/
947 /*!
948 * \brief Perform extra-operation related to Fb schemes when solving
949 * Navier-Stokes. Computation of the following quantities according to
950 * post-processing flags beeing activated.
951 * - The mass flux accross the boundaries.
952 * - The global mass in the computational domain
953 * - The norm of the velocity divergence
954 * - the cellwise mass flux balance
955 * - the kinetic energy
956 * - the velocity gradient
957 * - the pressure gradient
958 * - the vorticity
959 * - the helicity
960 * - the enstrophy
961 * - the stream function
962 *
963 * \param[in] nsp pointer to a \ref cs_navsto_param_t struct.
964 * \param[in] mesh pointer to a cs_mesh_t structure
965 * \param[in] quant pointer to a \ref cs_cdo_quantities_t struct.
966 * \param[in] connect pointer to a \ref cs_cdo_connect_t struct.
967 * \param[in] ts pointer to a \ref cs_time_step_t struct.
968 * \param[in,out] time_plotter pointer to a \ref cs_time_plot_t struct.
969 * \param[in] adv_field pointer to a \ref cs_adv_field_t struct.
970 * \param[in] mass_flux scalar-valued mass flux for each face
971 * \param[in] p_cell scalar-valued pressure in each cell
972 * \param[in] u_cell vector-valued velocity in each cell
973 * \param[in] u_face vector-valued velocity on each face
974 */
975 /*----------------------------------------------------------------------------*/
976
977 void
cs_cdofb_navsto_extra_op(const cs_navsto_param_t * nsp,const cs_mesh_t * mesh,const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect,const cs_time_step_t * ts,cs_time_plot_t * time_plotter,const cs_adv_field_t * adv_field,const cs_real_t * mass_flux,const cs_real_t * p_cell,const cs_real_t * u_cell,const cs_real_t * u_face)978 cs_cdofb_navsto_extra_op(const cs_navsto_param_t *nsp,
979 const cs_mesh_t *mesh,
980 const cs_cdo_quantities_t *quant,
981 const cs_cdo_connect_t *connect,
982 const cs_time_step_t *ts,
983 cs_time_plot_t *time_plotter,
984 const cs_adv_field_t *adv_field,
985 const cs_real_t *mass_flux,
986 const cs_real_t *p_cell,
987 const cs_real_t *u_cell,
988 const cs_real_t *u_face)
989 {
990 CS_UNUSED(adv_field);
991
992 const cs_boundary_t *boundaries = nsp->boundaries;
993 const cs_real_t *bmass_flux = mass_flux + quant->n_i_faces;
994
995 /* 1. Compute for each boundary the integrated mass flux to perform mass
996 * balance
997 */
998
999 bool *belong_to_default = NULL;
1000 BFT_MALLOC(belong_to_default, quant->n_b_faces, bool);
1001 # pragma omp parallel for if (quant->n_b_faces > CS_THR_MIN)
1002 for (cs_lnum_t i = 0; i < quant->n_b_faces; i++)
1003 belong_to_default[i] = true;
1004
1005 cs_real_t *boundary_fluxes = NULL;
1006 BFT_MALLOC(boundary_fluxes, boundaries->n_boundaries + 1, cs_real_t);
1007 memset(boundary_fluxes, 0, (boundaries->n_boundaries + 1)*sizeof(cs_real_t));
1008
1009 for (int b_id = 0; b_id < boundaries->n_boundaries; b_id++) {
1010
1011 const cs_zone_t *z = cs_boundary_zone_by_id(boundaries->zone_ids[b_id]);
1012
1013 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
1014 const cs_lnum_t bf_id = z->elt_ids[i];
1015 belong_to_default[bf_id] = false;
1016 boundary_fluxes[b_id] += bmass_flux[bf_id];
1017 }
1018
1019 } /* Loop on domain boundaries */
1020
1021 /* Update the flux through the default boundary */
1022
1023 cs_lnum_t default_case_count = 0;
1024 for (cs_lnum_t i = 0; i < quant->n_b_faces; i++) {
1025 if (belong_to_default[i]) {
1026 default_case_count += 1;
1027 boundary_fluxes[boundaries->n_boundaries] += bmass_flux[i];
1028 }
1029 }
1030
1031 /* Parallel synchronization if needed */
1032
1033 cs_parall_sum(boundaries->n_boundaries + 1, CS_REAL_TYPE, boundary_fluxes);
1034 cs_parall_counter_max(&default_case_count, 1);
1035
1036 /* Output result */
1037
1038 cs_log_printf(CS_LOG_DEFAULT,
1039 "\n- Balance of the mass flux across the boundaries:\n");
1040
1041 char descr[32];
1042 for (int b_id = 0; b_id < boundaries->n_boundaries; b_id++) {
1043
1044 const cs_zone_t *z = cs_boundary_zone_by_id(boundaries->zone_ids[b_id]);
1045
1046 cs_boundary_get_type_descr(boundaries, boundaries->types[b_id], 32, descr);
1047
1048 cs_log_printf(CS_LOG_DEFAULT, "b %-32s | %-32s |% -8.6e\n",
1049 descr, z->name, boundary_fluxes[b_id]);
1050
1051 } /* Loop on boundaries */
1052
1053 /* Default boundary (if something to do) */
1054
1055 if (default_case_count > 0) {
1056
1057 cs_boundary_get_type_descr(boundaries, boundaries->default_type, 32, descr);
1058 cs_log_printf(CS_LOG_DEFAULT, "b %-32s | %-32s |% -8.6e\n",
1059 descr, "default boundary",
1060 boundary_fluxes[boundaries->n_boundaries]);
1061
1062 }
1063
1064 /* Free temporary buffers */
1065
1066 BFT_FREE(belong_to_default);
1067 BFT_FREE(boundary_fluxes);
1068
1069 /* Predefined post-processing */
1070 /* ========================== */
1071
1072 /* There are five values if all flags are activated for the monitoring plot */
1073
1074 int n_cols = 0;
1075 cs_real_t col_vals[5] = {0, 0, 0, 0, 0};
1076
1077 if (nsp->post_flag & CS_NAVSTO_POST_VELOCITY_DIVERGENCE) {
1078
1079 double div_norm2 = 0.;
1080 cs_field_t *vel_div = cs_field_by_name("velocity_divergence");
1081 assert(vel_div != NULL);
1082
1083 if (nsp->coupling != CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY)
1084 cs_field_current_to_previous(vel_div);
1085
1086 /* Only the face velocity is used */
1087
1088 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN) \
1089 reduction(+:div_norm2)
1090 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1091
1092 double div_c = cs_cdofb_navsto_cell_divergence(c_id,
1093 quant,
1094 connect->c2f,
1095 u_face);
1096
1097 vel_div->val[c_id] = div_c;
1098 div_norm2 += quant->cell_vol[c_id] * div_c * div_c;
1099
1100 } /* Loop on cells */
1101
1102 cs_parall_sum(1, CS_DOUBLE, &div_norm2);
1103 col_vals[n_cols++] = sqrt(div_norm2);
1104
1105 } /* Velocity divergence */
1106
1107 if (nsp->post_flag & CS_NAVSTO_POST_MASS_DENSITY) {
1108
1109 double mass_integral = 0.;
1110 cs_field_t *rho = cs_field_by_name("mass_density");
1111 assert(rho != NULL);
1112
1113 cs_field_current_to_previous(rho);
1114
1115 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN) \
1116 reduction(+:mass_integral)
1117 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1118
1119 double boussi_coef = 1;
1120
1121 for (int i = 0; i < nsp->n_boussinesq_terms; i++) {
1122
1123 cs_navsto_param_boussinesq_t *bp = nsp->boussinesq_param + i;
1124 boussi_coef += -bp->beta*(bp->var[c_id] - bp->var0);
1125
1126 } /* Loop on Boussinesq terms */
1127
1128 double rho_c = nsp->mass_density->ref_value * boussi_coef;
1129 rho->val[c_id] = rho_c;
1130 mass_integral += quant->cell_vol[c_id] * rho_c;
1131
1132 } /* Loop on cells */
1133
1134 cs_parall_sum(1, CS_DOUBLE, &mass_integral);
1135 col_vals[n_cols++] = mass_integral;
1136
1137 } /* Mass density */
1138
1139 if (nsp->post_flag & CS_NAVSTO_POST_CELL_MASS_FLUX_BALANCE) {
1140
1141 cs_field_t *mf_balance = cs_field_by_name("mass_flux_balance");
1142 assert(mf_balance != NULL);
1143
1144 cs_field_current_to_previous(mf_balance);
1145
1146 const cs_adjacency_t *c2f = connect->c2f;
1147
1148 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
1149 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1150
1151 double balance = 0.;
1152 for (cs_lnum_t j = c2f->idx[c_id]; j < c2f->idx[c_id+1]; j++) {
1153 balance += c2f->sgn[j]*mass_flux[c2f->ids[j]];
1154
1155 }
1156 mf_balance->val[c_id] = balance;
1157
1158 } /* Loop on cells */
1159
1160 } /* Cell mass flux balance */
1161
1162 if (nsp->post_flag & CS_NAVSTO_POST_PRESSURE_GRADIENT) {
1163
1164 cs_field_t *pr_grd = cs_field_by_name("pressure_gradient");
1165 assert(pr_grd != NULL);
1166
1167 cs_field_current_to_previous(pr_grd);
1168
1169 /* Compute a face pressure */
1170
1171 cs_real_t *p_face = NULL;
1172 BFT_MALLOC(p_face, quant->n_faces, cs_real_t);
1173
1174 cs_cdofb_navsto_compute_face_pressure(mesh,
1175 connect,
1176 quant,
1177 ts,
1178 nsp,
1179 p_cell,
1180 p_face);
1181
1182 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
1183 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++)
1184 cs_reco_grad_cell_from_fb_dofs(c_id, connect, quant,
1185 p_cell, p_face,
1186 pr_grd->val + 3*c_id);
1187
1188 BFT_FREE(p_face);
1189
1190 } /* Pressure gradient */
1191
1192 if (nsp->post_flag & CS_NAVSTO_POST_KINETIC_ENERGY) {
1193
1194 double k_integral = 0.;
1195 cs_field_t *kinetic_energy = cs_field_by_name("kinetic_energy");
1196 assert(kinetic_energy != NULL);
1197
1198 cs_field_current_to_previous(kinetic_energy);
1199
1200 if (cs_property_is_uniform(nsp->mass_density)) {
1201
1202 /* This can be any cell but one assumes that there is at least one cell by
1203 MPI rank */
1204
1205 const cs_real_t rho = cs_property_get_cell_value(0, /* cell_id */
1206 ts->t_cur,
1207 nsp->mass_density);
1208
1209 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN) \
1210 reduction(+:k_integral)
1211 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1212
1213 double kc = 0.5*rho*cs_math_3_square_norm(u_cell + 3*c_id);
1214
1215 kinetic_energy->val[c_id] = kc;
1216 k_integral += quant->cell_vol[c_id]*kc;
1217
1218 }
1219
1220 }
1221 else { /* Mass density is not uniform in space */
1222
1223 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN) \
1224 reduction(+:k_integral)
1225 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1226
1227 cs_real_t rho_c = cs_property_get_cell_value(c_id,
1228 ts->t_cur,
1229 nsp->mass_density);
1230 double kc = 0.5*rho_c*cs_math_3_square_norm(u_cell + 3*c_id);
1231
1232 kinetic_energy->val[c_id] = kc;
1233 k_integral += quant->cell_vol[c_id] * kc;
1234
1235 }
1236
1237 }
1238
1239 cs_parall_sum(1, CS_DOUBLE, &k_integral); /* Sync. parallel computations */
1240 col_vals[n_cols++] = k_integral;
1241
1242 } /* Kinetic energy */
1243
1244 if (nsp->post_flag & CS_NAVSTO_POST_VELOCITY_GRADIENT) {
1245
1246 cs_field_t *velocity_gradient = cs_field_by_name("velocity_gradient");
1247 assert(velocity_gradient != NULL);
1248
1249 cs_field_current_to_previous(velocity_gradient);
1250
1251 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
1252 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++)
1253 cs_reco_grad_33_cell_from_fb_dofs(c_id, connect, quant,
1254 u_cell, u_face,
1255 velocity_gradient->val + 9*c_id);
1256
1257 } /* Velocity gradient */
1258
1259 cs_flag_t mask_velgrd[3] = { CS_NAVSTO_POST_VORTICITY,
1260 CS_NAVSTO_POST_HELICITY,
1261 CS_NAVSTO_POST_ENSTROPHY };
1262
1263 if (cs_flag_at_least(nsp->post_flag, 3, mask_velgrd)) {
1264
1265 cs_field_t *vorticity = cs_field_by_name("vorticity");
1266 assert(vorticity != NULL);
1267 cs_field_current_to_previous(vorticity);
1268
1269 double e_integral = 0.;
1270 cs_field_t *enstrophy = cs_field_by_name_try("enstrophy");
1271 if (nsp->post_flag & CS_NAVSTO_POST_ENSTROPHY)
1272 cs_field_current_to_previous(enstrophy);
1273
1274 double h_integral = 0.;
1275 cs_field_t *helicity = cs_field_by_name_try("helicity");
1276 if (nsp->post_flag & CS_NAVSTO_POST_HELICITY)
1277 cs_field_current_to_previous(helicity);
1278
1279 cs_field_t *velocity_gradient = cs_field_by_name_try("velocity_gradient");
1280
1281 if (velocity_gradient == NULL) {
1282
1283 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN) \
1284 reduction(+:e_integral) reduction(+:h_integral)
1285 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1286
1287 cs_real_t grd_uc[9];
1288
1289 /* Compute the velocity gradient */
1290
1291 cs_reco_grad_33_cell_from_fb_dofs(c_id, connect, quant,
1292 u_cell, u_face, grd_uc);
1293
1294 /* Compute the cell vorticity */
1295
1296 cs_real_t *w = vorticity->val + 3*c_id;
1297 w[0] = grd_uc[7] - grd_uc[5];
1298 w[1] = grd_uc[2] - grd_uc[6];
1299 w[2] = grd_uc[3] - grd_uc[1];
1300
1301 if (nsp->post_flag & CS_NAVSTO_POST_ENSTROPHY) {
1302
1303 double ec = cs_math_3_square_norm(w);
1304 enstrophy->val[c_id] = ec;
1305 e_integral += quant->cell_vol[c_id] * ec;
1306
1307 }
1308
1309 if (nsp->post_flag & CS_NAVSTO_POST_HELICITY) {
1310
1311 double hc = cs_math_3_dot_product(u_cell + 3*c_id, w);
1312 helicity->val[c_id] = hc;
1313 h_integral += quant->cell_vol[c_id] * hc;
1314
1315 }
1316
1317 } /* Loop on cells */
1318
1319 }
1320 else {
1321
1322 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN) \
1323 reduction(+:e_integral) reduction(+:h_integral)
1324 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1325
1326 cs_real_t *grd_uc = velocity_gradient->val + 9*c_id;
1327
1328 /* Compute the cell vorticity */
1329
1330 cs_real_t *w = vorticity->val + 3*c_id;
1331 w[0] = grd_uc[7] - grd_uc[5];
1332 w[1] = grd_uc[2] - grd_uc[6];
1333 w[2] = grd_uc[3] - grd_uc[1];
1334
1335 if (nsp->post_flag & CS_NAVSTO_POST_ENSTROPHY) {
1336
1337 double ec = cs_math_3_square_norm(w);
1338 enstrophy->val[c_id] = ec;
1339 e_integral += quant->cell_vol[c_id] * ec;
1340
1341 }
1342
1343 if (nsp->post_flag & CS_NAVSTO_POST_HELICITY) {
1344
1345 double hc = cs_math_3_dot_product(u_cell + 3*c_id, w);
1346 helicity->val[c_id] = hc;
1347 h_integral += quant->cell_vol[c_id] * hc;
1348
1349 }
1350
1351 } /* Loop on cells */
1352
1353 } /* velocity gradient has been computed previously */
1354
1355 if (nsp->post_flag & CS_NAVSTO_POST_ENSTROPHY) {
1356
1357 cs_parall_sum(1, CS_DOUBLE, &e_integral);
1358 col_vals[n_cols++] = e_integral;
1359
1360 }
1361
1362 if (nsp->post_flag & CS_NAVSTO_POST_HELICITY) {
1363
1364 cs_parall_sum(1, CS_DOUBLE, &h_integral);
1365 col_vals[n_cols++] = h_integral;
1366
1367 }
1368
1369 } /* vorticity, helicity or enstrophy computations */
1370
1371 if (cs_glob_rank_id < 1 && time_plotter != NULL)
1372 cs_time_plot_vals_write(time_plotter,
1373 ts->nt_cur,
1374 ts->t_cur,
1375 n_cols,
1376 col_vals);
1377
1378 /* Stream function */
1379 /* --------------- */
1380
1381 if (nsp->post_flag & CS_NAVSTO_POST_STREAM_FUNCTION) {
1382
1383 cs_equation_t *eq = cs_equation_by_name(CS_NAVSTO_STREAM_EQNAME);
1384 assert(eq != NULL);
1385 cs_equation_solve_steady_state(mesh, eq);
1386
1387 cs_equation_param_t *eqp = cs_equation_get_param(eq);
1388 if (eqp->n_bc_defs == 0) {
1389
1390 /* Since this is an equation solved with only homogeneous Neumann BCs, one
1391 * substracts the mean value to get a unique solution */
1392
1393 cs_real_t mean_value;
1394 cs_equation_integrate_variable(connect, quant, eq, &mean_value);
1395 mean_value /= quant->vol_tot;
1396
1397 cs_real_t *psi_v = cs_equation_get_vertex_values(eq, false);
1398 for (cs_lnum_t i = 0; i < quant->n_vertices; i++)
1399 psi_v[i] -= mean_value;
1400
1401 } /* If homogeneous Neumann everywhere */
1402
1403 } /* Computation of the stream function is requested */
1404 }
1405
1406 /*----------------------------------------------------------------------------*/
1407 /*!
1408 * \brief Take into account a Dirichlet BCs on the three velocity components.
1409 * For instance for a velocity inlet boundary or a wall
1410 * Handle the velocity-block in the global algebraic system in case of
1411 * an algebraic technique.
1412 * This prototype matches the function pointer cs_cdo_apply_boundary_t
1413 *
1414 * \param[in] f face id in the cell mesh numbering
1415 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1416 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1417 * \param[in] pty pointer to a \ref cs_property_data_t structure
1418 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1419 * \param[in, out] csys structure storing the cellwise system
1420 */
1421 /*----------------------------------------------------------------------------*/
1422
1423 void
cs_cdofb_block_dirichlet_alge(short int f,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * pty,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1424 cs_cdofb_block_dirichlet_alge(short int f,
1425 const cs_equation_param_t *eqp,
1426 const cs_cell_mesh_t *cm,
1427 const cs_property_data_t *pty,
1428 cs_cell_builder_t *cb,
1429 cs_cell_sys_t *csys)
1430 {
1431 CS_UNUSED(eqp);
1432 CS_UNUSED(cm);
1433 CS_UNUSED(pty);
1434
1435 double *x_dir = cb->values;
1436 double *ax_dir = cb->values + 3;
1437 cs_sdm_t *m = csys->mat;
1438 cs_sdm_block_t *bd = m->block_desc;
1439 assert(bd != NULL);
1440 assert(bd->n_row_blocks == cm->n_fc || bd->n_row_blocks == cm->n_fc + 1);
1441
1442 /* Build x_dir */
1443 bool is_non_homogeneous = false; /* Assume homogeneous by default */
1444
1445 memset(cb->values, 0, 6*sizeof(double));
1446
1447 for (int k = 0; k < 3; k++) {
1448 if (csys->dof_flag[3*f+k] & CS_CDO_BC_DIRICHLET) {
1449 x_dir[k] = csys->dir_values[3*f+k];
1450 is_non_homogeneous = true;
1451 }
1452 }
1453
1454 if (is_non_homogeneous) {
1455
1456 for (int bi = 0; bi < bd->n_row_blocks; bi++) {
1457
1458 if (bi == f)
1459 continue;
1460
1461 cs_real_t *_rhs = csys->rhs + 3*bi;
1462 cs_sdm_t *mIF = cs_sdm_get_block(m, bi, f);
1463
1464 cs_sdm_square_matvec(mIF, x_dir, ax_dir);
1465 for (int k = 0; k < 3; k++) _rhs[k] -= ax_dir[k];
1466
1467 }
1468
1469 } /* Non-homogeneous Dirichlet BC */
1470
1471 /* Set RHS to the Dirichlet value for the related face */
1472 for (int k = 0; k < 3; k++)
1473 csys->rhs[3*f+k] = x_dir[k];
1474
1475 /* Second pass: Replace the Dirichlet block by a diagonal block and fill with
1476 * zero the remaining row and column */
1477 for (int bi = 0; bi < bd->n_row_blocks; bi++) {
1478
1479 if (bi != f) {
1480
1481 /* Reset block (I,F) which is a 3x3 block */
1482 cs_sdm_t *mIF = cs_sdm_get_block(m, bi, f);
1483 memset(mIF->val, 0, 9*sizeof(double));
1484
1485 }
1486 else { /* bi == f */
1487
1488 /* Reset block (I==F,J) which is a 3x3 block */
1489 for (int bj = 0; bj < bd->n_col_blocks; bj++) {
1490 cs_sdm_t *mFJ = cs_sdm_get_block(m, f, bj);
1491 memset(mFJ->val, 0, 9*sizeof(double));
1492 }
1493
1494 cs_sdm_t *mFF = cs_sdm_get_block(m, f, f);
1495 assert((mFF->n_cols == 3) && (mFF->n_rows == 3));
1496 for (int k = 0; k < 3; k++)
1497 mFF->val[4*k] = 1; /* 4 == mFF->n_rows + 1 */
1498
1499 }
1500
1501 } /* Block bi */
1502 }
1503
1504 /*----------------------------------------------------------------------------*/
1505 /*!
1506 * \brief Take into account a Dirichlet BCs on the three velocity components.
1507 * For instance for a velocity inlet boundary or a wall
1508 * Handle the velocity-block in the global algebraic system in case of
1509 * a penalization technique (with a large coefficient).
1510 * One assumes that static condensation has been performed and that
1511 * the velocity-block has size 3*n_fc
1512 * This prototype matches the function pointer cs_cdo_apply_boundary_t
1513 *
1514 * \param[in] f face id in the cell mesh numbering
1515 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1516 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1517 * \param[in] pty pointer to a \ref cs_property_data_t structure
1518 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1519 * \param[in, out] csys structure storing the cellwise system
1520 */
1521 /*----------------------------------------------------------------------------*/
1522
1523 void
cs_cdofb_block_dirichlet_pena(short int f,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * pty,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1524 cs_cdofb_block_dirichlet_pena(short int f,
1525 const cs_equation_param_t *eqp,
1526 const cs_cell_mesh_t *cm,
1527 const cs_property_data_t *pty,
1528 cs_cell_builder_t *cb,
1529 cs_cell_sys_t *csys)
1530 {
1531 CS_UNUSED(cb);
1532 CS_UNUSED(cm);
1533 CS_UNUSED(pty);
1534
1535 assert(csys != NULL);
1536
1537 cs_sdm_t *m = csys->mat;
1538 assert(m->block_desc != NULL);
1539 assert(m->block_desc->n_row_blocks == cm->n_fc);
1540
1541 const cs_flag_t *_flag = csys->dof_flag + 3*f;
1542 const cs_real_t *_dir_val = csys->dir_values + 3*f;
1543
1544 bool is_non_homogeneous = true;
1545 for (int k = 0; k < 3; k++) {
1546 if (_flag[k] & CS_CDO_BC_DIRICHLET)
1547 is_non_homogeneous = true;
1548 }
1549
1550 /* Penalize diagonal entry (and its rhs if needed) */
1551 cs_sdm_t *mFF = cs_sdm_get_block(m, f, f);
1552 assert((mFF->n_rows == 3) && (mFF->n_cols == 3));
1553
1554 if (is_non_homogeneous) {
1555
1556 cs_real_t *_rhs = csys->rhs + 3*f;
1557 for (int k = 0; k < 3; k++) {
1558 mFF->val[4*k] += eqp->strong_pena_bc_coeff; /* 4 == mFF->n_rows + 1 */
1559 _rhs[k] += _dir_val[k] * eqp->strong_pena_bc_coeff;
1560 }
1561
1562 }
1563 else {
1564
1565 for (int k = 0; k < 3; k++)
1566 mFF->val[4*k] += eqp->strong_pena_bc_coeff; /* 4 == mFF->n_rows + 1 */
1567
1568 } /* Homogeneous BC */
1569 }
1570
1571 /*----------------------------------------------------------------------------*/
1572 /*!
1573 * \brief Take into account a Dirichlet BCs on the three velocity components.
1574 * For instance for a velocity inlet boundary or a wall
1575 * Handle the velocity-block in the global algebraic system in case of
1576 * a weak penalization technique (Nitsche).
1577 * One assumes that static condensation has not been performed yet and
1578 * that the velocity-block has size 3*(n_fc + 1)
1579 * This prototype matches the function pointer cs_cdo_apply_boundary_t
1580 *
1581 * \param[in] fb face id in the cell mesh numbering
1582 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1583 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1584 * \param[in] pty pointer to a \ref cs_property_data_t structure
1585 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1586 * \param[in, out] csys structure storing the cellwise system
1587 */
1588 /*----------------------------------------------------------------------------*/
1589
1590 void
cs_cdofb_block_dirichlet_weak(short int fb,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * pty,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1591 cs_cdofb_block_dirichlet_weak(short int fb,
1592 const cs_equation_param_t *eqp,
1593 const cs_cell_mesh_t *cm,
1594 const cs_property_data_t *pty,
1595 cs_cell_builder_t *cb,
1596 cs_cell_sys_t *csys)
1597 {
1598 /* Sanity checks */
1599 assert(cm != NULL && cb != NULL && csys != NULL && pty != NULL);
1600 assert(pty->is_iso == true);
1601
1602 /* 0) Pre-compute the product between diffusion property and the
1603 face vector areas */
1604
1605 cs_real_3_t *kappa_f = cb->vectors;
1606 for (short int f = 0; f < cm->n_fc; f++) {
1607 const double coef = cm->face[f].meas*pty->value;
1608 for (short int k = 0; k < 3; k++)
1609 kappa_f[f][k] = coef * cm->face[f].unitv[k];
1610 }
1611
1612 /* 1) Build the bc_op matrix (scalar-valued version) */
1613
1614 /* Initialize the matrix related to the flux reconstruction operator */
1615 const short int n_dofs = cm->n_fc + 1; /* n_blocks or n_scalar_dofs */
1616 cs_sdm_t *bc_op = cb->loc;
1617 cs_sdm_square_init(n_dofs, bc_op);
1618
1619 /* Compute \int_f du/dn v and update the matrix */
1620 _normal_flux_reco(fb, eqp->diffusion_hodgep.coef, cm,
1621 (const cs_real_t (*)[3])kappa_f,
1622 bc_op);
1623
1624 /* 2) Update the bc_op matrix and the RHS with the Dirichlet values. */
1625
1626 /* coeff * \meas{f} / h_f */
1627 const cs_real_t pcoef = eqp->weak_pena_bc_coeff * sqrt(cm->face[fb].meas);
1628
1629 bc_op->val[fb*(n_dofs + 1)] += pcoef; /* Diagonal term */
1630
1631 for (short int k = 0; k < 3; k++)
1632 csys->rhs[3*fb + k] += pcoef * csys->dir_values[3*fb + k];
1633
1634 /* 3) Update the local system matrix */
1635
1636 for (int bi = 0; bi < n_dofs; bi++) { /* n_(scalar)_dofs == n_blocks */
1637 for (int bj = 0; bj < n_dofs; bj++) {
1638
1639 /* Retrieve the 3x3 matrix */
1640 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bi, bj);
1641 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
1642
1643 const cs_real_t _val = bc_op->val[n_dofs*bi + bj];
1644 /* Update diagonal terms only */
1645 bij->val[0] += _val;
1646 bij->val[4] += _val;
1647 bij->val[8] += _val;
1648
1649 }
1650 }
1651
1652 }
1653
1654 /*----------------------------------------------------------------------------*/
1655 /*!
1656 * \brief Take into account a Dirichlet BCs on the three velocity components.
1657 * For instance for a velocity inlet boundary or a wall
1658 * Handle the velocity-block in the global algebraic system in case of
1659 * a weak penalization technique (symmetrized Nitsche).
1660 * One assumes that static condensation has not been performed yet and
1661 * that the velocity-block has size 3*(n_fc + 1)
1662 * This prototype matches the function pointer cs_cdo_apply_boundary_t
1663 *
1664 * \param[in] fb face id in the cell mesh numbering
1665 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1666 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1667 * \param[in] pty pointer to a \ref cs_property_data_t structure
1668 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1669 * \param[in, out] csys structure storing the cellwise system
1670 */
1671 /*----------------------------------------------------------------------------*/
1672
1673 void
cs_cdofb_block_dirichlet_wsym(short int fb,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * pty,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1674 cs_cdofb_block_dirichlet_wsym(short int fb,
1675 const cs_equation_param_t *eqp,
1676 const cs_cell_mesh_t *cm,
1677 const cs_property_data_t *pty,
1678 cs_cell_builder_t *cb,
1679 cs_cell_sys_t *csys)
1680 {
1681 /* Sanity checks */
1682 assert(cm != NULL && cb != NULL && csys != NULL && pty != NULL);
1683 assert(cs_equation_param_has_diffusion(eqp));
1684 assert(pty->is_iso == true);
1685
1686 /* 0) Pre-compute the product between diffusion property and the
1687 face vector areas */
1688
1689 cs_real_3_t *kappa_f = cb->vectors;
1690 for (short int f = 0; f < cm->n_fc; f++) {
1691 const double coef = cm->face[f].meas*pty->value;
1692 for (int k = 0; k < 3; k++)
1693 kappa_f[f][k] = coef * cm->face[f].unitv[k];
1694 }
1695
1696 /* 1) Build the bc_op matrix (scalar-valued version) */
1697
1698 /* Initialize the matrix related to the flux reconstruction operator */
1699 const short int n_dofs = cm->n_fc + 1; /* n_blocks or n_scalar_dofs */
1700 cs_sdm_t *bc_op = cb->loc, *bc_op_t = cb->aux;
1701 cs_sdm_square_init(n_dofs, bc_op);
1702
1703 /* Compute \int_f du/dn v and update the matrix. Only the line associated to
1704 fb has non-zero values */
1705 _normal_flux_reco(fb, eqp->diffusion_hodgep.coef, cm,
1706 (const cs_real_t (*)[3])kappa_f,
1707 bc_op);
1708
1709 /* 2) Update the bc_op matrix and the RHS with the Dirichlet values */
1710
1711 /* Update bc_op = bc_op + transp and transp = transpose(bc_op)
1712 cb->loc plays the role of the flux operator */
1713 cs_sdm_square_add_transpose(bc_op, bc_op_t);
1714
1715 /* Update the RHS with bc_op_t * dir_fb */
1716 for (int k = 0; k < 3; k++) {
1717
1718 /* Only the fb column has non-zero values in bc_op_t
1719 One thus simplifies the matrix-vector operation */
1720 const cs_real_t dir_fb = csys->dir_values[3*fb+k];
1721 for (short int i = 0; i < n_dofs; i++)
1722 csys->rhs[3*i+k] += bc_op_t->val[i*n_dofs+fb] * dir_fb;
1723
1724 } /* Loop on components */
1725
1726 /* 3) Update the bc_op matrix and the RHS with the penalization */
1727
1728 /* coeff * \meas{f} / h_f \equiv coeff * h_f */
1729 const double pcoef = eqp->weak_pena_bc_coeff * sqrt(cm->face[fb].meas);
1730
1731 bc_op->val[fb*(n_dofs + 1)] += pcoef; /* Diagonal term */
1732
1733 for (short int k = 0; k < 3; k++)
1734 csys->rhs[3*fb + k] += pcoef * csys->dir_values[3*fb + k];
1735
1736 /* 4) Update the local system matrix */
1737
1738 for (int bi = 0; bi < n_dofs; bi++) { /* n_(scalar)_dofs == n_blocks */
1739 for (int bj = 0; bj < n_dofs; bj++) {
1740
1741 /* Retrieve the 3x3 matrix */
1742 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bi, bj);
1743 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
1744
1745 const cs_real_t _val = bc_op->val[n_dofs*bi + bj];
1746 /* Update diagonal terms only */
1747 bij->val[0] += _val;
1748 bij->val[4] += _val;
1749 bij->val[8] += _val;
1750
1751 }
1752 }
1753
1754 }
1755
1756 /*----------------------------------------------------------------------------*/
1757 /*!
1758 * \brief Take into account a boundary defined as 'symmetry' (treated as a
1759 * sliding BCs on the three velocity components.)
1760 * A weak penalization technique (symmetrized Nitsche) is used. One
1761 * assumes that static condensation has not been performed yet and that
1762 * the velocity-block has (n_fc + 1) blocks of size 3x3.
1763 * This prototype matches the function pointer cs_cdo_apply_boundary_t
1764 *
1765 * \param[in] fb face id in the cell mesh numbering
1766 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1767 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1768 * \param[in] pty pointer to a \ref cs_property_data_t structure
1769 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1770 * \param[in, out] csys structure storing the cellwise system
1771 */
1772 /*----------------------------------------------------------------------------*/
1773
1774 void
cs_cdofb_symmetry(short int fb,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * pty,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1775 cs_cdofb_symmetry(short int fb,
1776 const cs_equation_param_t *eqp,
1777 const cs_cell_mesh_t *cm,
1778 const cs_property_data_t *pty,
1779 cs_cell_builder_t *cb,
1780 cs_cell_sys_t *csys)
1781 {
1782 /* Sanity checks */
1783 assert(cm != NULL && cb != NULL && csys != NULL && pty != NULL);
1784 assert(pty->is_iso == true); /* if not the case something else TODO ? */
1785 assert(cs_equation_param_has_diffusion(eqp));
1786
1787 /* 0) Pre-compute the product between diffusion property and the
1788 face vector areas */
1789
1790 cs_real_3_t *kappa_f = cb->vectors;
1791 for (short int f = 0; f < cm->n_fc; f++) {
1792 const double coef = cm->face[f].meas*pty->value;
1793 for (int k = 0; k < 3; k++)
1794 kappa_f[f][k] = coef * cm->face[f].unitv[k];
1795 }
1796
1797 /* 1) Build the bc_op matrix (scalar-valued version) */
1798
1799 /* Initialize the matrix related this flux reconstruction operator */
1800 const short int n_dofs = cm->n_fc + 1; /* n_blocks or n_scalar_dofs */
1801 cs_sdm_t *bc_op = cb->aux;
1802 cs_sdm_square_init(n_dofs, bc_op);
1803
1804 /* Compute \int_f du/dn v and update the matrix. Only the line associated to
1805 fb has non-zero values */
1806 _normal_flux_reco(fb, eqp->diffusion_hodgep.coef, cm,
1807 (const cs_real_t (*)[3])kappa_f,
1808 bc_op);
1809
1810 /* 2) Update the bc_op matrix and nothing is added to the RHS since a sliding
1811 means homogeneous Dirichlet values on the normal component and hommogeneous
1812 Neumann on the tangential flux */
1813
1814 const cs_quant_t pfq = cm->face[fb];
1815 const cs_real_t *nf = pfq.unitv;
1816 const cs_real_33_t nf_nf = { {nf[0]*nf[0], nf[0]*nf[1], nf[0]*nf[2]},
1817 {nf[1]*nf[0], nf[1]*nf[1], nf[1]*nf[2]},
1818 {nf[2]*nf[0], nf[2]*nf[1], nf[2]*nf[2]} };
1819
1820 /* coeff * \meas{f} / h_f \equiv coeff * h_f */
1821 const double pcoef = eqp->weak_pena_bc_coeff * sqrt(pfq.meas);
1822
1823 /* Handle the diagonal block: Retrieve the 3x3 matrix */
1824 cs_sdm_t *bFF = cs_sdm_get_block(csys->mat, fb, fb);
1825 assert(bFF->n_rows == bFF->n_cols && bFF->n_rows == 3);
1826
1827 const cs_real_t _val = pcoef + 2*bc_op->val[fb*(n_dofs+1)];
1828 for (short int k = 0; k < 3; k++) {
1829 bFF->val[3*k ] += nf_nf[0][k] * _val;
1830 bFF->val[3*k+1] += nf_nf[1][k] * _val;
1831 bFF->val[3*k+2] += nf_nf[2][k] * _val;
1832 }
1833
1834 for (short int xj = 0; xj < n_dofs; xj++) {
1835
1836 if (xj == fb)
1837 continue;
1838
1839 /* It should be done both for face- and cell-defined DoFs */
1840 /* Retrieve the 3x3 matrix */
1841 cs_sdm_t *bFJ = cs_sdm_get_block(csys->mat, fb, xj);
1842 assert(bFJ->n_rows == bFJ->n_cols && bFJ->n_rows == 3);
1843 cs_sdm_t *bJF = cs_sdm_get_block(csys->mat, xj, fb);
1844 assert(bJF->n_rows == bJF->n_cols && bJF->n_rows == 3);
1845
1846 const cs_real_t op_fj = bc_op->val[n_dofs*fb + xj];
1847 const cs_real_t op_jf = bc_op->val[n_dofs*xj + fb];
1848 const cs_real_t _val_fj_jf = op_fj + op_jf;
1849
1850 for (int k = 0; k < 3; k++) {
1851
1852 bFJ->val[3*k ] += nf_nf[0][k] * _val_fj_jf;
1853 bFJ->val[3*k+1] += nf_nf[1][k] * _val_fj_jf;
1854 bFJ->val[3*k+2] += nf_nf[2][k] * _val_fj_jf;
1855
1856 /* nf_nf is symmetric */
1857 bJF->val[3*k ] += nf_nf[0][k] * _val_fj_jf;
1858 bJF->val[3*k+1] += nf_nf[1][k] * _val_fj_jf;
1859 bJF->val[3*k+2] += nf_nf[2][k] * _val_fj_jf;
1860
1861 }
1862
1863 } /* Loop on xj */
1864 }
1865
1866 /*----------------------------------------------------------------------------*/
1867 /*!
1868 * \brief Take into account a wall BCs by a weak enforcement using Nitsche
1869 * technique plus a symmetric treatment.
1870 * This prototype matches the function pointer cs_cdo_apply_boundary_t
1871 *
1872 * \param[in] fb face id in the cell mesh numbering
1873 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1874 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1875 * \param[in] pty pointer to a \ref cs_property_data_t structure
1876 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1877 * \param[in, out] csys structure storing the cellwise system
1878 */
1879 /*----------------------------------------------------------------------------*/
1880
1881 void
cs_cdofb_fixed_wall(short int fb,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * pty,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1882 cs_cdofb_fixed_wall(short int fb,
1883 const cs_equation_param_t *eqp,
1884 const cs_cell_mesh_t *cm,
1885 const cs_property_data_t *pty,
1886 cs_cell_builder_t *cb,
1887 cs_cell_sys_t *csys)
1888 {
1889 CS_UNUSED(cb);
1890 CS_UNUSED(pty);
1891
1892 assert(cm != NULL && csys != NULL); /* Sanity checks */
1893
1894 const cs_quant_t pfq = cm->face[fb];
1895 const cs_real_t *ni = pfq.unitv;
1896 const cs_real_t ni_ni[9] = { ni[0]*ni[0], ni[0]*ni[1], ni[0]*ni[2],
1897 ni[1]*ni[0], ni[1]*ni[1], ni[1]*ni[2],
1898 ni[2]*ni[0], ni[2]*ni[1], ni[2]*ni[2] };
1899
1900 /* coeff * \meas{fb} / h_fb \equiv coeff * h_fb */
1901 const cs_real_t pcoef = eqp->weak_pena_bc_coeff * sqrt(pfq.meas);
1902
1903 cs_sdm_t *bii = cs_sdm_get_block(csys->mat, fb, fb);
1904 assert(bii->n_rows == bii->n_cols && bii->n_rows == 3);
1905
1906 for (short int k = 0; k < 9; k++)
1907 bii->val[k] += pcoef * ni_ni[k];
1908 }
1909
1910 /*----------------------------------------------------------------------------*/
1911 /*!
1912 * \brief Test if one has to do one more non-linear iteration.
1913 * Test if performed on the relative norm on the increment between
1914 * two iterations
1915 *
1916 * \param[in] nl_algo_type type of non-linear algorithm
1917 * \param[in] pre_iterate previous state of the mass flux iterate
1918 * \param[in] cur_iterate current state of the mass flux iterate
1919 * \param[in, out] algo pointer to a cs_iter_algo_t structure
1920 *
1921 * \return the convergence state
1922 */
1923 /*----------------------------------------------------------------------------*/
1924
1925 cs_sles_convergence_state_t
cs_cdofb_navsto_nl_algo_cvg(cs_param_nl_algo_t nl_algo_type,const cs_real_t * pre_iterate,cs_real_t * cur_iterate,cs_iter_algo_t * algo)1926 cs_cdofb_navsto_nl_algo_cvg(cs_param_nl_algo_t nl_algo_type,
1927 const cs_real_t *pre_iterate,
1928 cs_real_t *cur_iterate,
1929 cs_iter_algo_t *algo)
1930 {
1931 assert(algo != NULL);
1932
1933 if (nl_algo_type == CS_PARAM_NL_ALGO_ANDERSON && algo->n_algo_iter > 0) {
1934
1935 cs_iter_algo_param_aa_t aap = cs_iter_algo_get_anderson_param(algo);
1936
1937 switch (aap.dp_type) {
1938
1939 case CS_PARAM_DOTPROD_EUCLIDEAN:
1940 cs_iter_algo_aa_update(algo,
1941 cur_iterate,
1942 pre_iterate,
1943 cs_cdo_blas_dotprod_face,
1944 cs_cdo_blas_square_norm_face);
1945 break;
1946
1947 case CS_PARAM_DOTPROD_CDO:
1948 cs_iter_algo_aa_update(algo,
1949 cur_iterate,
1950 pre_iterate,
1951 cs_cdo_blas_dotprod_pfsf,
1952 cs_cdo_blas_square_norm_pfsf);
1953 break;
1954
1955 default:
1956 bft_error(__FILE__, __LINE__, 0, "%s: Invalid case.\n", __func__);
1957
1958 }
1959
1960 } /* Anderson acceleration */
1961
1962 /* Update the residual values. Compute the norm of the difference between the
1963 two mass fluxes (the current one and the previous one) */
1964
1965 algo->prev_res = algo->res;
1966 algo->res = cs_cdo_blas_square_norm_pfsf_diff(pre_iterate, cur_iterate);
1967 assert(algo->res > -DBL_MIN);
1968 algo->res = sqrt(algo->res);
1969
1970 if (algo->n_algo_iter < 1) /* Store the first residual to detect a
1971 divergence */
1972 algo->res0 = algo->res;
1973
1974 /* Update the convergence members */
1975
1976 cs_iter_algo_update_cvg(algo);
1977
1978 if (algo->param.verbosity > 0) {
1979
1980 switch (nl_algo_type) {
1981
1982 case CS_PARAM_NL_ALGO_ANDERSON:
1983 if (algo->n_algo_iter == 1)
1984 _nl_algo_print_header("## Anderson");
1985 _nl_algo_print_entry("## Anderson", algo);
1986 break;
1987
1988 case CS_PARAM_NL_ALGO_PICARD:
1989 if (algo->n_algo_iter == 1)
1990 _nl_algo_print_header("## Picard");
1991 _nl_algo_print_entry("## Picard", algo);
1992 break;
1993
1994 default:
1995 if (algo->n_algo_iter == 1)
1996 _nl_algo_print_header("## ");
1997 _nl_algo_print_entry("## ", algo);
1998 break;
1999
2000 }
2001
2002 } /* verbosity > 0 */
2003
2004 return algo->cvg;
2005 }
2006
2007 /*----------------------------------------------------------------------------*/
2008 /*!
2009 * \brief Set the function pointer computing the source term in the momentum
2010 * equation related to the gravity effect (hydrostatic pressure or the
2011 * Boussinesq approximation)
2012 *
2013 * \param[in] nsp set of parameters for the Navier-Stokes system
2014 * \param[out] p_func way to compute the gravity effect
2015 */
2016 /*----------------------------------------------------------------------------*/
2017
2018 void
cs_cdofb_navsto_set_gravity_func(const cs_navsto_param_t * nsp,cs_cdofb_navsto_source_t ** p_func)2019 cs_cdofb_navsto_set_gravity_func(const cs_navsto_param_t *nsp,
2020 cs_cdofb_navsto_source_t **p_func)
2021 {
2022 if (nsp->model_flag & CS_NAVSTO_MODEL_BOUSSINESQ) {
2023
2024 switch (cs_cdofb_navsto_boussinesq_type) {
2025
2026 case CS_CDOFB_NAVSTO_BOUSSINESQ_FACE_DOF:
2027 *p_func = cs_cdofb_navsto_boussinesq_at_face;
2028 break;
2029
2030 case CS_CDOFB_NAVSTO_BOUSSINESQ_CELL_DOF:
2031 *p_func = cs_cdofb_navsto_boussinesq_at_cell;
2032 break;
2033
2034 default:
2035 bft_error(__FILE__, __LINE__, 0,
2036 "%s: Invalid type of algorithm to compute the Boussinesq"
2037 " approximation.\n", __func__);
2038 }
2039
2040 }
2041 else if (nsp->model_flag & CS_NAVSTO_MODEL_GRAVITY_EFFECTS)
2042 *p_func = cs_cdofb_navsto_gravity_term;
2043
2044 else
2045 *p_func = NULL;
2046 }
2047
2048 /*----------------------------------------------------------------------------*/
2049 /*!
2050 * \brief Take into account the gravity effects.
2051 * Compute and add the source term to the local RHS.
2052 * This is a special treatment since of face DoFs are involved
2053 * contrary to the standard case where only the cell DoFs is involved.
2054 *
2055 * \param[in] nsp set of parameters to handle the Navier-Stokes system
2056 * \param[in] cm pointer to a cs_cell_mesh_t structure
2057 * \param[in] nsb pointer to a builder structure for the NavSto system
2058 * \param[in, out] csys pointer to a cs_cell_sys_t structure
2059 */
2060 /*----------------------------------------------------------------------------*/
2061
2062 void
cs_cdofb_navsto_gravity_term(const cs_navsto_param_t * nsp,const cs_cell_mesh_t * cm,const cs_cdofb_navsto_builder_t * nsb,cs_cell_sys_t * csys)2063 cs_cdofb_navsto_gravity_term(const cs_navsto_param_t *nsp,
2064 const cs_cell_mesh_t *cm,
2065 const cs_cdofb_navsto_builder_t *nsb,
2066 cs_cell_sys_t *csys)
2067 {
2068 assert(nsp->model_flag & CS_NAVSTO_MODEL_GRAVITY_EFFECTS);
2069
2070 const cs_real_t *gravity_vector = nsp->phys_constants->gravity;
2071 const cs_real_t cell_coef = _dp3(gravity_vector, cm->xc);
2072
2073 for (int f = 0; f < cm->n_fc; f++) {
2074
2075 const cs_real_t *_div_f = nsb->div_op + 3*f;
2076 const cs_quant_t pfq = cm->face[f];
2077 const cs_real_t face_coef = _dp3(gravity_vector, pfq.center);
2078
2079 /* div_op is built such that _div_f[k] = -i_{f,c} * |f| * n_f[k] */
2080
2081 for (int k = 0; k < 3; k++)
2082 csys->rhs[3*f+k] += _div_f[k] * nsb->rho_c * (cell_coef - face_coef);
2083
2084 } /* Loop on cell faces */
2085
2086 }
2087
2088 /*----------------------------------------------------------------------------*/
2089 /*!
2090 * \brief Take into account the buoyancy force with the Boussinesq approx.
2091 * Compute and add the source term to the local RHS.
2092 * This is the standard case where the face DoFs are used for the
2093 * constant part rho0 . g[] and only the cell DoFs are involved for the
2094 * remaining part (the Boussinesq approximation).
2095 *
2096 * \param[in] nsp set of parameters to handle the Navier-Stokes system
2097 * \param[in] cm pointer to a cs_cell_mesh_t structure
2098 * \param[in] nsb pointer to a builder structure for the NavSto system
2099 * \param[in, out] csys pointer to a cs_cell_sys_t structure
2100 */
2101 /*----------------------------------------------------------------------------*/
2102
2103 void
cs_cdofb_navsto_boussinesq_at_cell(const cs_navsto_param_t * nsp,const cs_cell_mesh_t * cm,const cs_cdofb_navsto_builder_t * nsb,cs_cell_sys_t * csys)2104 cs_cdofb_navsto_boussinesq_at_cell(const cs_navsto_param_t *nsp,
2105 const cs_cell_mesh_t *cm,
2106 const cs_cdofb_navsto_builder_t *nsb,
2107 cs_cell_sys_t *csys)
2108 {
2109 CS_UNUSED(nsb);
2110 assert(nsp->model_flag & CS_NAVSTO_MODEL_BOUSSINESQ);
2111
2112 /* Boussinesq term: rho0 * g[] * ( 1 - beta * (var[c] - var0) ). The
2113 * remaining part rho0 * g[] * ( -beta * (var - var_c) has a zero mean-value
2114 * if one considers the reconstruction var = var_c + grad(vard)|_c * ( x -
2115 * x_c) which has a mean value equal to var_c */
2116
2117 const cs_real_t rho0 = nsp->mass_density->ref_value;
2118 const cs_real_t *gravity_vector = nsp->phys_constants->gravity;
2119
2120 cs_real_t rho0g[3] = { rho0 * gravity_vector[0],
2121 rho0 * gravity_vector[1],
2122 rho0 * gravity_vector[2] };
2123
2124 /* Constant part: rho_ref * g[] => Should be in balance with the pressure
2125 gradient in order to retrieve the hydrostatic case */
2126
2127 for (int f = 0; f < cm->n_fc; f++) {
2128 const cs_real_t *_div_f = nsb->div_op + 3*f;
2129 for (int k = 0; k < 3; k++)
2130 csys->rhs[3*f+k] += rho0g[k] * _div_f[k] * cm->xc[k];
2131 }
2132
2133 /* Volume part */
2134
2135 double boussi_coef = 0;
2136
2137 for (int i = 0; i < nsp->n_boussinesq_terms; i++) {
2138
2139 cs_navsto_param_boussinesq_t *bp = nsp->boussinesq_param + i;
2140 boussi_coef += -bp->beta*(bp->var[cm->c_id] - bp->var0);
2141
2142 }
2143
2144 for (int k = 0; k < 3; k++)
2145 csys->rhs[3*cm->n_fc+k] += rho0g[k] * boussi_coef * cm->vol_c;
2146 }
2147
2148 /*----------------------------------------------------------------------------*/
2149 /*!
2150 * \brief Take into account the buoyancy force with the Boussinesq approx.
2151 * Compute and add the source term to the local RHS.
2152 * This way to compute the Boussinesq approximation applies only to
2153 * face DoFs. This should enable to keep a stable (no velocity) in
2154 * case of a stratified configuration.
2155 *
2156 * \param[in] nsp set of parameters to handle the Navier-Stokes system
2157 * \param[in] cm pointer to a cs_cell_mesh_t structure
2158 * \param[in] nsb pointer to a builder structure for the NavSto system
2159 * \param[in, out] csys pointer to a cs_cell_sys_t structure
2160 */
2161 /*----------------------------------------------------------------------------*/
2162
2163 void
cs_cdofb_navsto_boussinesq_at_face(const cs_navsto_param_t * nsp,const cs_cell_mesh_t * cm,const cs_cdofb_navsto_builder_t * nsb,cs_cell_sys_t * csys)2164 cs_cdofb_navsto_boussinesq_at_face(const cs_navsto_param_t *nsp,
2165 const cs_cell_mesh_t *cm,
2166 const cs_cdofb_navsto_builder_t *nsb,
2167 cs_cell_sys_t *csys)
2168 {
2169 CS_UNUSED(nsb);
2170 assert(nsp->model_flag & CS_NAVSTO_MODEL_BOUSSINESQ);
2171
2172 /* Boussinesq term: rho0 * g[] * ( 1 - beta * (var[c] - var0) ). */
2173
2174 /* 1. Compute the mass density for this cell taking into account the
2175 Boussinesq approximation */
2176
2177 double boussi_coef = 1;
2178
2179 for (int i = 0; i < nsp->n_boussinesq_terms; i++) {
2180
2181 cs_navsto_param_boussinesq_t *bp = nsp->boussinesq_param + i;
2182 boussi_coef += -bp->beta*(bp->var[cm->c_id] - bp->var0);
2183
2184 } /* Loop on Boussinesq terms */
2185
2186 const cs_real_t rho_c = nsp->mass_density->ref_value * boussi_coef;
2187 const cs_real_t *gravity_vector = nsp->phys_constants->gravity;
2188 const double cell_coef = _dp3(gravity_vector, cm->xc);
2189
2190 for (int f = 0; f < cm->n_fc; f++) {
2191
2192 /* div_op is built such that _div_f[k] = -i_{f,c} * |f| * n_f[k] */
2193
2194 const cs_real_t *_div_f = nsb->div_op + 3*f;
2195 const cs_quant_t pfq = cm->face[f];
2196 const double face_coef = _dp3(gravity_vector, pfq.center);
2197 const double rhs_coef = rho_c * (cell_coef - face_coef);
2198
2199 for (int k = 0; k < 3; k++)
2200 csys->rhs[3*f+k] += rhs_coef * _div_f[k];
2201
2202 } /* Loop on cell faces */
2203
2204 }
2205
2206 /*----------------------------------------------------------------------------*/
2207 /*!
2208 * \brief Get the source term for computing the stream function.
2209 * This relies on the prototype associated to the generic function
2210 * pointer \ref cs_dof_function_t
2211 *
2212 * \param[in] n_elts number of elements to consider
2213 * \param[in] elt_ids list of elements ids
2214 * \param[in] dense_output perform an indirection in retval or not
2215 * \param[in] input NULL or pointer to a structure cast on-the-fly
2216 * \param[in, out] retval result of the function. Must be allocated.
2217 */
2218 /*----------------------------------------------------------------------------*/
2219
2220 void
cs_cdofb_navsto_stream_source_term(cs_lnum_t n_elts,const cs_lnum_t * elt_ids,bool dense_output,void * input,cs_real_t * retval)2221 cs_cdofb_navsto_stream_source_term(cs_lnum_t n_elts,
2222 const cs_lnum_t *elt_ids,
2223 bool dense_output,
2224 void *input,
2225 cs_real_t *retval)
2226 {
2227 assert(input != NULL);
2228 assert(retval != NULL);
2229
2230 /* input is a pointer to the vorticity field */
2231 const cs_real_t *w = (cs_real_t *)input;
2232
2233 for (cs_lnum_t i = 0; i < n_elts; i++) {
2234
2235 cs_lnum_t id = (elt_ids == NULL) ? i : elt_ids[i];
2236 cs_lnum_t r_id = dense_output ? i : id;
2237
2238 retval[r_id] = w[3*id+2]; /* Extract the z component */
2239
2240 }
2241 }
2242
2243 /*----------------------------------------------------------------------------*/
2244
2245 #undef _dp3
2246
2247 END_C_DECLS
2248