TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Calcul_Production_K_VEF.cpp
1/****************************************************************************
2* Copyright (c) 2015 - 2016, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <Convection_Diffusion_Concentration.h>
17#include <Calcul_Production_K_VEF.h>
18#include <TRUSTTab.h>
19#include <Domaine_VEF.h>
20#include <Domaine_Cl_VEF.h>
21#include <TRUSTVect.h>
22#include <TRUSTTrav.h>
23#include <Champ_P1NC.h>
24#include <Periodique.h>
25#include <Champ_Uniforme.h>
26#include <Champ_Don_base.h>
27#include <K_Omega_Eps_constants.h>
28
29/*! @brief Compute the production term for the turbulent kinetic energy.
30 *
31 * The total production term writes
32 * \f$ P = -\frac{2}{3}k\text{div}(U) + 2\nu_t S_{ij}S_{ij}\f$
33 * Being in a incompressible flow, the first part is not computed.
34 *
35 * Like the TKE and epsilon, P is discretised on face centers.
36 *
37 * @param[in] const Domaine_VEF& domaine_VEF,
38 * @param[in] const Domaine_Cl_VEF& zcl_VEF,
39 * @param[in] DoubleTab& prodK,
40 * @param[in] const DoubleTab& K_eps,
41 * @param[in] const DoubleTab& vit,
42 * @param[in] const DoubleTab& visco_turb,
43 * @param[in] const int& interpol_visco,
44 * @param[in] const double& limiteur
45 * @param[out]
46 * @return DoubleTab& prodK
47 */
48DoubleTab& Calcul_Production_K_VEF::calculer_terme_production_K(const Domaine_VEF& domaine_VEF,
49 const Domaine_Cl_VEF& zcl_VEF,
50 DoubleTab& prodK,
51 const DoubleTab& K_eps,
52 const DoubleTab& vit,
53 const DoubleTab& visco_turb,
54 const int& interpol_visco,
55 const double& limiteur,
56 const bool& deactivate_production_limiter,
57 const double& cst_production_limiter) const
58{
59 prodK = 0;
60
61 // Compute the velocity gradient
62 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
63 const int dimension = Objet_U::dimension;
64 DoubleTrav gradient_elem(nb_elem_tot, dimension, dimension);
65 Champ_P1NC::calcul_gradient(vit, gradient_elem, zcl_VEF);
66
67 const IntTab& face_voisins = domaine_VEF.face_voisins();
68 const DoubleVect& volumes = domaine_VEF.volumes();
69
70 // Boundary conditions
71 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
72 {
73 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
74 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
75 const int first_boundary_face = le_bord.num_premiere_face();
76 const int last_boundary_face = first_boundary_face + le_bord.nb_faces();
77
78 if (sub_type(Periodique, la_cl.valeur()))
79 loop_for_internal_or_periodic_faces(prodK, gradient_elem, visco_turb, volumes,
80 face_voisins, first_boundary_face, last_boundary_face,
81 interpol_visco, limiteur, K_eps,deactivate_production_limiter,cst_production_limiter);
82 else
83 loop_for_non_periodic_boundaries(prodK, gradient_elem, visco_turb, volumes,
84 face_voisins, first_boundary_face, last_boundary_face,
85 interpol_visco, limiteur, K_eps,deactivate_production_limiter,cst_production_limiter);
86 }
87
88 // Internal faces
89 const int premiere_face_int = domaine_VEF.premiere_face_int();
90 const int nb_faces_ = domaine_VEF.nb_faces();
91 loop_for_internal_or_periodic_faces(prodK, gradient_elem, visco_turb, volumes,
92 face_voisins, premiere_face_int, nb_faces_,
93 interpol_visco, limiteur,K_eps,deactivate_production_limiter,cst_production_limiter);
94 return prodK;
95}
96
97/*! @brief Compute production term on non periodic boundary faces
98 *
99 * Using the velocity gradient, the loop computes the production term written
100 * \f$P = 2 \nu_t S_{ij}S_{ij}\f$
101 * with \f$S_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})\f$.
102 *
103 * @param[in] DoubleTab& prodK
104 * @param[in] const DoubleTab& gradient_elem,
105 * @param[in] const DoubleTab& visco_turb,
106 * @param[in] const DoubleVect& volumes,
107 * @param[in] const IntTab& face_voisins,
108 * @param[in] const int nfaceinit,
109 * @param[in] const int nfaceend,
110 * @param[in] const int interpol_visco,
111 * @param[in] const double limiteur
112 * @return
113 */
114void Calcul_Production_K_VEF::loop_for_non_periodic_boundaries(DoubleTab& tab_prodK,
115 const DoubleTab& tab_gradient_elem,
116 const DoubleTab& tab_visco_turb,
117 const DoubleVect& volumes,
118 const IntTab& tab_face_voisins,
119 const int nfaceinit,
120 const int nfaceend,
121 const int interpol_visco,
122 const double limiteur,
123 const DoubleTab& tab_K_Omega,
124 const bool& deactivate_production_limiter,
125 const double& cst_production_limiter
126 ) const
127{
128 int dim = Objet_U::dimension;
129 CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
130 CDoubleTabView K_Omega = tab_K_Omega.view_ro();
131 CDoubleArrView visco_turb = static_cast<const ArrOfDouble&>(tab_visco_turb).view_ro();
132 CIntTabView face_voisins = tab_face_voisins.view_ro();
133 DoubleArrView prodK = static_cast<ArrOfDouble&>(tab_prodK).view_wo();
134 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(nfaceinit, nfaceend), KOKKOS_LAMBDA(const int fac)
135 {
136 const int poly1 = face_voisins(fac, 0);
137 const double visco_face = visco_turb(poly1);
138
139 const double du_dx = gradient_elem(poly1, 0, 0);
140 const double du_dy = gradient_elem(poly1, 0, 1);
141 const double dv_dx = gradient_elem(poly1, 1, 0);
142 const double dv_dy = gradient_elem(poly1, 1, 1);
143
144 // Determination du terme de production
145 prodK(fac) = (2*(du_dx*du_dx + dv_dy*dv_dy)
146 + ((du_dy + dv_dx)*(du_dy + dv_dx)))*visco_face;
147
148 if (dim==3)
149 {
150 const double du_dz = gradient_elem(poly1, 0, 2);
151 const double dv_dz = gradient_elem(poly1, 1, 2);
152 const double dw_dx = gradient_elem(poly1, 2, 0);
153 const double dw_dy = gradient_elem(poly1, 2, 1);
154 const double dw_dz = gradient_elem(poly1, 2, 2);
155
156 // Determination du terme de production
157 prodK(fac) = (2*(du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz)
158 + ((du_dy + dv_dx) * (du_dy + dv_dx)
159 + (du_dz + dw_dx) * (du_dz + dw_dx)
160 + (dw_dy + dv_dz) * (dw_dy + dv_dz)))*visco_face;
161 }
162 if (!deactivate_production_limiter)
163 prodK(fac)=std::min(prodK(fac),cst_production_limiter*BETA_K*K_Omega(fac,0)*K_Omega(fac,1));
164 });
165 end_gpu_timer(__KERNEL_NAME__);
166}
167
168/*! @brief Compute production term on internal and periodic boundary faces
169 *
170 * Using the velocity gradient, the loop computes the production term written
171 * \f$P = 2 \nu_t S_{ij}S_{ij}\f$
172 * with \f$S_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})\f$.
173 *
174 * @param[in] DoubleTab& prodK
175 * @param[in] const DoubleTab& gradient_elem,
176 * @param[in] const DoubleTab& visco_turb,
177 * @param[in] const DoubleVect& volumes,
178 * @param[in] const IntTab& face_voisins,
179 * @param[in] const int nfaceinit,
180 * @param[in] const int nfaceend,
181 * @param[in] const int interpol_visco,
182 * @param[in] const double limiteur
183 * @return
184 */
185void Calcul_Production_K_VEF::loop_for_internal_or_periodic_faces(DoubleTab& tab_prodK,
186 const DoubleTab& tab_gradient_elem,
187 const DoubleTab& tab_visco_turb,
188 const DoubleVect& tab_volumes,
189 const IntTab& tab_face_voisins,
190 const int nfaceinit,
191 const int nfaceend,
192 const int interpol_visco,
193 const double limiteur,
194 const DoubleTab& tab_K_Omega,
195 const bool& deactivate_production_limiter,
196 const double& cst_production_limiter
197 ) const
198{
199 int dim = Objet_U::dimension;
200 CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
201 CDoubleTabView K_Omega = tab_K_Omega.view_ro();
202 CDoubleArrView visco_turb = static_cast<const ArrOfDouble&>(tab_visco_turb).view_ro();
203 CDoubleArrView volumes = tab_volumes.view_ro();
204 CIntTabView face_voisins = tab_face_voisins.view_ro();
205 DoubleArrView prodK = static_cast<ArrOfDouble&>(tab_prodK).view_wo();
206 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(nfaceinit, nfaceend), KOKKOS_LAMBDA(const int fac)
207 {
208 const int poly1 = face_voisins(fac, 0);
209 const int poly2 = face_voisins(fac, 1);
210 const double a = volumes(poly1)/(volumes(poly1) + volumes(poly2));
211 const double b = volumes(poly2)/(volumes(poly1) + volumes(poly2));
212
213 const double visco_face = get_turbulent_viscosity(visco_turb, volumes,
214 interpol_visco, poly1, poly2,
215 limiteur);
216
217 const double du_dx = a*gradient_elem(poly1, 0, 0) + b*gradient_elem(poly2, 0, 0);
218 const double du_dy = a*gradient_elem(poly1, 0, 1) + b*gradient_elem(poly2, 0, 1);
219 const double dv_dx = a*gradient_elem(poly1, 1, 0) + b*gradient_elem(poly2, 1, 0);
220 const double dv_dy = a*gradient_elem(poly1, 1, 1) + b*gradient_elem(poly2, 1, 1);
221
222 // Determination du terme de production
223
224 prodK(fac) = (2*(du_dx*du_dx + dv_dy*dv_dy)
225 + ((du_dy + dv_dx)*(du_dy + dv_dx)))*visco_face;
226
227 if (dim==3)
228 {
229 const double du_dz = a*gradient_elem(poly1, 0, 2) + b*gradient_elem(poly2, 0, 2);
230 const double dv_dz = a*gradient_elem(poly1, 1, 2) + b*gradient_elem(poly2, 1, 2);
231 const double dw_dx = a*gradient_elem(poly1, 2, 0) + b*gradient_elem(poly2, 2, 0);
232 const double dw_dy = a*gradient_elem(poly1, 2, 1) + b*gradient_elem(poly2, 2, 1);
233 const double dw_dz = a*gradient_elem(poly1, 2, 2) + b*gradient_elem(poly2, 2, 2);
234
235 // Determination du terme de production
236 prodK(fac) = (2*(du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz)
237 + ((du_dy + dv_dx)*(du_dy + dv_dx)
238 + (du_dz + dw_dx)*(du_dz + dw_dx)
239 + (dw_dy + dv_dz)*(dw_dy + dv_dz)))*visco_face;
240 }
241 if (!deactivate_production_limiter)
242 prodK(fac)=std::min(prodK(fac),cst_production_limiter*BETA_K*K_Omega(fac,0)*K_Omega(fac,1));
243 });
244 end_gpu_timer(__KERNEL_NAME__);
245}
246
247DoubleTab& Calcul_Production_K_VEF::
248calculer_terme_production_K_BiK(const Domaine_VEF& domaine_VEF,const Domaine_Cl_VEF& zcl_VEF,
249 DoubleTab& P,const DoubleTab& K,const DoubleTab& eps,
250 const DoubleTab& vit,const DoubleTab& visco_turb, const int& interpol_visco, const double& limiteur) const
251{
252 // P est discretise comme K et Eps i.e au centre des faces
253 //
254 // P(elem) = -(2/3)*k(i)*div_U(i) + nu_t(i) * F(u,v,w)
255 //
256 // 2 2 2
257 // avec F(u,v,w) = 2[(du/dx) + (dv/dy) + (dw/dz) ] +
258 //
259 // 2 2 2
260 // (du/dy+dv/dx) + (du/dz+dw/dx) + (dw/dy+dv/dz)
261 //
262 // Rqs: On se place dans le cadre incompressible donc on neglige
263 // le terme (2/3)*k(i)*div_U(i)
264
265 P= 0;
266
267 // Calcul de F(u,v,w):
268 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
269 const IntTab& face_voisins = domaine_VEF.face_voisins();
270 const DoubleVect& volumes = domaine_VEF.volumes();
271 const int dimension = Objet_U::dimension;
272
273 DoubleTab gradient_elem(nb_elem_tot, dimension, dimension);
274 gradient_elem = 0.;
275
276 ///////////////////////////////////////////////////////////////////////////////////////////////
277 // <
278 // calcul des gradients; < [ Ujp*np/vol(j) ]
279 // j
280 ////////////////////////////////////////////////////////////////////////////////////////////////
281
282 Champ_P1NC::calcul_gradient(vit, gradient_elem, zcl_VEF);
283
284 // Calcul des du/dx dv/dy et des derivees croisees sur les faces de chaque elements dans le cas 2D
285
286 // Boucle sur les bords pour traiter les faces de bord
287 // en distinguant le cas periodicite
288 ToDo_Kokkos("critical");
289 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
290 {
291 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
292 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
293 const int ndeb = le_bord.num_premiere_face();
294 const int nfin = ndeb + le_bord.nb_faces();
295
296 if (sub_type(Periodique, la_cl.valeur()))
297 {
298 for (int fac = ndeb; fac < nfin; fac++)
299 {
300 const int poly1 = face_voisins(fac, 0);
301 const int poly2 = face_voisins(fac, 1);
302 const double a = volumes(poly1)/(volumes(poly1) + volumes(poly2));
303 const double b = volumes(poly2)/(volumes(poly1) + volumes(poly2));
304
305 const double visco_face = get_turbulent_viscosity(visco_turb, volumes,
306 interpol_visco, poly1, poly2,
307 limiteur);
308
309 const double du_dx = a*gradient_elem(poly1, 0, 0) + b*gradient_elem(poly2, 0, 0);
310 const double du_dy = a*gradient_elem(poly1, 0, 1) + b*gradient_elem(poly2, 0, 1);
311 const double dv_dx = a*gradient_elem(poly1, 1, 0) + b*gradient_elem(poly2, 1, 0);
312 const double dv_dy = a*gradient_elem(poly1, 1, 1) + b*gradient_elem(poly2, 1, 1);
313
314 P(fac) = visco_face*( 2*(du_dx*du_dx + dv_dy*dv_dy)
315 + ((du_dy+dv_dx)*(du_dy+dv_dx) ) );
316
317 if (dimension == 3)
318 {
319 const double du_dz = a*gradient_elem(poly1, 0, 2) + b*gradient_elem(poly2, 0, 2);
320 const double dv_dz = a*gradient_elem(poly1, 1, 2) + b*gradient_elem(poly2, 1, 2);
321 const double dw_dx = a*gradient_elem(poly1, 2, 0) + b*gradient_elem(poly2, 2, 0);
322 const double dw_dy = a*gradient_elem(poly1, 2, 1) + b*gradient_elem(poly2, 2, 1);
323 const double dw_dz = a*gradient_elem(poly1, 2, 2) + b*gradient_elem(poly2, 2, 2);
324
325 // Determination du terme de production
326 P(fac) = (2*(du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz)
327 + ((du_dy + dv_dx)*(du_dy + dv_dx)
328 + (du_dz + dw_dx)*(du_dz + dw_dx)
329 + (dw_dy + dv_dz)*(dw_dy + dv_dz)))*visco_face;
330 }
331 }
332 }
333 else
334 {
335 for (int fac = ndeb; fac < nfin; fac++)
336 {
337 const int poly1 = face_voisins(fac,0);
338 const double visco_face = visco_turb(poly1);
339
340 const double du_dx = gradient_elem(poly1, 0, 0);
341 const double du_dy = gradient_elem(poly1, 0, 1);
342 const double dv_dx = gradient_elem(poly1, 1, 0);
343 const double dv_dy = gradient_elem(poly1, 1, 1);
344
345 P(fac) = visco_face*( 2*(du_dx*du_dx + dv_dy*dv_dy) + ((du_dy+dv_dx)*(du_dy+dv_dx)));
346
347 if (dimension == 3)
348 {
349 const double du_dz = gradient_elem(poly1, 0, 2);
350 const double dv_dz = gradient_elem(poly1, 1, 2);
351 const double dw_dx = gradient_elem(poly1, 2, 0);
352 const double dw_dy = gradient_elem(poly1, 2, 1);
353 const double dw_dz = gradient_elem(poly1, 2, 2);
354
355 P(fac) = (2*(du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz)
356 + ( (du_dy + dv_dx)*(du_dy + dv_dx)
357 + (du_dz + dw_dx)*(du_dz + dw_dx)
358 + (dw_dy + dv_dz)*(dw_dy + dv_dz)))*visco_face;
359
360 }
361 }
362 }
363 }
364
365 // Traitement des faces internes
366 const int premiere_face_int = domaine_VEF.premiere_face_int();
367 const int nb_faces_ = domaine_VEF.nb_faces();
368
369 for (int fac = premiere_face_int; fac < nb_faces_; fac++)
370 {
371 const int poly1 = face_voisins(fac, 0);
372 const int poly2 = face_voisins(fac, 1);
373 const double a = volumes(poly1)/(volumes(poly1) + volumes(poly2));
374 const double b = volumes(poly2)/(volumes(poly1) + volumes(poly2));
375
376 const double visco_face = get_turbulent_viscosity(visco_turb, volumes,
377 interpol_visco, poly1, poly2,
378 limiteur);
379
380 const double du_dx = a*gradient_elem(poly1, 0, 0) + b*gradient_elem(poly2, 0, 0);
381 const double du_dy = a*gradient_elem(poly1, 0, 1) + b*gradient_elem(poly2, 0, 1);
382 const double dv_dx = a*gradient_elem(poly1, 1, 0) + b*gradient_elem(poly2, 1, 0);
383 const double dv_dy = a*gradient_elem(poly1, 1, 1) + b*gradient_elem(poly2, 1, 1);
384
385 P(fac) = (2*(du_dx*du_dx + dv_dy*dv_dy)
386 + ((du_dy + dv_dx)*(du_dy + dv_dx)))*visco_face;
387 if (dimension == 3)
388 {
389 const double du_dz = a*gradient_elem(poly1, 0, 2) + b*gradient_elem(poly2, 0, 2);
390 const double dv_dz = a*gradient_elem(poly1, 1, 2) + b*gradient_elem(poly2, 1, 2);
391 const double dw_dx = a*gradient_elem(poly1, 2, 0) + b*gradient_elem(poly2, 2, 0);
392 const double dw_dy = a*gradient_elem(poly1, 2, 1) + b*gradient_elem(poly2, 2, 1);
393 const double dw_dz = a*gradient_elem(poly1, 2, 2) + b*gradient_elem(poly2, 2, 2);
394
395 P(fac) = (2*(du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz)
396 + ((du_dy + dv_dx)*(du_dy + dv_dx)
397 + (du_dz + dw_dx)*(du_dz + dw_dx)
398 + (dw_dy + dv_dz)*(dw_dy + dv_dz)))*visco_face;
399 }
400 }
401 return P;
402}
403
404
405
406DoubleTab& Calcul_Production_K_VEF::calculer_terme_production_K_EASM(const Domaine_VEF& domaine_VEF,
407 const Domaine_Cl_VEF& zcl_VEF,
408 DoubleTab& P,
409 const DoubleTab& K_eps,
410 const DoubleTab& gradient_elem,
411 const DoubleTab& visco_turb,
412 const DoubleTab& Re,
413 const int& interpol_visco,
414 const double& limiteur) const
415{
416 // We treat all exception at the beginning
417 if (interpol_visco != 0 && interpol_visco != 1 && interpol_visco != 2)
418 {
419 Cerr << "interpol_visco must be equal to [0, 1, 2] and its current value is " << interpol_visco << ".\n";
420 Process::exit();
421 }
422
423 // P : Production
424 P = 0;
425
426 // Compute velocity tensor on faces for
427 const int nb_faces_ = domaine_VEF.nb_faces();
428 const int dimension = Objet_U::dimension;
429
430 DoubleTab gradient_face(nb_faces_, dimension, dimension);
431 calcul_tenseur_face(gradient_face, gradient_elem, domaine_VEF, zcl_VEF);
432
433 DoubleTab Re_face(nb_faces_, dimension, dimension);
434 calcul_tenseur_face(Re_face, Re, domaine_VEF, zcl_VEF);
435
436 const IntTab& face_voisins = domaine_VEF.face_voisins();
437 const DoubleVect& volumes = domaine_VEF.volumes();
438
439 // Loop on boundary conditions
440 ToDo_Kokkos("critical");
441 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
442 {
443 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
444 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
445 const int ndeb = le_bord.num_premiere_face();
446 const int nfin = ndeb + le_bord.nb_faces();
447
448 if (sub_type(Periodique, la_cl.valeur()))
449 {
450 for (int fac = ndeb; fac < nfin; fac++)
451 {
452 const int poly1 = face_voisins(fac, 0);
453 const int poly2 = face_voisins(fac, 1);
454 const double visco_face = get_turbulent_viscosity(visco_turb, volumes,
455 interpol_visco, poly1, poly2,
456 limiteur);
457 compute_production_term_EASM(fac, visco_face, Re_face, gradient_face, P);
458 }
459 }
460 else
461 {
462 for (int fac = ndeb; fac < nfin; fac++)
463 {
464 const int poly1 = face_voisins(fac, 0);
465 const double visco_face = visco_turb(poly1);
466 compute_production_term_EASM(fac, visco_face, Re_face, gradient_face, P);
467 }
468 }
469 }
470
471 // Loop on internal faces
472 const int first_internal_face = domaine_VEF.premiere_face_int();
473 for (int fac = first_internal_face; fac < nb_faces_; fac++)
474 {
475 const int poly1 = face_voisins(fac, 0);
476 const int poly2 = face_voisins(fac, 1);
477 const double visco_face = get_turbulent_viscosity(visco_turb, volumes,
478 interpol_visco, poly1, poly2,
479 limiteur);
480 compute_production_term_EASM(fac, visco_face, Re_face, gradient_face, P);
481 }
482 return P;
483}
484
485void Calcul_Production_K_VEF::compute_production_term_EASM(const int face,
486 const double visco_face,
487 const DoubleTab& Re_face,
488 const DoubleTab& gradient_face,
489 DoubleTab& P) const
490{
491 for (int i = 0; i < Objet_U::dimension; i++)
492 for (int j = 0; j < Objet_U::dimension; j++)
493 P(face) += Re_face(face, i, j) * gradient_face(face, i, j);
494 P(face) *= visco_face;
495}
496
497/*! @brief Get the turbulent viscosity depending on the interpolation choice.
498 *
499 * if type_interpo == 0: arithmetic interpolation (instable)
500 * if type_interpo == 1: harmonic mean (used for the realisable k-epsilon)
501 * if type_interpo == 2: weighted harmonic mean (used for the realisable k-epsilon)
502 *
503 * @param[in] const DoubleTab& visco_turb
504 * @param[in] const DoubleVect& volumes
505 * @param[in] const int type_interpo
506 * @param[in] const int poly1
507 * @param[in] const int poly2
508 * @param[in] const double limiteur
509 * @param[out]
510 * @return double
511 */
512KOKKOS_INLINE_FUNCTION
513double Calcul_Production_K_VEF::get_turbulent_viscosity(CDoubleArrView visco_turb,
514 CDoubleArrView volumes,
515 const int type_interpo,
516 const int poly1,
517 const int poly2,
518 const double limiteur) const
519{
520 switch(type_interpo)
521 {
522 case 0: // cas initial, toutes les interpolation arithmetiques ponderees sont fortement instables
523 return 0.5*(visco_turb(poly1) + visco_turb(poly2));
524
525 case 1: //Moyenne harmonique (uniquement utilisee dans le cas du keps realisable)
526 if (visco_turb(poly1) > 1.e-10 && visco_turb(poly2) > 1.e-10)
527 return limiteur/(1./visco_turb(poly1) + 1./visco_turb(poly2));
528 else
529 return limiteur*0.5*(visco_turb(poly1) + visco_turb(poly2));
530
531 case 2: //Moyenne harmonique ponderee pour garantir la continuite du tenseur des contraintes a la face (uniquement utilisee dans le cas du keps realisable)
532 if (visco_turb(poly1) > 1.e-10 && visco_turb(poly2) > 1.e-10)
533 {
534 const double a = volumes(poly1)/(volumes(poly1) + volumes(poly2));
535 const double b = volumes(poly2)/(volumes(poly1) + volumes(poly2));
536 return limiteur*(visco_turb(poly1)*visco_turb(poly2))/(b*visco_turb(poly1) + a*visco_turb(poly2));
537 }
538 else
539 return limiteur*0.5*(visco_turb(poly1) + visco_turb(poly2));
540
541 default:
542 printf("Interpolation type: %d\n",type_interpo);
543 Process::Kokkos_exit("It should be 0, 1 or 2");
544 return 1;
545 }
546}
547// Soon obsolete, once all is ported on Kokkos
548double Calcul_Production_K_VEF::get_turbulent_viscosity(const DoubleTab& visco_turb,
549 const DoubleVect& volumes,
550 const int type_interpo,
551 const int poly1,
552 const int poly2,
553 const double limiteur) const
554{
555 switch(type_interpo)
556 {
557 case 0: // cas initial, toutes les interpolation arithmetiques ponderees sont fortement instables
558 return 0.5*(visco_turb(poly1) + visco_turb(poly2));
559
560 case 1: //Moyenne harmonique (uniquement utilisee dans le cas du keps realisable)
561 if (visco_turb(poly1) > 1.e-10 && visco_turb(poly2) > 1.e-10)
562 return limiteur/(1./visco_turb(poly1) + 1./visco_turb(poly2));
563 else
564 return limiteur*0.5*(visco_turb(poly1) + visco_turb(poly2));
565
566 case 2: //Moyenne harmonique ponderee pour garantir la continuite du tenseur des contraintes a la face (uniquement utilisee dans le cas du keps realisable)
567 if (visco_turb(poly1) > 1.e-10 && visco_turb(poly2) > 1.e-10)
568 {
569 const double a = volumes(poly1)/(volumes(poly1) + volumes(poly2));
570 const double b = volumes(poly2)/(volumes(poly1) + volumes(poly2));
571 return limiteur*(visco_turb(poly1)*visco_turb(poly2))/(b*visco_turb(poly1) + a*visco_turb(poly2));
572 }
573 else
574 return limiteur*0.5*(visco_turb(poly1) + visco_turb(poly2));
575
576 default:
577 Cerr << "Interpolation type should be 0, 1 or 2, not " << type_interpo << "\n";
578 Process::exit();
579 return 1;
580 }
581}
582
583// Utilisée dans le bas Reynolds anisotherme
584DoubleTab& Calcul_Production_K_VEF::calculer_terme_destruction_K_gen(
585 const Domaine_VEF& domaine_VEF,
586 const Domaine_Cl_VEF& zcl_VEF,
587 DoubleTab& G,
588 const DoubleTab& inconnue, // scalar, concentration, temperature, etc.
589 const DoubleTab& alpha_turb,
590 const Champ_Don_base& ch_beta,
591 const DoubleVect& gravite,
592 int nb_consti) const
593{
594 // G est discretise comme K et Eps i.e au centre des faces
595 // G(face) = beta alpha_t(face) G . gradT(face)
596
597 if (nb_consti < 0)
598 Process::exit("nb_consti is negative");
599
600 const int nb_compo = sub_type(Champ_Uniforme, ch_beta) ? 0 : ch_beta.nb_comp();
601
602 if (nb_compo < 0)
603 Process::exit("nb_compo is negative");
604
605 const int dimension = Objet_U::dimension;
606
607 const IntTab& face_voisins = domaine_VEF.face_voisins();
608 const DoubleVect& volumes = domaine_VEF.volumes();
609 const DoubleTab& tab_beta = ch_beta.valeurs();
610
611 G = 0;
612
613 if (nb_consti == 0 || nb_consti == 1)
614 {
615 // Compute the gradient of the equation unknwon
616 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
617 DoubleTrav gradient_elem(nb_elem_tot, dimension);
618 gradient_elem = 0;
619 Champ_P1NC::calcul_gradient(inconnue, gradient_elem, zcl_VEF);
620
621 // Compute u_theta
622 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
623 DoubleTrav u_theta(nb_faces_tot, dimension);
624 u_theta = 0;
625
626 if (nb_compo == 0)
627 compute_utheta_nbConsti_le_1_nbCompo_eq_0(domaine_VEF, zcl_VEF, face_voisins,
628 volumes, tab_beta, alpha_turb, gradient_elem,
629 u_theta);
630 else if (nb_compo == 1)
631 compute_utheta_nbConsti_le_1_nbCompo_eq_1(domaine_VEF, zcl_VEF, face_voisins,
632 volumes, tab_beta, alpha_turb, gradient_elem,
633 u_theta);
634
635 else if (nb_compo > 1)
636 compute_utheta_nbConsti_le_1_nbCompo_gt_1(domaine_VEF, zcl_VEF, face_voisins,
637 volumes, tab_beta, alpha_turb, gradient_elem,
638 u_theta);
639
640
641 // Calcul de gravite . u_teta
642 const int nb_faces_ = domaine_VEF.nb_faces();
643 ToDo_Kokkos("critical");
644 for (int fac = 0; fac < nb_faces_; fac++)
645 {
646 G[fac] = gravite(0)*u_theta(fac, 0) + gravite(1)*u_theta(fac, 1);
647 if (dimension == 3)
648 G[fac] += gravite(2)*u_theta(fac, 2);
649 }
650 }
651 else if (nb_consti > 1)
652 {
653 // Compute the gradient of the equation unknwon
654 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
655 DoubleTrav gradient_elem(nb_elem_tot, nb_consti, dimension);
656 gradient_elem = 0;
657 Champ_P1NC::calcul_gradient(inconnue, gradient_elem, zcl_VEF);
658
659 // Compute u_theta
660 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
661 DoubleTrav u_theta(nb_faces_tot, nb_consti, dimension);
662 u_theta = 0;
663
664 if (nb_compo == 0)
665 compute_utheta_nbConsti_gt_1_nbCompo_eq_0(domaine_VEF, zcl_VEF, face_voisins,
666 volumes, tab_beta, alpha_turb, gradient_elem,
667 nb_consti, u_theta);
668 else if (nb_compo == 1)
669 compute_utheta_nbConsti_gt_1_nbCompo_eq_1(domaine_VEF, zcl_VEF, face_voisins,
670 volumes, tab_beta, alpha_turb, gradient_elem,
671 nb_consti, u_theta);
672 else if (nb_compo > 1)
673 compute_utheta_nbConsti_gt_1_nbCompo_gt_1(domaine_VEF, zcl_VEF, face_voisins,
674 volumes, tab_beta, alpha_turb, gradient_elem,
675 nb_consti, u_theta);
676
677 // Calcul de gravite . u_teta
678 const int nb_faces_ = domaine_VEF.nb_faces();
679 ToDo_Kokkos("critical");
680 for (int fac = 0; fac < nb_faces_; fac++)
681 for (int k = 0; k < nb_consti; k++)
682 {
683 G[fac] = gravite(0)*u_theta(fac, k, 0) + gravite(1)*u_theta(fac, k, 1);
684 if (dimension == 3)
685 G[fac] += gravite(2)*u_theta(fac, k, 2);
686 }
687 }
688
689 return G;
690}
691
692void Calcul_Production_K_VEF::compute_utheta_nbConsti_le_1_nbCompo_eq_0(const Domaine_VEF& domaine_VEF,
693 const Domaine_Cl_VEF& zcl_VEF,
694 const IntTab& face_voisins,
695 const DoubleVect& volumes,
696 const DoubleTab& tab_beta,
697 const DoubleTab& alpha_turb,
698 const DoubleTrav& gradient_elem,
699 DoubleTrav& u_theta) const
700{
701 // we treat the boundaries
702 ToDo_Kokkos("critical");
703 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
704 {
705 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
706 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
707 const int ndeb = le_bord.num_premiere_face();
708 const int nfin = ndeb + le_bord.nb_faces();
709
710 if (sub_type(Periodique, la_cl.valeur()))
711 {
712 for (int fac = ndeb ; fac < nfin; fac++)
713 {
714 const int elem1 = face_voisins(fac, 0);
715 const int elem2 = face_voisins(fac, 1);
716 const double invVol = 1/(volumes(elem1) + volumes(elem2));
717 const double a = volumes(elem1)*invVol;
718 const double b = volumes(elem2)*invVol;
719 for (int i = 0; i < Objet_U::dimension; i++)
720 u_theta(fac, i) = a*tab_beta(0, 0)*alpha_turb(elem1)*gradient_elem(elem1, i)
721 + b*tab_beta(0, 0)*alpha_turb(elem2)*gradient_elem(elem2, i);
722 }
723 }
724 else
725 {
726 for (int fac = ndeb; fac < nfin; fac++)
727 {
728 const int elem1 = face_voisins(fac, 0);
729 for (int i = 0; i < Objet_U::dimension; i++)
730 u_theta(fac, i) = tab_beta(0, 0)*alpha_turb(elem1)*gradient_elem(elem1, i);
731 }
732 }
733 }
734
735 // we treat the internal faces
736 ToDo_Kokkos("critical");
737 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
738 for (int fac = 0; fac < nb_faces_tot; fac++)
739 {
740 const int elem1 = face_voisins(fac, 0);
741 const int elem2 = face_voisins(fac, 1);
742
743 if ((elem1 >= 0) && (elem2 >= 0))
744 {
745 const double a = volumes(elem1)/(volumes(elem1) + volumes(elem2));
746 const double b = volumes(elem2)/(volumes(elem1) + volumes(elem2));
747
748 for (int i = 0; i < Objet_U::dimension; i++)
749 u_theta(fac, i) = a*tab_beta(0, 0)*alpha_turb(elem1)*gradient_elem(elem1, i)
750 + b*tab_beta(0, 0)*alpha_turb(elem2)*gradient_elem(elem2, i);
751 }
752 }
753}
754
755void Calcul_Production_K_VEF::compute_utheta_nbConsti_le_1_nbCompo_eq_1(const Domaine_VEF& domaine_VEF,
756 const Domaine_Cl_VEF& zcl_VEF,
757 const IntTab& face_voisins,
758 const DoubleVect& volumes,
759 const DoubleTab& tab_beta,
760 const DoubleTab& alpha_turb,
761 const DoubleTrav& gradient_elem,
762 DoubleTrav& u_theta) const
763{
764 // we treat the boundaries
765 ToDo_Kokkos("critical");
766 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
767 {
768 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
769 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
770 const int ndeb = le_bord.num_premiere_face();
771 const int nfin = ndeb + le_bord.nb_faces();
772
773 if (sub_type(Periodique, la_cl.valeur()))
774 {
775 for (int fac = ndeb; fac < nfin; fac++)
776 {
777 const int elem1 = face_voisins(fac, 0);
778 const int elem2 = face_voisins(fac, 1);
779 const double a = volumes(elem1)/(volumes(elem1) + volumes(elem2));
780 const double b = volumes(elem2)/(volumes(elem1) + volumes(elem2));
781
782 for (int i = 0; i < Objet_U::dimension; i++)
783 u_theta(fac, i) = a*tab_beta(elem1)*alpha_turb(elem1)*gradient_elem(elem1, i)
784 + b*tab_beta(elem2)*alpha_turb(elem2)*gradient_elem(elem2, i);
785 }
786 }
787 else
788 {
789 for (int fac = ndeb; fac < nfin; fac++)
790 {
791 const int elem1 = face_voisins(fac, 0);
792
793 for (int i = 0; i < Objet_U::dimension; i++)
794 u_theta(fac, i) = tab_beta(elem1)*alpha_turb(elem1)*gradient_elem(elem1, i);
795 }
796 }
797 }
798 // we treat the internal faces
799 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
800 for (int fac = 0; fac < nb_faces_tot; fac++)
801 {
802 const int elem1 = face_voisins(fac, 0);
803 const int elem2 = face_voisins(fac, 1);
804
805 if ((elem1 >= 0) && (elem2 >= 0))
806 {
807 const double a = volumes(elem1)/(volumes(elem1) + volumes(elem2));
808 const double b = volumes(elem2)/(volumes(elem1) + volumes(elem2));
809
810 for (int i = 0; i < Objet_U::dimension; i++)
811 u_theta(fac, i) = a*tab_beta(elem1)*alpha_turb(elem1)*gradient_elem(elem1, i)
812 + b*tab_beta(elem2)*alpha_turb(elem2)*gradient_elem(elem2, i);
813 }
814 }
815}
816
817void Calcul_Production_K_VEF::compute_utheta_nbConsti_le_1_nbCompo_gt_1(const Domaine_VEF& domaine_VEF,
818 const Domaine_Cl_VEF& zcl_VEF,
819 const IntTab& face_voisins,
820 const DoubleVect& volumes,
821 const DoubleTab& tab_beta,
822 const DoubleTab& alpha_turb,
823 const DoubleTrav& gradient_elem,
824 DoubleTrav& u_theta) const
825{
826 // we treat the boundaries
827 ToDo_Kokkos("critical");
828 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
829 {
830 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
831 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
832 const int ndeb = le_bord.num_premiere_face();
833 const int nfin = ndeb + le_bord.nb_faces();
834
835 if (sub_type(Periodique, la_cl.valeur()))
836 {
837 for (int fac = ndeb; fac < nfin; fac++)
838 {
839 const int elem1 = face_voisins(fac, 0);
840 const int elem2 = face_voisins(fac, 1);
841 const double a = volumes(elem1)/(volumes(elem1) + volumes(elem2));
842 const double b = volumes(elem2)/(volumes(elem1) + volumes(elem2));
843
844 for (int i = 0; i < Objet_U::dimension; i++)
845 u_theta(fac, i) = a*tab_beta(elem1, 0)*alpha_turb(elem1)*gradient_elem(elem1, i)
846 + b*tab_beta(elem2, 0)*alpha_turb(elem2)*gradient_elem(elem2, i);
847
848 }
849 }
850 else
851 {
852 for (int fac = ndeb; fac < nfin; fac++)
853 {
854 const int elem1 = face_voisins(fac, 0);
855
856 for (int i = 0; i < Objet_U::dimension; i++)
857 u_theta(fac, i) = tab_beta(elem1, 0)*alpha_turb(elem1)*gradient_elem(elem1, i);
858 }
859 }
860 }
861 // we treat the internal faces
862 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
863 ToDo_Kokkos("critical");
864 for (int fac = 0; fac < nb_faces_tot; fac++)
865 {
866 const int elem1 = face_voisins(fac, 0);
867 const int elem2 = face_voisins(fac, 1);
868
869 if ((elem1 >= 0) && (elem2 >= 0))
870 {
871 const double a = volumes(elem1)/(volumes(elem1) + volumes(elem2));
872 const double b = volumes(elem2)/(volumes(elem1) + volumes(elem2));
873
874 for (int i = 0; i < Objet_U::dimension; i++)
875 u_theta(fac, i) = a*tab_beta(elem1, 0)*alpha_turb(elem1)*gradient_elem(elem1, i)
876 + b*tab_beta(elem2, 0)*alpha_turb(elem2)*gradient_elem(elem2, i);
877 }
878 }
879}
880
881void Calcul_Production_K_VEF::compute_utheta_nbConsti_gt_1_nbCompo_eq_0(const Domaine_VEF& domaine_VEF,
882 const Domaine_Cl_VEF& zcl_VEF,
883 const IntTab& face_voisins,
884 const DoubleVect& volumes,
885 const DoubleTab& tab_beta,
886 const DoubleTab& alpha_turb,
887 const DoubleTrav& gradient_elem,
888 const int nb_consti,
889 DoubleTrav& u_theta) const
890{
891 // we treat the boundaries
892 ToDo_Kokkos("critical");
893 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
894 {
895 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
896 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
897 const int ndeb = le_bord.num_premiere_face();
898 const int nfin = ndeb + le_bord.nb_faces();
899
900 if (sub_type(Periodique, la_cl.valeur()))
901 {
902 for (int fac = ndeb; fac < nfin; fac++)
903 {
904 const int elem1 = face_voisins(fac, 0);
905 const int elem2 = face_voisins(fac, 1);
906 const double invVol = 1/(volumes(elem1) + volumes(elem2));
907 const double a = volumes(elem1)*invVol;
908 const double b = volumes(elem2)*invVol;
909
910 for (int i = 0; i < nb_consti; i++)
911 for (int k = 0; k < Objet_U::dimension; k++)
912 u_theta(fac, i, k) = a*tab_beta(0, 0)*alpha_turb(elem1)*gradient_elem(elem1, i, k)
913 + b*tab_beta(0, 0)*alpha_turb(elem2)*gradient_elem(elem2, i, k);
914 }
915 }
916 else
917 {
918 for (int fac = ndeb; fac < nfin; fac++)
919 {
920 const int elem1 = face_voisins(fac, 0);
921
922 for (int i = 0; i < nb_consti; i++)
923 for (int k = 0; k < Objet_U::dimension; k++)
924 u_theta(fac, i, k) = tab_beta(0, 0)*alpha_turb(elem1)*gradient_elem(elem1, i, k);
925 }
926 }
927 }
928
929 // we treat the internal faces
930 ToDo_Kokkos("critical");
931 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
932 for (int fac = 0; fac < nb_faces_tot; fac++)
933 {
934 const int elem1 = face_voisins(fac, 0);
935 const int elem2 = face_voisins(fac, 1);
936
937 if ((elem1 >= 0) && (elem2 >= 0))
938 {
939 const double invVol = 1/(volumes(elem1) + volumes(elem2));
940 const double a = volumes(elem1)*invVol;
941 const double b = volumes(elem2)*invVol;
942
943 for (int i = 0; i < nb_consti; i++)
944 for (int k = 0; k < Objet_U::dimension; k++)
945 u_theta(fac, i, k) = a*tab_beta(0, 0)*alpha_turb(elem1)*gradient_elem(elem1, i, k)
946 + b*tab_beta(0, 0)*alpha_turb(elem2)*gradient_elem(elem2, i, k);
947 }
948 }
949}
950
951void Calcul_Production_K_VEF::compute_utheta_nbConsti_gt_1_nbCompo_eq_1(const Domaine_VEF& domaine_VEF,
952 const Domaine_Cl_VEF& zcl_VEF,
953 const IntTab& face_voisins,
954 const DoubleVect& volumes,
955 const DoubleTab& tab_beta,
956 const DoubleTab& alpha_turb,
957 const DoubleTrav& gradient_elem,
958 const int nb_consti,
959 DoubleTrav& u_theta) const
960{
961 // we treat the boundaries
962 ToDo_Kokkos("critical");
963 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
964 {
965 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
966 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
967 const int ndeb = le_bord.num_premiere_face();
968 const int nfin = ndeb + le_bord.nb_faces();
969 if (sub_type(Periodique, la_cl.valeur()))
970 {
971 for (int fac = ndeb; fac < nfin; fac++)
972 {
973 const int elem1 = face_voisins(fac, 0);
974 const int elem2 = face_voisins(fac, 1);
975 const double invVol = 1/(volumes(elem1) + volumes(elem2));
976 const double a = volumes(elem1)*invVol;
977 const double b = volumes(elem2)*invVol;
978
979 for (int i = 0; i < nb_consti; i++)
980 for (int k = 0; k < Objet_U::dimension; k++)
981 u_theta(fac, i, k) = a*tab_beta(elem1)*alpha_turb(elem1)*gradient_elem(elem1, i, k)
982 + b*tab_beta(elem2)*alpha_turb(elem2)*gradient_elem(elem2, i, k);
983 }
984 }
985 else
986 {
987 for (int fac = ndeb; fac < nfin; fac++)
988 {
989 const int elem1 = face_voisins(fac, 0);
990
991 for (int i = 0; i < nb_consti; i++)
992 for (int k = 0; k < Objet_U::dimension; k++)
993 u_theta(fac, i, k) = tab_beta(elem1)*alpha_turb(elem1)*gradient_elem(elem1, i, k);
994 }
995 }
996 }
997 // we treat the internal faces
998 ToDo_Kokkos("critical");
999 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
1000 for (int fac = 0; fac < nb_faces_tot; fac++)
1001 {
1002 const int elem1 = face_voisins(fac, 0);
1003 const int elem2 = face_voisins(fac, 1);
1004
1005 if ((elem1 >= 0) && (elem2 >= 0))
1006 {
1007 const double invVol = 1/(volumes(elem1) + volumes(elem2));
1008 const double a = volumes(elem1)*invVol;
1009 const double b = volumes(elem2)*invVol;
1010 for (int i = 0; i < nb_consti; i++)
1011 for (int k = 0; k < Objet_U::dimension; k++)
1012 u_theta(fac, i, k) = a*tab_beta(elem1)*alpha_turb(elem1)*gradient_elem(elem1, i, k)
1013 + b*tab_beta(elem2)*alpha_turb(elem2)*gradient_elem(elem2, i, k);
1014 }
1015 }
1016
1017}
1018
1019void Calcul_Production_K_VEF::compute_utheta_nbConsti_gt_1_nbCompo_gt_1(const Domaine_VEF& domaine_VEF,
1020 const Domaine_Cl_VEF& zcl_VEF,
1021 const IntTab& face_voisins,
1022 const DoubleVect& volumes,
1023 const DoubleTab& tab_beta,
1024 const DoubleTab& alpha_turb,
1025 const DoubleTrav& gradient_elem,
1026 const int nb_consti,
1027 DoubleTrav& u_theta) const
1028{
1029 // we treat the boundaries
1030 ToDo_Kokkos("critical");
1031 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
1032 {
1033 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
1034 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1035 const int ndeb = le_bord.num_premiere_face();
1036 const int nfin = ndeb + le_bord.nb_faces();
1037 if (sub_type(Periodique, la_cl.valeur()))
1038 {
1039 for (int fac = ndeb; fac < nfin; fac++)
1040 {
1041 const int elem1 = face_voisins(fac, 0);
1042 const int elem2 = face_voisins(fac, 1);
1043 const double invVol = 1/(volumes(elem1) + volumes(elem2));
1044 const double a = volumes(elem1)*invVol;
1045 const double b = volumes(elem2)*invVol;
1046
1047 for (int i = 0; i < nb_consti; i++)
1048 for (int k = 0; k < Objet_U::dimension; k++)
1049 u_theta(fac, i, k) = a*tab_beta(elem1, 0)*alpha_turb(elem1)*gradient_elem(elem1, i, k)
1050 + b*tab_beta(elem2, 0)*alpha_turb(elem2)*gradient_elem(elem2, i, k);
1051 }
1052 }
1053 else
1054 {
1055 for (int fac = ndeb; fac < nfin; fac++)
1056 {
1057 const int elem1 = face_voisins(fac, 0);
1058
1059 for (int i = 0; i < nb_consti; i++)
1060 for (int k = 0; k < Objet_U::dimension; k++)
1061 u_theta(fac, i, k) = tab_beta(elem1, 0)*alpha_turb(elem1)*gradient_elem(elem1, i, k);
1062 }
1063 }
1064 }
1065 // we treat the internal faces
1066 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
1067 ToDo_Kokkos("critical");
1068 for (int fac = 0; fac < nb_faces_tot; fac++)
1069 {
1070 const int elem1 = face_voisins(fac, 0);
1071 const int elem2 = face_voisins(fac, 1);
1072 if ((elem1 >= 0) && (elem2 >= 0))
1073 {
1074 const double invVol = 1/(volumes(elem1) + volumes(elem2));
1075 const double a = volumes(elem1)*invVol;
1076 const double b = volumes(elem2)*invVol;
1077
1078 for (int i = 0; i < nb_consti; i++)
1079 for (int k = 0; k < Objet_U::dimension; k++)
1080 u_theta(fac, i, k) = a*tab_beta(elem1, 0)*alpha_turb(elem1)*gradient_elem(elem1, i, k)
1081 + b*tab_beta(elem2, 0)*alpha_turb(elem2)*gradient_elem(elem2, i, k);
1082 }
1083 }
1084}
1085
1086// Calcul d'un tenseur aux faces a partir d'un tenseur aux elements
1087DoubleTab& Calcul_Production_K_VEF::calcul_tenseur_face(DoubleTab& Tenseur_face,
1088 const DoubleTab& Tenseur_elem,
1089 const Domaine_VEF& domaine_VEF,
1090 const Domaine_Cl_VEF& domaine_Cl_VEF) const
1091{
1092
1093 const int dimension = Objet_U::dimension;
1094 const IntTab& face_voisins = domaine_VEF.face_voisins();
1095
1096 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
1097 const int nb_cl = les_cl.size();
1098 const DoubleVect& volumes = domaine_VEF.volumes();
1099 ToDo_Kokkos("critical");
1100 for (int n_bord = 0; n_bord < nb_cl; n_bord++)
1101 {
1102 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1103 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
1104 const int ndeb = le_bord.num_premiere_face();
1105 const int nfin = ndeb + le_bord.nb_faces();
1106
1107 if (sub_type(Periodique, la_cl.valeur()))
1108 {
1109 for (int fac = ndeb; fac < nfin; fac++)
1110 {
1111 const int poly1 = face_voisins(fac, 0);
1112 const int poly2 = face_voisins(fac, 1);
1113 const double invVol = 1/(volumes(poly1) + volumes(poly2));
1114 const double a = volumes(poly1)*invVol;
1115 const double b = volumes(poly2)*invVol;
1116
1117 for (int i = 0; i < dimension; i++)
1118 for (int j = 0; j < dimension; j++)
1119 Tenseur_face(fac, i, j) = a*Tenseur_elem(poly1, i, j)
1120 + b*Tenseur_elem(poly2, i, j);
1121 }
1122 }
1123 else
1124 {
1125 for (int fac = ndeb; fac < nfin; fac++)
1126 {
1127 const int poly1 = face_voisins(fac, 0);
1128
1129 for (int i = 0; i < dimension; i++)
1130 for (int j = 0; j < dimension; j++)
1131 Tenseur_face(fac, i, j) = Tenseur_elem(poly1, i, j);
1132 }
1133 }
1134 }
1135
1136 const int n0 = domaine_VEF.premiere_face_int();
1137 const int nb_faces = domaine_VEF.nb_faces();
1138 ToDo_Kokkos("critical");
1139 for (int fac = n0; fac < nb_faces; fac++)
1140 {
1141 const int poly1 = face_voisins(fac, 0);
1142 const int poly2 = face_voisins(fac, 1);
1143 const double invVol = 1/(volumes(poly1) + volumes(poly2));
1144 const double a = volumes(poly1)*invVol;
1145 const double b = volumes(poly2)*invVol;
1146
1147 for (int i = 0; i < dimension; i++)
1148 for (int j = 0; j < dimension; j++)
1149 Tenseur_face(fac, i, j) = a*Tenseur_elem(poly1, i, j) + b*Tenseur_elem(poly2, i, j);
1150 }
1151
1152 return Tenseur_face;
1153}
Calcul_Production_K_VEF::compute_utheta_nbConsti_le_1_nbCompo_eq_0
void compute_utheta_nbConsti_le_1_nbCompo_eq_0(const Domaine_VEF &domaine_VEF, const Domaine_Cl_VEF &zcl_VEF, const IntTab &face_voisins, const DoubleVect &volumes, const DoubleTab &tab_beta, const DoubleTab &alpha_turb, const DoubleTrav &gradient_elem, DoubleTrav &u_theta) const
Definition Calcul_Production_K_VEF.cpp:692
Calcul_Production_K_VEF::compute_utheta_nbConsti_gt_1_nbCompo_eq_1
void compute_utheta_nbConsti_gt_1_nbCompo_eq_1(const Domaine_VEF &domaine_VEF, const Domaine_Cl_VEF &zcl_VEF, const IntTab &face_voisins, const DoubleVect &volumes, const DoubleTab &tab_beta, const DoubleTab &alpha_turb, const DoubleTrav &gradient_elem, const int nb_consti, DoubleTrav &u_theta) const
Definition Calcul_Production_K_VEF.cpp:951
Calcul_Production_K_VEF::compute_utheta_nbConsti_le_1_nbCompo_eq_1
void compute_utheta_nbConsti_le_1_nbCompo_eq_1(const Domaine_VEF &domaine_VEF, const Domaine_Cl_VEF &zcl_VEF, const IntTab &face_voisins, const DoubleVect &volumes, const DoubleTab &tab_beta, const DoubleTab &alpha_turb, const DoubleTrav &gradient_elem, DoubleTrav &u_theta) const
Definition Calcul_Production_K_VEF.cpp:755
Calcul_Production_K_VEF::loop_for_non_periodic_boundaries
void loop_for_non_periodic_boundaries(DoubleTab &prodK, const DoubleTab &gradient_elem, const DoubleTab &visco_turb, const DoubleVect &volumes, const IntTab &face_voisins, const int nfaceinit, const int nfaceend, const int interpol_visco, const double limiteur, const DoubleTab &K_Omega, const bool &activate_production_limiter=false, const double &cst_production_limiter=0.) const
Compute production term on non periodic boundary faces.
Definition Calcul_Production_K_VEF.cpp:114
Calcul_Production_K_VEF::compute_production_term_EASM
void compute_production_term_EASM(const int face, const double visco_face, const DoubleTab &Re_face, const DoubleTab &gradient_face, DoubleTab &P) const
Definition Calcul_Production_K_VEF.cpp:485
Calcul_Production_K_VEF::compute_utheta_nbConsti_gt_1_nbCompo_gt_1
void compute_utheta_nbConsti_gt_1_nbCompo_gt_1(const Domaine_VEF &domaine_VEF, const Domaine_Cl_VEF &zcl_VEF, const IntTab &face_voisins, const DoubleVect &volumes, const DoubleTab &tab_beta, const DoubleTab &alpha_turb, const DoubleTrav &gradient_elem, const int nb_consti, DoubleTrav &u_theta) const
Definition Calcul_Production_K_VEF.cpp:1019
Calcul_Production_K_VEF::calculer_terme_destruction_K_gen
DoubleTab & calculer_terme_destruction_K_gen(const Domaine_VEF &, const Domaine_Cl_VEF &, DoubleTab &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &, const DoubleVect &, int) const
Definition Calcul_Production_K_VEF.cpp:584
Calcul_Production_K_VEF::calculer_terme_production_K
DoubleTab & calculer_terme_production_K(const Domaine_VEF &, const Domaine_Cl_VEF &, DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const int &interpol_visco, const double &limiteur, const bool &deactivate_production_limiter=false, const double &cst_production_limiter=0.) const
Compute the production term for the turbulent kinetic energy.
Definition Calcul_Production_K_VEF.cpp:48
Calcul_Production_K_VEF::compute_utheta_nbConsti_le_1_nbCompo_gt_1
void compute_utheta_nbConsti_le_1_nbCompo_gt_1(const Domaine_VEF &domaine_VEF, const Domaine_Cl_VEF &zcl_VEF, const IntTab &face_voisins, const DoubleVect &volumes, const DoubleTab &tab_beta, const DoubleTab &alpha_turb, const DoubleTrav &gradient_elem, DoubleTrav &u_theta) const
Definition Calcul_Production_K_VEF.cpp:817
Calcul_Production_K_VEF::calcul_tenseur_face
DoubleTab & calcul_tenseur_face(DoubleTab &, const DoubleTab &, const Domaine_VEF &, const Domaine_Cl_VEF &) const
Definition Calcul_Production_K_VEF.cpp:1087
Calcul_Production_K_VEF::calculer_terme_production_K_EASM
DoubleTab & calculer_terme_production_K_EASM(const Domaine_VEF &, const Domaine_Cl_VEF &, DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const int &interpol_visco, const double &limiteur) const
Definition Calcul_Production_K_VEF.cpp:406
Calcul_Production_K_VEF::get_turbulent_viscosity
KOKKOS_INLINE_FUNCTION double get_turbulent_viscosity(CDoubleArrView visco_turb, CDoubleArrView volumes, const int type_interpo, const int poly1, const int poly2, const double limiteur) const
Get the turbulent viscosity depending on the interpolation choice.
Definition Calcul_Production_K_VEF.cpp:513
Calcul_Production_K_VEF::calculer_terme_production_K_BiK
DoubleTab & calculer_terme_production_K_BiK(const Domaine_VEF &, const Domaine_Cl_VEF &, DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const int &interpol_visco, const double &limiteur) const
Definition Calcul_Production_K_VEF.cpp:248
Calcul_Production_K_VEF::loop_for_internal_or_periodic_faces
void loop_for_internal_or_periodic_faces(DoubleTab &prodK, const DoubleTab &gradient_elem, const DoubleTab &visco_turb, const DoubleVect &volumes, const IntTab &face_voisins, const int nfaceinit, const int nfaceend, const int interpol_visco, const double limiteur, const DoubleTab &K_Omega, const bool &activate_production_limiter=false, const double &cst_production_limiter=0.) const
Compute production term on internal and periodic boundary faces.
Definition Calcul_Production_K_VEF.cpp:185
Calcul_Production_K_VEF::compute_utheta_nbConsti_gt_1_nbCompo_eq_0
void compute_utheta_nbConsti_gt_1_nbCompo_eq_0(const Domaine_VEF &domaine_VEF, const Domaine_Cl_VEF &zcl_VEF, const IntTab &face_voisins, const DoubleVect &volumes, const DoubleTab &tab_beta, const DoubleTab &alpha_turb, const DoubleTrav &gradient_elem, const int nb_consti, DoubleTrav &u_theta) const
Definition Calcul_Production_K_VEF.cpp:881
Champ_Don_base
classe Champ_Don_base classe de base des Champs donnes (non calcules)
Definition Champ_Don_base.h:32
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_P1NC::calcul_gradient
static DoubleTab & calcul_gradient(const DoubleTab &, DoubleTab &, const Domaine_Cl_VEF &)
Definition Champ_P1NC.cpp:915
Champ_Uniforme
classe Champ_Uniforme Represente un champ constant dans l'espace et dans le temps.
Definition Champ_Uniforme.h:26
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Conds_lim
classe Conds_lim Cette classe represente un vecteur de conditions aux limites.
Definition Conds_lim.h:32
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::nb_faces_tot
int nb_faces_tot() const
renvoie le nombre total de faces.
Definition Domaine_VF.h:481
Domaine_VF::volumes
double volumes(int i) const
Definition Domaine_VF.h:113
Domaine_VF::premiere_face_int
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_dis_base::nb_elem_tot
int nb_elem_tot() const
Definition Domaine_dis_base.h:50
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Field_base::nb_comp
virtual int nb_comp() const
Definition Field_base.h:56
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
Process::Kokkos_exit
static KOKKOS_INLINE_FUNCTION void Kokkos_exit(const char *)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.h:173
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
TRUSTTab::view_ro
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, ConstView< _TYPE_, _SHAPE_ > > view_ro() const
Definition TRUSTTab.h:261