1 /*============================================================================
2 * Base turbulence boundary conditions.
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------*/
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39 #include <math.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include "bft_mem.h"
46 #include "bft_error.h"
47 #include "bft_printf.h"
48
49 #include "cs_log.h"
50 #include "cs_field.h"
51 #include "cs_field_pointer.h"
52 #include "cs_map.h"
53 #include "cs_math.h"
54 #include "cs_parall.h"
55 #include "cs_mesh_location.h"
56
57 #include "cs_turbulence_model.h"
58
59 /*----------------------------------------------------------------------------
60 * Header for the current file
61 *----------------------------------------------------------------------------*/
62
63 #include "cs_turbulence_bc.h"
64
65 /*----------------------------------------------------------------------------*/
66
67 BEGIN_C_DECLS
68
69 /*=============================================================================
70 * Additional doxygen documentation
71 *============================================================================*/
72
73 /*!
74 \file cs_turbulence_bc.c
75 Base turbulence boundary conditions.
76 */
77
78 /*----------------------------------------------------------------------------*/
79
80 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
81
82 /*=============================================================================
83 * Macro definitions
84 *============================================================================*/
85
86 /*============================================================================
87 * Type definitions
88 *============================================================================*/
89
90 /* Variable ids for main turbulence variables
91
92 CAUTION:
93
94 A redesign of the icodcl/rcodcl boundary condition logic
95 and simpler mapping to fields should remove the need for this
96 function in the future.
97
98 So this structure allows us an efficient access to current "rcodcl"
99 values for turbulent variables in the meantime, but should not
100 be a starting point for more general developments.
101
102 */
103
104 typedef struct {
105
106 int k; /* variable id for k */
107 int eps; /* variable id for epsilon */
108
109 int rij; /* variable id for R_ij tensor */
110
111 int phi; /* variable id for phi */
112 int f_bar; /* variable id for f_bar */
113 int alp_bl; /* variable id for blending alpha (dynamic) */
114
115 int omg; /* variable id for omega */
116 int nusa; /* variable id for nu_t (SA model) */
117
118 int size_ut; /* size of array of variable ids for turbulent fluxes */
119 int size_alp_bl_t; /* size of array of variable ids for blending alpha */
120 int *ut; /* variable ids for turbulent fluxes */
121 int *alp_bl_t; /* variable ids for blending alpha */
122
123 } cs_turb_bc_id_t;
124
125 /*============================================================================
126 * Static global variables
127 *============================================================================*/
128
129 /* Variable ids for main turbulence variables */
130
131 static cs_turb_bc_id_t
132 _turb_bc_id =
133 {
134 -1, /* k */
135 -1, /* eps */
136
137 -1, /* rij */
138
139 -1, /* phi */
140 -1, /* f_bar */
141 -1, /* alp_bl */
142
143 -1, /* omg */
144 -1, /* nusa */
145
146 0, /* size of ut */
147 0, /* size of alp_bl_t */
148 NULL, /* ut */
149 NULL /* alp_bl_t */
150 };
151
152 /*============================================================================
153 * Prototypes for functions intended for use only by Fortran wrappers.
154 * (descriptions follow, with function bodies).
155 *============================================================================*/
156
157 void
158 cs_f_turbulence_bc_inlet_hyd_diam(cs_lnum_t face_num,
159 double uref2,
160 double dh,
161 double rho,
162 double mu,
163 double *rcodcl);
164
165 void
166 cs_f_turbulence_bc_inlet_turb_intensity(cs_lnum_t face_num,
167 double uref2,
168 double t_intensity,
169 double dh,
170 double *rcodcl);
171
172 void
173 cs_f_turbulence_bc_inlet_k_eps(cs_lnum_t face_num,
174 double k,
175 double eps,
176 double vel_dir[],
177 double shear_dir[],
178 double *rcodcl);
179
180 void
181 cs_f_turbulence_bc_init_inlet_k_eps(cs_lnum_t face_num,
182 double k,
183 double eps,
184 double *rcodcl);
185
186 void
187 cs_f_turbulence_bc_set_uninit_inlet_k_eps(cs_lnum_t face_num,
188 double k,
189 double eps,
190 double vel_dir[],
191 double shear_dir[],
192 double *rcodcl);
193
194 /*============================================================================
195 * Private function definitions
196 *============================================================================*/
197
198 /*----------------------------------------------------------------------------
199 * Calculation of \f$ u^\star \f$, \f$ k \f$ and \f$\varepsilon \f$
200 * from a diameter \f$ D_H \f$ and the reference velocity \f$ U_{ref} \f$
201 * for a circular duct flow with smooth wall
202 * (for inlet boundary conditions).
203 *
204 * Both \f$ u^\star \f$ and\f$ (k,\varepsilon )\f$ are returned, so that
205 * the user may compute other values of \f$ k \f$ and \f$ \varepsilon \f$
206 * with \f$ u^\star \f$.
207 *
208 * We use the laws from Idel'Cik, i.e.
209 * the head loss coefficient \f$ \lambda \f$ is defined by:
210 * \f[ |\dfrac{\Delta P}{\Delta x}| =
211 * \dfrac{\lambda}{D_H} \frac{1}{2} \rho U_{ref}^2 \f]
212 *
213 * then the relation reads \f$u^\star = U_{ref} \sqrt{\dfrac{\lambda}{8}}\f$.
214 * \f$\lambda \f$ depends on the hydraulic Reynolds number
215 * \f$ Re = \dfrac{U_{ref} D_H}{ \nu} \f$ and is given by:
216 * - for \f$ Re < 2000 \f$
217 * \f[ \lambda = \dfrac{64}{Re} \f]
218 *
219 * - for \f$ Re > 4000 \f$
220 * \f[ \lambda = \dfrac{1}{( 1.8 \log_{10}(Re)-1.64 )^2} \f]
221 *
222 * - for \f$ 2000 < Re < 4000 \f$, we complete by a straight line
223 * \f[ \lambda = 0.021377 + 5.3115. 10^{-6} Re \f]
224 *
225 * From \f$ u^\star \f$, we can estimate \f$ k \f$ and \f$ \varepsilon\f$
226 * from the well known formulae of developped turbulence
227 *
228 * \f[ k = \dfrac{u^{\star 2}}{\sqrt{C_\mu}} \f]
229 * \f[ \varepsilon = \dfrac{ u^{\star 3}}{(\kappa D_H /10)} \f]
230 *
231 * parameters:
232 * uref2 <-- square of the reference flow velocity
233 * dh <-- hydraulic diameter \f$ D_H \f$
234 * rho <-- mass density \f$ \rho \f$
235 * mu <-- dynamic viscosity \f$ \nu \f$
236 * ustar2 --> square of friction speed
237 * k --> calculated turbulent intensity \f$ k \f$
238 * eps --> calculated turbulent dissipation
239 * \f$ \varepsilon \f$
240 *----------------------------------------------------------------------------*/
241
242 static inline void
_ke_hyd_diam(double uref2,double dh,double rho,double mu,double * ustar2,double * k,double * eps)243 _ke_hyd_diam(double uref2,
244 double dh,
245 double rho,
246 double mu,
247 double *ustar2,
248 double *k,
249 double *eps)
250 {
251 double re = sqrt(uref2)*dh*rho/mu;
252
253 if (re < 2000) {
254
255 /* in this case we calculate directly \f$u*^2\f$ to avoid an issue with
256 \f$ xlmbda= \dfrac{64}{Re} \f$ when Re->0 */
257
258 *ustar2 = 8.0*mu*sqrt(uref2)/rho/dh;
259
260 }
261 else if (re < 4000) {
262
263 double xlmbda = 0.021377 + 5.3115e-6*re;
264 *ustar2 = uref2*xlmbda/8.0;
265
266 }
267 else {
268
269 /* xlmbda = 1.0/(1.8*log(re) / log(10.0)-1.64)^2;
270 ustar2 = uref2*xlmbda/8.0; */
271
272 double a = 1.8*log(re) / log(10.0)-1.64;
273 *ustar2 = uref2*0.125/(a*a);
274
275 }
276
277 *k = *ustar2 / sqrt(cs_turb_cmu);
278 *eps = pow(*ustar2, 1.5) / (cs_turb_xkappa*dh*0.1);
279 }
280
281 /*----------------------------------------------------------------------------
282 * Calculation of \f$ k \f$ and \f$\varepsilon\f$
283 * from a diameter \f$ D_H \f$, a turbulent intensity \f$ I \f$
284 * and the reference velocity \f$ U_{ref} \f$
285 * for a circular duct flow with smooth wall
286 * (for inlet boundary conditions).
287 *
288 * \f[ k = 1.5 I {U_{ref}}^2 \f]
289 * \f[ \varepsilon = 10 \dfrac{{C_\mu}^{0.75} k^{1.5}}{ \kappa D_H} \f]
290 *
291 * parameters:
292 * uref2 <-- square of the reference flow velocity
293 * t_intensity <-- turbulent intensity \f$ I \f$
294 * dh <-- hydraulic diameter \f$ D_H \f$
295 * k --> calculated turbulent intensity \f$ k \f$
296 * eps --> calculated turbulent dissipation
297 * \f$ \varepsilon \f$
298 *----------------------------------------------------------------------------*/
299
300 static inline void
_ke_turb_intensity(double uref2,double t_intensity,double dh,double * k,double * eps)301 _ke_turb_intensity(double uref2,
302 double t_intensity,
303 double dh,
304 double *k,
305 double *eps)
306 {
307 *k = 1.5*uref2*t_intensity*t_intensity;
308 *eps = 10.0*pow(cs_turb_cmu, 0.75)*pow(*k, 1.5)/(cs_turb_xkappa*dh);
309 }
310
311 /*----------------------------------------------------------------------------*
312 * Assign turbulent boundary condition to a given face
313 *
314 * parameters:
315 * face_id <-- face id
316 * k <-- k
317 * eps <-- epsilon
318 * vel_dir <-- velocity direction
319 * shear_dir <-- shear direction, it also contains the level of
320 * anisotropy (Rnt = a_nt k)
321 * rcodcl <-> boundary condition values
322 *----------------------------------------------------------------------------*/
323
324 static inline void
_inlet_bc(cs_lnum_t face_id,double k,double eps,double vel_dir[],double shear_dir[],double * rcodcl)325 _inlet_bc(cs_lnum_t face_id,
326 double k,
327 double eps,
328 double vel_dir[],
329 double shear_dir[],
330 double *rcodcl)
331 {
332 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
333 const cs_turb_model_t *turb_model = cs_get_glob_turb_model();
334 assert(turb_model != NULL);
335
336 if (turb_model->itytur == 2) {
337
338 rcodcl[_turb_bc_id.k *n_b_faces + face_id] = k;
339 rcodcl[_turb_bc_id.eps*n_b_faces + face_id] = eps;
340
341 }
342
343 else if (turb_model->order == CS_TURB_SECOND_ORDER) {
344
345 cs_lnum_t s_r_11 = _turb_bc_id.rij * n_b_faces;
346 cs_lnum_t s_r_22 = (_turb_bc_id.rij + 1) * n_b_faces;
347 cs_lnum_t s_r_33 = (_turb_bc_id.rij + 2) * n_b_faces;
348 cs_lnum_t s_r_12 = (_turb_bc_id.rij + 3) * n_b_faces;
349 cs_lnum_t s_r_23 = (_turb_bc_id.rij + 4) * n_b_faces;
350 cs_lnum_t s_r_13 = (_turb_bc_id.rij + 5) * n_b_faces;
351
352 double d2s3 = 2./3.;
353
354 rcodcl[s_r_11 + face_id] = d2s3 * k;
355 rcodcl[s_r_22 + face_id] = d2s3 * k;
356 rcodcl[s_r_33 + face_id] = d2s3 * k;
357 if (vel_dir != NULL) {
358 cs_math_3_normalize(vel_dir, vel_dir);
359 /* Note: do not normalize shear_dir because it contains the level of anisotropy */
360 /* Rxy */
361 rcodcl[s_r_12 + face_id]
362 = k * (vel_dir[0]*shear_dir[1] + vel_dir[1]*shear_dir[0]);
363 /* Ryz */
364 rcodcl[s_r_23 + face_id]
365 = k * (vel_dir[1]*shear_dir[2] + vel_dir[2]*shear_dir[1]);
366 /* Rxz */
367 rcodcl[s_r_13 + face_id]
368 = k * (vel_dir[0]*shear_dir[2] + vel_dir[2]*shear_dir[0]);
369 }
370 else {
371 rcodcl[s_r_12 + face_id] = 0.;
372 rcodcl[s_r_23 + face_id] = 0.;
373 rcodcl[s_r_13 + face_id] = 0.;
374 }
375 rcodcl[_turb_bc_id.eps*n_b_faces + face_id] = eps;
376
377 if (turb_model->iturb == CS_TURB_RIJ_EPSILON_EBRSM)
378 rcodcl[_turb_bc_id.alp_bl*n_b_faces + face_id] = 1.;
379
380 /* Initialization of the turbulent fluxes to 0 if DFM or
381 * EB-DFM are used for scalars (turbulence_flux_model = 30 or 31)
382 * Alpha_theta for EB-DFM / EB-AFM / EB-GGDH is
383 * initialize to 1 */
384
385 if (_turb_bc_id.size_ut > 0) {
386 for (int var_id = 0; var_id < _turb_bc_id.size_ut; var_id++) {
387 rcodcl[_turb_bc_id.ut[var_id]*n_b_faces + face_id] = 0.;
388 rcodcl[(_turb_bc_id.ut[var_id]+1)*n_b_faces + face_id] = 0.;
389 rcodcl[(_turb_bc_id.ut[var_id]+2)*n_b_faces + face_id] = 0.;
390 }
391 }
392
393 if (_turb_bc_id.size_alp_bl_t > 0) {
394 for (int var_id = 0; var_id < _turb_bc_id.size_alp_bl_t; var_id++) {
395 rcodcl[_turb_bc_id.alp_bl_t[var_id]*n_b_faces + face_id] = 1.;
396 }
397 }
398
399 }
400 else if (turb_model->itytur == 5) {
401
402 rcodcl[_turb_bc_id.k *n_b_faces + face_id] = k;
403 rcodcl[_turb_bc_id.eps*n_b_faces + face_id] = eps;
404
405 rcodcl[_turb_bc_id.phi*n_b_faces + face_id] = 2./3.;
406 if (turb_model->iturb == 50) {
407 rcodcl[_turb_bc_id.f_bar*n_b_faces + face_id] = 0.;
408 }
409 else if (turb_model->iturb == 51) {
410 rcodcl[_turb_bc_id.alp_bl*n_b_faces + face_id] = 0.;
411 }
412
413 }
414 else if (turb_model->iturb == CS_TURB_K_OMEGA) {
415
416 rcodcl[_turb_bc_id.k *n_b_faces + face_id] = k;
417 rcodcl[_turb_bc_id.omg*n_b_faces + face_id] = eps/cs_turb_cmu/k;
418
419 }
420 else if (turb_model->iturb == CS_TURB_SPALART_ALLMARAS) {
421
422 rcodcl[_turb_bc_id.nusa*n_b_faces + face_id] = cs_turb_cmu*k*k/eps;
423
424 }
425 }
426
427 /*----------------------------------------------------------------------------*
428 * Assign turbulent boundary condition to a given face only if unitialized.
429 *
430 * parameters:
431 * face_id <-- face id
432 * k <-- k
433 * eps <-- epsilon
434 * vel_dir <-- velocity direction
435 * shear_dir <-- shear direction, it also contains the level of
436 * anisotropy (Rnt = a_nt k)
437 * rcodcl <-> boundary condition values
438 *----------------------------------------------------------------------------*/
439
440 static inline void
_set_uninit_inlet_bc(cs_lnum_t face_id,double k,double eps,double vel_dir[],double shear_dir[],double * rcodcl)441 _set_uninit_inlet_bc(cs_lnum_t face_id,
442 double k,
443 double eps,
444 double vel_dir[],
445 double shear_dir[],
446 double *rcodcl)
447 {
448 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
449 const cs_turb_model_t *turb_model = cs_get_glob_turb_model();
450
451 if (turb_model->itytur == 2) {
452
453 if (rcodcl[_turb_bc_id.k *n_b_faces + face_id] > 0.5*cs_math_infinite_r)
454 rcodcl[_turb_bc_id.k *n_b_faces + face_id] = k;
455
456 if (rcodcl[_turb_bc_id.eps*n_b_faces + face_id] > 0.5*cs_math_infinite_r)
457 rcodcl[_turb_bc_id.eps*n_b_faces + face_id] = eps;
458
459 }
460
461 else if (turb_model->order == CS_TURB_SECOND_ORDER) {
462
463 cs_lnum_t s_r_11 = _turb_bc_id.rij * n_b_faces;
464 cs_lnum_t s_r_22 = (_turb_bc_id.rij + 1) * n_b_faces;
465 cs_lnum_t s_r_33 = (_turb_bc_id.rij + 2) * n_b_faces;
466 cs_lnum_t s_r_12 = (_turb_bc_id.rij + 3) * n_b_faces;
467 cs_lnum_t s_r_23 = (_turb_bc_id.rij + 4) * n_b_faces;
468 cs_lnum_t s_r_13 = (_turb_bc_id.rij + 5) * n_b_faces;
469
470 double d2s3 = 2./3.;
471 if (rcodcl[s_r_11 + face_id] > 0.5*cs_math_infinite_r)
472 rcodcl[s_r_11 + face_id] = d2s3 * k;
473 if (rcodcl[s_r_22 + face_id] > 0.5*cs_math_infinite_r)
474 rcodcl[s_r_22 + face_id] = d2s3 * k;
475 if (rcodcl[s_r_33 + face_id] > 0.5*cs_math_infinite_r)
476 rcodcl[s_r_33 + face_id] = d2s3 * k;
477 if (vel_dir != NULL) {
478 cs_math_3_normalize(vel_dir, vel_dir);
479 if (rcodcl[s_r_12 + face_id] > 0.5*cs_math_infinite_r)
480 rcodcl[s_r_12 + face_id]
481 = k * (vel_dir[0]*shear_dir[1] + vel_dir[1]*shear_dir[0]);
482 if (rcodcl[s_r_23 + face_id] > 0.5*cs_math_infinite_r)
483 rcodcl[s_r_23 + face_id]
484 = k * (vel_dir[1]*shear_dir[2] + vel_dir[2]*shear_dir[1]);
485 if (rcodcl[s_r_13 + face_id] > 0.5*cs_math_infinite_r)
486 rcodcl[s_r_13 + face_id]
487 = k * (vel_dir[0]*shear_dir[2] + vel_dir[2]*shear_dir[0]);
488 }
489 else {
490 if (rcodcl[s_r_12 + face_id] > 0.5*cs_math_infinite_r)
491 rcodcl[s_r_12 + face_id] = 0.;
492 if (rcodcl[s_r_23 + face_id] > 0.5*cs_math_infinite_r)
493 rcodcl[s_r_23 + face_id] = 0.;
494 if (rcodcl[s_r_13 + face_id] > 0.5*cs_math_infinite_r)
495 rcodcl[s_r_13 + face_id] = 0.;
496 }
497 if (rcodcl[_turb_bc_id.eps*n_b_faces + face_id] > 0.5*cs_math_infinite_r)
498 rcodcl[_turb_bc_id.eps*n_b_faces + face_id] = eps;
499
500 if (turb_model->iturb == 32)
501 if ( rcodcl[_turb_bc_id.alp_bl*n_b_faces + face_id]
502 > 0.5*cs_math_infinite_r)
503 rcodcl[_turb_bc_id.alp_bl*n_b_faces + face_id] = 1.;
504
505 /* Initialization of the turbulent fluxes to 0 if DFM or
506 * EB-DFM are used for scalars (turbulence_flux_model = 30 or 31)
507 * Alpha_theta for EB-DFM / EB-AFM / EB-GGDH is
508 * initialize to 1 */
509
510 if (_turb_bc_id.size_ut > 0) {
511 for (int var_id = 0; var_id < _turb_bc_id.size_ut; var_id++) {
512 if ( rcodcl[_turb_bc_id.ut[var_id]*n_b_faces + face_id]
513 > 0.5*cs_math_infinite_r)
514 rcodcl[_turb_bc_id.ut[var_id]*n_b_faces + face_id] = 0.;
515 if ( rcodcl[(_turb_bc_id.ut[var_id]+1)*n_b_faces + face_id]
516 > 0.5*cs_math_infinite_r)
517 rcodcl[(_turb_bc_id.ut[var_id]+1)*n_b_faces + face_id] = 0.;
518 if ( rcodcl[(_turb_bc_id.ut[var_id]+2)*n_b_faces + face_id]
519 > 0.5*cs_math_infinite_r)
520 rcodcl[(_turb_bc_id.ut[var_id]+2)*n_b_faces + face_id] = 0.;
521 }
522 }
523
524 if (_turb_bc_id.size_alp_bl_t > 0) {
525 for (int var_id = 0; var_id < _turb_bc_id.size_alp_bl_t; var_id++) {
526 if ( rcodcl[_turb_bc_id.alp_bl_t[var_id]*n_b_faces + face_id]
527 > 0.5*cs_math_infinite_r)
528 rcodcl[_turb_bc_id.alp_bl_t[var_id]*n_b_faces + face_id] = 1.;
529 }
530 }
531
532 }
533 else if (turb_model->itytur == 5) {
534
535 if (rcodcl[_turb_bc_id.k *n_b_faces + face_id] > 0.5*cs_math_infinite_r)
536 rcodcl[_turb_bc_id.k *n_b_faces + face_id] = k;
537 if (rcodcl[_turb_bc_id.eps*n_b_faces + face_id] > 0.5*cs_math_infinite_r)
538 rcodcl[_turb_bc_id.eps*n_b_faces + face_id] = eps;
539
540 if (rcodcl[_turb_bc_id.phi*n_b_faces + face_id] > 0.5*cs_math_infinite_r)
541 rcodcl[_turb_bc_id.phi*n_b_faces + face_id] = 2./3.;
542 if (turb_model->iturb == 50) {
543 if ( rcodcl[_turb_bc_id.f_bar*n_b_faces + face_id]
544 > 0.5*cs_math_infinite_r)
545 rcodcl[_turb_bc_id.f_bar*n_b_faces + face_id] = 0.;
546 }
547 else if (turb_model->iturb == 51) {
548 if ( rcodcl[_turb_bc_id.alp_bl*n_b_faces + face_id]
549 > 0.5*cs_math_infinite_r)
550 rcodcl[_turb_bc_id.alp_bl*n_b_faces + face_id] = 0.;
551 }
552
553 }
554 else if (turb_model->iturb == CS_TURB_K_OMEGA) {
555
556 if (rcodcl[_turb_bc_id.k *n_b_faces + face_id] > 0.5*cs_math_infinite_r)
557 rcodcl[_turb_bc_id.k *n_b_faces + face_id] = k;
558 if (rcodcl[_turb_bc_id.omg*n_b_faces + face_id] > 0.5*cs_math_infinite_r)
559 rcodcl[_turb_bc_id.omg*n_b_faces + face_id] = eps/cs_turb_cmu/k;
560
561 }
562 else if (turb_model->iturb == CS_TURB_SPALART_ALLMARAS) {
563 if (rcodcl[_turb_bc_id.nusa*n_b_faces + face_id] > 0.5*cs_math_infinite_r)
564 rcodcl[_turb_bc_id.nusa*n_b_faces + face_id] = cs_turb_cmu*k*k/eps;
565
566 }
567 }
568
569 /*============================================================================
570 * Fortran wrapper function definitions
571 *============================================================================*/
572
573 /*----------------------------------------------------------------------------
574 * Equivalent of cs_turbulence_bc_inlet_ke_hyd_diam for Fortran calls
575 * (using 1-based face number instead of id).
576 *
577 * parameters:
578 * face_num <-- face number
579 * uref2 <-- square of the reference flow velocity
580 * dh <-- hydraulic diameter \f$ D_H \f$
581 * rho <-- mass density \f$ \rho \f$
582 * mu <-- dynamic viscosity \f$ \nu \f$
583 * rcodcl <-> boundary condition values
584 *----------------------------------------------------------------------------*/
585
586 void
cs_f_turbulence_bc_inlet_hyd_diam(cs_lnum_t face_num,double uref2,double dh,double rho,double mu,double * rcodcl)587 cs_f_turbulence_bc_inlet_hyd_diam(cs_lnum_t face_num,
588 double uref2,
589 double dh,
590 double rho,
591 double mu,
592 double *rcodcl)
593 {
594 double ustar2, k, eps;
595
596 _ke_hyd_diam(uref2, dh, rho, mu, &ustar2, &k, &eps);
597
598 _inlet_bc(face_num - 1, k, eps, NULL, NULL, rcodcl);
599 }
600
601 /*----------------------------------------------------------------------------
602 * Equivalent of cs_turbulence_bc_inlet_ke_hyd_diam for Fortran calls
603 * (using 1-based face number instead of id).
604 *
605 * parameters:
606 * face_num <-- face number
607 * uref2 <-- square of the reference flow velocity
608 * t_intensity <-- turbulence intensity
609 * dh <-- hydraulic diameter \f$ D_H \f$
610 * rcodcl <-> boundary condition values
611 *----------------------------------------------------------------------------*/
612
613 void
cs_f_turbulence_bc_inlet_turb_intensity(cs_lnum_t face_num,double uref2,double t_intensity,double dh,double * rcodcl)614 cs_f_turbulence_bc_inlet_turb_intensity(cs_lnum_t face_num,
615 double uref2,
616 double t_intensity,
617 double dh,
618 double *rcodcl)
619 {
620 double k, eps;
621
622 _ke_turb_intensity(uref2, t_intensity, dh, &k, &eps);
623
624 _inlet_bc(face_num - 1, k, eps, NULL, NULL, rcodcl);
625 }
626
627 /*----------------------------------------------------------------------------
628 * Equivalent of cs_turbulence_bc_inlet_ke for Fortran calls
629 * (using 1-based face number instead of id).
630 *
631 * parameters:
632 * face_num <-- face number
633 * k <-- turbulent kinetic energy
634 * eps <-- turbulent dissipation
635 * vel_dir <-- velocity direction
636 * shear_dir <-- shear direction
637 * rcodcl <-> boundary condition values
638 *----------------------------------------------------------------------------*/
639
640 void
cs_f_turbulence_bc_inlet_k_eps(cs_lnum_t face_num,double k,double eps,double vel_dir[],double shear_dir[],double * rcodcl)641 cs_f_turbulence_bc_inlet_k_eps(cs_lnum_t face_num,
642 double k,
643 double eps,
644 double vel_dir[],
645 double shear_dir[],
646 double *rcodcl)
647 {
648 _inlet_bc(face_num - 1, k, eps, vel_dir, shear_dir, rcodcl);
649 }
650
651 /*----------------------------------------------------------------------------
652 * Equivalent of cs_turbulence_bc_set_uninit_inlet_ke for Fortran calls
653 * (using 1-based face number instead of id).
654 *
655 * parameters:
656 * face_num <-- face number
657 * k <-- turbulent kinetic energy
658 * eps <-- turbulent dissipation
659 * vel_dir <-- velocity direction
660 * shear_dir <-- shear direction
661 * rcodcl <-> boundary condition values
662 *----------------------------------------------------------------------------*/
663
664 void
cs_f_turbulence_bc_set_uninit_inlet_k_eps(cs_lnum_t face_num,double k,double eps,double vel_dir[],double shear_dir[],double * rcodcl)665 cs_f_turbulence_bc_set_uninit_inlet_k_eps(cs_lnum_t face_num,
666 double k,
667 double eps,
668 double vel_dir[],
669 double shear_dir[],
670 double *rcodcl)
671 {
672 _set_uninit_inlet_bc(face_num - 1, k, eps, vel_dir, shear_dir, rcodcl);
673 }
674
675 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
676
677 /*=============================================================================
678 * Public function definitions
679 *============================================================================*/
680
681 /*----------------------------------------------------------------------------*/
682 /*!
683 * \brief Initialize turbulence model boundary condition ids.
684 */
685 /*----------------------------------------------------------------------------*/
686
687 void
cs_turbulence_model_init_bc_ids(void)688 cs_turbulence_model_init_bc_ids(void)
689 {
690 /* CAUTION: see note about the cs_turb_bc_id structure above. */
691
692 const int var_key_id = cs_field_key_id("variable_id");
693 const int k_turbt = cs_field_key_id("turbulent_flux_model");
694 const int k_f_turbt = cs_field_key_id("turbulent_flux_id");
695 const int k_f_turbt_alp_bl = cs_field_key_id("alpha_turbulent_flux_id");
696 const int k_sca = cs_field_key_id("scalar_id");
697
698 if (CS_F_(k) != NULL)
699 _turb_bc_id.k = cs_field_get_key_int(CS_F_(k), var_key_id) -1;
700 if (CS_F_(eps) != NULL)
701 _turb_bc_id.eps = cs_field_get_key_int(CS_F_(eps), var_key_id) -1;
702
703 if (CS_F_(rij) != NULL)
704 _turb_bc_id.rij = cs_field_get_key_int(CS_F_(rij), var_key_id) -1;
705
706 if (CS_F_(phi) != NULL)
707 _turb_bc_id.phi = cs_field_get_key_int(CS_F_(phi), var_key_id) -1;
708 if (CS_F_(f_bar) != NULL)
709 _turb_bc_id.f_bar = cs_field_get_key_int(CS_F_(f_bar), var_key_id) -1;
710 if (CS_F_(alp_bl) != NULL)
711 _turb_bc_id.alp_bl = cs_field_get_key_int(CS_F_(alp_bl), var_key_id) -1;
712
713 if (CS_F_(omg) != NULL)
714 _turb_bc_id.omg = cs_field_get_key_int(CS_F_(omg), var_key_id) -1;
715 if (CS_F_(nusa) != NULL)
716 _turb_bc_id.nusa = cs_field_get_key_int(CS_F_(nusa), var_key_id) -1;
717
718 int n_fields = cs_field_n_fields();
719 int n_sca_ut = 0;
720 int n_sca_alp_bl = 0;
721
722 /* For scalar turbulent fluxes, loop over all scalars to determine:
723 * - number of scalars with (EB)DFM (turbulence_flux_model = 30 or 31)
724 * - number of scalars using an elliptic blending model
725 * (turbulence_flux_model = 11 or 21 or 31) */
726
727 for (int f_id = 0; f_id < n_fields; f_id++) {
728 const cs_field_t *f = cs_field_by_id(f_id);
729 if (f->type & CS_FIELD_VARIABLE) {
730 int s_num = cs_field_get_key_int(f, k_sca);
731 if (s_num > 0) {
732 int f_turbt = cs_field_get_key_int(f, k_turbt) ;
733 if (f_turbt / 10 == 3)
734 n_sca_ut ++;
735 if (f_turbt == 11 || f_turbt == 21 || f_turbt == 31)
736 n_sca_alp_bl ++;
737 }
738 }
739 }
740
741 _turb_bc_id.size_ut = n_sca_ut;
742 _turb_bc_id.size_alp_bl_t = n_sca_alp_bl;
743
744 if (_turb_bc_id.size_ut > 0)
745 BFT_MALLOC(_turb_bc_id.ut , n_sca_ut , int);
746 if (_turb_bc_id.size_alp_bl_t > 0)
747 BFT_MALLOC( _turb_bc_id.alp_bl_t, n_sca_alp_bl, int);
748
749 n_sca_ut = 0;
750 n_sca_alp_bl = 0;
751
752 for (int f_id = 0; f_id < n_fields; f_id++) {
753 const cs_field_t *f = cs_field_by_id(f_id);
754 if (f->type & CS_FIELD_VARIABLE) {
755 int s_num = cs_field_get_key_int(f, k_sca);
756 if (s_num > 0) {
757 int f_turbt = cs_field_get_key_int(f, k_turbt) ;
758 if (f_turbt / 10 == 3) {
759 int fid_turbt = cs_field_get_key_int(f, k_f_turbt);
760 _turb_bc_id.ut[n_sca_ut] =
761 cs_field_get_key_int(cs_field_by_id(fid_turbt), var_key_id) -1;
762 n_sca_ut ++;
763 }
764 if (f_turbt == 11 || f_turbt == 21 || f_turbt == 31) {
765 int fid_turbt = cs_field_get_key_int(f, k_f_turbt_alp_bl);
766 _turb_bc_id.alp_bl_t[n_sca_alp_bl] =
767 cs_field_get_key_int(cs_field_by_id(fid_turbt), var_key_id) -1;
768 n_sca_alp_bl ++;
769 }
770 }
771 }
772 }
773 }
774
775 /*----------------------------------------------------------------------------*/
776 /*!
777 * \brief Free memory allocations for turbulence boundary conditions ids.
778 */
779 /*----------------------------------------------------------------------------*/
780
781 void
cs_turbulence_model_free_bc_ids(void)782 cs_turbulence_model_free_bc_ids(void)
783 {
784 if (_turb_bc_id.size_ut > 0)
785 BFT_FREE(_turb_bc_id.ut);
786 if (_turb_bc_id.size_alp_bl_t > 0)
787 BFT_FREE( _turb_bc_id.alp_bl_t);
788 }
789
790 /*----------------------------------------------------------------------------*/
791 /*!
792 * \brief Calculation of \f$ u^\star \f$, \f$ k \f$ and \f$\varepsilon \f$
793 * from a diameter \f$ D_H \f$ and the reference velocity \f$ U_{ref} \f$
794 * for a circular duct flow with smooth wall
795 * (use for inlet boundary conditions).
796 *
797 * Both \f$ u^\star \f$ and\f$ (k,\varepsilon )\f$ are returned, so that
798 * the user may compute other values of \f$ k \f$ and \f$ \varepsilon \f$
799 * with \f$ u^\star \f$.
800 *
801 * We use the laws from Idel'Cik, i.e.
802 * the head loss coefficient \f$ \lambda \f$ is defined by:
803 * \f[ |\dfrac{\Delta P}{\Delta x}| =
804 * \dfrac{\lambda}{D_H} \frac{1}{2} \rho U_{ref}^2 \f]
805 *
806 * then the relation reads \f$u^\star = U_{ref} \sqrt{\dfrac{\lambda}{8}}\f$.
807 * \f$\lambda \f$ depends on the hydraulic Reynolds number
808 * \f$ Re = \dfrac{U_{ref} D_H}{ \nu} \f$ and is given by:
809 * - for \f$ Re < 2000 \f$
810 * \f[ \lambda = \dfrac{64}{Re} \f]
811 *
812 * - for \f$ Re > 4000 \f$
813 * \f[ \lambda = \dfrac{1}{( 1.8 \log_{10}(Re)-1.64 )^2} \f]
814 *
815 * - for \f$ 2000 < Re < 4000 \f$, we complete by a straight line
816 * \f[ \lambda = 0.021377 + 5.3115. 10^{-6} Re \f]
817 *
818 * From \f$ u^\star \f$, we can estimate \f$ k \f$ and \f$ \varepsilon\f$
819 * from the well known formulae of developped turbulence
820 *
821 * \f[ k = \dfrac{u^{\star 2}}{\sqrt{C_\mu}} \f]
822 * \f[ \varepsilon = \dfrac{ u^{\star 3}}{(\kappa D_H /10)} \f]
823 *
824 * \param[in] uref2 square of the reference flow velocity
825 * \param[in] dh hydraulic diameter \f$ D_H \f$
826 * \param[in] rho mass density \f$ \rho \f$
827 * \param[in] mu dynamic viscosity \f$ \nu \f$
828 * \param[out] ustar2 square of friction speed
829 * \param[out] k calculated turbulent intensity \f$ k \f$
830 * \param[out] eps calculated turbulent dissipation
831 * \f$ \varepsilon \f$
832 */
833 /*----------------------------------------------------------------------------*/
834
835 void
cs_turbulence_bc_ke_hyd_diam(double uref2,double dh,double rho,double mu,double * ustar2,double * k,double * eps)836 cs_turbulence_bc_ke_hyd_diam(double uref2,
837 double dh,
838 double rho,
839 double mu,
840 double *ustar2,
841 double *k,
842 double *eps)
843 {
844 _ke_hyd_diam(uref2, dh, rho, mu, ustar2, k, eps);
845 }
846
847 /*----------------------------------------------------------------------------*/
848 /*!
849 * \brief Calculation of \f$ k \f$ and \f$\varepsilon\f$
850 * from a diameter \f$ D_H \f$, a turbulent intensity \f$ I \f$
851 * and the reference velocity \f$ U_{ref} \f$
852 * for a circular duct flow with smooth wall
853 * (for inlet boundary conditions).
854 *
855 * \f[ k = 1.5 I {U_{ref}}^2 \f]
856 * \f[ \varepsilon = 10 \dfrac{{C_\mu}^{0.75} k^{1.5}}{ \kappa D_H} \f]
857 *
858 * \param[in] uref2 square of the reference flow velocity
859 * \param[in] t_intensity turbulent intensity \f$ I \f$
860 * \param[in] dh hydraulic diameter \f$ D_H \f$
861 * \param[out] k calculated turbulent intensity \f$ k \f$
862 * \param[out] eps calculated turbulent dissipation
863 * \f$ \varepsilon \f$
864 */
865 /*----------------------------------------------------------------------------*/
866
867 void
cs_turbulence_bc_ke_turb_intensity(double uref2,double t_intensity,double dh,double * k,double * eps)868 cs_turbulence_bc_ke_turb_intensity(double uref2,
869 double t_intensity,
870 double dh,
871 double *k,
872 double *eps)
873 {
874 _ke_turb_intensity(uref2, t_intensity, dh, k, eps);
875 }
876
877 /*----------------------------------------------------------------------------*/
878 /*!
879 * \brief Set inlet boundary condition values for turbulence variables based
880 * on a diameter \f$ D_H \f$ and the reference velocity \f$ U_{ref} \f$
881 * for a circular duct flow with smooth wall.
882 *
883 * We use the laws from Idel'Cik, i.e.
884 * the head loss coefficient \f$ \lambda \f$ is defined by:
885 * \f[ |\dfrac{\Delta P}{\Delta x}| =
886 * \dfrac{\lambda}{D_H} \frac{1}{2} \rho U_{ref}^2 \f]
887 *
888 * then the relation reads \f$u^\star = U_{ref} \sqrt{\dfrac{\lambda}{8}}\f$.
889 * \f$\lambda \f$ depends on the hydraulic Reynolds number
890 * \f$ Re = \dfrac{U_{ref} D_H}{ \nu} \f$ and is given by:
891 * - for \f$ Re < 2000 \f$
892 * \f[ \lambda = \dfrac{64}{Re} \f]
893 *
894 * - for \f$ Re > 4000 \f$
895 * \f[ \lambda = \dfrac{1}{( 1.8 \log_{10}(Re)-1.64 )^2} \f]
896 *
897 * - for \f$ 2000 < Re < 4000 \f$, we complete by a straight line
898 * \f[ \lambda = 0.021377 + 5.3115. 10^{-6} Re \f]
899 *
900 * From \f$ u^\star \f$, we can estimate \f$ k \f$ and \f$ \varepsilon\f$
901 * from the well known formulae of developped turbulence
902 *
903 * \param[in] face_id boundary face id
904 * \param[in] uref2 square of the reference flow velocity
905 * \param[in] dh hydraulic diameter \f$ D_H \f$
906 * \param[in] rho mass density \f$ \rho \f$
907 * \param[in] mu dynamic viscosity \f$ \nu \f$
908 * \param[out] rcodcl boundary condition values
909 */
910 /*----------------------------------------------------------------------------*/
911
912 void
cs_turbulence_bc_inlet_hyd_diam(cs_lnum_t face_id,double uref2,double dh,double rho,double mu,double * rcodcl)913 cs_turbulence_bc_inlet_hyd_diam(cs_lnum_t face_id,
914 double uref2,
915 double dh,
916 double rho,
917 double mu,
918 double *rcodcl)
919 {
920 double ustar2, k, eps;
921
922 _ke_hyd_diam(uref2, dh, rho, mu, &ustar2, &k, &eps);
923
924 _inlet_bc(face_id, k, eps, NULL, NULL, rcodcl);
925 }
926
927 /*----------------------------------------------------------------------------*/
928 /*!
929 * \brief Set inlet boundary condition values for turbulence variables based
930 * on a diameter \f$ D_H \f$, a turbulent intensity \f$ I \f$
931 * and the reference velocity \f$ U_{ref} \f$
932 * for a circular duct flow with smooth wall.
933 *
934 * \param[in] face_id boundary face id
935 * \param[in] uref2 square of the reference flow velocity
936 * \param[in] t_intensity turbulent intensity \f$ I \f$
937 * \param[in] dh hydraulic diameter \f$ D_H \f$
938 * \param[out] rcodcl boundary condition values
939 */
940 /*----------------------------------------------------------------------------*/
941
942 void
cs_turbulence_bc_inlet_turb_intensity(cs_lnum_t face_id,double uref2,double t_intensity,double dh,double * rcodcl)943 cs_turbulence_bc_inlet_turb_intensity(cs_lnum_t face_id,
944 double uref2,
945 double t_intensity,
946 double dh,
947 double *rcodcl)
948 {
949 double k, eps;
950
951 _ke_turb_intensity(uref2, t_intensity, dh, &k, &eps);
952
953 _inlet_bc(face_id, k, eps, NULL, NULL, rcodcl);
954 }
955
956 /*----------------------------------------------------------------------------*/
957 /*!
958 * \brief Set inlet boundary condition values for turbulence variables based
959 * on given k and epsilon values.
960 *
961 * \param[in] face_id boundary face id
962 * \param[in] k turbulent kinetic energy
963 * \param[in] eps turbulent dissipation
964 * \param[out] rcodcl boundary condition values
965 */
966 /*----------------------------------------------------------------------------*/
967
968 void
cs_turbulence_bc_inlet_k_eps(cs_lnum_t face_id,double k,double eps,double * rcodcl)969 cs_turbulence_bc_inlet_k_eps(cs_lnum_t face_id,
970 double k,
971 double eps,
972 double *rcodcl)
973 {
974 _inlet_bc(face_id, k, eps, NULL, NULL, rcodcl);
975 }
976
977 /*----------------------------------------------------------------------------*/
978 /*!
979 * \brief Set inlet boundary condition values for turbulence variables based
980 * on given k and epsilon values only if not already initialized.
981 *
982 * \param[in] face_id boundary face id
983 * \param[in] k turbulent kinetic energy
984 * \param[in] eps turbulent dissipation
985 * \param[in] vel_dir velocity direction
986 * \param[in] shear_dir shear direction, it also contains the level of
987 * anisotropy (Rnt = a_nt k)
988 * \param[out] rcodcl boundary condition values
989 */
990 /*----------------------------------------------------------------------------*/
991
992 void
cs_turbulence_bc_set_uninit_inlet_k_eps(cs_lnum_t face_id,double k,double eps,double vel_dir[],double shear_dir[],double * rcodcl)993 cs_turbulence_bc_set_uninit_inlet_k_eps(cs_lnum_t face_id,
994 double k,
995 double eps,
996 double vel_dir[],
997 double shear_dir[],
998 double *rcodcl)
999 {
1000 _set_uninit_inlet_bc(face_id, k, eps, vel_dir, shear_dir, rcodcl);
1001 }
1002
1003 /*----------------------------------------------------------------------------*/
1004 /*!
1005 * \brief Compute matrix \f$\tens{\alpha}\f$ used in the computation of the
1006 * Reynolds stress tensor boundary conditions.
1007 *
1008 * We note \f$\tens{R}_g\f$ the Reynolds Stress tensor in the global reference
1009 * frame (mesh reference frame) and \f$\tens{R}_l\f$ the Reynolds stress
1010 * tensor in the local reference frame (reference frame associated to the
1011 * boundary face).
1012 *
1013 * \f$\tens{P}_{lg}\f$ is the change of basis orthogonal matrix from local
1014 * to global reference frame.
1015
1016 * \f$\tens{\alpha}\f$ is a 6 by 6 matrix defined such that:
1017 * \f[
1018 * \vect{R}_{g,\fib} = \tens{\alpha} \vect{R}_{g,\centip} + \vect{R}_{g}^*
1019 * \f]
1020 * where symetric tensors \f$\tens{R}_g\f$ have been unfolded as follows:
1021 * \f[
1022 * \vect{R}_g = \transpose{\left(R_{g,11},R_{g,22},R_{g,33},
1023 * R_{g,12},R_{g,13},R_{g,23}\right)}
1024 * \f].
1025 *
1026 * \f$ \tens{R}_{g,\fib} \f$ should be computed
1027 * as a function of \f$\tens{R}_{g,\centip}\f$ as follows:
1028 * \f[
1029 * \tens{R}_{g,\fib}=\tens{P}_{lg}\tens{R}_{l,\fib}\transpose{\tens{P}_{lg}}
1030 * \f]
1031 *
1032 * with
1033 * \f[
1034 * \tens{R}_{l,\fib} =
1035 * \begin{bmatrix}
1036 * R_{l,11,\centip} & 0 & c R_{l,13,\centip}\\
1037 * 0 & R_{l,22,\centip} & 0 \\
1038 * c R_{l,13,\centip} & 0 & R_{l,33,\centip}
1039 * \end{bmatrix} +
1040 * \underbrace{\begin{bmatrix}
1041 * 0 & (1-c) u^* u_k & 0 \\
1042 * (1-c) u^* u_k & 0 & 0 \\
1043 * 0 & 0 & 0
1044 * \end{bmatrix}}_{\vect{R}_l^*}
1045 * \f]
1046 *
1047 * and
1048 * \f$\tens{R}_{l,\centip}=\transpose{\tens{P}_{lg}}\tens{R}_{g,\centip}
1049 * \tens{P}_{lg}\f$.
1050 *
1051 * Constant c is chosen depending on the type of the boundary face:
1052 * \f$c = 0\f$ at a wall face, \f$c = 1\f$ at a symmetry face.
1053 *
1054 *
1055 *
1056 * \param[in] is_sym Constant c in description above
1057 * (1 at a symmetry face, 0 at a wall face)
1058 * \param[in] p_lg change of basis matrix (local to global)
1059 * \param[out] alpha transformation matrix
1060 */
1061 /*----------------------------------------------------------------------------*/
1062
1063 void
cs_turbulence_bc_rij_transform(int is_sym,cs_real_t p_lg[3][3],cs_real_t alpha[6][6])1064 cs_turbulence_bc_rij_transform(int is_sym,
1065 cs_real_t p_lg[3][3],
1066 cs_real_t alpha[6][6])
1067 {
1068 cs_real_t p_lg2[3][3];
1069 for (int ii = 0; ii < 3; ii++)
1070 for (int jj = 0; jj < 3; jj++)
1071 p_lg2[ii][jj] = cs_math_pow2(p_lg[ii][jj]);
1072
1073 /* alpha(i,j) for i in [1,3] and j in [1,3]: 9 terms */
1074 for (int ii = 0; ii < 3; ii++) {
1075 for (int jj = 0; jj < 3; jj++) {
1076 alpha[jj][ii] = p_lg2[0][ii] * p_lg2[0][jj]
1077 + p_lg2[1][ii] * p_lg2[1][jj]
1078 + p_lg2[2][ii] * p_lg2[2][jj]
1079 + 2. * is_sym * p_lg[0][ii] * p_lg[2][ii]
1080 * p_lg[0][jj] * p_lg[2][jj];
1081 }
1082 }
1083
1084 /* alpha(i,j) for i in [1,3] and j in [4,6]: 9 terms */
1085
1086 /* Correcponding line to j */
1087 const int _jj_to_kk[3] = {0, 1, 0};
1088 /* Corresponding column to j*/
1089 const int _jj_to_pp[3] = {1, 2, 2};
1090
1091 for (int ii = 0; ii < 3; ii++) {
1092 for (int jj = 0; jj < 3; jj++) {
1093
1094 int kk = _jj_to_kk[jj];
1095 int pp = _jj_to_pp[jj];
1096
1097 alpha[jj + 3][ii] =
1098 2. * ( p_lg2[0][ii] * p_lg[0][kk] * p_lg[0][pp]
1099 + p_lg2[1][ii] * p_lg[1][kk] * p_lg[1][pp]
1100 + p_lg2[2][ii] * p_lg[2][kk] * p_lg[2][pp]
1101 + is_sym * p_lg[2][ii] * p_lg[0][ii]
1102 * ( p_lg[0][kk]*p_lg[2][pp]
1103 + p_lg[2][kk]*p_lg[0][pp]));
1104 }
1105 }
1106
1107 /* alpha(i,j) for i in [4,6] and j in [1,3]: 9 terms */
1108
1109 for (int ii = 0; ii < 3; ii++) {
1110 for (int jj = 0; jj < 3; jj++) {
1111
1112 int kk = _jj_to_kk[ii];
1113 int pp = _jj_to_pp[ii];
1114
1115 /* Note: could be simplified because it is 0.5*alpha[jj+3][ii] */
1116 alpha[jj][ii + 3] =
1117 p_lg[0][kk] * p_lg[0][pp] * p_lg2[0][jj]
1118 + p_lg[1][kk] * p_lg[1][pp] * p_lg2[1][jj]
1119 + p_lg[2][kk] * p_lg[2][pp] * p_lg2[2][jj]
1120 + is_sym * p_lg[2][jj] * p_lg[0][jj]
1121 * ( p_lg[0][kk]*p_lg[2][pp]
1122 + p_lg[2][kk]*p_lg[0][pp]);
1123 }
1124 }
1125
1126 /* alpha(i,j) for i in [4,6] and j in [4,6]: 9 terms */
1127 for (int ii = 0; ii < 3; ii++) {
1128 for (int jj = 0; jj < 3; jj++) {
1129
1130 int kk = _jj_to_kk[ii];
1131 int pp = _jj_to_pp[ii];
1132
1133 int jj1 = _jj_to_kk[jj];
1134 int jj2 = _jj_to_pp[jj];
1135
1136 alpha[jj + 3][ii + 3] =
1137 2. * ( p_lg[0][kk] * p_lg[0][pp] * p_lg[0][jj1] * p_lg[0][jj2]
1138 + p_lg[1][kk] * p_lg[1][pp] * p_lg[1][jj1] * p_lg[1][jj2]
1139 + p_lg[2][kk] * p_lg[2][pp] * p_lg[2][jj1] * p_lg[2][jj2])
1140 + is_sym * ( p_lg[0][kk]*p_lg[2][pp]
1141 + p_lg[2][kk]*p_lg[0][pp])
1142 * ( p_lg[2][jj1]*p_lg[0][jj2]
1143 + p_lg[0][jj1]*p_lg[2][jj2]);
1144 }
1145 }
1146 }
1147
1148 /*----------------------------------------------------------------------------*/
1149
1150 END_C_DECLS
1151