17#include <ParoiVDF_TBLE.h>
18#include <Domaine_Cl_VDF.h>
19#include <Dirichlet_paroi_fixe.h>
20#include <Dirichlet_paroi_defilante.h>
21#include <Champ_Uniforme.h>
22#include <Fluide_base.h>
25#include <Schema_Temps_base.h>
26#include <Operateur_Grad.h>
27#include <Probleme_base.h>
28#include <Modele_turbulence_scal_base.h>
29#include <Paroi_TBLE_scal_VDF.h>
30#include <Terme_Boussinesq_VDF_Face.h>
32#include <EcrFicPartage.h>
33#include <Navier_Stokes_std.h>
34#include <Modele_turbulence_hyd_base.h>
47 Cerr <<
"Reading of data for a " <<
que_suis_je() <<
" wall law" << finl;
50 param.lire_avec_accolades_depuis(is);
88 Cerr <<
"debut init_lois_paroi" << finl;
91 const IntVect& orientation = domaine_VDF.
orientation();
93 const IntTab& elem_faces = domaine_VDF.
elem_faces();
94 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
97 int compteur_faces_paroi = 0;
106 corresp.resize(le_dom_dis_->nb_faces_bord());
108 compteur_faces_paroi = 0;
110 for (
int n_bord = 0; n_bord < domaine_VDF.
nb_front_Cl(); n_bord++)
112 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
118 const int nfin = ndeb + le_bord.
nb_faces();
120 for (
int num_face = ndeb; num_face < nfin; num_face++)
122 corresp[num_face]=compteur_faces_paroi;
123 compteur_faces_paroi++;
128 compteur_faces_paroi = 0;
130 for (
int n_bord = 0; n_bord < domaine_VDF.
nb_front_Cl(); n_bord++)
136 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
142 const int nfin = ndeb + le_bord.
nb_faces();
147 for (
int num_face = ndeb; num_face < nfin; num_face++)
150 int ori = orientation(num_face);
152 if ((elem = face_voisins(num_face, 0)) == -1)
153 elem = face_voisins(num_face, 1);
164 eq_vit[corresp[num_face]].set_y0(0.);
165 eq_vit[corresp[num_face]].set_yn(dist);
168 eq_vit[corresp[num_face]].set_u_y0(0,0.);
170 vmoy = 0.5*(vit(elem_faces(elem, (ori + 1) % 4)) + vit(elem_faces(elem, (ori + 3) % 4)));
172 eq_vit[corresp[num_face]].set_u_yn(0, vmoy);
176 eq_vit[corresp[num_face]].set_Uinit_lin(0, 0., vmoy);
178 for (
int itble = 0; itble <
nb_pts + 1; itble++)
179 eq_vit[corresp[num_face]].set_Uinit(0, itble,
valeurs_reprises(corresp[num_face], 0, itble));
181 compteur_faces_paroi++;
187 for (
int num_face = ndeb; num_face < nfin; num_face++)
190 const int ori = orientation(num_face);
192 if ((elem = face_voisins(num_face, 0)) == -1)
193 elem = face_voisins(num_face, 1);
201 eq_vit[corresp[num_face]].set_y0(0.);
202 eq_vit[corresp[num_face]].set_yn(dist);
206 eq_vit[corresp[num_face]].set_u_y0(0, 0.);
207 eq_vit[corresp[num_face]].set_u_y0(1, 0.);
209 vmoy = 0.5*(vit(elem_faces(elem, (ori + 1) % 6)) + vit(elem_faces(elem, (ori + 4) % 6)));
211 eq_vit[corresp[num_face]].set_u_yn(0, vmoy);
215 eq_vit[corresp[num_face]].set_Uinit_lin(0, 0., vmoy);
217 for (
int itble = 0; itble <
nb_pts + 1; itble++)
218 eq_vit[corresp[num_face]].set_Uinit(0, itble,
valeurs_reprises(corresp[num_face], 0, itble));
220 vmoy = 0.5*(vit(elem_faces(elem, (ori + 2) % 6)) + vit(elem_faces(elem, (ori + 5) % 6)));
221 eq_vit[corresp[num_face]].set_u_yn(1, vmoy);
223 eq_vit[corresp[num_face]].set_Uinit_lin(1, 0., vmoy);
225 for (
int itble = 0; itble <
nb_pts + 1; itble++)
226 eq_vit[corresp[num_face]].set_Uinit(1, itble,
valeurs_reprises(corresp[num_face], 1, itble));
228 compteur_faces_paroi++;
241 int nb_sources=les_sources.size();
243 for (
int j=0; j<nb_sources; j++)
252 const DoubleTab& tab_champ_beta_t = ch_beta_t.
valeurs();
255 beta_t = std::max(tab_champ_beta_t(0,0),DMINFLOAT);
259 Cerr <<
" On ne sait pas traiter un beta_t non uniforme dans TBLE !!!"<< finl;
270 Cerr <<
"fin init_lois_paroi" << finl;
298 const IntVect& orientation = domaine_VDF.
orientation();
299 const IntTab& face_voisins = domaine_VDF.
face_voisins();
300 const IntTab& elem_faces = domaine_VDF.
elem_faces();
302 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
306 const DoubleTab& tab_visco = ch_visco_cin.
valeurs();
314 Cerr <<
" On ne sait pas traiter la gravite non uniforme dans TBLE !!!"<< finl;
323 visco = std::max(tab_visco(0,0),DMINFLOAT);
333 Cerr <<
"In ParoiVDF_TBLE::calculer_hyd_BiK : visco = " << tab_visco.
local_min_vect() <<
" <= 0 ? " << finl;
345 double gradient_de_pression0 = 0., gradient_de_pression1 = 0., vmoy = 0., ts0 =0., ts1=0.;
350 DoubleTab grad_p(vit);
353 gradient.calculer(p, grad_p);
354 const DoubleTab& visco_turb = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
357 DoubleTab termes_sources;
362 int nb_sources=les_sources.size();
363 for (
int j=0; j<nb_sources; j++)
377 int compteur_faces_paroi = 0;
382 compteur_faces_paroi = 0;
384 for (
int n_bord=0; n_bord<domaine_VDF.
nb_front_Cl(); n_bord++)
391 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
411 for (
int num_face=ndeb; num_face<nfin; num_face++)
413 ori = orientation(num_face);
418 if ((elem = face_voisins(num_face,0)) == -1)
420 elem = face_voisins(num_face,1);
427 eq_vit[corresp[num_face]].set_dt(1.e12);
429 eq_vit[corresp[num_face]].set_dt(dt);
434 int face1 = elem_faces(elem,(ori+1));
435 int face2 = elem_faces(elem,(ori+3)%4);
442 DoubleVect ts_boussi(
nb_pts+1);
447 const DoubleVect& gravite = ch_gravite.
valeurs();
461 gn+=gravite(idim)*domaine_VDF.
face_normales(num_face,idim)/norm_n;
468 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.
face_normales(num_face,idim)/norm_n;
471 double g_t=gt_vect(1-ori);
473 for(
int i=0; i<
nb_pts+1; i++)
484 vmoy = 0.5*(vit(face1) + vit(face2));
488 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
489 + termes_sources(face2)/volumes_entrelaces(face2));
491 eq_vit[corresp[num_face]].set_u_y0(0,0.);
493 eq_vit[corresp[num_face]].set_u_yn(0,vmoy);
499 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
501 for(
int i=0 ; i<
nb_pts+1; i++)
502 eq_vit[compteur_faces_paroi].set_F(0, i, -gradient_de_pression0+ts0+ts_boussi(i));
509 eq_vit[compteur_faces_paroi].aller_au_temps(tps);
514 if(itmax<eq_vit[corresp[num_face]].get_it())
516 itmax = eq_vit[compteur_faces_paroi].get_it();
518 Cerr <<
"WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
553 d_visco = tab_visco[elem];
555 tab_u_star_(num_face) = pow(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0),0.25);
566 else if (y_plus>8 && y_plus<30)
570 double deriv = Fdypar_direct(lm_plus);
571 double x = lm_plus*u_star*deriv;
572 tab_k(elem) = x*x/sqrt(
Cmu_) ;
573 tab_eps(elem) = (tab_k(elem)*u_star*u_star*deriv)*sqrt(
Cmu_)/d_visco;
578 tab_k(elem) =us2/sqrt(
Cmu_);
583 compteur_faces_paroi++;
599 for (
int num_face=ndeb; num_face<nfin; num_face++)
602 ori = orientation(num_face);
609 if ((elem = face_voisins(num_face,0)) == -1)
611 elem = face_voisins(num_face,1);
618 eq_vit[corresp[num_face]].set_dt(1.e12);
620 eq_vit[corresp[num_face]].set_dt(dt);
624 int face1 = elem_faces(elem,(ori+1));
625 int face2 = elem_faces(elem,(ori+4)%6);
627 int face3 = elem_faces(elem,(ori+2));
628 int face4 = elem_faces(elem,(ori+5)%6);
635 vmoy = 0.5*(vit(face1) + vit(face2));
639 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
640 + termes_sources(face2)/volumes_entrelaces(face2));
642 eq_vit[corresp[num_face]].set_u_y0(0,0.);
644 eq_vit[corresp[num_face]].set_u_yn(0,vmoy);
650 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
656 vmoy = 0.5*(vit(face3) + vit(face4));
658 ts1 = 0.5*(termes_sources(face3)/volumes_entrelaces(face3)
659 + termes_sources(face4)/volumes_entrelaces(face4));
661 eq_vit[corresp[num_face]].set_u_y0(1,0.);
663 eq_vit[corresp[num_face]].set_u_yn(1,vmoy);
671 gradient_de_pression1 = 0.5*(grad_p(face3)/volumes_entrelaces(face3)+grad_p(face4)/volumes_entrelaces(face4));
674 DoubleTab ts_boussi(
nb_pts+1,2);
679 const DoubleVect& gravite = ch_gravite.
valeurs();
691 gn+=gravite(idim)*domaine_VDF.
face_normales(num_face,idim)/norm_n;
698 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.
face_normales(num_face,idim)/norm_n;
701 for(
int i=0; i<
nb_pts+1; i++)
712 calculer_convection(num_face,face1, face2, face3, face4, elem, ndeb, nfin, ori, gradient_de_pression0, ts0, gradient_de_pression1, ts1);
715 for(
int i=0 ; i<
nb_pts+1 ; i++)
717 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 + ts_boussi(i,0));
718 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 + ts_boussi(i,1));
734 eq_vit[corresp[num_face]].aller_au_temps(tps);
738 if(itmax<eq_vit[corresp[num_face]].get_it())
740 itmax = eq_vit[corresp[num_face]].get_it();
742 Cerr <<
"WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
781 d_visco = tab_visco[elem];
783 tab_u_star_(num_face) = sqrt(sqrt(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0)
784 +eq_vit[corresp[num_face]].get_cis(1)*eq_vit[corresp[num_face]].get_cis(1)));
796 else if (y_plus>8 && y_plus<30)
800 double deriv = Fdypar_direct(lm_plus);
801 double x = lm_plus*u_star*deriv;
802 tab_k(elem) = x*x/sqrt(
Cmu_);
803 tab_eps(elem) = (tab_k(elem)*u_star*u_star*deriv)*sqrt(
Cmu_)/d_visco;
808 tab_k(elem)=us2/sqrt(
Cmu_);
812 compteur_faces_paroi++;
822 fic << tps <<
" " << itmax << finl;
846 const IntVect& orientation = domaine_VDF.
orientation();
847 const IntTab& face_voisins = domaine_VDF.
face_voisins();
848 const IntTab& elem_faces = domaine_VDF.
elem_faces();
850 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
854 const DoubleTab& tab_visco = ch_visco_cin.
valeurs();
862 Cerr <<
" On ne sait pas traiter la gravite non uniforme dans TBLE !!!"<< finl;
871 visco = std::max(tab_visco(0,0),DMINFLOAT);
881 Cerr <<
"In ParoiVDF_TBLE::calculer_hyd : visco = " << tab_visco.
local_min_vect() <<
" <= 0 ? " << finl;
893 double gradient_de_pression0 = 0., gradient_de_pression1 = 0., vmoy = 0., ts0 =0., ts1=0.;
898 DoubleTab grad_p(vit);
901 gradient.calculer(p, grad_p);
902 const DoubleTab& visco_turb = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
905 DoubleTab termes_sources;
910 int nb_sources=les_sources.size();
911 for (
int j=0; j<nb_sources; j++)
925 int compteur_faces_paroi = 0;
930 compteur_faces_paroi = 0;
932 for (
int n_bord=0; n_bord<domaine_VDF.
nb_front_Cl(); n_bord++)
939 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
959 for (
int num_face=ndeb; num_face<nfin; num_face++)
961 ori = orientation(num_face);
966 if ((elem = face_voisins(num_face,0)) == -1)
968 elem = face_voisins(num_face,1);
975 eq_vit[corresp[num_face]].set_dt(1.e12);
977 eq_vit[corresp[num_face]].set_dt(dt);
982 int face1 = elem_faces(elem,(ori+1));
983 int face2 = elem_faces(elem,(ori+3)%4);
990 DoubleVect ts_boussi(
nb_pts+1);
995 const DoubleVect& gravite = ch_gravite.
valeurs();
1007 for (
int idim=0; idim<
dimension; idim++)
1009 gn+=gravite(idim)*domaine_VDF.
face_normales(num_face,idim)/norm_n;
1014 for (
int idim=0; idim<
dimension; idim++)
1016 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.
face_normales(num_face,idim)/norm_n;
1019 double g_t=gt_vect(1-ori);
1021 for(
int i=0; i<
nb_pts+1; i++)
1035 vmoy = 0.5*(vit(face1) + vit(face2));
1039 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1040 + termes_sources(face2)/volumes_entrelaces(face2));
1042 eq_vit[corresp[num_face]].set_u_y0(0,0.);
1044 eq_vit[corresp[num_face]].set_u_yn(0,vmoy);
1050 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
1052 for(
int i=0 ; i<
nb_pts+1; i++)
1053 eq_vit[compteur_faces_paroi].set_F(0, i, -gradient_de_pression0+ts0+ts_boussi(i));
1060 eq_vit[compteur_faces_paroi].aller_au_temps(tps);
1065 if(itmax<eq_vit[corresp[num_face]].get_it())
1067 itmax = eq_vit[compteur_faces_paroi].get_it();
1069 Cerr <<
"WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
1104 d_visco = tab_visco[elem];
1106 tab_u_star_(num_face) = pow(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0),0.25);
1116 tab1(elem)=eq_vit[corresp[num_face]].get_nu_t_yn();
1127 else if (y_plus>8 && y_plus<30)
1131 double deriv = Fdypar_direct(lm_plus);
1132 double x = lm_plus*u_star*deriv;
1133 tab1(elem,0) = x*x/sqrt(
Cmu_) ;
1134 tab1(elem,1) = (tab1(elem,0)*u_star*u_star*deriv)*sqrt(
Cmu_)/d_visco;
1139 tab1(elem,0)=us2/sqrt(
Cmu_);
1144 compteur_faces_paroi++;
1160 for (
int num_face=ndeb; num_face<nfin; num_face++)
1163 ori = orientation(num_face);
1170 if ((elem = face_voisins(num_face,0)) == -1)
1172 elem = face_voisins(num_face,1);
1179 eq_vit[corresp[num_face]].set_dt(1.e12);
1181 eq_vit[corresp[num_face]].set_dt(dt);
1185 int face1 = elem_faces(elem,(ori+1));
1186 int face2 = elem_faces(elem,(ori+4)%6);
1188 int face3 = elem_faces(elem,(ori+2));
1189 int face4 = elem_faces(elem,(ori+5)%6);
1196 vmoy = 0.5*(vit(face1) + vit(face2));
1200 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1201 + termes_sources(face2)/volumes_entrelaces(face2));
1203 eq_vit[corresp[num_face]].set_u_y0(0,0.);
1205 eq_vit[corresp[num_face]].set_u_yn(0,vmoy);
1211 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
1217 vmoy = 0.5*(vit(face3) + vit(face4));
1219 ts1 = 0.5*(termes_sources(face3)/volumes_entrelaces(face3)
1220 + termes_sources(face4)/volumes_entrelaces(face4));
1222 eq_vit[corresp[num_face]].set_u_y0(1,0.);
1224 eq_vit[corresp[num_face]].set_u_yn(1,vmoy);
1232 gradient_de_pression1 = 0.5*(grad_p(face3)/volumes_entrelaces(face3)+grad_p(face4)/volumes_entrelaces(face4));
1235 DoubleTab ts_boussi(
nb_pts+1,2);
1240 const DoubleVect& gravite = ch_gravite.
valeurs();
1250 for (
int idim=0; idim<
dimension; idim++)
1252 gn+=gravite(idim)*domaine_VDF.
face_normales(num_face,idim)/norm_n;
1257 for (
int idim=0; idim<
dimension; idim++)
1259 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.
face_normales(num_face,idim)/norm_n;
1262 for(
int i=0; i<
nb_pts+1; i++)
1279 calculer_convection(num_face,face1, face2, face3, face4, elem, ndeb, nfin, ori, gradient_de_pression0, ts0, gradient_de_pression1, ts1);
1282 for(
int i=0 ; i<
nb_pts+1 ; i++)
1284 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 + ts_boussi(i,0));
1285 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 + ts_boussi(i,1));
1301 eq_vit[corresp[num_face]].aller_au_temps(tps);
1305 if(itmax<eq_vit[corresp[num_face]].get_it())
1307 itmax = eq_vit[corresp[num_face]].get_it();
1309 Cerr <<
"WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
1348 d_visco = tab_visco[elem];
1350 tab_u_star_(num_face) = sqrt(sqrt(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0)
1351 +eq_vit[corresp[num_face]].get_cis(1)*eq_vit[corresp[num_face]].get_cis(1)));
1359 tab1(elem)=eq_vit[corresp[num_face]].get_nu_t_yn();
1373 else if (y_plus>8 && y_plus<30)
1377 double deriv = Fdypar_direct(lm_plus);
1378 double x = lm_plus*u_star*deriv;
1379 tab1(elem,0) = x*x/sqrt(
Cmu_);
1380 tab1(elem,1) = (tab1(elem,0)*u_star*u_star*deriv)*sqrt(
Cmu_)/d_visco;
1385 tab1(elem,0)=us2/sqrt(
Cmu_);
1390 compteur_faces_paroi++;
1399 SFichier fic(
"iter.dat", ios::app);
1400 fic << tps <<
" " << itmax << finl;
1415int ParoiVDF_TBLE::calculer_stats()
1418 const IntVect& orientation = domaine_VDF.
orientation();
1420 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1434 for(
int i=0 ; i<
nb_pts+1 ; i++)
1449 for(
int i=0 ; i<
nb_pts+1 ; i++)
1461 for(
int i=0 ; i<
nb_pts+1 ; i++)
1473 for(
int i=0 ; i<
nb_pts+1 ; i++)
1492 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1496 const IntVect& orientation = domaine_VDF.
orientation();
1511 for(
int i=0; i<
nb_pts+1; i++)
1517 for(
int i=0; i<
nb_pts+1; i++)
1521 for(
int i=0; i<
nb_pts+1; i++)
1525 for(
int i=0; i<
nb_pts+1; i++)
1540 double tps = mon_modele_turb_hyd->equation().inconnue().temps();
1550 double tps_reprise = mon_modele_turb_hyd->equation().schema_temps().temps_courant();
1562void ParoiVDF_TBLE::calculer_convection(
int num_face,
int face1,
int face2,
int face3,
int face4,
int elem,
int ndeb,
int nfin,
int ori,
double gradient_de_pression0,
double ts0,
double gradient_de_pression1,
double ts1)
1565 const IntTab& face_voisins = domaine_VDF.
face_voisins();
1566 const IntTab& elem_faces = domaine_VDF.
elem_faces();
1569 double d00d0 = 0., d0vdy= 0., d01d1= 0., d10d0= 0., d1vdy= 0., d11d1=0.;
1570 int elem_gauche0, elem_gauche1, elem_droite0, elem_droite1;
1571 int face_gauche0, face_gauche1, face_droite0, face_droite1;
1580 elem_gauche0 = face_voisins(face1,0);
1581 elem_droite0 = face_voisins(face2,1);
1582 if (elem_gauche0 == elem)
1584 elem_droite0 = face_voisins(face1,1);
1585 elem_gauche0 = face_voisins(face2,0);
1588 elem_gauche1 = face_voisins(face3,0);
1589 elem_droite1 = face_voisins(face4,1);
1590 if (elem_gauche1 == elem)
1592 elem_droite1 = face_voisins(face3,1);
1593 elem_gauche1 = face_voisins(face4,0);
1597 if((elem_droite1!=-1)&&(elem_droite0!=-1)&&(elem_gauche1!=-1)&&(elem_gauche0!=-1))
1600 face_gauche0 = elem_faces(elem_gauche0, ori);
1601 if ((face_voisins(face_gauche0,0) != -1)&&(face_voisins(face_gauche0,1) != -1))
1603 face_gauche0 = elem_faces(elem_gauche0, ori+3);
1605 face_droite0 = elem_faces(elem_droite0, ori);
1606 if ((face_voisins(face_droite0,0) != -1)&&(face_voisins(face_droite0,1) != -1))
1608 face_droite0 = elem_faces(elem_droite0, ori+3);
1611 face_gauche1 = elem_faces(elem_gauche1, ori);
1612 if ((face_voisins(face_gauche1,0) != -1)&&(face_voisins(face_gauche1,1) != -1))
1614 face_gauche1 = elem_faces(elem_gauche1, ori+3);
1616 face_droite1 = elem_faces(elem_droite1, ori);
1617 if ((face_voisins(face_droite1,0) != -1)&&(face_voisins(face_droite1,1) != -1))
1619 face_droite1 = elem_faces(elem_droite1, ori+3);
1623 if((face_droite1<nfin)&&(face_droite0<nfin)&&(face_gauche1<nfin)&&(face_gauche0<nfin)&&
1624 (ndeb<=face_droite1)&&(ndeb<=face_droite0)&&(ndeb<=face_gauche1)&&(ndeb<=face_gauche0))
1633 if(eq_vit[corresp[num_face]].get_Unp1(0,1)>0)
1634 d0d0_1 = (eq_vit[corresp[face_droite0]].get_Unp1(0,1)-eq_vit[corresp[num_face]].get_Unp1(0,1))
1637 d0d0_1 = (eq_vit[corresp[num_face]].get_Unp1(0,1)-eq_vit[corresp[face_gauche0]].get_Unp1(0,1))
1640 if(eq_vit[corresp[num_face]].get_Unp1(0,1)>0)
1641 d1d1_1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,1)-eq_vit[corresp[num_face]].get_Unp1(1,1))
1644 d1d1_1 = (eq_vit[corresp[num_face]].get_Unp1(1,1)-eq_vit[corresp[face_gauche1]].get_Unp1(1,1))
1647 v += -(d0d0_1+d1d1_1)
1648 *0.5*(eq_vit[corresp[num_face]].get_yc(1)-eq_vit[corresp[num_face]].get_yc(0));
1650 eq_vit[corresp[num_face]].set_v(i,v);
1656 for(i=2 ; i<
nb_pts ; i++)
1659 (eq_vit[corresp[face_droite0]].get_Unp1(0,i)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1662 (eq_vit[corresp[face_droite1]].get_Unp1(1,i)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1665 *0.5*(eq_vit[corresp[num_face]].get_yc(i)-eq_vit[corresp[num_face]].get_yc(i-1))
1667 (eq_vit[corresp[face_droite0]].get_Unp1(0,i-1)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i-1))
1670 (eq_vit[corresp[face_droite1]].get_Unp1(1,i-1)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i-1))
1673 *0.5*(eq_vit[corresp[num_face]].get_yc(i-1)-eq_vit[corresp[num_face]].get_yc(i-2));
1675 eq_vit[corresp[num_face]].set_v(i,v);
1686 (eq_vit[corresp[face_droite0]].get_Unp1(0,i)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1689 (eq_vit[corresp[face_droite1]].get_Unp1(1,i)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1692 *0.5*(eq_vit[corresp[num_face]].get_yc(
nb_pts-1)-eq_vit[corresp[num_face]].get_yc(
nb_pts-2));
1694 eq_vit[corresp[num_face]].set_v(
nb_pts,v);
1705 eq_vit[corresp[num_face]].set_F(0, 0, - gradient_de_pression0 + ts0);
1706 eq_vit[corresp[num_face]].set_F(1, 0, - gradient_de_pression1 + ts1);
1708 for(i=1 ; i<
nb_pts ; i++)
1710 d00d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(0,i)
1711 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1715 d01d1 = (eq_vit[corresp[face_droite1]].get_Unp1(0,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1716 -eq_vit[corresp[face_gauche1]].get_Unp1(0,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1720 d10d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(1,i)
1721 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(1,i))
1725 d11d1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1726 -eq_vit[corresp[face_gauche1]].get_Unp1(1,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1730 d0vdy = (eq_vit[corresp[num_face]].get_Unp1(0,i+1)*eq_vit[corresp[num_face]].get_v(i+1)
1731 -eq_vit[corresp[num_face]].get_Unp1(0,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1732 /(eq_vit[corresp[num_face]].get_y(i+1)-eq_vit[corresp[num_face]].get_y(i-1));
1734 d1vdy = (eq_vit[corresp[num_face]].get_Unp1(1,i+1)*eq_vit[corresp[num_face]].get_v(i+1)
1735 -eq_vit[corresp[num_face]].get_Unp1(1,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1736 /(eq_vit[corresp[num_face]].get_y(i+1)-eq_vit[corresp[num_face]].get_y(i-1));
1738 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 - (d00d0+d01d1+d0vdy));
1739 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 - (d11d1+d10d0+d1vdy));
1745 d00d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(0,i)
1746 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1750 d01d1 = (eq_vit[corresp[face_droite1]].get_Unp1(0,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1751 -eq_vit[corresp[face_gauche1]].get_Unp1(0,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1755 d10d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(1,i)
1756 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(1,i))
1760 d11d1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1761 -eq_vit[corresp[face_gauche1]].get_Unp1(1,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1765 d0vdy = (eq_vit[corresp[num_face]].get_Unp1(0,i)*eq_vit[corresp[num_face]].get_v(i)
1766 -eq_vit[corresp[num_face]].get_Unp1(0,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1767 /(eq_vit[corresp[num_face]].get_y(i)-eq_vit[corresp[num_face]].get_y(i-1));
1769 d1vdy = (eq_vit[corresp[num_face]].get_Unp1(1,i)*eq_vit[corresp[num_face]].get_v(i)
1770 -eq_vit[corresp[num_face]].get_Unp1(1,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1771 /(eq_vit[corresp[num_face]].get_y(i)-eq_vit[corresp[num_face]].get_y(i-1));
1777 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 - (d00d0+d01d1+d0vdy));
1778 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 - (d11d1+d10d0+d1vdy));
1786 for(
int i=0 ; i<
nb_pts+1 ; i++)
1788 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0);
1789 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1);
classe Champ_Don_base classe de base des Champs donnes (non calcules)
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
double temps() const
Renvoie le temps du champ.
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Classe Diffu_totale_hyd_base Classe abstraite calculant la diffusivite totale (somme diffusivite.
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
double dist_norm_bord_axi(int num_face) const
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
double dist_elem_period(int, int, int) const
double face_normales(int, int) const override
double dist_norm_bord(int num_face) const override
virtual const DoubleVect & face_surfaces() const
int nb_faces() const
renvoie le nombre global de faces.
DoubleVect & volumes_entrelaces()
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Sortie & syncfile() override
Provoque l'ecriture sur disque des donnees accumulees sur les differents processeurs depuis le dernie...
Class defining operators and methods for all reading operation in an input flow (file,...
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
virtual const Milieu_base & milieu() const =0
virtual const RefObjU & get_modele(Type_modele type) const
Sources & sources()
Renvoie les termes sources asssocies a l'equation.
virtual const Champ_Inc_base & inconnue() const =0
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
classe Fluide_base Cette classe represente un d'un fluide incompressible ainsi que
const Champ_Don_base & viscosite_cinematique() const
int num_premiere_face() const
virtual const Champ_Don_base & beta_t() const
Renvoie beta_t du milieu.
virtual const Champ_Don_base & gravite() const
Renvoie la gravite du milieu si elle a ete associe provoque une erreur sinon.
Classe Modele_turbulence_scal_base Cette classe represente un modele de turbulence pour une equation ...
const Turbulence_paroi_scal_base & loi_paroi() const
Renvoie la loi de turbulence sur la paroi (version const).
Une chaine de caractere (Nom) en majuscules.
classe Navier_Stokes_std Cette classe porte les termes de l'equation de la dynamique
Operateur_Grad & operateur_gradient()
Renvoie l'operateur de calcul du gradient associe a l'equation.
Champ_Inc_base & pression()
class Nom Une chaine de caractere pour nommer les objets de TRUST
virtual void set_param(Param &) const
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
virtual const Nom & le_nom() const
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Classe Operateur_Grad Classe generique de la hierarchie des operateurs calculant le gradient.
Helper class to factorize the readOn method of Objet_U classes.
void ajouter_flag(const char *keyword, const bool *value)
Register a boolean flag whose mere presence switches it to true.
void ajouter_non_std(const char *keyword, const Objet_U *value, Param::Nature nat=Param::OPTIONAL)
Register a keyword handled by Objet_U::lire_motcle_non_standard.
void set_param(Param ¶m) const override
int reprendre(Entree &) override
Reprise d'un Objet_U sur un flot d'entree Methode a surcharger.
int calculer_hyd_BiK(DoubleTab &, DoubleTab &) override
void imprimer_ustar(Sortie &) const override
int lire_motcle_non_standard(const Motcle &, Entree &) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
int init_lois_paroi() override
int sauvegarder(Sortie &) const override
Sauvegarde d'un Objet_U sur un flot de sortie Methode a surcharger.
int calculer_hyd(DoubleTab &) override
const Probleme_base & getPbBase() const override
int est_initialise() const
Eq_couch_lim & get_eq_couche_lim_T(int i)
double tps_start_stat_nu_t_dyn
void calculer_stat(int indice_post, int indice_maillage, double Fx, double Fy, double Fz, double u, double v, double w, double dt)
void imprimer_stat(Sortie &, double) const
void set_param(Param ¶m) const
int reprendre(Entree &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, double tps)
IntTab num_global_faces_post
int sauvegarder(Sortie &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, double tps) const
DoubleTab valeurs_reprises
int lire_motcle_non_standard(const Motcle &, Entree &)
int init_lois_paroi(const Domaine_VF &, const Domaine_Cl_dis_base &)
CLASS: Paroi_ODVM_scal_VDF.
void set_param(Param ¶m) const
double calcul_lm_plus(double d_plus)
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
virtual const Equation_base & equation(int) const =0
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Cette classe est a la classe C++ ofstream ce que la classe Sortie est a la classe C++ ostream Elle re...
double temps_courant() const
Renvoie le temps courant.
double pas_temps_min() const
Renvoie le pas de temps minimum.
double pas_de_temps() const
Renvoie le pas de temps (delta_t) courant.
Classe de base des flux de sortie.
classe Source_base Un objet Source_base est un terme apparaissant au second membre d'une
virtual DoubleTab & ajouter(DoubleTab &) const
class Sources Sources represente une liste de Source.
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
_TYPE_ local_min_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
const Objet_U & valeur() const
class Terme_Boussinesq_scalaire_VDF_Face
Classe Terme_Boussinesq_base Cette classe represente le terme de gravite qui figure dans l'equation.
double Scalaire0(int i) const
const DoubleVect & tab_d_plus() const
const DoubleVect & tab_u_star() const
DoubleTab Cisaillement_paroi_
Classe Turbulence_paroi_scal_base Classe de base pour la hierarchie des classes representant les mode...