TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Source_PDF_VEF.cpp
1/****************************************************************************
2* Copyright (c) 2026, 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 <Source_PDF_VEF.h>
17#include <Domaine.h>
18#include <Domaine_VEF.h>
19#include <Domaine_Cl_VEF.h>
20#include <Probleme_base.h>
21#include <Equation_base.h>
22#include <Schema_Temps_base.h>
23#include <Discretisation_base.h>
24#include <Matrice_Morse.h>
25#include <Matrice_Bloc.h>
26#include <TRUSTTrav.h>
27#include <Interpolation_IBM_elem_fluid.h>
28#include <Interpolation_IBM_mean_gradient.h>
29#include <Interpolation_IBM_hybrid.h>
30#include <Interpolation_IBM_power_law_tbl.h>
31#include <Interpolation_IBM_power_law_tbl_u_star.h>
32#include <Interpolation_IBM_thermal_wall_law.h>
33#include <Dirichlet.h>
34#include <SFichier.h>
35#include <Champ_P1NC.h>
36#include <Op_Diff_VEF_base.h>
37
38Implemente_instanciable(Source_PDF_VEF,"Source_PDF_VEF_P1NC",Source_PDF_base);
39
40/*##################################################################################################
41####################################################################################################
42################################# READON AND PRINTON ###############################################
43####################################################################################################
44##################################################################################################*/
45
46
47
48Entree& Source_PDF_VEF::readOn(Entree& s)
49{
50 Source_PDF_base::readOn(s);
51 return s;
52}
53
54Sortie& Source_PDF_VEF::printOn(Sortie& s ) const
55{
56 Source_PDF_base::printOn(s);
57 return s;
58}
59
60/*##################################################################################################
61####################################################################################################
62################################# PROBLEM ASSOCIATION AND ZONES ####################################
63####################################################################################################
64##################################################################################################*/
65
66void Source_PDF_VEF::associer_pb(const Probleme_base& pb)
67{
68 Source_PDF_base::associer_pb(pb);
69
70 int nb_som=le_dom_VEF->nb_som();
71 tab_u_star_ibm_.resize(nb_som);
72 tab_y_plus_ibm_.resize(nb_som);
73}
74
79
80void Source_PDF_VEF::associer_domaines(const Domaine_dis_base& domaine_dis,
81 const Domaine_Cl_dis_base& domaine_Cl_dis)
82{
83 le_dom_VEF = ref_cast(Domaine_VEF, domaine_dis);
84 le_dom_Cl_VEF = ref_cast(Domaine_Cl_VEF, domaine_Cl_dis);
85}
86
87void Source_PDF_VEF::multiply_coeff_volume(DoubleTab& coeff) const
88{
89 const DoubleVect& vol_diam=ref_cast(Domaine_VEF, le_dom_VEF.valeur()).volumes_entrelaces();
90 int n = vol_diam.size_totale();
91 int nc = coeff.dimension_tot(0) ;
92 if (n != nc)
93 {
94 Cerr<<" dimensions differentes n nc = "<<n<<" "<<nc<<finl;
95 exit();
96 }
97 int ndim = coeff.dimension(1) ;
98 for (int i=0; i<n; i++)
99 {
100 for (int comp=0; comp<ndim; comp++) coeff(i,comp) = coeff(i,comp)*vol_diam(i);
101 }
102}
103
104/*##################################################################################################
105####################################################################################################
106################################# AJOUTER_ #########################################################
107####################################################################################################
108##################################################################################################*/
109
110/*
111This function redirects toward the ajouter_ which correspond to the model chosen.
112*/
113
114DoubleTab& Source_PDF_VEF::ajouter_(const DoubleTab& variable, DoubleTab& resu, const int i_traitement_special) const
115{
116 // calcul de : coefficient rho/dt * coeff_relax * (volume_thilde / 8) * variable
117 /* i_traitement_special = 0 => traitement classique; coefficient (1/eta) */
118 /* i_traitement_special = 1 => coefficient (1/eta -> 1) */
119 /* i_traitement_special = 2 => coefficient (1/eta -> 1 + 1/eta) */
120
121 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
122 const IntTab& elems= domaine_VEF.elem_faces() ;
123 int nb_face_elem=domaine_VEF.domaine().nb_faces_elem();
124 int nb_elems=domaine_VEF.domaine().nb_elem_tot();
125 const DoubleVect& volume=domaine_VEF.volumes_entrelaces();
126 int dim_esp = Objet_U::dimension ;
127 int dim_var = variable.dimension(1);
128 // Cerr << "erreur ici" << finl;
129 if (dim_var != 1 && dim_var != dim_esp)
130 {
131 Cerr<<"ajouter_: dimension variable differente 1 ou "<<dim_esp<<"! "<<finl;
132 exit();
133 }
134 ArrOfDouble tuvw(dim_var);
135 const DoubleTab& rotation=champ_rotation_->valeurs();
136 const DoubleTab& aire=champ_aire_->valeurs();
137 //OWN_PTR(Champ_Don_base) rho_test;
138 //champ_rho_->affecter(equation().probleme().get_champ("masse_volumique"));
139 const DoubleTab& rho_m=champ_rho_->valeurs();
140 DoubleTab pond;
141
142 // return = rho/dt * coeff_relax * volume
143 int coef1 = nb_face_elem;
144 pond = compute_pond(rho_m, aire, volume, coef1, nb_elems);
145
146 for (int num_elem=0; num_elem<nb_elems; num_elem++)
147 {
148 if (aire(num_elem)>0.)
149 {
150 if (i_traitement_special == 0)
151 {
152 for (int k=0; k<dim_var; k++) tuvw[k] = 1.0 / modele_lu_.eta_;
153 }
154 else if (i_traitement_special == 1) //terme temps en rho v
155 {
156 for (int k=0; k<dim_var; k++) tuvw[k] = 1.0;
157 }
158 else if (i_traitement_special == 101) //terme temps en v
159 {
160 for (int k=0; k<dim_var; k++) tuvw[k] = 1.0 / rho_m(num_elem);
161 }
162 else if (i_traitement_special == 2)
163 {
164 for (int k=0; k<dim_var; k++) tuvw[k] = 1.0 + 1.0 / modele_lu_.eta_;
165 }
166 else if (i_traitement_special == 102)
167 {
168 for (int k=0; k<dim_var; k++) tuvw[k] = 1.0 / rho_m(num_elem) + 1.0 / modele_lu_.eta_;
169 }
170 else
171 {
172 Cerr << "Source_PDF_VEF::ajouter_ : i_traitement_special should be 0, 1, 2, 101 or 102" << finl;
173 exit();
174 }
175 for (int comp=0; comp<dim_var; comp++)
176 {
177 double tijvj=0;
178 if (dim_var == 1)
179 {
180 tijvj = tuvw[0];
181 }
182 else
183 // dim_var = dim_esp
184 {
185 for (int c=0; c<dim_var; c++)
186 {
187 double coeff_im=0;
188 for(int k=0; k<dim_esp; k++)
189 {
190 coeff_im+=rotation(num_elem,3*k+comp)*tuvw[k]*rotation(num_elem,3*k+c);
191 }
192 tijvj+=coeff_im ;
193 }
194 }
195 // Cerr<< "Source_PDF_VEF::ajouter_ tijvj=" <<tijvj<<" pond("<<num_elem<<") = "<<pond(num_elem)<< finl;
196 // cas generique : 1/eta rho/dt * coeff_relax * volume * variable
197 for (int i=0; i<nb_face_elem; i++)
198 {
199 resu(elems(num_elem,i),comp)-=pond(num_elem)*coef1*tijvj*variable(elems(num_elem,i),comp);
200 // Cerr<< "Source_PDF_VEF::ajouter_ noeud: "<<i <<" "<<pond(num_elem)*coef1*tijvj*variable(elems(num_elem,i),comp)<<" variable = "<<variable(elems(num_elem,i),comp)<< finl;
201 }
202 }
203 }
204 }
205 return resu;
206}
207
208DoubleTab& Source_PDF_VEF::ajouter_(const DoubleTab& variable, DoubleTab& resu) const
209{
210 const int i_traitement_special = 0 ;
211 ajouter_(variable, resu, i_traitement_special) ;
212 return resu;
213}
214
215/*##################################################################################################
216####################################################################################################
217################################# CONTRIBUER_A_AVEC ################################################
218####################################################################################################
219##################################################################################################*/
220
221/*
222This function redirects toward the contribuer_avec_ which correspond to the model chosen.
223*/
224
225void Source_PDF_VEF::contribuer_a_avec(const DoubleTab& inco, Matrice_Morse& matrice) const
226{
227 // calcul : 1/eta rho/dt * coeff_relax * volume
228 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
229 const IntTab& elems= domaine_VEF.elem_faces();
230 int nb_face_elem=domaine_VEF.domaine().nb_faces_elem();
231 int nb_elems=domaine_VEF.domaine().nb_elem_tot();
232 const DoubleVect& volume=domaine_VEF.volumes_entrelaces();
233 const DoubleTab& variable=equation().inconnue().valeurs();
234 int dim_var = variable.dimension(1);
235 int dim_esp = Objet_U::dimension ;
236 if (dim_var != 1 && dim_var != dim_esp)
237 {
238 Cerr<<"contribuer_a_avec: dimension variable differente 1 ou "<<dim_esp<<"! "<<finl;
239 exit();
240 }
241 ArrOfDouble tuvw(dim_var);
242 const DoubleTab& rotation=champ_rotation_->valeurs();
243 const DoubleTab& aire=champ_aire_->valeurs();
244 //champ_rho_->affecter(equation().probleme().get_champ("masse_volumique"));
245 const DoubleTab& rho_m=champ_rho_->valeurs();
246 DoubleTab pond ;
247
248 // return = rho/dt * coeff_relax * volume pour elements interceptes
249 int coef1 = nb_face_elem;
250 pond = compute_pond(rho_m, aire, volume, coef1, nb_elems);
251
252 for (int num_elem=0; num_elem<nb_elems; num_elem++) //elements loop
253 {
254 if (aire(num_elem)>0.)
255 {
256 for (int k=0; k<dim_var; k++) tuvw[k] = 1.0 / modele_lu_.eta_;
257 // Cerr<< "Source_PDF_VEF::contribuer_a_avec: tuvw=" <<tuvw<<" pond("<<num_elem<<") = "<<pond(num_elem)<< finl;
258
259 for (int comp=0; comp<dim_var; comp++)
260 {
261 for (int comp2=0; comp2<dim_var; comp2++)
262 {
263 for (int i=0; i<nb_face_elem; i++)
264 {
265 for (int j=0; j<nb_face_elem; j++)
266 {
267 int s1=elems(num_elem,i);
268
269 double coeff_im=0;
270 if (dim_var == 1)
271 {
272 coeff_im = tuvw[0];
273 }
274 else
275 // dim_var = dim_esp
276 {
277 for (int k=0; k<dim_esp; k++)
278 {
279 coeff_im+=rotation(num_elem,3*k+comp)*tuvw[k]*rotation(num_elem,3*k+comp2);
280 }
281 }
282 int c1=s1*dim_var+comp;
283 // Cerr << c1 << finl;
284
285 // 1.0/eta * rho/dt * coeff_relax * volume
286 matrice.coef(c1,c1)+=pond(num_elem)*coeff_im;
287 // Cerr<< "Source_PDF_VEF::contribuer_a_avec: contribution matrice c1c1 : "<<c1<<" = "<<pond(num_elem)*coeff_im<< finl;
288 }
289 }
290 }
291 }
292 }
293 }
294 // verif_ajouter_contrib(variable, matrice) ;
295}
296
297void Source_PDF_VEF::verif_ajouter_contrib(const DoubleTab& variable, Matrice_Morse& matrice) const
298{
299 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
300 const IntTab& elems= domaine_VEF.elem_faces();
301 int nb_face_elem=domaine_VEF.domaine().nb_faces_elem();
302 int nb_elems=domaine_VEF.domaine().nb_elem_tot();
303 const DoubleTab& aire=champ_aire_->valeurs();
304
305 DoubleTrav force(variable);
306 int dim_var = variable.dimension(1);
307 ajouter_(variable,force);
308 DoubleTrav force2(variable);
309 matrice.ajouter_multvect_(variable, force2) ;
310 DoubleTab diforce(force);
311 diforce += force2 ;
312 double difmax = 0.;
313 for (int num_elem=0; num_elem<nb_elems; num_elem++)
314 {
315 if (aire(num_elem)>0.)
316 {
317 for (int i=0; i<nb_face_elem; i++)
318 {
319 int s1=elems(num_elem,i);
320 double difmod2 = 0.;
321 for (int k=0; k<dim_var; k++) difmod2 += diforce(s1,k)*diforce(s1,k);
322 if (difmod2 > 1.0e-5)
323 {
324 if (difmod2 > difmax) difmax = difmod2 ;
325 // for (int k=0; k<dim_var; k++) Cerr << " ajouter multvect diforce x = "<< force(s1,k) << " "<< force2(s1,k) << " " << diforce(s1,k)<<finl;
326 }
327 }
328 }
329 }
330 if (difmax > 0.)
331 {
332 Cerr<< "Source_PDF_VEF: Max norme caree diff. force absolue = "<<difmax<<finl;
333 }
334}
335
336/*##################################################################################################
337####################################################################################################
338################################# calculer_variable_imposee ########################################
339####################################################################################################
340##################################################################################################*/
341
342void Source_PDF_VEF::calculer_variable_imposee_elem_fluid()
343{
344 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
345 int nb_faces=domaine_VEF.nb_faces();
346 int dim_esp = Objet_U::dimension;
347 const Domaine& dom = domaine_VEF.domaine();
348 const IntTab& faces_sommets=domaine_VEF.face_sommets();
349 int nb_som_face=faces_sommets.dimension(1);
350 int nb_som = domaine_VEF.nb_som();
351 DoubleTab& variable_imposee_mod = modele_lu_.variable_imposee_->valeurs();
352 DoubleTab& variable_imposee_calculee = variable_imposee_;
353 Interpolation_IBM_elem_fluid& interp = ref_cast(Interpolation_IBM_elem_fluid,interpolation_lue_.valeur());
354 DoubleTab& fluid_points = interp.fluid_points_->valeurs();
355 DoubleTab& solid_points = interp.solid_points_->valeurs();
356 DoubleTab& fluid_elems = interp.fluid_elems_->valeurs();
357 Champ_P1NC& champ_variable_inconnue = ref_cast(Champ_P1NC,equation().inconnue());
358 DoubleTab& val_variable_inconnue = champ_variable_inconnue.valeurs();
359 int nb_comp = val_variable_inconnue.dimension(1);
360 DoubleTab variable_imposee_sommet(nb_som,nb_comp);
361 IntTab sommet_interp(1, nb_som);
362 DoubleTab xf(1, dim_esp);
363 DoubleTab vf(1, nb_comp);
364 ArrOfInt cells(1);
365 double eps = 1e-12;
366 variable_imposee_calculee = 0.0;
367 variable_imposee_sommet = 0.0;
368 sommet_interp = 0;
369 for (int i = 0; i < nb_som; i++)
370 {
371 if ((fluid_elems(i) >= 0.0))
372 {
373 sommet_interp(0, i) = 1;
374 double d1 = 0.0;
375 double d2 = 0.0;
376 for(int j = 0; j < dim_esp; j++)
377 {
378 double xj = dom.coord(i,j);
379 double xjf = fluid_points(i,j);
380 xf(0, j) = xjf;
381 double xjs = solid_points(i,j);
382 d1 += (xj-xjs)*(xj-xjs);
383 d2 += (xjf-xj)*(xjf-xj);
384 }
385 d1 = sqrt(d1);
386 d2 = sqrt(d2);
387 double inv_d ;
388 if ( (d1 + d2) > eps ) inv_d = 1.0 / (d1 + d2);
389 else inv_d = 0. ;
390 cells[0] = int(fluid_elems(i));
391 champ_variable_inconnue.valeur_a_elem(xf, vf, cells[0]);
392
393 for (int j = 0; j < nb_comp; j++)
394 {
395 double vjs = variable_imposee_mod(i,j);
396 variable_imposee_sommet(i,j) = vjs + (vf(0, j)-vjs)*inv_d*d1;
397 }
398 }
399 else
400 {
401 for(int j = 0; j < nb_comp; j++)
402 {
403 variable_imposee_sommet(i,j) = variable_imposee_mod(i,j);
404 }
405 }
406 }
407 for (int f = 0; f < nb_faces; f++)
408 {
409 int compteur = 0;
410 int som_f = 0;
411 while ((som_f < nb_som_face) and (sommet_interp(0, faces_sommets(f,som_f))==1))
412 {
413 compteur++;
414 som_f++;
415 }
416 if (compteur == nb_som_face)
417 {
418 for (int som_=0; som_ < nb_som_face; som_++)
419 {
420 int i = faces_sommets(f,som_);
421 for(int j = 0; j < nb_comp; j++)
422 {
423 variable_imposee_calculee(f,j) += variable_imposee_sommet(i,j)/nb_som_face;
424 }
425 }
426 }
427 }
428}
429
430void Source_PDF_VEF::calculer_variable_imposee_mean_grad()
431{
432 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
433 int nb_faces=domaine_VEF.nb_faces();
434 int dim_esp = Objet_U::dimension;
435 const Domaine& dom = domaine_VEF.domaine();
436 const IntTab& faces_sommets=domaine_VEF.face_sommets();
437 int nb_som_face=faces_sommets.dimension(1);
438 int nb_som = domaine_VEF.nb_som();
439 DoubleTab& variable_imposee_mod = modele_lu_.variable_imposee_->valeurs();
440 DoubleTab& variable_imposee_calculee = variable_imposee_;
441 Interpolation_IBM_mean_gradient& interp = ref_cast(Interpolation_IBM_mean_gradient,interpolation_lue_.valeur());
442 DoubleTab& solid_points = interp.solid_points_->valeurs();
443 DoubleTab& is_dirichlet = interp.is_dirichlet_->valeurs();
444 double eps = 1e-12;
445 DoubleTab& variable_inconnue = equation().inconnue().valeurs();
446 int nb_comp = variable_inconnue.dimension(1);
447 variable_imposee_mod.echange_espace_virtuel();
448 DoubleTab variable_imposee_sommet(nb_som,nb_comp);
449 IntTab sommet_interp(1, nb_som);
450 variable_imposee_calculee = 0.0;
451 variable_imposee_sommet = 0.0;
452 sommet_interp = 0;
453 for (int i = 0; i < nb_som; i++)
454 {
455 if (is_dirichlet(i) > 0.0)
456 {
457 sommet_interp(0, i) = 1;
458 double d2 = 0.0;
459 for(int j = 0; j < dim_esp; j++)
460 {
461 double x = dom.coord(i,j);
462 double xp = solid_points(i,j);
463 d2 += (x - xp)*(x - xp);
464 }
465 for(int j = 0; j <nb_comp ; j++) variable_imposee_sommet(i,j) = variable_imposee_mod(i,j);
466 d2 = sqrt(d2);
467 if (d2 > eps)
468 {
469 IntList& voisins = interp.getSommetsVoisinsOf(i);
470 if (!(voisins.est_vide()))
471 {
472 ArrOfDouble mean_grad(nb_comp);
473 mean_grad = 0.0;
474 int taille = voisins.size();
475 int nb_contrib = 0;
476 for (int k = 0; k < taille; k++)
477 {
478 int num_som = voisins[k];
479 double d1 = 0.0;
480 for (int j = 0; j < dim_esp; j++)
481 {
482 double xf = dom.coord(num_som,j);
483 double xpf = solid_points(num_som,j);
484 d1 += (xf - xpf)*(xf - xpf);
485 }
486 d1 = sqrt(d1);
487 if (d1 > eps)
488 {
489 for (int j = 0; j < nb_comp; j++)
490 {
491 double vpf = variable_imposee_mod(num_som,j);
492 double vf = variable_inconnue(num_som,j);
493 mean_grad[j] += (vf - vpf)/d1;
494 }
495 nb_contrib += 1;
496 }
497 }
498 for (int j = 0; j < nb_comp; j++)
499 {
500 mean_grad[j] /= nb_contrib;
501 variable_imposee_sommet(i,j) += mean_grad[j]*d2;
502 }
503 }
504 }
505 }
506 else
507 {
508 for(int j = 0; j < nb_comp; j++)
509 {
510 variable_imposee_sommet(i,j) = variable_imposee_mod(i,j);
511 }
512 }
513 }
514 for (int f = 0; f < nb_faces; f++)
515 {
516 int compteur = 0;
517 int som_f = 0;
518 while ((som_f < nb_som_face) and (sommet_interp(0, faces_sommets(f,som_f))==1))
519 {
520 compteur++;
521 som_f++;
522 }
523 if (compteur == nb_som_face)
524 {
525 for (int som_=0; som_ < nb_som_face; som_++)
526 {
527 int i = faces_sommets(f,som_);
528 for(int j = 0; j < nb_comp; j++)
529 {
530 variable_imposee_calculee(f,j) += variable_imposee_sommet(i,j)/nb_som_face;
531 }
532 }
533 }
534 }
535}
536
537void Source_PDF_VEF::calculer_variable_imposee_hybrid()
538{
539 Cerr << "on passe ici" << finl;
540 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
541 int nb_faces=domaine_VEF.nb_faces();
542 int dim_esp = Objet_U::dimension;
543 const Domaine& dom = domaine_VEF.domaine();
544 const IntTab& faces_sommets=domaine_VEF.face_sommets();
545 int nb_som_face=faces_sommets.dimension(1);
546 int nb_som = domaine_VEF.nb_som();
547 DoubleTab& variable_imposee_mod = modele_lu_.variable_imposee_->valeurs();
548 DoubleTab& variable_imposee_calculee = variable_imposee_;
549 Interpolation_IBM_hybrid& interp = ref_cast(Interpolation_IBM_hybrid,interpolation_lue_.valeur());
550 DoubleTab& solid_points = interp.solid_points_->valeurs();
551 DoubleTab& fluid_elems = interp.fluid_elems_->valeurs();
552 DoubleTab& fluid_points = interp.fluid_points_->valeurs();
553 double eps = 1e-12;
554 Champ_P1NC& champ_variable_inconnue = ref_cast(Champ_P1NC,equation().inconnue());
555 DoubleTab& variable_inconnue = champ_variable_inconnue.valeurs();
556 int nb_comp = variable_inconnue.dimension(1);
557 DoubleTab xf(1, dim_esp);
558 DoubleTab vf(1, nb_comp);
559 ArrOfInt cells(1);
560 DoubleTab variable_imposee_sommet(nb_som,nb_comp);
561 IntTab sommet_interp(1, nb_som);
562 variable_imposee_calculee = 0.0;
563 variable_imposee_sommet = 0.0;
564 sommet_interp = 0;
565 Cerr << "on passe ici" << finl;
566 for (int i = 0; i < nb_som; i++)
567 {
568 Cerr << "on passe la" << finl;
569 if (fluid_elems(i) >= 0.0)
570 {
571 sommet_interp(0, i) = 1;
572 double d1 = 0.0;
573 double d2 = 0.0;
574 for(int j = 0; j < dim_esp; j++)
575 {
576 double xj = dom.coord(i,j);
577 double xjf = fluid_points(i,j);
578 xf(0, j) = xjf;
579 double xjs = solid_points(i,j);
580 d1 += (xj-xjs)*(xj-xjs);
581 d2 += (xjf-xj)*(xjf-xj);
582 }
583 d1 = sqrt(d1);
584 d2 = sqrt(d2);
585 double inv_d = 1.0 / (d1 + d2);
586 cells[0] = int(fluid_elems(i));
587 champ_variable_inconnue.valeur_a_elem(xf, vf, cells[0]);
588 for (int j = 0; j < nb_comp; j++)
589 {
590 double vjs = variable_imposee_mod(i,j);
591 variable_imposee_sommet(i,j) = vjs + (vf(0, j)-vjs)*inv_d*d1;
592 }
593 }
594 else if (fluid_elems(i) > -1.5 && fluid_elems(i) < -0.5)
595 {
596 sommet_interp(0, i) = 1;
597 double d2 = 0.0;
598 for(int j = 0; j < dim_esp; j++)
599 {
600 double x = dom.coord(i,j);
601 double xp = solid_points(i,j);
602 d2 += (x - xp)*(x - xp);
603 }
604 d2 = sqrt(d2);
605 for(int j = 0; j <nb_comp ; j++) variable_imposee_sommet(i,j) = variable_imposee_mod(i,j);
606 if (d2 > eps)
607 {
608 IntList& voisins = interp.getSommetsVoisinsOf(i);
609 if (!(voisins.est_vide()))
610 {
611 ArrOfDouble mean_grad(nb_comp);
612 mean_grad = 0.0;
613 int taille = voisins.size();
614 int nb_contrib = 0;
615 for (int k = 0; k < taille; k++)
616 {
617 int num_som = voisins[k];
618 double d1 = 0.0;
619 for (int j = 0; j < dim_esp; j++)
620 {
621 double xjf = dom.coord(num_som,j);
622 double xpf = solid_points(num_som,j);
623 d1 += (xjf - xpf)*(xjf - xpf);
624 }
625 d1 = sqrt(d1);
626 if (d1 > eps)
627 {
628 for (int j = 0; j < nb_comp; j++)
629 {
630 double vpf = variable_imposee_mod(num_som,j);
631 double vjf = variable_inconnue(num_som,j);
632 mean_grad[j] += (vjf - vpf)/d1;
633 }
634 nb_contrib += 1;
635 }
636 }
637 for (int j = 0; j < nb_comp; j++)
638 {
639 mean_grad[j] /= nb_contrib;
640 variable_imposee_sommet(i,j) += mean_grad[j]*d2;
641 }
642 }
643 }
644 }
645 else
646 {
647 for(int j = 0; j < nb_comp; j++)
648 {
649 variable_imposee_sommet(i,j) = variable_imposee_mod(i,j);
650 }
651 }
652 }
653 for (int f = 0; f < nb_faces; f++)
654 {
655 int compteur = 0;
656 int som_f = 0;
657 while ((som_f < nb_som_face) and (sommet_interp(0, faces_sommets(f,som_f))==1))
658 {
659 compteur++;
660 som_f++;
661 }
662 if (compteur == nb_som_face)
663 {
664 for (int som_=0; som_ < nb_som_face; som_++)
665 {
666 int i = faces_sommets(f,som_);
667 for(int j = 0; j < nb_comp; j++)
668 {
669 variable_imposee_calculee(f,j) += variable_imposee_sommet(i,j)/nb_som_face;
670 }
671 }
672 }
673 }
674}
675
676/*##################################################################################################
677####################################################################################################
678################################# calculer_vitesse_imposee #########################################
679####################################################################################################
680##################################################################################################*/
681
682void Source_PDF_VEF::calculer_vitesse_imposee_power_law_tbl()
683{
684 assert(Objet_U::dimension==3);
685 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
686 int nb_faces=domaine_VEF.nb_faces();
687 int nb_comp = Objet_U::dimension;
688 const Domaine& dom = domaine_VEF.domaine();
689 const IntTab& faces_sommets=domaine_VEF.face_sommets();
690 int nb_som_face=faces_sommets.dimension(1);
691 int nb_som=domaine_VEF.nb_som();
692 DoubleTab& vitesse_imposee_mod = modele_lu_.variable_imposee_->valeurs();
693 DoubleTab& vitesse_imposee_calculee = variable_imposee_;
694 Interpolation_IBM_power_law_tbl& interp = ref_cast(Interpolation_IBM_power_law_tbl,interpolation_lue_.valeur());
695 int form_WJSP = interp.get_formulation_WJSP();
696 DoubleTab& fluid_points = interp.fluid_points_->valeurs();
697 DoubleTab& solid_points = interp.solid_points_->valeurs();
698 DoubleTab& fluid_elems = interp.fluid_elems_->valeurs();
699 double A_pwl = interp.get_A_pwl(form_WJSP);
700 double B_pwl = interp.get_B_pwl();
701 double y_c_p_pwl = 0;
702 double C_pwl_WJSP = 0;
703 double D_pwl_WJSP = 0;
704 double p_pwl_WJSP = 0;
705 double y_c1_p_pwl_WJSP = 0;
706 double y_c2_p_pwl_WJSP = 0;
707 if (!form_WJSP)
708 {
709 y_c_p_pwl = interp.get_y_c_p_pwl();
710 }
711 else
712 {
713 C_pwl_WJSP = interp.get_C_pwl_WJSP();
714 D_pwl_WJSP = interp.get_D_pwl_WJSP();
715 p_pwl_WJSP = interp.get_p_pwl_WJSP();
716 y_c1_p_pwl_WJSP = interp.get_y_c1_p_pwl_WJSP();
717 y_c2_p_pwl_WJSP = interp.get_y_c2_p_pwl_WJSP();
718 }
719 int impr_yplus = interp.get_impr() ;
720 Champ_P1NC& champ_vitesse_inconnue = ref_cast(Champ_P1NC,equation().inconnue());
721 DoubleTab& val_vitesse_inconnue = champ_vitesse_inconnue.valeurs();
722 if (nb_comp != val_vitesse_inconnue.dimension(1))
723 {
724 Cerr<<"Source_PDF_VEF::calculer_vitesse_imposee_power_law_tbl: variable must be a vector of dimension "<<Objet_U::dimension<<" ! "<<finl;
725 abort();
726 }
727 double form_lin_pwl = interp.get_formulation_linear_pwl();
728 DoubleTab xf(1, nb_comp);
729 DoubleTab vf(1, nb_comp);
730 ArrOfInt cells(1);
731 DoubleTab vitesse_imposee_sommet(nb_som,nb_comp);
732 IntTab sommet_interp(1, nb_som);
733 vitesse_imposee_calculee = 0.0;
734 vitesse_imposee_sommet = 0.0;
735 sommet_interp = 0;
736 // operateur(0) : diffusivite
737 if (equation().nombre_d_operateurs()<1)
738 {
739 Cerr << "Source_PDF_VEF : nombre_d_operateurs = "<<equation().nombre_d_operateurs()<<" < 1"<<finl;
740 exit();
741 }
742
743 int flag = 0;
744 const Op_Diff_VEF_base& opdiffu = ref_cast(Op_Diff_VEF_base,equation().operateur(0).l_op_base());
745 flag = opdiffu.diffusivite().valeurs().dimension(0)>1 ? 1 : 0;
746
747 double eps = 1e-12;
748 double d1_min = 1.0e+10;
749 double d1_max = 0.0;
750 double d1_mean = 0.0;
751 double yplus_min = 1.0e+10;
752 double yplus_max = 0.0;
753 double yplus_mean = 0.0;
754 double u_tau_min = 1.0e+10;
755 double u_tau_max = 0.0;
756 double u_tau_mean = 0.0;
757 double yplus_ref_min = 1.0e+10;
758 double yplus_ref_max = 0.0;
759 double yplus_ref_mean = 0.0;
760 double u_tau_ref_min = 1.0e+10;
761 double u_tau_ref_max = 0.0;
762 double u_tau_ref_mean = 0.0;
763 double h_yplus_min = 1.0e+10;
764 double h_yplus_max = 0.0;
765 double h_yplus_mean = 0.0;
766 int yplus_count=0 ;
767 double y_plus=0;
768
769 int N = interp.get_N_histo();
770 DoubleTab tab_h(1, nb_som);
771 DoubleTab abs_h(1, N+1);
772 DoubleTab compteur_h(1, N);
773 compteur_h = 0.;
774 tab_h = 0.;
775
776 tab_u_star_ibm_ = 0.;
777 tab_y_plus_ibm_ = 0.;
778
779 for (int i = 0; i < nb_som; i++)
780 {
781 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = vitesse_imposee_mod(i,j);
782 int itisok = 1;
783 if ((fluid_elems(i) >= 0.0))// and (is_dirichlet(i) == 1))
784 {
785 sommet_interp(0, i) = 1;
786 double d1 = 0.0;
787 double d2 = 0.0;
788 DoubleTab normale(1, nb_comp);
789 double norme_de_la_normale= 0.0;
790 for(int j = 0; j < nb_comp; j++)
791 {
792 double xj = dom.coord(i,j);
793 double xjf = fluid_points(i,j);
794 xf(0, j) = xjf;
795 double xjs = solid_points(i,j);
796 d1 += (xj-xjs)*(xj-xjs);
797 d2 += (xjf-xj)*(xjf-xj);
798 normale(0,j) = xjf - xjs;
799 norme_de_la_normale += normale(0,j)*normale(0,j);
800 }
801 d1 = sqrt(d1);
802 d2 = sqrt(d2);
803 double y_ref= d1+d2;
804
805 // traitement des exceptions
806 if (d1 < eps) itisok = 0; //le point P est sur la frontiere immergee : affectation
807 norme_de_la_normale = sqrt(norme_de_la_normale);
808 if ( norme_de_la_normale > eps )
809 for(int j = 0; j < nb_comp; j++) normale(0,j) /= norme_de_la_normale; // la on met la norme unite tout en gardant la direction, on normalise quoi
810 else
811 {
812 for(int j = 0; j < nb_comp; j++) normale(0,j) = 0.; // le point fluide est sur la frontiere immergee : affectation
813 itisok = 0;
814 }
815 if ( y_ref < eps ) // le point fluide est sur la frontiere immergee : affectation
816 {
817 y_ref = 1.;
818 itisok = 0;
819 }
820 // calcul vitesse power law
821 if (itisok)
822 {
823 cells(0) = int(fluid_elems(i));
824 champ_vitesse_inconnue.valeur_a_elem(xf, vf, cells[0]); // vf la vitesse totale interpolee au pt fluide
825 double Vn = 0.;
826 for(int j = 0; j < nb_comp; j++) Vn += vf(0, j) * normale(0,j);
827 DoubleTab v_ref_t(1, nb_comp);
828
829 for(int j = 0; j < nb_comp; j++) v_ref_t(0, j) =vf(0, j) - Vn*normale(0,j);
830
831 double norme_v_ref_t = 0.;
832 for(int j = 0; j < nb_comp; j++) norme_v_ref_t += v_ref_t(0, j)*v_ref_t(0,j) ;
833 norme_v_ref_t = sqrt ( norme_v_ref_t );
834 if (norme_v_ref_t < 1.0e-6) norme_v_ref_t = 1.e-6;
835
836 double nu = (flag ? opdiffu.diffusivite().valeurs()(cells) : opdiffu.diffusivite().valeurs()(0,0));
837 // On calcule U_tau_ref pour pouvoir calculer les quantites adimensionees ( y+ ...)
838 // Hypothese : sous-couche puissance
839 double u_tau_ref = pow ( norme_v_ref_t , (1/(1+B_pwl)) ) * pow ( A_pwl, (-1/(1+B_pwl)) ) * pow ( y_ref , (-B_pwl/(1+B_pwl)) ) * pow ( nu,(B_pwl/(1+B_pwl)) ) ;
840 // Cerr<<"u_tau_ref y_ref nu ="<<u_tau_ref<<" "<<y_ref<<" "<<nu<<finl;;
841 double y_ref_p = y_ref * u_tau_ref / nu; //la on a enfin tout ce qu'il faut pour etablir la loi polynomiale pour y r+
842 double test_ref;
843
844 if (!form_WJSP)
845 {
846 if ( y_ref_p > y_c_p_pwl) // a partir de la commence l'expression de la loi de paroi polynomiale turbulente
847 {
848 if (!form_lin_pwl)
849 {
850 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) * pow ( d1 / y_ref , B_pwl ) ;
851 }
852 else // Formulation lineaire de la loi polynomilale ------------------------------
853 {
854 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) *(1. - B_pwl*(1-d1/y_ref)) ;
855 }
856 test_ref = 1.;
857 // Cerr << "zone log/sous-couche inertielle en coherence avec le calcul de u tau" << finl;
858 }
859 else
860 {
861 // En fait, on est en sous-couche visqueuse
862 if (!form_lin_pwl)
863 {
864 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) * ( d1 / y_ref ) ;
865 }
866 else
867 {
868 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) * (1. - pow(u_tau_ref, 2.)*(y_ref-d1)/(nu *norme_v_ref_t));
869 }
870 test_ref = -1. ;
871 u_tau_ref = pow ( (nu * norme_v_ref_t / y_ref) , 0.5 ) ; // on recalcule u_tau et y_plus en lineaire au pt reference
872 y_ref_p = y_ref * u_tau_ref / nu;
873 if ( y_ref_p > y_c_p_pwl )
874 {
875 // Incoherence : on n utilise pas de lois de paroi
876 itisok = 0;
877 }
878 // Cerr << "zone lineaire/sous-couche visqueuse" << finl;
879 } // Fin LdPturb
880
881 double norme_v_imp_c = 0;
882 for(int j=0 ; j < nb_comp; j++) norme_v_imp_c += vitesse_imposee_sommet(i ,j) * vitesse_imposee_sommet(i ,j);
883 norme_v_imp_c = sqrt( norme_v_imp_c );
884
885 double u_tau = pow ( norme_v_imp_c , (1/(1+B_pwl)) ) * pow ( A_pwl, (-1/(1+B_pwl)) ) * pow ( d1 , (-B_pwl/(1+B_pwl)) ) * pow ( nu,(B_pwl/(1+B_pwl)) ) ; // Hypothese : sous-couche puissance
886 y_plus = u_tau * d1 / nu;
887 double test_;
888 if (y_plus >y_c_p_pwl)
889 {
890 test_ = 1.;
891 // Cerr<<"u_tau ="<<u_tau<<" y+ = "<<y_plus<<" y+ref = "<<y_ref_p<<finl;
892 }
893 else
894 {
895 test_ = -1.;
896 u_tau = pow ( (nu * norme_v_imp_c/ d1) , 0.5 ) ; // on recalcule u_tau et y_plus en lineaire au pt force
897 y_plus = u_tau * d1 / nu;
898 if ( y_plus > y_c_p_pwl )
899 {
900 // Incoherence : on n utilise pas de lois de paroi
901 Cerr<<"INCOHERENCE: u_tau ="<<u_tau<<" y+ = "<<y_plus<<" y+ref = "<<y_ref_p<<finl;
902 itisok = 0;
903 }
904 }
905 // Cerr<<"u_tau d1 nu ="<<u_tau<<" "<<d1<<" "<<nu<<finl;;
906
907 // ici on effectue le test pour savoir si le noeud de frontiere et le point fluide se trouvent dans la meme zone: si oui ok, si non on impose la vitesse
908
909 // Traitement des exceptions
910 if ( (test_ * test_ref < 0) || (itisok == 0) ) //si non on affecte la vitesse paroi
911 {
912 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = vitesse_imposee_mod(i,j);
913
914 // Cerr << "erreur" << finl;
915 itisok = 0;
916 }
917 else
918 {
919 tab_u_star_ibm_(i) = u_tau;
920 tab_y_plus_ibm_(i) = y_plus;
921 }
922 if (impr_yplus && itisok && (indicateur_nodal_champ_aire_(i)==1.))
923 {
924 if (d1 > d1_max) d1_max = d1;
925 if (d1 < d1_min) d1_min = d1;
926 d1_mean += d1;
927 if (y_plus > yplus_max) yplus_max = y_plus;
928 if (y_plus < yplus_min) yplus_min = y_plus;
929 yplus_mean += y_plus;
930 if (u_tau > u_tau_max) u_tau_max = u_tau;
931 if (u_tau < u_tau_min) u_tau_min = u_tau;
932 u_tau_mean += u_tau;
933 if (y_ref_p > yplus_ref_max) yplus_ref_max = y_ref_p;
934 if (y_ref_p < yplus_ref_min) yplus_ref_min = y_ref_p;
935 yplus_ref_mean += y_ref_p;
936 if (u_tau_ref > u_tau_ref_max) u_tau_ref_max = u_tau_ref;
937 if (u_tau_ref < u_tau_ref_min) u_tau_ref_min = u_tau_ref;
938 u_tau_ref_mean += u_tau_ref;
939 // Distribution des y+ au voisinage des parois
940 if (y_plus > h_yplus_max) h_yplus_max = y_plus;
941 if (y_plus < h_yplus_min) h_yplus_min = y_plus;
942 h_yplus_mean += y_plus;
943 tab_h(0,yplus_count) = y_plus;
944 yplus_count += 1;
945 }
946 }
947 else
948 {
949 //double u_tau;
950 if ( y_ref_p > y_c2_p_pwl_WJSP) // a partir de la commence l'expression de la loi de paroi polynomiale turbulente
951 {
952 y_plus = u_tau_ref * d1/nu;
953 if (y_plus > y_c2_p_pwl_WJSP)
954 {
955 if (!form_lin_pwl)
956 {
957 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) * pow ( d1 / y_ref , B_pwl ) ;
958 }
959 else // Formulation lineaire de la loi polynomilale ------------------------------
960 {
961 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) * (1. - B_pwl*(1-d1/y_ref)) ;
962 }
963 }
964 else if (y_plus < y_c2_p_pwl_WJSP and y_plus > y_c1_p_pwl_WJSP)
965 {
966 if (!form_lin_pwl)
967 {
968 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = (C_pwl_WJSP*pow(y_plus,2*p_pwl_WJSP-1) *u_tau_ref + D_pwl_WJSP*pow(y_plus,p_pwl_WJSP-1) *u_tau_ref)* v_ref_t(0,j)/norme_v_ref_t;
969 }
970 else // Formulation lineaire de la loi polynomilale ------------------------------
971 {
972 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = 0 ;
973 }
974 }
975 else
976 {
977 if (!form_lin_pwl)
978 {
979 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = y_plus * u_tau_ref * v_ref_t(0,j)/norme_v_ref_t;
980 }
981 else // Formulation lineaire de la loi polynomilale ------------------------------
982 {
983 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) * (1. - B_pwl*(1-d1/y_ref)) ;
984 }
985 }
986 // Cerr << "zone log/sous-couche inertielle en coherence avec le calcul de u tau" << finl;
987 }
988 else // hypothèse buffer layer
989 {
990 u_tau_ref = (nu/y_ref) * pow(-(D_pwl_WJSP/(2*C_pwl_WJSP)) * (1 + pow(1 + 4 * C_pwl_WJSP/pow(D_pwl_WJSP,2) * norme_v_ref_t * y_ref /nu ,1/2)),1/p_pwl_WJSP) ;
991 y_ref_p = y_ref * u_tau_ref / nu;
992 if (y_ref_p < y_c2_p_pwl_WJSP and y_ref_p > y_c1_p_pwl_WJSP)
993 {
994 y_plus = u_tau_ref * d1/nu;
995 if (y_plus < y_c2_p_pwl_WJSP and y_plus > y_c1_p_pwl_WJSP)
996 {
997 if (!form_lin_pwl)
998 {
999 double phi = -pow(1 + 4 * C_pwl_WJSP/pow(D_pwl_WJSP,2) * norme_v_ref_t * y_ref /nu,1/2) - 1;
1000 double u_norm = pow(D_pwl_WJSP,2) * nu / (4 * C_pwl_WJSP * y_ref) * (pow ( d1 / y_ref , 2*p_pwl_WJSP - 1 ) * pow(phi,2) + 2 * pow ( d1 / y_ref , p_pwl_WJSP - 1 ) * phi);
1001 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = u_norm * v_ref_t(0,j)/norme_v_ref_t;
1002 }
1003 else // Formulation lineaire de la loi polynomilale ------------------------------
1004 {
1005 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = 0 ; //à faire
1006 }
1007 }
1008 else
1009 {
1010 if (!form_lin_pwl)
1011 {
1012 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = y_plus * u_tau_ref * v_ref_t(0,j)/norme_v_ref_t ;
1013 }
1014 else // Formulation lineaire de la loi polynomilale ------------------------------
1015 {
1016 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = 0 ;
1017 }
1018 }
1019 }
1020 else
1021 {
1022 // En fait, on est en sous-couche visqueuse
1023 if (!form_lin_pwl)
1024 {
1025 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) * ( d1 / y_ref ) ;
1026 }
1027 else
1028 {
1029 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = v_ref_t(0,j) * (1. - pow(u_tau_ref, 2.)*(y_ref-d1)/(nu *norme_v_ref_t));
1030 }
1031 u_tau_ref = pow ( (nu * norme_v_ref_t / y_ref) , 0.5 ) ; // on recalcule u_tau et y_plus en lineaire au pt reference
1032 y_ref_p = y_ref * u_tau_ref / nu;
1033 if ( y_ref_p > y_c1_p_pwl_WJSP)
1034 {
1035 // Incoherence : on n utilise pas de lois de paroi
1036 itisok = 0;
1037 Cerr << "Incoherence" << finl;
1038 }
1039 }
1040 }
1041 // Cerr<<"u_tau d1 nu ="<<u_tau<<" "<<d1<<" "<<nu<<finl;;
1042
1043 // ici on effectue le test pour savoir si le noeud de frontiere et le point fluide se trouvent dans la meme zone: si oui ok, si non on impose la vitesse
1044
1045 // Traitement des exceptions
1046 if (itisok == 0) //si non on affecte la vitesse paroi
1047 {
1048 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = vitesse_imposee_mod(i,j);
1049
1050 // Cerr << "erreur" << finl;
1051 itisok = 0;
1052 }
1053 else
1054 {
1055 tab_u_star_ibm_(i) = u_tau_ref;
1056 tab_y_plus_ibm_(i) = y_plus;
1057 }
1058 if (impr_yplus && itisok && (indicateur_nodal_champ_aire_(i)==1.))
1059 {
1060 if (d1 > d1_max) d1_max = d1;
1061 if (d1 < d1_min) d1_min = d1;
1062 d1_mean += d1;
1063 if (y_plus > yplus_max) yplus_max = y_plus;
1064 if (y_plus < yplus_min) yplus_min = y_plus;
1065 yplus_mean += y_plus;
1066 if (u_tau_ref > u_tau_max) u_tau_max = u_tau_ref;
1067 if (u_tau_ref < u_tau_min) u_tau_min = u_tau_ref;
1068 u_tau_mean += u_tau_ref;
1069 if (y_ref_p > yplus_ref_max) yplus_ref_max = y_ref_p;
1070 if (y_ref_p < yplus_ref_min) yplus_ref_min = y_ref_p;
1071 yplus_ref_mean += y_ref_p;
1072 if (u_tau_ref > u_tau_ref_max) u_tau_ref_max = u_tau_ref;
1073 if (u_tau_ref < u_tau_ref_min) u_tau_ref_min = u_tau_ref;
1074 u_tau_ref_mean += u_tau_ref;
1075 // Distribution des y+ au voisinage des parois
1076 if (y_plus > h_yplus_max) h_yplus_max = y_plus;
1077 if (y_plus < h_yplus_min) h_yplus_min = y_plus;
1078 h_yplus_mean += y_plus;
1079 tab_h(0,yplus_count) = y_plus;
1080 yplus_count += 1;
1081 }
1082 }
1083 }
1084 }
1085 }
1086 get_champ("u_star_ibm");
1087 get_champ("y_plus_ibm");
1088 for (int f = 0; f < nb_faces; f++)
1089 {
1090 int compteur = 0;
1091 int som_f = 0;
1092 while ((som_f < nb_som_face) and (sommet_interp(0, faces_sommets(f,som_f))==1))
1093 {
1094 compteur++;
1095 som_f++;
1096 }
1097 if (compteur == nb_som_face)
1098 {
1099 for (int som_=0; som_ < nb_som_face; som_++)
1100 {
1101 int i = faces_sommets(f,som_);
1102 for(int j = 0; j < nb_comp; j++)
1103 {
1104 vitesse_imposee_calculee(f,j) += vitesse_imposee_sommet(i,j)/nb_som_face;
1105 }
1106 }
1107 }
1108 }
1109
1110 if (impr_yplus && (yplus_count >= 1) )
1111 {
1112 d1_mean /= yplus_count;
1113 yplus_mean /= yplus_count;
1114 yplus_ref_mean /= yplus_count;
1115 u_tau_mean /= yplus_count;
1116 u_tau_ref_mean /= yplus_count;
1117 Cerr<<"min mean max y = "<<d1_min<<" "<<d1_mean<<" "<<d1_max<<finl;
1118 Cerr<<"min mean max y+ = "<<yplus_min<<" "<<yplus_mean<<" "<<yplus_max<<finl;
1119 Cerr<<"min mean u_tau = "<<u_tau_min<<" "<<u_tau_mean<<" "<<u_tau_max<<finl;
1120 Cerr<<"min mean max y+_ref = "<<yplus_ref_min<<" "<<yplus_ref_mean<<" "<<yplus_ref_max<<finl;
1121 Cerr<<"min mean u_tau_ref = "<<u_tau_ref_min<<" "<<u_tau_ref_mean<<" "<<u_tau_ref_max<<finl;
1122
1123 h_yplus_mean /= yplus_count;
1124 double range = ( h_yplus_max - h_yplus_min)/N;
1125 for (int k = 0; k < N+1; k++) abs_h(0,k) = (h_yplus_min + k*range);
1126 for (int i=0; i < yplus_count; i++)
1127 {
1128 for (int j = 0; j < N; j++)
1129 {
1130 if ( abs_h(0,j) < tab_h(0,i) && tab_h(0,i) <= abs_h(0,j+1) ) compteur_h(0,j) +=1;
1131 }
1132 }
1133 Cerr<<"histogramme y+ compteur =";
1134 for (int j = 0; j < N; j++) Cerr<<" "<<compteur_h(0,j)<<" ";
1135 Cerr<<finl;
1136 Cerr<<"histogramme y+ abscisse =";
1137 for (int j = 0; j < N+1; j++) Cerr<<" "<<abs_h(0,j)<<" ";
1138 Cerr<<finl;
1139 Cerr<<"histogramme y+ : min = "<<h_yplus_min<<" - mean = "<<h_yplus_mean<<" - max = "<<h_yplus_max<<finl;
1140
1141 Cerr<<"Moyenne sur "<<yplus_count<<" points"<<finl;
1142 }
1143}
1144
1145void Source_PDF_VEF::calculer_vitesse_imposee_power_law_tbl_u_star()
1146{
1147 assert(Objet_U::dimension==3);
1148 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
1149 int nb_faces=domaine_VEF.nb_faces();
1150 int nb_comp = Objet_U::dimension;
1151 const Domaine& dom = domaine_VEF.domaine();
1152 const IntTab& faces_sommets=domaine_VEF.face_sommets();
1153 int nb_som_face=faces_sommets.dimension(1);
1154 int nb_som = domaine_VEF.nb_som();
1155 DoubleTab& vitesse_imposee_mod = modele_lu_.variable_imposee_->valeurs();
1156 DoubleTab& vitesse_imposee_calculee = variable_imposee_;
1157 Interpolation_IBM_power_law_tbl_u_star& interp = ref_cast(Interpolation_IBM_power_law_tbl_u_star,interpolation_lue_.valeur());
1158 DoubleTab& is_dirichlet = interp.is_dirichlet_->valeurs();
1159 DoubleTab& solid_points = interp.solid_points_->valeurs();
1160 DoubleTab& solid_elems = interp.solid_elems_->valeurs();
1161 double A_pwl = interp.get_A_pwl(0);
1162 double B_pwl = interp.get_B_pwl();
1163 double y_c_p_pwl = interp.get_y_c_p_pwl();
1164 int impr_yplus = interp.get_impr() ;
1165 DoubleTab& vitesse_inconnue = equation().inconnue().valeurs();
1166 DoubleTab vitesse_imposee_sommet(nb_som,nb_comp);
1167 IntTab sommet_interp(1, nb_som);
1168 vitesse_imposee_calculee = 0.0;
1169 vitesse_imposee_sommet = 0.0;
1170 sommet_interp = 0;
1171 if (nb_comp != vitesse_inconnue.dimension(1))
1172 {
1173 Cerr<<"Source_PDF_EF::calculer_vitesse_imposee_power_law_tbl_u_star: variable must be a vector of dimension "<<Objet_U::dimension<<" ! "<<finl;
1174 abort();
1175 }
1176
1177 ArrOfInt cells(1);
1178
1179 // operateur(0) : diffusivite
1180 if (equation().nombre_d_operateurs()<1)
1181 {
1182 Cerr << "Source_PDF_VEF : nombre_d_operateurs = "<<equation().nombre_d_operateurs()<<" < 1"<<finl;
1183 exit();
1184 }
1185
1186 int flag = 0;
1187 const Op_Diff_VEF_base& opdiffu = ref_cast(Op_Diff_VEF_base,equation().operateur(0).l_op_base());
1188 flag = opdiffu.diffusivite().valeurs().dimension(0)>1 ? 1 : 0;
1189
1190 double eps = 1e-12;
1191 double d1_min = 1.0e+10;
1192 double d1_max = 0.0;
1193 double d1_mean = 0.0;
1194 double yplus_min = 1.0e+10;
1195 double yplus_max = 0.0;
1196 double yplus_mean = 0.0;
1197 double u_tau_min = 1.0e+10;
1198 double u_tau_max = 0.0;
1199 double u_tau_mean = 0.0;
1200 double yplus_ref_min = 1.0e+10;
1201 double yplus_ref_max = 0.0;
1202 double yplus_ref_mean = 0.0;
1203 double u_tau_ref_min = 1.0e+10;
1204 double u_tau_ref_max = 0.0;
1205 double u_tau_ref_mean = 0.0;
1206 double h_yplus_min = 1.0e+10;
1207 double h_yplus_max = 0.0;
1208 double h_yplus_mean = 0.0;
1209 int yplus_count=0 ;
1210 int yplus_ref_count=0 ;
1211
1212 int N = interp.get_N_histo();
1213 DoubleTab tab_h(1, nb_som);
1214 DoubleTab abs_h(1, N+1);
1215 DoubleTab compteur_h(1, N);
1216 compteur_h = 0.;
1217 tab_h = 0.;
1218
1219 vitesse_imposee_mod.echange_espace_virtuel();
1220 tab_u_star_ibm_ = 0.;
1221 tab_y_plus_ibm_ = 0.;
1222 // DoubleTrav nb_vois(is_dirichlet);
1223
1224 for (int i = 0; i < nb_som; i++)
1225 {
1226 double x = 0.;
1227 if (is_dirichlet(i) > 0.0)
1228 {
1229 sommet_interp(0, i) = 1;
1230 int itisok_P = 1;
1231 double d1P = 0.0;
1232 DoubleTab normaleP(1, nb_comp);
1233 double norme_de_la_normaleP= 0.0;
1234 // normal et distance normale au point P
1235 for(int j = 0; j < nb_comp; j++)
1236 {
1237 vitesse_imposee_sommet(i,j) = vitesse_imposee_mod(i,j);
1238 x = dom.coord(i,j);
1239 double xp = solid_points(i,j);
1240 normaleP(0,j) = x - xp;
1241 norme_de_la_normaleP += normaleP(0,j)*normaleP(0,j);
1242 }
1243 norme_de_la_normaleP = sqrt(norme_de_la_normaleP);
1244 d1P = norme_de_la_normaleP;
1245 if (d1P > eps)
1246 for(int j = 0; j < nb_comp; j++) normaleP(0,j) /= norme_de_la_normaleP;
1247 else normaleP = 0.;
1248
1249 cells(0) = int(solid_elems(i)); //pour definir la viscosite laminaire
1250
1251 double nu = (flag ? opdiffu.diffusivite().valeurs()(cells) : opdiffu.diffusivite().valeurs()(0,0));
1252
1253 double u_tau = 0.;
1254 double y_plus = 0.;
1255 if (d1P > eps)
1256 //le point P n'est pas sur la frontiere immergee
1257 {
1258 IntList& voisins = interp.getSommetsVoisinsOf(i);
1259 if (!(voisins.est_vide()))
1260 // il y a une liste de dof voisins
1261 {
1262 double u_tau_ref;
1263 double y_plus_ref;
1264 DoubleTab mean_u_tau_neighbour(1, nb_comp); // moyenne de u_tau vectoriel sur les dof voisins
1265 mean_u_tau_neighbour = 0.0;
1266 double mean_y_plus_neighbour = 0.; // moyenne de y_plus sur les dof voisins
1267 int taille = voisins.size();
1268 // nb_vois(i) = taille * 1.0 ;
1269 double pond_tot = 0.;
1270 double ponderation_k;
1271 for (int k = 0; k < taille; k++)
1272 // boucle sur les dof voisins
1273 {
1274 int num_som = voisins[k];
1275 int itisok = 1;
1276 double d1 = 0.0;
1277 DoubleTab a(1, nb_comp);
1278 DoubleTab b(1, nb_comp);
1279 double norme_a = 0.0;
1280 double norme_b = 0.0;
1281
1282 for (int j = 0; j < nb_comp; j++)
1283 {
1284 double xf = dom.coord(num_som,j);
1285 double xpf = solid_points(num_som,j);
1286 d1 += (xf - xpf)*(xf - xpf);
1287 a(0,j) = (xf - x);
1288 b(0,j) = (xf - x)*normaleP(0,j);
1289 norme_a += a(0,j)*a(0,j);
1290 norme_b += b(0,j)*b(0,j);
1291 }
1292 d1 = sqrt(d1);
1293 if ( d1 < eps ) itisok = 0;
1294
1295 if ( (norme_a - norme_b) > eps )
1296 ponderation_k = 1/(sqrt(norme_a - norme_b));
1297 else
1298 ponderation_k = 1./eps;
1299
1300 u_tau_ref = 0.;
1301 y_plus_ref = 0.;
1302 if (itisok)
1303 {
1304 // normale et tangente vitesse vis a vis de normalP
1305 double Vn = 0.;
1306 for(int j = 0; j < nb_comp; j++) Vn += vitesse_inconnue(num_som,j) * normaleP(0,j);
1307 DoubleTab vtf_k(1, nb_comp);
1308 for(int j = 0; j < nb_comp; j++) vtf_k(0, j) = vitesse_inconnue(num_som, j) - Vn * normaleP(0,j);
1309 double norme_vtf_k = 0.;
1310 for(int j = 0; j < nb_comp; j++) norme_vtf_k += vtf_k(0, j)*vtf_k(0,j);
1311 norme_vtf_k = sqrt(norme_vtf_k);
1312 u_tau_ref = pow ( norme_vtf_k , (1/(1+B_pwl)) ) * pow ( A_pwl, (-1/(1+B_pwl)) ) * pow ( d1 , (-B_pwl/(1+B_pwl)) ) * pow ( nu,(B_pwl/(1+B_pwl)) ) ; // vitesse de frottement power law au point fluide k
1313 y_plus_ref = d1 * u_tau_ref / nu;
1314 if (y_plus_ref < y_c_p_pwl ) // En fait, on est en lois lineaire !
1315 {
1316 u_tau_ref = pow ( (nu * norme_vtf_k / d1) , 0.5 ) ; // on recalcule u_tau_ref et y_plus_ref en lineaire au point fluide k
1317 y_plus_ref = d1 * u_tau_ref / nu;
1318 if ( y_plus_ref > y_c_p_pwl )
1319 {
1320 // Incoherence : on n utilise pas de lois de paroi
1321 itisok = 0;
1322 }
1323 }
1324
1325 // On moyenne les vecteurs vitesses de frottement tangents (vis a vis de P) par l'inverse des distances
1326 if (norme_vtf_k < eps) itisok = 0;
1327 if (itisok)
1328 {
1329 for(int j = 0; j < nb_comp; j++) mean_u_tau_neighbour(0,j) += (u_tau_ref * vtf_k(0, j) / norme_vtf_k) * ponderation_k;
1330 mean_y_plus_neighbour += y_plus_ref * ponderation_k;
1331 pond_tot += ponderation_k;
1332 }
1333 }
1334
1335 if (impr_yplus && itisok)
1336 {
1337 if (y_plus_ref > yplus_ref_max) yplus_ref_max = y_plus_ref;
1338 if (y_plus_ref < yplus_ref_min) yplus_ref_min = y_plus_ref;
1339 yplus_ref_mean += y_plus_ref;
1340 if (u_tau_ref > u_tau_ref_max) u_tau_ref_max = u_tau_ref;
1341 if (u_tau_ref < u_tau_ref_min) u_tau_ref_min = u_tau_ref;
1342 u_tau_ref_mean += u_tau_ref;
1343 yplus_ref_count += 1;
1344 }
1345
1346 }
1347
1348 // On calcul u_tau moyen, y_plus et la vitesse tangente au point force P
1349 if (pond_tot > 0.)
1350 {
1351 mean_u_tau_neighbour /= pond_tot; // u_tau moyen vectoriel
1352 mean_y_plus_neighbour /= pond_tot; // y_plus moyen
1353 u_tau = 0. ;
1354 for(int j = 0; j < nb_comp; j++) u_tau += mean_u_tau_neighbour(0,j) * mean_u_tau_neighbour(0,j);
1355 u_tau = sqrt(u_tau); // u_tau = norme de u_tau moyen vectoriel
1356 y_plus = u_tau * d1P / nu;
1357
1358 if (u_tau > eps) mean_u_tau_neighbour /= u_tau;
1359 else mean_u_tau_neighbour = 0. ;// tangente en P
1360 if (y_plus > y_c_p_pwl ) // loi polynomiale
1361 {
1362 if (mean_y_plus_neighbour > y_c_p_pwl ) // coherence loi polynomiale
1363 {
1364 double vit_coeff = A_pwl * pow (d1P, B_pwl) / pow (nu, B_pwl);
1365 for (int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = pow (u_tau, (1+B_pwl)) * vit_coeff * mean_u_tau_neighbour(0,j);
1366 }
1367 else itisok_P = 0; // non coherence loi polynomiale
1368 }
1369 else // loi lineaire
1370 {
1371 if (mean_y_plus_neighbour < y_c_p_pwl ) // coherence loi lineaire
1372 {
1373 for (int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = (d1P * u_tau * u_tau / nu) * mean_u_tau_neighbour(0,j);
1374 }
1375 else itisok_P = 0; // non coherence loi lineaire
1376 }
1377 }
1378 // il n y a pas de contribution des dof voisins
1379 else itisok_P = 0;
1380 }
1381 else
1382 // il n y a pas de liste des dof voisins
1383 itisok_P = 0;
1384 }
1385 //le point P est sur la frontiere immergee : affectation
1386 else itisok_P = 0;
1387
1388 // pour post-pro
1389 if(itisok_P)
1390 {
1391 // pour post-pro & nécessité pour T wall law
1392 tab_u_star_ibm_(i) = u_tau;
1393 tab_y_plus_ibm_(i) = y_plus;
1394 }
1395 if (impr_yplus && itisok_P)
1396 {
1397 if (d1P > d1_max) d1_max = d1P;
1398 if (d1P < d1_min) d1_min = d1P;
1399 d1_mean += d1P;
1400 if (y_plus > yplus_max) yplus_max = y_plus;
1401 if (y_plus < yplus_min) yplus_min = y_plus;
1402 yplus_mean += y_plus;
1403 if (u_tau > u_tau_max) u_tau_max = u_tau;
1404 if (u_tau < u_tau_min) u_tau_min = u_tau;
1405 u_tau_mean += u_tau;
1406 // Distribution des y+ au voisinage des parois
1407 if (y_plus > h_yplus_max) h_yplus_max = y_plus;
1408 if (y_plus < h_yplus_min) h_yplus_min = y_plus;
1409 h_yplus_mean += y_plus;
1410 tab_h(0,yplus_count) = y_plus;
1411 yplus_count += 1;
1412 }
1413 }
1414 else
1415 {
1416 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = vitesse_imposee_mod(i,j);
1417 }
1418 }
1419 //vitesse_imposee_calculee.echange_espace_virtuel();
1420 for (int f = 0; f < nb_faces; f++)
1421 {
1422 int compteur = 0;
1423 int som_f = 0;
1424 while ((som_f < nb_som_face) and (sommet_interp(0, faces_sommets(f,som_f))==1))
1425 {
1426 compteur++;
1427 som_f++;
1428 }
1429 if (compteur == nb_som_face)
1430 {
1431 for (int som_=0; som_ < nb_som_face; som_++)
1432 {
1433 int i = faces_sommets(f,som_);
1434 for(int j = 0; j < nb_comp; j++)
1435 {
1436 vitesse_imposee_sommet(f,j) += vitesse_imposee_sommet(i,j)/nb_som_face;
1437 }
1438 }
1439 }
1440 }
1441
1442 if (impr_yplus && (yplus_count >= 1) && (yplus_ref_count >= 1) )
1443 {
1444 d1_mean /= yplus_count;
1445 yplus_mean /= yplus_count;
1446 yplus_ref_mean /= yplus_ref_count;
1447 u_tau_mean /= yplus_count;
1448 u_tau_ref_mean /= yplus_ref_count;
1449 Cerr<<"min mean max y = "<<d1_min<<" "<<d1_mean<<" "<<d1_max<<finl;
1450 Cerr<<"min mean max y+ = "<<yplus_min<<" "<<yplus_mean<<" "<<yplus_max<<finl;
1451 Cerr<<"min mean u_tau = "<<u_tau_min<<" "<<u_tau_mean<<" "<<u_tau_max<<finl;
1452 Cerr<<"min mean max y+_ref = "<<yplus_ref_min<<" "<<yplus_ref_mean<<" "<<yplus_ref_max<<finl;
1453 Cerr<<"min mean u_tau_ref = "<<u_tau_ref_min<<" "<<u_tau_ref_mean<<" "<<u_tau_ref_max<<finl;
1454
1455 h_yplus_mean /= yplus_count;
1456 double range = ( h_yplus_max - h_yplus_min)/N;
1457 for (int k = 0; k < N+1; k++) abs_h(0,k) = (h_yplus_min + k*range);
1458 for (int i=0; i < yplus_count; i++)
1459 {
1460 for (int j = 0; j < N; j++)
1461 {
1462 if ( abs_h(0,j) < tab_h(0,i) && tab_h(0,i) <= abs_h(0,j+1) ) compteur_h(0,j) +=1;
1463 }
1464 }
1465 Cerr<<"histogramme y+ compteur =";
1466 for (int j = 0; j < N; j++) Cerr<<" "<<compteur_h(0,j)<<" ";
1467 Cerr<<finl;
1468 Cerr<<"histogramme y+ abscisse =";
1469 for (int j = 0; j < N+1; j++) Cerr<<" "<<abs_h(0,j)<<" ";
1470 Cerr<<finl;
1471 Cerr<<"histogramme y+ : min = "<<h_yplus_min<<" - mean = "<<h_yplus_mean<<" - max = "<<h_yplus_max<<finl;
1472
1473 Cerr<<"Moyenne sur "<<yplus_count<<" points forces et "<<yplus_ref_count<<" points references"<<finl;
1474 }
1475
1476}
1477
1478
1479/*##################################################################################################
1480####################################################################################################
1481################################# calculer_temperature_imposee #########################################
1482####################################################################################################
1483##################################################################################################*/
1484
1485
1486void Source_PDF_VEF::calculer_temperature_imposee_wall_law()
1487{
1488 int dim_esp = Objet_U::dimension;
1489 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
1490 const Domaine& dom = domaine_VEF.domaine();
1491// Variables du domaine
1492 const IntTab& faces_sommets=domaine_VEF.face_sommets();
1493 int nb_faces=domaine_VEF.nb_faces();
1494 int nb_som_face=faces_sommets.dimension(1);
1495 int nb_som = domaine_VEF.nb_som();
1496// Tableaux valeurs imposees et calc
1497 DoubleTab& variable_imposee_mod = modele_lu_.variable_imposee_->valeurs();
1498 DoubleTab& temperature_imposee_calculee = variable_imposee_;
1499// Tableaux informations de la paroi immergée
1500 Interpolation_IBM_thermal_wall_law& interp = ref_cast(Interpolation_IBM_thermal_wall_law,interpolation_lue_.valeur());
1501 DoubleTab& solid_points = interp.solid_points_->valeurs();
1502 DoubleTab& fluid_points = interp.fluid_points_->valeurs();
1503 DoubleTab& fluid_elems = interp.fluid_elems_->valeurs();
1504// Discrétisation de la grandeur calculée
1505 Champ_P1NC& champ_variable_inconnue = ref_cast(Champ_P1NC,equation().inconnue());
1506 DoubleTab& val_variable_inconnue = champ_variable_inconnue.valeurs();
1507// Récupération des données d'une loi de paroi en vitesse
1508 const Champ_base& champ_ustar_som = equation().probleme().get_champ("u_star_ibm");
1509 const DoubleTab& val_ustar = champ_ustar_som.valeurs();
1510 const Champ_base& champ_yplus = equation().probleme().get_champ("y_plus_ibm");
1511 const DoubleTab& val_yplus = champ_yplus.valeurs();
1512// Composante température
1513 int nb_comp = val_variable_inconnue.dimension(1);
1514// Variables pour calculs
1515 DoubleTab temperature_imposee_sommet(nb_som,nb_comp); // Tableau des valeurs aux sommets calc.
1516 IntTab sommet_interp(nb_som); // Tableau disant si un noeud est forcé
1517 DoubleTab x(1, dim_esp);
1518 DoubleTab xf(1, dim_esp);
1519 DoubleTab Tref(1, nb_comp);
1520 ArrOfInt cells(1);
1521 double thetap = 0.0;
1522 double thetapref = 0.0;
1523// Choix de modélisation
1524 int form_Tp = interp.get_formulation_Tp(); // Choix de la loi de paroi utilisée
1525 int boundary_type = interp.get_boundary_type(); // Choix du type de condition aux limites
1526 double T_inlet = interp.get_T_inlet();
1527 double Pr = interp.get_Prandlt_mol();
1528// Initialisation tableaux
1529 temperature_imposee_calculee = 0.0;
1530 temperature_imposee_sommet = 0.0;
1531 sommet_interp = 0;
1532// Précision
1533 double eps = 1e-6;
1534
1535// Début du calcul de T_imp aux noeuds
1536 for (int i = 0; i < nb_som; i++) // On parcoure l'ensemble des noeuds du domaine
1537 {
1538 if ((fluid_elems(i) >= 0.0)) // Noeuds forcés uniquement
1539 {
1540 int itisok = 1;
1541 sommet_interp(i) = 1; //On note que le sommet i est forcé
1542 temperature_imposee_sommet(i,0) = variable_imposee_mod(i,0);
1543 // Obtention du y+ et u* au noeud considéré
1544 double utau = val_ustar(i);
1545 double yplus = val_yplus(i);
1546 // Calcul si utau != 0 => loi classique
1547 if (utau != 0)
1548 {
1549 // Obtention des y et yref, ainsi que la normale à la paroi
1550 double d1 = 0.0;
1551 double d2 = 0.0;
1552 DoubleTab normale(1, nb_comp);
1553 double norme_de_la_normale= 0.0;
1554 for(int j = 0; j < nb_comp; j++)
1555 {
1556 double xj = dom.coord(i,j);
1557 double xjf = fluid_points(i,j);
1558 xf(0, j) = xjf;
1559 double xjs = solid_points(i,j);
1560 d1 += (xj-xjs)*(xj-xjs);
1561 d2 += (xjf-xj)*(xjf-xj);
1562 normale(0,j) = xjf - xjs;
1563 norme_de_la_normale += normale(0,j)*normale(0,j);
1564 }
1565 d1 = sqrt(d1);
1566 d2 = sqrt(d2);
1567 double y_ref= d1+d2;
1568 // traitement des exceptions
1569 if (d1 < eps) itisok = 0; // le point P est sur la frontiere immergee : affectation
1570 else
1571 {
1572 norme_de_la_normale = sqrt(norme_de_la_normale);
1573 for(int j = 0; j < nb_comp; j++) normale(0,j) /= norme_de_la_normale;
1574 }
1575 // Obtention des caractéristiques du fluide utiles au calcul
1576 double nu = d1 * utau / yplus;
1577 double rho = equation().milieu().masse_volumique().valeurs()(0,0);
1578 double Cp = equation().milieu().capacite_calorifique().valeurs()(0, 0);
1579 //double alpha = equation().milieu().diffusivite().valeurs()(0, 0);
1580 // -------
1581 cells[0] = int(fluid_elems(i)); // N° élément fluide
1582 champ_variable_inconnue.valeur_a_elem(xf, Tref, cells[0]); // Détermination de la température au point fluide de référence
1583 double yplusref = y_ref * utau / nu;
1584 if (!form_Tp) //Loi de paroi de Kader
1585 {
1586 thetap = interp.Kader(yplus,Pr); // Calcul du theta plus avec Kader pour le point forcé (connaissance de y+)
1587 thetapref = interp.Kader(yplusref,Pr); // Calcul du theta plus avec Kader pour le point fluide
1588 }
1589 else
1590 {
1591 Cerr << "Loi non implementee, choix d'une autre loi obligatoire" << finl;
1592 exit();
1593 }
1594 if ((itisok == 1) and (!boundary_type))
1595 {
1596 double Tw = variable_imposee_mod(i,0);
1597 if (Tw > T_inlet)
1598 {
1599 temperature_imposee_sommet(i,0) = Tw - abs((thetap / thetapref) * (Tref(0,0) - Tw));
1600 }
1601 else
1602 {
1603 temperature_imposee_sommet(i,0) = Tw + abs((thetap / thetapref) * (Tref(0,0) - Tw));
1604 }
1605 }
1606 else if ((itisok == 1) and (boundary_type == 1))
1607 {
1608 double qw = variable_imposee_mod(i,0);
1609 double T_tau = qw/(rho * Cp * utau);
1610 temperature_imposee_sommet(i,0) = Tref(0,0) + T_tau * (thetapref - thetap);
1611 }
1612 else
1613 {
1614 Cerr << "Erreur dans le choix du type de condition aux limites" << finl;
1615 exit();
1616 }
1617 }
1618 }
1619 }
1620
1621// Calcul de la température pour les faces forcées
1622 for (int f = 0; f < nb_faces; f++)
1623 {
1624 int compteur = 0;
1625 int som_f = 0;
1626 while ((som_f < nb_som_face) and (sommet_interp(faces_sommets(f,som_f))==1)) //On récupère les faces dont l'ensemble des noeuds sont forcées
1627 {
1628 compteur++;
1629 som_f++;
1630 }
1631 if (compteur == nb_som_face)
1632 {
1633 for (int som_faces =0; som_faces < nb_som_face; som_faces++)
1634 {
1635 int i = faces_sommets(f,som_faces);
1636 temperature_imposee_calculee(f,0) += temperature_imposee_sommet(i,0)/nb_som_face; //Calcul de la température à la face par moyenne des valeurs aux noeuds
1637 }
1638 }
1639 }
1640}
1641
1642/*void Source_PDF_VEF::calculer_temperature_imposee_mean_wall_law()
1643{
1644//à faire
1645 int dim_esp = Objet_U::dimension;
1646 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
1647 const Domaine& dom = domaine_VEF.domaine();
1648// Variables du domaine
1649 const IntTab& faces_sommets=domaine_VEF.face_sommets();
1650 int nb_faces=domaine_VEF.nb_faces();
1651 int nb_som_face=faces_sommets.dimension(1);
1652 int nb_som = domaine_VEF.nb_som();
1653// Tableaux valeurs imposees et calc
1654 DoubleTab& variable_imposee_mod = modele_lu_.variable_imposee_->valeurs();
1655 DoubleTab& temperature_imposee_calculee = variable_imposee_;
1656// Tableaux informations de la paroi immergée
1657 Interpolation_IBM_thermal_wall_law& interp = ref_cast(Interpolation_IBM_thermal_wall_law,interpolation_lue_.valeur());
1658 DoubleTab& solid_points = interp.solid_points_->valeurs();
1659 DoubleTab& fluid_points = interp.fluid_points_->valeurs();
1660 DoubleTab& fluid_elems = interp.fluid_elems_->valeurs();
1661// Discrétisation de la grandeur calculée
1662 Champ_P1NC& champ_variable_inconnue = ref_cast(Champ_P1NC,equation().inconnue());
1663 DoubleTab& val_variable_inconnue = champ_variable_inconnue.valeurs();
1664// Récupération des données d'une loi de paroi en vitesse
1665 const Champ_base& champ_ustar_som = equation().probleme().get_champ("u_star_ibm");
1666 const DoubleTab& val_ustar = champ_ustar_som.valeurs();
1667 const Champ_base& champ_yplus = equation().probleme().get_champ("y_plus_ibm");
1668 const DoubleTab& val_yplus = champ_yplus.valeurs();
1669// Composante température
1670 int nb_comp = val_variable_inconnue.dimension(1);
1671// Variables pour calculs
1672 DoubleTab temperature_imposee_sommet(nb_som,nb_comp); // Tableau des valeurs aux sommets calc.
1673 IntTab sommet_interp(nb_som); // Tableau disant si un noeud est forcé
1674 DoubleTab x(1, dim_esp);
1675 DoubleTab xf(1, dim_esp);
1676 DoubleTab Tref(1, nb_comp);
1677 ArrOfInt cells(1);
1678 double thetap = 0.0;
1679 double thetapref = 0.0;
1680// Choix de modélisation
1681 int form_Tp = interp.get_formulation_Tp(); // Choix de la loi de paroi utilisée
1682 int boundary_type = interp.get_boundary_type(); // Choix du type de condition aux limites
1683 double T_inlet = interp.get_T_inlet();
1684 double Pr = interp.get_Prandlt_mol();
1685// Initialisation tableaux
1686 temperature_imposee_calculee = 0.0;
1687 temperature_imposee_sommet = 0.0;
1688 sommet_interp = 0;
1689// Précision
1690 double eps = 1e-6;
1691
1692// Début du calcul de T_imp aux noeuds
1693 for (int i = 0; i < nb_som; i++) // On parcoure l'ensemble des noeuds du domaine
1694 {
1695 if ((fluid_elems(i) >= 0.0)) // Noeuds forcés uniquement
1696 {
1697 int itisok = 1;
1698 sommet_interp(i) = 1; //On note que le sommet i est forcé
1699 temperature_imposee_sommet(i,0) = variable_imposee_mod(i,0);
1700 // Obtention du y+ et u* au noeud considéré
1701 double utau = val_ustar(i);
1702 double yplus = val_yplus(i);
1703 // Calcul si utau != 0 => loi classique
1704 if (utau != 0)
1705 {
1706 // Obtention des y et yref, ainsi que la normale à la paroi
1707 double d1 = 0.0;
1708 double d2 = 0.0;
1709 DoubleTab normale(1, nb_comp);
1710 double norme_de_la_normale= 0.0;
1711 for(int j = 0; j < nb_comp; j++)
1712 {
1713 double xj = dom.coord(i,j);
1714 double xjf = fluid_points(i,j);
1715 xf(0, j) = xjf;
1716 double xjs = solid_points(i,j);
1717 d1 += (xj-xjs)*(xj-xjs);
1718 d2 += (xjf-xj)*(xjf-xj);
1719 normale(0,j) = xjf - xjs;
1720 norme_de_la_normale += normale(0,j)*normale(0,j);
1721 }
1722 d1 = sqrt(d1);
1723 d2 = sqrt(d2);
1724 double y_ref= d1+d2;
1725 // traitement des exceptions
1726 if (d1 < eps) itisok = 0; // le point P est sur la frontiere immergee : affectation
1727 else
1728 {
1729 norme_de_la_normale = sqrt(norme_de_la_normale);
1730 for(int j = 0; j < nb_comp; j++) normale(0,j) /= norme_de_la_normale;
1731 }
1732 // Obtention des caractéristiques du fluide utiles au calcul
1733 double nu = d1 * utau / yplus;
1734 double rho = equation().milieu().masse_volumique().valeurs()(0,0);
1735 double Cp = equation().milieu().capacite_calorifique().valeurs()(0, 0);
1736 //double alpha = equation().milieu().diffusivite().valeurs()(0, 0);
1737 double Pr = 1; //nu/alpha;
1738 // -------
1739 cells[0] = int(fluid_elems(i)); // N° élément fluide
1740 champ_variable_inconnue.valeur_a_elem(xf, Tref, cells[0]); // Détermination de la température au point fluide de référence
1741 double yplusref = y_ref * utau / nu;
1742 if (!form_Tp) //Loi de paroi de Kader
1743 {
1744 thetap = interp.Kader(yplus,Pr); // Calcul du theta plus avec Kader pour le point forcé (connaissance de y+)
1745 thetapref = interp.Kader(yplusref,Pr); // Calcul du theta plus avec Kader pour le point fluide
1746 }
1747 else
1748 {
1749 Cerr << "Loi non implementee, choix d'une autre loi obligatoire" << finl;
1750 exit();
1751 }
1752 if ((itisok == 1) and (!boundary_type))
1753 {
1754 double Tw = variable_imposee_mod(i,0);
1755 if (Tw > T_inlet)
1756 {
1757 temperature_imposee_sommet(i,0) = Tw - abs((thetap / thetapref) * (Tref(0,0) - Tw));
1758 }
1759 else
1760 {
1761 temperature_imposee_sommet(i,0) = Tw + abs((thetap / thetapref) * (Tref(0,0) - Tw));
1762 }
1763 }
1764 else if ((itisok == 1) and (boundary_type == 1))
1765 {
1766 double qw = variable_imposee_mod(i,0);
1767 double T_tau = qw/(rho * Cp * utau);
1768 temperature_imposee_sommet(i,0) = Tref(0,0) + T_tau * (thetapref - thetap);
1769 }
1770 else
1771 {
1772 Cerr << "Erreur dans le choix du type de condition aux limites" << finl;
1773 exit();
1774 }
1775 }
1776 }
1777 }
1778
1779// Calcul de la température pour les faces forcées
1780 for (int f = 0; f < nb_faces; f++)
1781 {
1782 int compteur = 0;
1783 int som_f = 0;
1784 while ((som_f < nb_som_face) and (sommet_interp(faces_sommets(f,som_f))==1)) //On récupère les faces dont l'ensemble des noeuds sont forcées
1785 {
1786 compteur++;
1787 som_f++;
1788 }
1789 if (compteur == nb_som_face)
1790 {
1791 for (int som_faces =0; som_faces < nb_som_face; som_faces++)
1792 {
1793 int i = faces_sommets(f,som_faces);
1794 temperature_imposee_calculee(f,0) += temperature_imposee_sommet(i,0)/nb_som_face; //Calcul de la température à la face par moyenne des valeurs aux noeuds
1795 }
1796 }
1797 }
1798}*/
1799
1800/*##################################################################################################
1801####################################################################################################
1802################################# Pression #########################################################
1803####################################################################################################
1804##################################################################################################*/
1805
1806void Source_PDF_VEF::correct_incr_pressure(const DoubleTab& coeff_node, DoubleTab& correction_en_pression) const
1807{
1808 const DoubleTab& aire=champ_aire_->valeurs();
1809 // int nb_elem = correction_en_pression.size();
1810 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
1811 const IntTab& elems= domaine_VEF.elem_faces();
1812 int nb_face_elem=domaine_VEF.domaine().nb_faces_elem();
1813 int nb_elems=domaine_VEF.domaine().nb_elem_tot();
1814 int ncomp = coeff_node.dimension(1) ;
1815 DoubleTrav coef_elem(nb_elems);
1816
1817 // Filtre par face
1818 DoubleTrav flag_cl(coeff_node);
1819 filtre_CLD(flag_cl);
1820 // int nb_face=domaine_VEF.domaine().nb_faces();
1821 // for(int face=0; face<nb_face; face++)Cerr<<"face flag_c = "<<face<<" "<<flag_cl(face,0)<<" "<<flag_cl(face,1)<<" "<<flag_cl(face,2)<<finl;
1822
1823 for(int num_elem=0; num_elem<nb_elems; num_elem++)
1824 {
1825 if (aire(num_elem)>0.)
1826 {
1827 correction_en_pression(num_elem) = 0.;
1828 }
1829 else
1830 {
1831 coef_elem(num_elem) = 0. ;
1832 for (int i=0; i<nb_face_elem; i++)
1833 {
1834 int sii=elems(num_elem,i);
1835 for (int comp=0; comp<ncomp; comp++)
1836 {
1837 if ((coeff_node(sii,comp)!= 1.0) || (flag_cl(sii,comp) != 1.0)) coef_elem(num_elem) += 1. ;
1838 }
1839 }
1840 if (coef_elem(num_elem) == (ncomp*nb_face_elem)) correction_en_pression(num_elem) = 0.;
1841 }
1842 }
1843}
1844
1845void Source_PDF_VEF::correct_pressure(const DoubleTab& coeff_node, DoubleTab& pression, const DoubleTab& correction_en_pression) const
1846{
1847 const DoubleTab& aire=champ_aire_->valeurs();
1848 //int nb = pression.size();
1849 //Cerr << pression << finl;
1850 //Cerr << nb << finl;
1851 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
1852 const IntTab& elems= domaine_VEF.elem_faces();
1853 int nb_face_elem=domaine_VEF.domaine().nb_faces_elem();
1854 int nb_elem= domaine_VEF.domaine().nb_elem_tot();
1855 int ncomp = coeff_node.dimension(1) ;
1856 DoubleTrav coef_elem(nb_elem);
1857
1858 DoubleTrav flag_cl(coeff_node);
1859 filtre_CLD(flag_cl);
1860
1861 for(int num_elem=0; num_elem<nb_elem; num_elem++)
1862 {
1863 if (aire(num_elem)<=0.)
1864 {
1865 coef_elem(num_elem) = 0. ;
1866 for (int i=0; i<nb_face_elem; i++)
1867 {
1868 int sii=elems(num_elem,i);
1869 for (int comp=0; comp<ncomp; comp++)
1870 {
1871 if ((coeff_node(sii,comp)!= 1.0) || (flag_cl(sii,comp) != 1.0)) coef_elem(num_elem) += 1. ;
1872 }
1873 }
1874 if (coef_elem(num_elem) != (ncomp*nb_face_elem)) pression(num_elem)+=correction_en_pression(num_elem);
1875 }
1876 }
1877}
1878
1879void Source_PDF_VEF::filtre_CLD(DoubleTab& flag_cl) const
1880{
1881 // Filtre et dof par face
1882 const Domaine_Cl_VEF& le_dom_cl = le_dom_Cl_VEF.valeur();
1883 int nb_cond_lim = le_dom_cl.nb_cond_lim();
1884 int nb_comp = flag_cl.dimension(1);
1885 flag_cl = 1.;
1886 for (int ij=0; ij<nb_cond_lim; ij++)
1887 {
1888 const Cond_lim_base& la_cl_base = le_dom_cl.les_conditions_limites(ij).valeur();
1889 const Front_VF& le_bord = ref_cast(Front_VF,le_dom_cl.les_conditions_limites(ij)->frontiere_dis());
1890 int nfaces = le_bord.nb_faces_tot();
1891 if (sub_type(Dirichlet,la_cl_base))
1892 {
1893 for (int ind_face=0; ind_face < nfaces; ind_face++)
1894 {
1895 int face=le_bord.num_face(ind_face);
1896 for (int comp=0; comp<nb_comp; comp++) flag_cl(face,comp) = 0. ;
1897 }
1898 }
1899 }
1900 //vector; NS Attention seulement defini en EF
1901 // if (nb_comp == Objet_U::dimension) le_dom_cl.imposer_symetrie(flag_cl,0);
1902 flag_cl.echange_espace_virtuel();
1903 return;
1904}
1905
1906/*##################################################################################################
1907####################################################################################################
1908################################# Impression #######################################################
1909####################################################################################################
1910##################################################################################################*/
1911
1912int Source_PDF_VEF::impr(Sortie& os) const
1913{
1914 if (out_=="??")
1915 {
1916 Cerr << __FILE__ << (int)__LINE__ << " ERROR / Source_PDF_VEF::impr" << finl;
1917 Cerr << "No balance printed for " << que_suis_je() << finl;
1918 Cerr << "Because output file name is not specified." << finl;
1919 return 0;
1920 }
1921 else
1922 {
1923 int nb_compo=bilan_.size();
1924 if (nb_compo==0)
1925 {
1926 Cerr << __FILE__ << (int)__LINE__ << " ERROR / Source_PDF_VEF::impr" << finl;
1927 Cerr << "No balance printed for " << que_suis_je() << finl;
1928 Cerr << "Because bilan_ array is not filled." << finl;
1929 return 0;
1930 }
1931 else
1932 {
1933 const Schema_Temps_base& sch=equation().probleme().schema_temps();
1934 double temps=sch.temps_courant();
1935 double pdtps = sch.pas_de_temps();
1936 if (temps == pdtps)
1937 {
1938 int flag = Process::je_suis_maitre();
1939 if(flag) ouvrir_fichier(Flux,"",flag);
1940 return 0;
1941 }
1942 const DoubleTab& variable=equation().inconnue().valeurs();
1943 Nom espace=" \t";
1944
1945 int nb_comp = variable.dimension(1);
1946 assert(nb_comp == modele_lu_.dim_variable_);
1947 int pdf_dt_conv = modele_lu_.pdf_bilan();
1948
1949 DoubleTrav secmem_conv(variable);
1950 // defaut pdf_dt_conv = 0
1951 secmem_conv = 0.;
1952 int i_traitement_special = 0;
1953
1954 if (pdf_dt_conv == 1 || pdf_dt_conv == 2)
1955 // pdf_dt_conv = 1 ou 2
1956 {
1957 if (nb_comp == Objet_U::dimension)
1958 {
1959 //vector; NS
1960 Navier_Stokes_std& eq_NS = ref_cast_non_const(Navier_Stokes_std, equation());
1961 int transport_rhou=0;
1962 if (eq_NS.nombre_d_operateurs() > 1)
1963 {
1964 transport_rhou= (eq_NS.vitesse_pour_transport().le_nom()=="rho_u"?1:0);
1965 if (pdf_dt_conv == 2) eq_NS.derivee_en_temps_conv(secmem_conv , variable);
1966 }
1967 i_traitement_special = (transport_rhou==1?2:102);
1968 }
1969 else if (nb_comp == 1)
1970 {
1971 //scalar; temps en rho
1972 Equation_base& eq_Ba = ref_cast_non_const(Equation_base, equation());
1973 if (equation().nombre_d_operateurs() > 1)
1974 {
1975 if (pdf_dt_conv == 2) eq_Ba.derivee_en_temps_conv(secmem_conv , variable);
1976 }
1977 i_traitement_special = 2;
1978 }
1979 else
1980 {
1981 Cerr << "Source_PDF_VEF: for scalar or vector only; dim = " <<nb_comp<< finl;
1982 Process::exit();
1983 }
1984 }
1985 else if(pdf_dt_conv != 0 )
1986 {
1987 Cerr<<"Source_PDF_VEF: Modele pdf_bilan must be 0; 1 or 2 only"<<finl;
1988 Process::exit();
1989 }
1990
1991 Source_PDF_base& my_src = ref_cast_non_const(Source_PDF_base,*this);
1992 my_src.compute_source_term_PDF(i_traitement_special, secmem_conv, bilan_);
1993
1994 int flag = Process::je_suis_maitre();
1995 if(flag)
1996 {
1997 ouvrir_fichier(Flux,"",flag);
1998 if(nb_comp == Objet_U::dimension)
1999 {
2000 // vector
2001 assert(Objet_U::dimension==3);
2002 Flux << temps << espace << bilan_(0) << espace << bilan_(1) << espace << bilan_(2) << finl;
2003 }
2004 else
2005 {
2006 // scalar
2007 Flux << temps << espace << bilan_(0) << finl;
2008 }
2009 }
2010 }
2011 return 1;
2012 }
2013}
2014// void imprimer_ustar_yplus__mean_only(Sortie&, const Nom& ) const override;
2015
2016
2017/*##################################################################################################
2018####################################################################################################
2019################################# POSTRAITEMENT ####################################################
2020####################################################################################################
2021##################################################################################################*/
2022
2023void Source_PDF_VEF::creer_champ(const Motcle& motlu)
2024{
2025 Source_PDF_base::creer_champ(motlu);
2026
2027 if (motlu=="u_star_ibm" && !champ_u_star_ibm_ && imm_wall_law_)
2028 {
2029 int nb_comp = 1;
2030 Noms noms(1);
2031 noms[0]="u_star_ibm";
2032 Noms unites(1);
2033 unites[0] = "m/s";
2034 double temps=0.;
2035 const Discretisation_base& discr = equation().probleme().discretisation();
2036 discr.discretiser_champ("champ_sommets",equation().domaine_dis(),scalaire,noms,unites,nb_comp,temps,champ_u_star_ibm_);
2037 champs_compris_.ajoute_champ(champ_u_star_ibm_);
2038 }
2039 else if (motlu=="y_plus_ibm" && !champ_y_plus_ibm_ && imm_wall_law_)
2040 {
2041 int nb_comp = 1;
2042 Noms noms(1);
2043 noms[0]="y_plus_ibm";
2044 Noms unites(1);
2045 unites[0] = "-";
2046 double temps=0.;
2047 const Discretisation_base& discr = equation().probleme().discretisation();
2048 discr.discretiser_champ("champ_sommets",equation().domaine_dis(),scalaire,noms,unites,nb_comp,temps,champ_y_plus_ibm_);
2049 champs_compris_.ajoute_champ(champ_y_plus_ibm_);
2050 }
2051}
2052
2053bool Source_PDF_VEF::has_champ(const Motcle& nom, OBS_PTR(Champ_base) &ref_champ) const
2054{
2055 if (Source_PDF_base::has_champ(nom)) return Source_PDF_base::has_champ(nom, ref_champ);
2056
2057 if (nom == "u_star_ibm" && champ_u_star_ibm_)
2058 {
2059 ref_champ = Source_PDF_VEF::get_champ(nom);
2060 return true;
2061 }
2062 else if (nom == "y_plus_ibm" && champ_y_plus_ibm_)
2063 {
2064 ref_champ = Source_PDF_VEF::get_champ(nom);
2065 return true;
2066 }
2067 else
2068 return false; /* rien trouve */
2069}
2070
2071bool Source_PDF_VEF::has_champ(const Motcle& nom) const
2072{
2073 if (Source_PDF_base::has_champ(nom)) return true;
2074
2075 if (nom == "u_star_ibm" && champ_u_star_ibm_)
2076 return true;
2077 else if (nom == "y_plus_ibm" && champ_y_plus_ibm_)
2078 return true;
2079 else
2080 return false; /* rien trouve */
2081}
2082
2083const Champ_base& Source_PDF_VEF::get_champ(const Motcle& nom) const
2084{
2085 if (Source_PDF_base::has_champ(nom)) return Source_PDF_base::get_champ(nom);
2086 if (nom=="u_star_ibm")
2087 {
2088 if (!champ_u_star_ibm_)
2089 throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
2090 // Initialisation a 0 du champ volumique u_star
2091 DoubleTab& valeurs = champ_u_star_ibm_->valeurs();
2092 valeurs=0;
2093 if (tab_u_star_ibm_.size_array()>0)
2094 {
2095 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
2096 int nb_som=domaine_VEF.nb_som();
2097 for (int num_node=0; num_node<nb_som; num_node++)
2098 valeurs(num_node)=tab_u_star_ibm_(num_node);
2099 }
2100 valeurs.echange_espace_virtuel();
2101 champ_u_star_ibm_->mettre_a_jour(equation().probleme().schema_temps().temps_courant());
2102 return champs_compris_.get_champ(nom);
2103 }
2104 else if (nom=="y_plus_ibm")
2105 {
2106 if (!champ_y_plus_ibm_)
2107 throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
2108
2109 // Initialisation a 0 du champ volumique u_star
2110 DoubleTab& valeurs = champ_y_plus_ibm_->valeurs();
2111 valeurs=0;
2112 if (tab_y_plus_ibm_.size_array()>0)
2113 {
2114 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
2115 int nb_som=domaine_VEF.nb_som();
2116 for (int num_node=0; num_node<nb_som; num_node++)
2117 valeurs(num_node)=tab_y_plus_ibm_(num_node);
2118 }
2119 valeurs.echange_espace_virtuel();
2120 champ_y_plus_ibm_->mettre_a_jour(equation().probleme().schema_temps().temps_courant());
2121 return champs_compris_.get_champ(nom);
2122 }
2123 else
2124 return champs_compris_.get_champ(nom);
2125}
2126
2127void Source_PDF_VEF::get_noms_champs_postraitables(Noms& nom,Option opt) const
2128{
2129 Source_PDF_base::get_noms_champs_postraitables(nom,opt);
2130
2131 Noms noms_compris = champs_compris_.liste_noms_compris();
2132 if(imm_wall_law_)
2133 {
2134 noms_compris.add("u_star_ibm");
2135 noms_compris.add("y_plus_ibm");
2136 }
2137 if (opt==DESCRIPTION)
2138 Cerr<<" Source_PDF_VEF : "<< noms_compris <<finl;
2139 else
2140 nom.add(noms_compris);
2141}
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_P1NC
Definition Champ_P1NC.h:35
Champ_P1NC::valeur_a_elem
DoubleVect & valeur_a_elem(const DoubleVect &position, DoubleVect &val, int le_poly) const override
provoque une erreur ! doit etre surchargee par les classes derivees
Definition Champ_P1NC.h:67
Champ_Proto::valeurs
virtual DoubleTab & valeurs()=0
Champ_base
classe Champ_base Cette classe est la base de la hierarchie des champs.
Definition Champ_base.h:43
Cond_lim_base
classe Cond_lim_base Classe de base pour la hierarchie des classes qui representent les differentes c...
Definition Cond_lim_base.h:40
Dirichlet
classe Dirichlet Cette classe est la classe de base de la hierarchie des conditions aux limites de ty...
Definition Dirichlet.h:31
Discretisation_base
classe Discretisation_base Cette classe represente un schema de discretisation en espace,...
Definition Discretisation_base.h:45
Discretisation_base::discretiser_champ
void discretiser_champ(const Motcle &directive, const Domaine_dis_base &z, const Nom &nom, const Nom &unite, int nb_comp, int nb_pas_dt, double temps, OWN_PTR(Champ_Inc_base)&champ, const Nom &sous_type=NOM_VIDE) const
Definition Discretisation_base.cpp:59
Domaine_32_64::nb_elem_tot
int_t nb_elem_tot() const
Definition Domaine.h:132
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Renvoie le nombre de face de type i des elements geometriques constituants le domaine.
Definition Domaine.h:484
Domaine_32_64::coord
double coord(int_t i, int j) const
Definition Domaine.h:110
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_dis_base
classe Domaine_Cl_dis_base Les objets Domaine_Cl_dis_base representent les conditions aux limites
Definition Domaine_Cl_dis_base.h:37
Domaine_Cl_dis_base::nb_cond_lim
int nb_cond_lim() const
Renvoie le nombre de conditions aux limites.
Definition Domaine_Cl_dis_base.cpp:425
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::volumes_entrelaces
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
Domaine_VF::face_sommets
int face_sommets(int i, int j) const
renvoie le numero du ieme sommet de la face num_face.
Definition Domaine_VF.h:583
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
Domaine_dis_base
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
Definition Domaine_dis_base.h:39
Domaine_dis_base::nb_som
int nb_som() const
Definition Domaine_dis_base.h:51
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Equation_base
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
Definition Equation_base.h:76
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Equation_base::derivee_en_temps_conv
DoubleTab & derivee_en_temps_conv(DoubleTab &, const DoubleTab &)
Add convection term In: solution is the unknown I.
Definition Equation_base.cpp:1481
Equation_base::inconnue
virtual const Champ_Inc_base & inconnue() const =0
Equation_base::nombre_d_operateurs
virtual int nombre_d_operateurs() const =0
Equation_base::probleme
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
Definition Equation_base.cpp:747
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces_tot
int nb_faces_tot() const
Definition Front_VF.h:58
Front_VF::num_face
int num_face(const int) const
Definition Front_VF.h:68
Interpolation_IBM_base::get_N_histo
int get_N_histo()
Definition Interpolation_IBM_base.h:37
Interpolation_IBM_base::get_impr
bool get_impr()
Definition Interpolation_IBM_base.h:36
Interpolation_IBM_elem_fluid
Definition Interpolation_IBM_elem_fluid.h:22
Interpolation_IBM_hybrid
Definition Interpolation_IBM_hybrid.h:23
Interpolation_IBM_hybrid::getSommetsVoisinsOf
IntList & getSommetsVoisinsOf(int i)
Definition Interpolation_IBM_hybrid.h:27
Interpolation_IBM_mean_gradient
Definition Interpolation_IBM_mean_gradient.h:24
Interpolation_IBM_mean_gradient::getSommetsVoisinsOf
IntList & getSommetsVoisinsOf(int i)
Definition Interpolation_IBM_mean_gradient.h:30
Interpolation_IBM_power_law_tbl_proto::get_C_pwl_WJSP
double get_C_pwl_WJSP()
Definition Interpolation_IBM_power_law_tbl_proto.h:39
Interpolation_IBM_power_law_tbl_proto::get_y_c1_p_pwl_WJSP
double get_y_c1_p_pwl_WJSP()
Definition Interpolation_IBM_power_law_tbl_proto.h:49
Interpolation_IBM_power_law_tbl_proto::get_A_pwl
double get_A_pwl(int a)
Definition Interpolation_IBM_power_law_tbl_proto.h:31
Interpolation_IBM_power_law_tbl_proto::get_D_pwl_WJSP
double get_D_pwl_WJSP()
Definition Interpolation_IBM_power_law_tbl_proto.h:40
Interpolation_IBM_power_law_tbl_proto::get_y_c2_p_pwl_WJSP
double get_y_c2_p_pwl_WJSP()
Definition Interpolation_IBM_power_law_tbl_proto.h:50
Interpolation_IBM_power_law_tbl_proto::get_y_c_p_pwl
double get_y_c_p_pwl()
Definition Interpolation_IBM_power_law_tbl_proto.h:48
Interpolation_IBM_power_law_tbl_proto::get_p_pwl_WJSP
double get_p_pwl_WJSP()
Definition Interpolation_IBM_power_law_tbl_proto.h:45
Interpolation_IBM_power_law_tbl_proto::get_B_pwl
double get_B_pwl()
Definition Interpolation_IBM_power_law_tbl_proto.h:44
Interpolation_IBM_power_law_tbl_u_star
Definition Interpolation_IBM_power_law_tbl_u_star.h:24
Interpolation_IBM_power_law_tbl_u_star::getSommetsVoisinsOf
IntList & getSommetsVoisinsOf(int i)
Definition Interpolation_IBM_power_law_tbl_u_star.h:31
Interpolation_IBM_power_law_tbl
Definition Interpolation_IBM_power_law_tbl.h:23
Interpolation_IBM_power_law_tbl::get_formulation_linear_pwl
int get_formulation_linear_pwl()
Definition Interpolation_IBM_power_law_tbl.h:26
Interpolation_IBM_power_law_tbl::get_formulation_WJSP
int get_formulation_WJSP()
Definition Interpolation_IBM_power_law_tbl.h:27
Interpolation_IBM_thermal_wall_law
Definition Interpolation_IBM_thermal_wall_law.h:22
Interpolation_IBM_thermal_wall_law::get_boundary_type
int get_boundary_type()
Definition Interpolation_IBM_thermal_wall_law.h:28
Interpolation_IBM_thermal_wall_law::get_T_inlet
double get_T_inlet()
Definition Interpolation_IBM_thermal_wall_law.h:29
Interpolation_IBM_thermal_wall_law::get_Prandlt_mol
double get_Prandlt_mol()
Definition Interpolation_IBM_thermal_wall_law.h:30
Interpolation_IBM_thermal_wall_law::Kader
double Kader(double, double)
Definition Interpolation_IBM_thermal_wall_law.cpp:48
Interpolation_IBM_thermal_wall_law::get_formulation_Tp
int get_formulation_Tp()
Definition Interpolation_IBM_thermal_wall_law.h:27
Matrice_Morse
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
Definition Matrice_Morse.h:50
Matrice_Morse::ajouter_multvect_
DoubleVect & ajouter_multvect_(const DoubleVect &, DoubleVect &) const override
Operation de multiplication-accumulation (saxpy) matrice vecteur.
Definition Matrice_Morse.cpp:877
Matrice_Morse::coef
double coef(int i, int j) const
Definition Matrice_Morse.h:129
Milieu_base::capacite_calorifique
virtual const Champ_Don_base & capacite_calorifique() const
Renvoie la capacite calorifique du milieu.
Definition Milieu_base.cpp:867
Milieu_base::masse_volumique
virtual const Champ_base & masse_volumique() const
Renvoie la masse volumique du milieu.
Definition Milieu_base.cpp:797
MorEqn::equation
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
Motcle
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
Navier_Stokes_std
classe Navier_Stokes_std Cette classe porte les termes de l'equation de la dynamique
Definition Navier_Stokes_std.h:56
Navier_Stokes_std::vitesse_pour_transport
virtual const Champ_base & vitesse_pour_transport() const
Definition Navier_Stokes_std.cpp:274
Navier_Stokes_std::nombre_d_operateurs
int nombre_d_operateurs() const override
Renvoie le nombre d'operateurs de l'equation: Pour Navier Stokes Standard c'est 2.
Definition Navier_Stokes_std.cpp:444
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
Nom::getString
const std::string & getString() const
Definition Nom.h:92
Noms
Un tableau de chaine de caracteres (VECT(Nom)).
Definition Noms.h:26
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::le_nom
virtual const Nom & le_nom() const
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
Definition Objet_U.cpp:319
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Op_Diff_VEF_base
class Op_Diff_VEF_base
Definition Op_Diff_VEF_base.h:39
Operateur_Diff_base::diffusivite
virtual const Champ_base & diffusivite() const =0
Probleme_base
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
Definition Probleme_base.h:55
Probleme_base::get_champ
const Champ_base & get_champ(const Motcle &nom) const override
Definition Probleme_base.cpp:841
Probleme_base::discretisation
const Discretisation_base & discretisation() const
Renvoie la discretisation associee au probleme.
Definition Probleme_base.h:241
Probleme_base::schema_temps
const Schema_Temps_base & schema_temps() const
Renvoie le schema en temps associe au probleme.
Definition Probleme_base.cpp:578
Process::abort
static void abort()
Routine de sortie de Trio-U sur une erreur abort().
Definition Process.cpp:570
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Process::je_suis_maitre
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Definition Process.cpp:86
Schema_Temps_base
class Schema_Temps_base
Definition Schema_Temps_base.h:67
Schema_Temps_base::temps_courant
double temps_courant() const
Renvoie le temps courant.
Definition Schema_Temps_base.h:445
Schema_Temps_base::pas_de_temps
double pas_de_temps() const
Renvoie le pas de temps (delta_t) courant.
Definition Schema_Temps_base.h:393
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
Source_PDF_VEF
Definition Source_PDF_VEF.h:31
Source_PDF_VEF::compute_indicateur_nodal_champ_aire
void compute_indicateur_nodal_champ_aire() override
Definition Source_PDF_VEF.cpp:75
Source_PDF_VEF::tab_u_star_ibm_
DoubleVect tab_u_star_ibm_
valeurs des u* IBM calculees localement
Definition Source_PDF_VEF.h:74
Source_PDF_VEF::calculer_vitesse_imposee_power_law_tbl
void calculer_vitesse_imposee_power_law_tbl() override
Definition Source_PDF_VEF.cpp:682
Source_PDF_VEF::associer_pb
void associer_pb(const Probleme_base &) override
Definition Source_PDF_VEF.cpp:66
Source_PDF_VEF::get_noms_champs_postraitables
void get_noms_champs_postraitables(Noms &nom, Option opt=NONE) const override
Definition Source_PDF_VEF.cpp:2127
Source_PDF_VEF::correct_pressure
void correct_pressure(const DoubleTab &, DoubleTab &, const DoubleTab &) const override
Definition Source_PDF_VEF.cpp:1845
Source_PDF_VEF::calculer_temperature_imposee_wall_law
void calculer_temperature_imposee_wall_law() override
Definition Source_PDF_VEF.cpp:1486
Source_PDF_VEF::has_champ
bool has_champ(const Motcle &nom, OBS_PTR(Champ_base) &ref_champ) const override
Definition Source_PDF_VEF.cpp:2053
Source_PDF_VEF::calculer_vitesse_imposee_power_law_tbl_u_star
void calculer_vitesse_imposee_power_law_tbl_u_star() override
Definition Source_PDF_VEF.cpp:1145
Source_PDF_VEF::contribuer_a_avec
void contribuer_a_avec(const DoubleTab &, Matrice_Morse &) const override
contribution a la matrice implicite des termes sources par defaut pas de contribution
Definition Source_PDF_VEF.cpp:225
Source_PDF_VEF::verif_ajouter_contrib
void verif_ajouter_contrib(const DoubleTab &variable, Matrice_Morse &matrice) const
Definition Source_PDF_VEF.cpp:297
Source_PDF_VEF::calculer_variable_imposee_hybrid
void calculer_variable_imposee_hybrid() override
Definition Source_PDF_VEF.cpp:537
Source_PDF_VEF::get_champ
const Champ_base & get_champ(const Motcle &nom) const override
Definition Source_PDF_VEF.cpp:2083
Source_PDF_VEF::filtre_CLD
void filtre_CLD(DoubleTab &) const override
Definition Source_PDF_VEF.cpp:1879
Source_PDF_VEF::calculer_variable_imposee_elem_fluid
void calculer_variable_imposee_elem_fluid() override
Definition Source_PDF_VEF.cpp:342
Source_PDF_VEF::creer_champ
void creer_champ(const Motcle &motlu) override
Definition Source_PDF_VEF.cpp:2023
Source_PDF_VEF::multiply_coeff_volume
void multiply_coeff_volume(DoubleTab &) const override
Definition Source_PDF_VEF.cpp:87
Source_PDF_VEF::correct_incr_pressure
void correct_incr_pressure(const DoubleTab &, DoubleTab &) const override
Definition Source_PDF_VEF.cpp:1806
Source_PDF_VEF::impr
int impr(Sortie &) const override
Definition Source_PDF_VEF.cpp:1912
Source_PDF_VEF::ajouter_
DoubleTab & ajouter_(const DoubleTab &, DoubleTab &) const override
Definition Source_PDF_VEF.cpp:208
Source_PDF_VEF::associer_domaines
void associer_domaines(const Domaine_dis_base &, const Domaine_Cl_dis_base &) override
Definition Source_PDF_VEF.cpp:80
Source_PDF_VEF::OBS_PTR
OBS_PTR(Domaine_VEF) le_dom_VEF
Source_PDF_VEF::calculer_variable_imposee_mean_grad
void calculer_variable_imposee_mean_grad() override
Definition Source_PDF_VEF.cpp:430
Source_PDF_VEF::tab_y_plus_ibm_
DoubleVect tab_y_plus_ibm_
valeurs des d+ IBM calculees localement
Definition Source_PDF_VEF.h:75
Source_PDF_base
class Source_PDF_base Base class for the source terms for the penalisation of the momentum in the Imm...
Definition Source_PDF_base.h:37
Source_PDF_base::champ_aire_
champ_aire_
Definition Source_PDF_base.h:117
Source_PDF_base::compute_pond
DoubleTab compute_pond(const DoubleTab &rho_m, const DoubleTab &aire, const DoubleVect &volume_thilde, int &nb_som_elem, int &nb_elems) const
Definition Source_PDF_base.cpp:1546
Source_PDF_base::creer_champ
void creer_champ(const Motcle &motlu) override
Definition Source_PDF_base.cpp:1852
Source_PDF_base::get_noms_champs_postraitables
void get_noms_champs_postraitables(Noms &nom, Option opt=NONE) const override
Definition Source_PDF_base.cpp:1935
Source_PDF_base::compute_source_term_PDF
DoubleVect & compute_source_term_PDF(int, DoubleTab &, DoubleVect &)
Definition Source_PDF_base.cpp:1610
Source_PDF_base::compute_indicateur_nodal_champ_aire
virtual void compute_indicateur_nodal_champ_aire()
Definition Source_PDF_base.cpp:1245
Source_PDF_base::ouvrir_fichier
void ouvrir_fichier(SFichier &, const Nom &, const int) const override
Ouverture/creation d'un fichier d'impression d'un terme source A surcharger dans les classes derivees...
Definition Source_PDF_base.cpp:1785
Source_PDF_base::champ_rho_
champ_rho_
Definition Source_PDF_base.h:117
Source_PDF_base::indicateur_nodal_champ_aire_
DoubleTab indicateur_nodal_champ_aire_
Definition Source_PDF_base.h:123
Source_PDF_base::imm_wall_law_
bool imm_wall_law_
Definition Source_PDF_base.h:147
Source_PDF_base::modele_lu_
PDF_model modele_lu_
Definition Source_PDF_base.h:126
Source_PDF_base::champ_rotation_
champ_rotation_
Definition Source_PDF_base.h:117
Source_PDF_base::get_champ
const Champ_base & get_champ(const Motcle &) const override
Definition Source_PDF_base.cpp:2004
Source_PDF_base::associer_pb
void associer_pb(const Probleme_base &) override
Definition Source_PDF_base.cpp:98
Source_PDF_base::has_champ
bool has_champ(const Motcle &nom, OBS_PTR(Champ_base) &ref_champ) const override
Definition Source_PDF_base.cpp:1950
Source_PDF_base::variable_imposee_
DoubleTab variable_imposee_
Definition Source_PDF_base.h:125
Source_base::champs_compris_
Champs_compris champs_compris_
Definition Source_base.h:105
Source_base::out_
Nom out_
Definition Source_base.h:100
Source_base::Flux
SFichier Flux
Definition Source_base.h:104
Source_base::bilan_
DoubleVect bilan_
Definition Source_base.h:103
TRUSTListElem::est_vide
int est_vide() const
Definition TRUSTListElem.h:40
TRUSTList::size
int size() const
Definition TRUSTList.h:68
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
TRUSTVect::size_totale
_SIZE_ size_totale() const
Definition TRUSTVect.tpp:61
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282