TrioCFD 1.9.9_beta
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<<" different dimensions 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 // compute: coefficient rho/dt * coeff_relax * (volume_tilde / 8) * variable
117 /* i_traitement_special = 0 => standard treatment; 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 << "error here" << finl;
129 if (dim_var != 1 && dim_var != dim_esp)
130 {
131 Cerr<<"ajouter_: variable dimension differs from 1 or "<<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) //time term in rho v
155 {
156 for (int k=0; k<dim_var; k++) tuvw[k] = 1.0;
157 }
158 else if (i_traitement_special == 101) //time term in 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 // generic case: 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_ node: "<<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 // compute: 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: variable dimension differs from 1 or "<<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 for intercepted elements
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 squared norm of absolute force diff. = "<<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 << "entering here" << 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 << "entering here" << finl;
566 for (int i = 0; i < nb_som; i++)
567 {
568 Cerr << "passing here" << 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 // operator(0) : diffusivity
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 // handling of exceptions
806 if (d1 < eps) itisok = 0; //point P is on the immersed boundary: assignment
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; // set unit norm while keeping direction, i.e. normalize
810 else
811 {
812 for(int j = 0; j < nb_comp; j++) normale(0,j) = 0.; // the fluid point is on the immersed boundary: assignment
813 itisok = 0;
814 }
815 if ( y_ref < eps ) // the fluid point is on the immersed boundary: assignment
816 {
817 y_ref = 1.;
818 itisok = 0;
819 }
820 // compute power law velocity
821 if (itisok)
822 {
823 cells(0) = int(fluid_elems(i));
824 champ_vitesse_inconnue.valeur_a_elem(xf, vf, cells[0]); // vf: total velocity interpolated at the fluid point
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 // Compute U_tau_ref to evaluate dimensionless quantities (y+ ...)
838 // Assumption: power-law sublayer
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; // now we have everything needed to establish the polynomial wall law for y r+
842 double test_ref;
843
844 if (!form_WJSP)
845 {
846 if ( y_ref_p > y_c_p_pwl) // from here begins the expression of the turbulent polynomial wall law
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 // Linear formulation of the polynomial wall law ------------------------------
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 << "log zone/inertial sublayer consistent with the u tau computation" << finl;
858 }
859 else
860 {
861 // Actually, we are in the viscous sublayer
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 ) ; // recompute u_tau and y_plus linearly at the reference point
872 y_ref_p = y_ref * u_tau_ref / nu;
873 if ( y_ref_p > y_c_p_pwl )
874 {
875 // Inconsistency: wall laws are not applied
876 itisok = 0;
877 }
878 // Cerr << "linear zone/viscous sublayer" << finl;
879 } // End 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)) ) ; // Assumption: power-law sublayer
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 ) ; // recompute u_tau and y_plus linearly at the forcing point
897 y_plus = u_tau * d1 / nu;
898 if ( y_plus > y_c_p_pwl )
899 {
900 // Inconsistency: wall laws are not applied
901 Cerr<<"INCONSISTENCY: 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 // here we test whether the boundary node and the fluid point are in the same zone: if yes ok, otherwise impose velocity
908
909 // Handling of exceptions
910 if ( (test_ * test_ref < 0) || (itisok == 0) ) //otherwise impose wall velocity
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 of y+ near the walls
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) // from here starts the expression of the turbulent polynomial wall law
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 // Linear formulation of the polynomial wall law ------------------------------
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 // Linear formulation of the polynomial wall law ------------------------------
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 // Linear formulation of the polynomial wall law ------------------------------
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 << "log zone/inertial sublayer consistent with the u tau computation" << finl;
987 }
988 else // buffer layer hypothesis
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 // Linear formulation of the polynomial wall law ------------------------------
1004 {
1005 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = 0 ; //TODO
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 // Linear formulation of the polynomial wall law ------------------------------
1015 {
1016 for(int j = 0; j < nb_comp; j++) vitesse_imposee_sommet(i,j) = 0 ;
1017 }
1018 }
1019 }
1020 else
1021 {
1022 // In fact, we are in the viscous sublayer
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 ) ; // recompute u_tau and y_plus linearly at the reference point
1032 y_ref_p = y_ref * u_tau_ref / nu;
1033 if ( y_ref_p > y_c1_p_pwl_WJSP)
1034 {
1035 // Inconsistency: wall laws are not used
1036 itisok = 0;
1037 Cerr << "Inconsistency" << finl;
1038 }
1039 }
1040 }
1041 // Cerr<<"u_tau d1 nu ="<<u_tau<<" "<<d1<<" "<<nu<<finl;;
1042
1043 // here we test whether the boundary node and the fluid point are in the same zone: if yes ok, otherwise impose velocity
1044
1045 // Handling of exceptions
1046 if (itisok == 0) //otherwise impose wall velocity
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 of y+ near the walls
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<<"y+ histogram count =";
1134 for (int j = 0; j < N; j++) Cerr<<" "<<compteur_h(0,j)<<" ";
1135 Cerr<<finl;
1136 Cerr<<"y+ histogram abscissa =";
1137 for (int j = 0; j < N+1; j++) Cerr<<" "<<abs_h(0,j)<<" ";
1138 Cerr<<finl;
1139 Cerr<<"y+ histogram: min = "<<h_yplus_min<<" - mean = "<<h_yplus_mean<<" - max = "<<h_yplus_max<<finl;
1140
1141 Cerr<<"Average over "<<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 // operator(0) : diffusivity
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 and normal distance to 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)); //to define the laminar viscosity
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 //point P is not on the immersed boundary
1257 {
1258 IntList& voisins = interp.getSommetsVoisinsOf(i);
1259 if (!(voisins.est_vide()))
1260 // there is a list of neighbouring dofs
1261 {
1262 double u_tau_ref;
1263 double y_plus_ref;
1264 DoubleTab mean_u_tau_neighbour(1, nb_comp); // mean vector u_tau over the neighbouring dofs
1265 mean_u_tau_neighbour = 0.0;
1266 double mean_y_plus_neighbour = 0.; // mean y_plus over the neighbouring dofs
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 // loop over the neighbouring dofs
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 // normal and tangential velocity components with respect to 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)) ) ; // power-law friction velocity at fluid point k
1313 y_plus_ref = d1 * u_tau_ref / nu;
1314 if (y_plus_ref < y_c_p_pwl ) // actually we are in the linear law region!
1315 {
1316 u_tau_ref = pow ( (nu * norme_vtf_k / d1) , 0.5 ) ; // recompute u_tau_ref and y_plus_ref linearly at fluid point k
1317 y_plus_ref = d1 * u_tau_ref / nu;
1318 if ( y_plus_ref > y_c_p_pwl )
1319 {
1320 // Inconsistency: wall laws are not applied
1321 itisok = 0;
1322 }
1323 }
1324
1325 // Average the tangential friction velocity vectors (with respect to P) weighted by the inverse of 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 // Compute mean u_tau, y_plus and tangential velocity at the forcing point P
1349 if (pond_tot > 0.)
1350 {
1351 mean_u_tau_neighbour /= pond_tot; // vector mean u_tau
1352 mean_y_plus_neighbour /= pond_tot; // mean y_plus
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 = norm of the mean vector u_tau
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. ;// tangent at P
1360 if (y_plus > y_c_p_pwl ) // polynomial wall law
1361 {
1362 if (mean_y_plus_neighbour > y_c_p_pwl ) // consistent with polynomial wall law
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; // inconsistent with polynomial wall law
1368 }
1369 else // linear law
1370 {
1371 if (mean_y_plus_neighbour < y_c_p_pwl ) // consistent with linear law
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; // inconsistent with linear law
1376 }
1377 }
1378 // there is no contribution from the neighbouring dofs
1379 else itisok_P = 0;
1380 }
1381 else
1382 // there is no list of neighbouring dofs
1383 itisok_P = 0;
1384 }
1385 //point P is on the immersed boundary: assignment
1386 else itisok_P = 0;
1387
1388 // for post-processing
1389 if(itisok_P)
1390 {
1391 // for post-processing & required for 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 of y+ near the walls
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<<"y+ histogram count =";
1466 for (int j = 0; j < N; j++) Cerr<<" "<<compteur_h(0,j)<<" ";
1467 Cerr<<finl;
1468 Cerr<<"y+ histogram abscissa =";
1469 for (int j = 0; j < N+1; j++) Cerr<<" "<<abs_h(0,j)<<" ";
1470 Cerr<<finl;
1471 Cerr<<"y+ histogram: min = "<<h_yplus_min<<" - mean = "<<h_yplus_mean<<" - max = "<<h_yplus_max<<finl;
1472
1473 Cerr<<"Average over "<<yplus_count<<" forced points and "<<yplus_ref_count<<" reference points"<<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// Domain variables
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// Arrays of imposed and computed values
1497 DoubleTab& variable_imposee_mod = modele_lu_.variable_imposee_->valeurs();
1498 DoubleTab& temperature_imposee_calculee = variable_imposee_;
1499// Information arrays for the immersed wall
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// Discretisation of the computed quantity
1505 Champ_P1NC& champ_variable_inconnue = ref_cast(Champ_P1NC,equation().inconnue());
1506 DoubleTab& val_variable_inconnue = champ_variable_inconnue.valeurs();
1507// Retrieval of data from a velocity wall law
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// Temperature component
1513 int nb_comp = val_variable_inconnue.dimension(1);
1514// Variables for computations
1515 DoubleTab temperature_imposee_sommet(nb_som,nb_comp); // Array of computed values at nodes.
1516 IntTab sommet_interp(nb_som); // Array indicating whether a node is forced
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// Modelling choice
1524 int form_Tp = interp.get_formulation_Tp(); // Choice of the wall law used
1525 int boundary_type = interp.get_boundary_type(); // Choice of the boundary condition type
1526 double T_inlet = interp.get_T_inlet();
1527 double Pr = interp.get_Prandlt_mol();
1528// Array initialisation
1529 temperature_imposee_calculee = 0.0;
1530 temperature_imposee_sommet = 0.0;
1531 sommet_interp = 0;
1532// Precision
1533 double eps = 1e-6;
1534
1535// Start of T_imp computation at nodes
1536 for (int i = 0; i < nb_som; i++) // Loop over all nodes of the domain
1537 {
1538 if ((fluid_elems(i) >= 0.0)) // Forced nodes only
1539 {
1540 int itisok = 1;
1541 sommet_interp(i) = 1; //Mark node i as forced
1542 temperature_imposee_sommet(i,0) = variable_imposee_mod(i,0);
1543 // Retrieve y+ and u* at the current node
1544 double utau = val_ustar(i);
1545 double yplus = val_yplus(i);
1546 // Compute only if utau != 0 => classical wall law
1547 if (utau != 0)
1548 {
1549 // Retrieve y and yref, as well as the wall normal
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 // exception handling
1569 if (d1 < eps) itisok = 0; // point P is on the immersed boundary: assignment
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 // Retrieve fluid properties needed for the computation
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)); // Fluid element index
1582 champ_variable_inconnue.valeur_a_elem(xf, Tref, cells[0]); // Determine temperature at the reference fluid point
1583 double yplusref = y_ref * utau / nu;
1584 if (!form_Tp) //Kader wall law
1585 {
1586 thetap = interp.Kader(yplus,Pr); // Compute theta plus with Kader for the forced point (known y+)
1587 thetapref = interp.Kader(yplusref,Pr); // Compute theta plus with Kader for the fluid point
1588 }
1589 else
1590 {
1591 Cerr << "Law not implemented, another law must be chosen" << 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 << "Error in the choice of boundary condition type" << finl;
1615 exit();
1616 }
1617 }
1618 }
1619 }
1620
1621// Compute temperature for forced faces
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)) //Retrieve faces whose all nodes are forced
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; //Compute face temperature as average of nodal values
1637 }
1638 }
1639 }
1640}
1641
1642/*void Source_PDF_VEF::calculer_temperature_imposee_mean_wall_law()
1643{
1644//TODO
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// Domain variables
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// Imposed and computed value arrays
1654 DoubleTab& variable_imposee_mod = modele_lu_.variable_imposee_->valeurs();
1655 DoubleTab& temperature_imposee_calculee = variable_imposee_;
1656// Information arrays for the immersed wall
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// Discretisation of the computed quantity
1662 Champ_P1NC& champ_variable_inconnue = ref_cast(Champ_P1NC,equation().inconnue());
1663 DoubleTab& val_variable_inconnue = champ_variable_inconnue.valeurs();
1664// Retrieval of data from a velocity wall law
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// Temperature component
1670 int nb_comp = val_variable_inconnue.dimension(1);
1671// Variables for computations
1672 DoubleTab temperature_imposee_sommet(nb_som,nb_comp); // Array of computed values at nodes.
1673 IntTab sommet_interp(nb_som); // Array indicating whether a node is forced
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// Modelling choice
1681 int form_Tp = interp.get_formulation_Tp(); // Choice of the wall law used
1682 int boundary_type = interp.get_boundary_type(); // Choice of the boundary condition type
1683 double T_inlet = interp.get_T_inlet();
1684 double Pr = interp.get_Prandlt_mol();
1685// Array initialisation
1686 temperature_imposee_calculee = 0.0;
1687 temperature_imposee_sommet = 0.0;
1688 sommet_interp = 0;
1689// Precision
1690 double eps = 1e-6;
1691
1692// Start of T_imp computation at nodes
1693 for (int i = 0; i < nb_som; i++) // Loop over all nodes of the domain
1694 {
1695 if ((fluid_elems(i) >= 0.0)) // Forced nodes only
1696 {
1697 int itisok = 1;
1698 sommet_interp(i) = 1; //Mark node i as forced
1699 temperature_imposee_sommet(i,0) = variable_imposee_mod(i,0);
1700 // Retrieve y+ and u* at the current node
1701 double utau = val_ustar(i);
1702 double yplus = val_yplus(i);
1703 // Compute only if utau != 0 => classical wall law
1704 if (utau != 0)
1705 {
1706 // Retrieve y and yref, as well as the wall normal
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 // exception handling
1726 if (d1 < eps) itisok = 0; // point P is on the immersed boundary: assignment
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 // Retrieve fluid properties needed for the computation
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)); // Fluid element index
1740 champ_variable_inconnue.valeur_a_elem(xf, Tref, cells[0]); // Determine temperature at the reference fluid point
1741 double yplusref = y_ref * utau / nu;
1742 if (!form_Tp) //Kader wall law
1743 {
1744 thetap = interp.Kader(yplus,Pr); // Compute theta plus with Kader for the forced point (known y+)
1745 thetapref = interp.Kader(yplusref,Pr); // Compute theta plus with Kader for the fluid point
1746 }
1747 else
1748 {
1749 Cerr << "Law not implemented, another law must be chosen" << 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 << "Error in the choice of boundary condition type" << finl;
1773 exit();
1774 }
1775 }
1776 }
1777 }
1778
1779// Compute temperature for forced faces
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)) //Retrieve faces whose all nodes are forced
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; //Compute face temperature as average of nodal values
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 // Filter by 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 // Filter and dof by 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 Warning: only defined in 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 or 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; time derivative in 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; /* nothing found */
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; /* nothing found */
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 // Initialize the volumetric u_star field to 0
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 // Initialize the volumetric y_plus field to 0
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
Overrides Champ_base::valeurs() Returns the array of values.
Definition Champ_Don_base.h:49
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Returns the array of field values at the current time.
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
Causes an error! Must be overridden by derived classes.
Definition Champ_P1NC.h:67
Champ_Proto::valeurs
virtual DoubleTab & valeurs()=0
Champ_base
class Champ_base This class is the base of the fields hierarchy.
Definition Champ_base.h:43
Cond_lim_base
class Cond_lim_base Base class for the hierarchy of classes that represent the different boundary con...
Definition Cond_lim_base.h:40
Dirichlet
Dirichlet This class is the base class of the hierarchy of Dirichlet-type boundary conditions.
Definition Dirichlet.h:31
Discretisation_base
class Discretisation_base This class represents a spatial discretization scheme, which
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
Returns the number of faces of type i of the geometric elements that make up the domain.
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:38
Domaine_Cl_dis_base
class Domaine_Cl_dis_base Domaine_Cl_dis_base objects represent discretized boundary conditions
Definition Domaine_Cl_dis_base.h:37
Domaine_Cl_dis_base::nb_cond_lim
int nb_cond_lim() const
Returns the number of boundary conditions.
Definition Domaine_Cl_dis_base.cpp:425
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Returns the i-th boundary condition.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_VF::nb_faces
int nb_faces() const
Returns the total number of 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
Returns the index of the i-th vertex of face num_face.
Definition Domaine_VF.h:582
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
Returns the index of the i-th face of element num_elem; the face numbering convention is.
Definition Domaine_VF.h:542
Domaine_dis_base
class Domaine_dis_base This class is the base of the hierarchy of discretized domains.
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
class Equation_base The role of an equation is the calculation of one or more fields....
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()
Returns the problem associated with the 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
Matrice_Morse class - Represents a (sparse) matrix M, not necessarily square,.
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:876
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
Returns the heat capacity of the medium (const version).
Definition Milieu_base.cpp:857
Milieu_base::masse_volumique
virtual const Champ_base & masse_volumique() const
Returns the mass density of the medium (const version).
Definition Milieu_base.cpp:793
MorEqn::equation
const Equation_base & equation() const
Returns the reference to the equation pointed to by MorEqn::mon_equation.
Definition MorEqn.h:62
Motcle
A character string (Nom) in uppercase.
Definition Motcle.h:26
Navier_Stokes_std
Navier_Stokes_std This class carries the terms of the momentum equation.
Definition Navier_Stokes_std.h:55
Navier_Stokes_std::vitesse_pour_transport
virtual const Champ_base & vitesse_pour_transport() const
Definition Navier_Stokes_std.cpp:272
Navier_Stokes_std::nombre_d_operateurs
int nombre_d_operateurs() const override
Returns the number of operators in the equation: 2 for the standard Navier-Stokes equations.
Definition Navier_Stokes_std.cpp:442
Nom
class Nom: a character string for naming TRUST objects.
Definition Nom.h:31
Nom::getString
const std::string & getString() const
Definition Nom.h:92
Noms
An array of character strings (VECT(Nom)).
Definition Noms.h:26
Objet_U::dimension
static int dimension
Definition Objet_U.h:94
Objet_U::que_suis_je
const Nom & que_suis_je() const
Returns the string identifying the class.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
Objet_U::le_nom
virtual const Nom & le_nom() const
Returns the name of the Objet_U. Virtual method to override: returns "neant" in this implementation.
Definition Objet_U.cpp:317
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
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
class Probleme_base It is a Probleme_U that is not a coupling.
Definition Probleme_base.h:55
Probleme_base::get_champ
const Champ_base & get_champ(const Motcle &nom) const override
Definition Probleme_base.cpp:836
Probleme_base::discretisation
const Discretisation_base & discretisation() const
Returns the discretization associated with the problem.
Definition Probleme_base.h:241
Probleme_base::schema_temps
const Schema_Temps_base & schema_temps() const
Returns the time scheme associated with the problem.
Definition Probleme_base.cpp:577
Process::abort
static void abort()
Abort routine for TRUST on a fatal error.
Definition Process.cpp:573
Process::exit
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Process::je_suis_maitre
static int je_suis_maitre()
Returns 1 if on the master processor of the current group (i.e. me() == 0), 0 otherwise.
Definition Process.cpp:82
Schema_Temps_base
class Schema_Temps_base
Definition Schema_Temps_base.h:67
Schema_Temps_base::temps_courant
double temps_courant() const
Returns the current time.
Definition Schema_Temps_base.h:446
Schema_Temps_base::pas_de_temps
double pas_de_temps() const
Returns the current time step (delta_t).
Definition Schema_Temps_base.h:394
Sortie
Base class for output streams.
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_
IBM u* values computed locally.
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 to the implicit matrix of source terms, by default no 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_
IBM d+ values computed locally.
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
Opening/creation of a file for printing a source term. To override in derived classes.
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