TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Op_Diff_VDF_Face_Axi_base.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 <Op_Diff_VDF_Face_Axi_base.h>
17
18Implemente_base(Op_Diff_VDF_Face_Axi_base,"Op_Diff_VDF_Face_Axi_base",Op_Diff_VDF_Face_base);
19
20Sortie& Op_Diff_VDF_Face_Axi_base::printOn(Sortie& s ) const { return s << que_suis_je() ; }
21Entree& Op_Diff_VDF_Face_Axi_base::readOn(Entree& s ) { return s ; }
22
23void Op_Diff_VDF_Face_Axi_base::associer(const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_cl_dis, const Champ_Inc_base& ch_transporte)
24{
25 const Domaine_VDF& zvdf = ref_cast(Domaine_VDF,domaine_dis);
26 const Domaine_Cl_VDF& zclvdf = ref_cast(Domaine_Cl_VDF,domaine_cl_dis);
27 const Champ_Face_VDF& inco = ref_cast(Champ_Face_VDF,ch_transporte);
28 le_dom_vdf = zvdf;
29 la_zcl_vdf = zclvdf;
30 inconnue = inco;
31 surface.ref(zvdf.face_surfaces());
32 volumes_entrelaces.ref(zvdf.volumes_entrelaces());
33 orientation.ref(zvdf.orientation());
34 porosite.ref(la_zcl_vdf->equation().milieu().porosite_face());
35 xp.ref(zvdf.xp());
36 xv.ref(zvdf.xv());
37 Qdm.ref(zvdf.Qdm());
38 face_voisins.ref(zvdf.face_voisins());
39 elem_faces.ref(zvdf.elem_faces());
40 type_arete_bord.ref(zclvdf.type_arete_bord());
41}
42
43double Op_Diff_VDF_Face_Axi_base::calculer_dt_stab() const
44{
45 return Op_Diff_VDF_base::calculer_dt_stab_(le_dom_vdf.valeur()) ;
46}
47
48void Op_Diff_VDF_Face_Axi_base::ajouter_elem(const DoubleTab& inco, DoubleTab& resu) const
49{
50 if (inco.line_size() > 1) not_implemented(__func__);
51 for (int num_elem = 0; num_elem < le_dom_vdf->nb_elem(); num_elem++)
52 {
53 const int fx0 = elem_faces(num_elem,0), fx1 = elem_faces(num_elem,dimension), fy0 = elem_faces(num_elem,1), fy1 = elem_faces(num_elem,1+dimension);
54 // Calcul de tau11
55 const double tau11 = (inco[fx1]-inco[fx0])/(xv(fx1,0) - xv(fx0,0));
56 // Calcul de tau22
57 double R = xp(num_elem,0), d_teta = xv(fy1,1) - xv(fy0,1);
58 if (d_teta < 0) d_teta += deux_pi;
59 double tau22 = (inco[fy1]-inco[fy0])/(R*d_teta);
60 tau22 += 0.5*(inco[fx0]+inco[fx1])/R; // termes supplementaires en axi
61
62 const double flux_X = tau11*nu_(num_elem)*0.5*(surface(fx0)+surface(fx1)), flux_Y = tau22*nu_(num_elem)*0.5*(surface(fy0)+surface(fy1));
63 resu[fx0] += flux_X;
64 resu[fx1] -= flux_X;
65 resu[fy0] += flux_Y;
66 resu[fy1] -= flux_Y;
67
68 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
69 const double coef_laplacien_axi = +0.5*tau22*nu_(num_elem);
70 resu[fx0] -= coef_laplacien_axi*volumes_entrelaces(fx0)*porosite(fx0)/xv(fx0,0);
71 resu[fx1] -= coef_laplacien_axi*volumes_entrelaces(fx1)*porosite(fx1)/xv(fx1,0);
72 }
73}
74
75void Op_Diff_VDF_Face_Axi_base::ajouter_elem_3D(const DoubleTab& inco, DoubleTab& resu) const
76{
77 for (int num_elem = 0; num_elem < le_dom_vdf->nb_elem(); num_elem++)
78 {
79 const int fz0 = elem_faces(num_elem,2), fz1 = elem_faces(num_elem,2+dimension);
80 // Calcul de tau33
81 const double tau33 = (inco[fz1]-inco[fz0])/(xv(fz1,2) - xv(fz0,2)), flux_Z = tau33*nu_(num_elem)*0.5*(surface(fz0)+surface(fz1));
82 resu[fz0] += flux_Z;
83 resu[fz1] -= flux_Z;
84 }
85}
86
87void Op_Diff_VDF_Face_Axi_base::ajouter_aretes_bords(const DoubleTab& inco, DoubleTab& resu) const
88{
89 int ndeb = le_dom_vdf->premiere_arete_bord(), nfin = ndeb + le_dom_vdf->nb_aretes_bord();
90 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
91 {
92 const int n_type = type_arete_bord(n_arete-ndeb);
93
94 switch(n_type)
95 {
96 case TypeAreteBordVDF::PAROI_PAROI: /* fall through */
97 case TypeAreteBordVDF::FLUIDE_FLUIDE:
98 case TypeAreteBordVDF::PAROI_FLUIDE:
99 {
100 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), signe = Qdm(n_arete,3), ori1 = orientation(fac1), ori3 = orientation(fac3);
101 const int rang1 = fac1 - le_dom_vdf->premiere_face_bord(), rang2 = fac2 - le_dom_vdf->premiere_face_bord();
102 double vit_imp, dist3, tps = inconnue->temps();
103
104 if (n_type == TypeAreteBordVDF::PAROI_FLUIDE) // arete paroi_fluide :il faut determiner qui est la face fluide
105 {
106 if (est_egal(inco[fac1],0)) vit_imp = Champ_Face_get_val_imp_face_bord(tps,rang2,ori3,la_zcl_vdf.valeur());
107 else vit_imp = Champ_Face_get_val_imp_face_bord(tps,rang1,ori3,la_zcl_vdf.valeur());
108 }
109 else vit_imp = 0.5*(Champ_Face_get_val_imp_face_bord(tps,rang1,ori3,la_zcl_vdf.valeur())+Champ_Face_get_val_imp_face_bord(tps,rang2,ori3,la_zcl_vdf.valeur()));
110
111 const double db_diffusivite = nu_mean_2_pts_(face_voisins(fac3,0),face_voisins(fac3,1));
112
113 if (ori1 == 0) // bord d'equation R = cte
114 {
115 double flux1;
116 if (ori3 == 1) // flux de tau12 a travers le bord
117 {
118 dist3 = xv(fac3,0)-xv(fac1,0);
119 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
120
121 const double tau12 = (inco[fac3]-vit_imp)/dist3;
122 flux1 = db_diffusivite*tau12*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
123 }
124 else //if (ori3 == 2) flux de tau13 a travers le bord
125 {
126 assert(ori3 == 2);
127 dist3 = xv(fac3,0)-xv(fac1,0);
128 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
129 const double tau13 = (inco[fac3]-vit_imp)/dist3;
130 flux1 = db_diffusivite*tau13*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
131 }
132 resu[fac3] += signe*flux1;
133 }
134 else if (ori1 == 1) // bord d'equation teta = cte
135 {
136 double R = xv(fac3,0), d_teta = xv(fac3,1) - xv(fac1,1);
137 if (d_teta < 0) d_teta += deux_pi;
138
139 dist3 = R*d_teta;
140 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
141 if (ori3 == 0) // flux de tau21 a travers le bord
142 {
143 double tau21 = (inco[fac3]-vit_imp)/dist3;
144 tau21 -= 0.5*(inco[fac1]+inco[fac2])/R; // Terme supplementaire en axi
145 const double flux2 = db_diffusivite*tau21*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
146 resu[fac3]+=signe*flux2;
147
148 // Termes supplementaires dans le laplacien en axi
149 // Ils sont integres comme des termes sources
150 const double coef_laplacien_axi = 0.5*db_diffusivite*tau21;
151 resu(fac1) += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
152 resu(fac2) += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
153 }
154 else if (ori3 == 2) // flux de tau23 a travers le bord
155 {
156 // XXX : attention si ecart : dans le cas constant c'etait tau23 = signe*(vit_imp-inco[fac3])/dist3 (normalement c'est la meme mais bon)
157 const double tau23 = (inco[fac3]-vit_imp)/dist3;
158 const double flux3 = db_diffusivite*tau23*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
159 resu[fac3] += signe*flux3;
160 }
161 }
162 else // (ori1 == 2) bord d'equation Z = cte
163 {
164 double flux4;
165 if (ori3 == 0) // flux de tau31 a travers le bord
166 {
167 dist3 = xv(fac3,2)-xv(fac1,2);
168 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
169
170 const double tau31 = (inco[fac3]-vit_imp)/dist3;
171 flux4 = db_diffusivite*tau31*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
172 }
173 else // if (ori3 == 1) flux de tau32 a travers le bord
174 {
175 assert(ori3 == 1) ;
176 dist3 = xv(fac3,2)-xv(fac1,2);
177 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
178
179 const double tau32 = (inco[fac3]-vit_imp)/dist3;
180 flux4 = db_diffusivite*tau32*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
181 }
182 resu[fac3] += signe*flux4;
183 }
184 break;
185 }
186 case TypeAreteBordVDF::NAVIER_NAVIER: // pas de flux diffusif calcule
187 break;
188 default :
189 {
190 Cerr << "On a rencontre un type d'arete non prevu : [ num arete : " << n_arete << " ], [ type : " << n_type << " ]" << finl;
191 Process::exit();
192 break;
193 }
194 }
195 }
196}
197
198void Op_Diff_VDF_Face_Axi_base::ajouter_aretes_mixtes_internes(const DoubleTab& inco, DoubleTab& resu) const
199{
200 const int ndeb = le_dom_vdf->premiere_arete_mixte(), nfin = le_dom_vdf->nb_aretes();
201 for (int n_arete=ndeb; n_arete<nfin; n_arete++)
202 {
203 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), fac4 = Qdm(n_arete,3), ori1 = orientation(fac1), ori3 = orientation(fac3);
204 const double db_diffusivite = nu_mean_4_pts_(fac3,fac4);
205
206 if (ori1 == 1) // (seule possibilite : ori3 =0) Arete XY
207 {
208 double flux1;
209 // Calcul de tau21
210 const double R = xv(fac3,0);
211 double d_teta = xv(fac4,1) - xv(fac3,1);
212 if (d_teta < 0) d_teta += deux_pi;
213 double tau21 = (inco(fac4)-inco(fac3))/(R*d_teta);
214
215 // Terme supplementaire en axi
216 tau21 -= 0.5*(inco[fac1]+inco[fac2])/R;
217
218 // flux de tau21 sur la facette a cheval sur les faces fac1 et fac2
219 flux1 = db_diffusivite*tau21*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
220 resu(fac3) += flux1;
221 resu(fac4) -= flux1;
222
223 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
224 const double coef_laplacien_axi = 0.5*db_diffusivite*tau21;
225 resu(fac1) += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
226 resu(fac2) += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
227
228 // Calcul de tau12
229 const double tau12 = (inco(fac2)-inco(fac1))/(xv(fac2,0) - xv(fac1,0));
230 // flux de tau12 sur la facette a cheval sur les faces fac3 et fac4
231 flux1 = db_diffusivite*tau12*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
232 resu(fac1) += flux1;
233 resu(fac2) -= flux1;
234 }
235 else if (ori3 == 1) // (seule possibilite ori1 = 2) arete YZ
236 {
237 double flux2;
238 // Calcul de tau32
239 const double tau32 = (inco(fac4)-inco(fac3))/(xv(fac4,2) - xv(fac3,2));
240
241 // flux de tau32 sur la facette a cheval sur les faces fac1 et fac2
242 flux2 = db_diffusivite*tau32*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
243 resu(fac3) += flux2;
244 resu(fac4) -= flux2;
245
246 // Calcul de tau23
247 const double R = xv(fac1,0);
248 double d_teta = xv(fac2,1) - xv(fac1,1);
249 if (d_teta < 0) d_teta += deux_pi;
250 const double tau23 = (inco(fac2)-inco(fac1))/(R*d_teta);
251 // flux de tau23 sur la facette a cheval sur les faces fac3 et fac4
252 flux2 = db_diffusivite*tau23*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
253 resu(fac1) += flux2;
254 resu(fac2) -= flux2;
255 }
256 else // seule possibilite ori1 = 2 et ori3 = 0: arete XZ
257 {
258 double flux3;
259 // Calcul de tau31
260 const double tau31 = (inco(fac4)-inco(fac3))/(xv(fac4,2) - xv(fac3,2));
261
262 // flux de tau31 sur la facette a cheval sur les faces fac1 et fac2
263 flux3 = db_diffusivite*tau31*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
264 resu(fac3) += flux3;
265 resu(fac4) -= flux3;
266
267 // Calcul de tau13
268 const double tau13 = (inco(fac2)-inco(fac1))/(xv(fac2,0) - xv(fac1,0));
269 // flux de tau13 sur la facette a cheval sur les faces fac3 et fac4
270 flux3 = db_diffusivite*tau13*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
271 resu(fac1) += flux3;
272 resu(fac2) -= flux3;
273 }
274 }
275}
276
277DoubleTab& Op_Diff_VDF_Face_Axi_base::ajouter(const DoubleTab& inco, DoubleTab& resu) const
278{
279 ajouter_elem(inco,resu); // Boucle sur les elements
280 if (dimension == 3) ajouter_elem_3D(inco,resu); // Boucle sur les elements supplementaires si 3D
281 ajouter_aretes_bords(inco, resu); // Boucle sur les aretes bord
282 ajouter_aretes_mixtes_internes(inco, resu); // Boucle sur les aretes mixtes et internes
283 return resu;
284}
285
286DoubleTab& Op_Diff_VDF_Face_Axi_base::calculer(const DoubleTab& inco, DoubleTab& resu) const
287{
288 resu = 0;
289 return ajouter(inco,resu);
290}
291
292void Op_Diff_VDF_Face_Axi_base::fill_coeff_matrice_morse(const int fac1, const int fac2, const double flux, Matrice_Morse& matrice) const
293{
294 const auto& tab1 = matrice.get_set_tab1();
295 const auto& tab2 = matrice.get_set_tab2();
296 auto& coeff = matrice.get_set_coeff();
297 for (auto k = tab1[fac1]-1; k < tab1[fac1+1]-1; k++)
298 {
299 if (tab2[k]-1 == fac1) coeff[k] += flux;
300 if (tab2[k]-1 == fac2) coeff[k] -= flux;
301 }
302 for (auto k = tab1[fac2]-1; k < tab1[fac2+1]-1; k++)
303 {
304 if (tab2[k]-1 == fac1) coeff[k] -= flux;
305 if (tab2[k]-1 == fac2) coeff[k] += flux;
306 }
307}
308
309void Op_Diff_VDF_Face_Axi_base::ajouter_contribution_elem(const DoubleTab& inco, Matrice_Morse& matrice) const
310{
311 if (inco.line_size() > 1) not_implemented(__func__);
312
313 const auto& tab1 = matrice.get_set_tab1();
314 const auto& tab2 = matrice.get_set_tab2();
315 auto& coeff = matrice.get_set_coeff();
316 for (int num_elem = 0; num_elem < le_dom_vdf->nb_elem(); num_elem++)
317 {
318 const int fx0 = elem_faces(num_elem,0), fx1 = elem_faces(num_elem,dimension), fy0 = elem_faces(num_elem,1), fy1 = elem_faces(num_elem,1+dimension);
319 // Calcul de tau11
320 const double tau11 = 1/(xv(fx1,0) - xv(fx0,0));
321
322 // Calcul de tau22
323 const double R = xp(num_elem,0);
324 double d_teta = xv(fy1,1) - xv(fy0,1);
325 if (d_teta < 0) d_teta += deux_pi;
326
327 double tau22 = 1/(R*d_teta);
328 // termes supplementaires en axi
329 tau22 += 0.5/R;
330 const double flux_X = tau11*nu_(num_elem)*0.5*(surface(fx0)+surface(fx1));
331 const double flux_Y = tau22*nu_(num_elem)*0.5*(surface(fy0)+surface(fy1));
332 fill_coeff_matrice_morse(fx0,fx1,flux_X,matrice);
333 fill_coeff_matrice_morse(fy0,fy1,flux_Y,matrice);
334
335 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
336 const double coef_laplacien_axi = +0.5*tau22*nu_(num_elem);
337
338 for (auto k = tab1[fx0]-1; k < tab1[fx0+1]-1; k++)
339 if (tab2[k]-1 == fx0) coeff[k] += coef_laplacien_axi*volumes_entrelaces(fx0)*porosite(fx0)/xv(fx0,0);
340
341 for (auto k=tab1[fx1]-1; k<tab1[fx1+1]-1; k++)
342 if (tab2[k]-1 == fx1) coeff[k] += coef_laplacien_axi*volumes_entrelaces(fx1)*porosite(fx1)/xv(fx1,0);
343 }
344}
345
346void Op_Diff_VDF_Face_Axi_base::ajouter_contribution_elem_3D(Matrice_Morse& matrice) const
347{
348 for (int num_elem = 0; num_elem < le_dom_vdf->nb_elem(); num_elem++)
349 {
350 const int fz0 = elem_faces(num_elem,2), fz1 = elem_faces(num_elem,2+dimension);
351 // Calcul de tau33
352 const double tau33 = 1/(xv(fz1,2) - xv(fz0,2));
353 const double flux_Z = tau33*nu_(num_elem)*0.5*(surface(fz0)+surface(fz1));
354 fill_coeff_matrice_morse(fz0,fz1,flux_Z,matrice);
355 }
356}
357
358void Op_Diff_VDF_Face_Axi_base::ajouter_contribution_aretes_bords(Matrice_Morse& matrice) const
359{
360 const auto& tab1 = matrice.get_set_tab1();
361 const auto& tab2 = matrice.get_set_tab2();
362 auto& coeff = matrice.get_set_coeff();
363 const int ndeb = le_dom_vdf->premiere_arete_bord(), nfin = ndeb + le_dom_vdf->nb_aretes_bord();
364 for (int n_arete=ndeb; n_arete<nfin; n_arete++)
365 {
366 const int n_type = type_arete_bord(n_arete-ndeb);
367 switch(n_type)
368 {
369 case TypeAreteBordVDF::PAROI_PAROI: /* fall through */
370 case TypeAreteBordVDF::FLUIDE_FLUIDE:
371 case TypeAreteBordVDF::PAROI_FLUIDE:
372 {
373 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), signe = Qdm(n_arete,3), ori1 = orientation(fac1), ori3 = orientation(fac3);
374 const double db_diffusivite = nu_mean_2_pts_(face_voisins(fac3,0),face_voisins(fac3,1));
375
376 if (ori1 == 0) // bord d'equation R = cte
377 {
378 double flux1;
379 if (ori3 == 1) // flux de tau12 a travers le bord
380 {
381 double dist3 = xv(fac3,0)-xv(fac1,0);
382 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
383 const double tau12 = 1/dist3;
384 flux1 = db_diffusivite*tau12*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
385 }
386 else //if (ori3 == 2) flux de tau13 a travers le bord
387 {
388 assert (ori3 == 2);
389 double dist3 = xv(fac3,0)-xv(fac1,0);
390 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
391 const double tau13 = 1/dist3;
392 flux1 = db_diffusivite*tau13*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
393 }
394
395 for (auto k = tab1[fac3]-1; k < tab1[fac3+1]-1; k++)
396 if (tab2[k]-1 == fac3) coeff[k] += signe*flux1;
397 }
398 else if (ori1 == 1) // bord d'equation teta = cte
399 {
400 const double R = xv(fac3,0);
401 double d_teta = xv(fac3,1) - xv(fac1,1);
402 if (d_teta < 0) d_teta += deux_pi;
403
404 double dist3 = R*d_teta;
405 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
406
407 if (ori3 == 0) // flux de tau21 a travers le bord
408 {
409 double tau21 = 1/dist3;
410
411 // Terme supplementaire en axi
412 tau21 -= 0.5/R;
413
414 const double flux2 = db_diffusivite*tau21*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
415
416 for (auto k = tab1[fac3]-1; k < tab1[fac3+1]-1; k++)
417 if (tab2[k]-1 == fac3) coeff[k] += signe*flux2;
418
419 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
420 const double coef_laplacien_axi = 0.5*db_diffusivite*tau21;
421
422 for (auto k = tab1[fac1]-1; k < tab1[fac1+1]-1; k++)
423 if (tab2[k]-1 == fac1) coeff[k] += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
424
425 for (auto k = tab1[fac2]-1; k < tab1[fac2+1]-1; k++)
426 if (tab2[k]-1 == fac2) coeff[k] += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
427 }
428 else // if (ori3 == 2) flux de tau23 a travers le bord
429 {
430 assert(ori3 == 2);
431 const double tau23 = 1/dist3;
432 const double flux3 = db_diffusivite*tau23*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
433 for (auto k = tab1[fac3]-1; k < tab1[fac3+1]-1; k++)
434 if (tab2[k]-1 == fac3) coeff[k] += signe*flux3;
435 }
436 }
437 else // (ori1 == 2) bord d'equation Z = cte
438 {
439 double flux4;
440 if (ori3 == 0) // flux de tau31 a travers le bord
441 {
442 double dist3 = xv(fac3,2)-xv(fac1,2);
443 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
444 const double tau31 = 1/dist3;
445 flux4 = db_diffusivite*tau31*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
446 }
447 else //if (ori3 == 1) flux de tau32 a travers le bord
448 {
449 assert(ori3 == 1);
450 double dist3 = xv(fac3,2)-xv(fac1,2);
451 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
452 const double tau32 = 1/dist3;
453 flux4 = db_diffusivite*tau32*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
454 }
455
456 for (auto k = tab1[fac3]-1; k < tab1[fac3+1]-1; k++)
457 if (tab2[k]-1 == fac3) coeff[k] += signe*flux4;
458 }
459 break;
460 }
461 case TypeAreteBordVDF::NAVIER_NAVIER: // pas de flux diffusif calcule
462 break;
463 default :
464 {
465 Cerr << "On a rencontre un type d'arete non prevu : [ num arete : " << n_arete << " ], [ type : " << n_type << " ]" << finl;
466 Process::exit();
467 break;
468 }
469 }
470 }
471}
472
473void Op_Diff_VDF_Face_Axi_base::ajouter_contribution_aretes_mixtes_internes(Matrice_Morse& matrice) const
474{
475 const auto& tab1 = matrice.get_set_tab1();
476 const auto& tab2 = matrice.get_set_tab2();
477 auto& coeff = matrice.get_set_coeff();
478 const int ndeb = le_dom_vdf->premiere_arete_mixte(), nfin = le_dom_vdf->nb_aretes();
479 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
480 {
481 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), fac4 = Qdm(n_arete,3), ori1 = orientation(fac1), ori3 = orientation(fac3);
482 const double db_diffusivite = nu_mean_4_pts_(fac3,fac4);
483
484 if (ori1 == 1) // (seule possibilite : ori3 =0) Arete XY
485 {
486 double flux1;
487 // Calcul de tau21
488 const double R = xv(fac3,0);
489 double d_teta = xv(fac4,1) - xv(fac3,1);
490 if (d_teta < 0) d_teta += deux_pi;
491 double tau21 = 1/(R*d_teta);
492 // Terme supplementaire en axi
493 tau21 -= 0.5/R;
494
495 // flux de tau21 sur la facette a cheval sur les faces fac1 et fac2
496 flux1 = db_diffusivite*tau21*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
497 fill_coeff_matrice_morse(fac3,fac4,flux1,matrice);
498
499 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
500 const double coef_laplacien_axi = 0.5*db_diffusivite*tau21;
501 for (auto k = tab1[fac1]-1; k < tab1[fac1+1]-1; k++)
502 if (tab2[k]-1 == fac1) coeff[k] += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
503
504 for (auto k = tab1[fac2]-1; k < tab1[fac2+1]-1; k++)
505 if (tab2[k]-1 == fac2) coeff[k] += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
506
507 // Calcul de tau12
508 const double tau12 = 1/(xv(fac2,0) - xv(fac1,0));
509
510 // flux de tau12 sur la facette a cheval sur les faces fac3 et fac4
511 flux1 = db_diffusivite*tau12*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
512 fill_coeff_matrice_morse(fac1,fac2,flux1,matrice);
513 }
514 else if (ori3 == 1) // (seule possibilite ori1 = 2) arete YZ
515 {
516 double flux2;
517 // Calcul de tau32
518 const double tau32 = 1/(xv(fac4,2) - xv(fac3,2));
519
520 // flux de tau32 sur la facette a cheval sur les faces fac1 et fac2
521 flux2 = db_diffusivite*tau32*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
522 fill_coeff_matrice_morse(fac3,fac4,flux2,matrice);
523
524 // Calcul de tau23
525 const double R = xv(fac1,0);
526 double d_teta = xv(fac2,1) - xv(fac1,1);
527 if (d_teta < 0) d_teta += deux_pi;
528 const double tau23 = 1/(R*d_teta);
529
530 // flux de tau23 sur la facette a cheval sur les faces fac3 et fac4
531 flux2 = db_diffusivite*tau23*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
532 fill_coeff_matrice_morse(fac1,fac2,flux2,matrice);
533 }
534 else // seule possibilite ori1 = 2 et ori3 = 0: arete XZ
535 {
536 double flux3;
537 // Calcul de tau31
538 const double tau31 = 1/(xv(fac4,2) - xv(fac3,2));
539
540 // flux de tau31 sur la facette a cheval sur les faces fac1 et fac2
541 flux3 = db_diffusivite*tau31*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
542 fill_coeff_matrice_morse(fac3,fac4,flux3,matrice);
543
544 // Calcul de tau13
545 const double tau13 = 1/(xv(fac2,0) - xv(fac1,0));
546
547 // flux de tau13 sur la facette a cheval sur les faces fac3 et fac4
548 flux3 = db_diffusivite*tau13*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
549 fill_coeff_matrice_morse(fac1,fac2,flux3,matrice);
550 }
551 }
552}
553
554void Op_Diff_VDF_Face_Axi_base::ajouter_contribution(const DoubleTab& inco, Matrice_Morse& matrice ) const
555{
556 ajouter_contribution_elem(inco,matrice); // Boucle sur les elements
557 if (dimension == 3) ajouter_contribution_elem_3D(matrice); // Boucle sur les elements supplementaires si 3D
558 ajouter_contribution_aretes_bords(matrice); // Boucle sur les aretes bord
559 ajouter_contribution_aretes_mixtes_internes(matrice); // Boucle sur les aretes mixtes et internes
560}
561
562void Op_Diff_VDF_Face_Axi_base::contribue_au_second_membre(DoubleTab& resu) const
563{
564 const int ndeb = le_dom_vdf->premiere_arete_bord(), nfin = ndeb + le_dom_vdf->nb_aretes_bord();
565 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
566 {
567 const int n_type = type_arete_bord(n_arete-ndeb);
568 switch(n_type)
569 {
570 case TypeAreteBordVDF::PAROI_PAROI: /* fall through */
571 case TypeAreteBordVDF::FLUIDE_FLUIDE:
572 case TypeAreteBordVDF::PAROI_FLUIDE:
573 {
574 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), signe = Qdm(n_arete,3);
575 const int ori1 = orientation(fac1), ori3 = orientation(fac3), rang1 = fac1 - le_dom_vdf->premiere_face_bord(), rang2 = fac2 - le_dom_vdf->premiere_face_bord();
576 double vit_imp, tps = inconnue->temps();
577
578 if (n_type == TypeAreteBordVDF::PAROI_FLUIDE) // arete paroi_fluide :il faut determiner qui est la face fluide
579 {
580 if (est_egal(inconnue->valeurs()(fac1), 0))
581 vit_imp = Champ_Face_get_val_imp_face_bord(tps, rang2, ori3, la_zcl_vdf.valeur());
582 else
583 vit_imp = Champ_Face_get_val_imp_face_bord(tps, rang1, ori3, la_zcl_vdf.valeur());
584 }
585 else vit_imp = 0.5*(Champ_Face_get_val_imp_face_bord(tps,rang1,ori3,la_zcl_vdf.valeur())+Champ_Face_get_val_imp_face_bord(tps,rang2,ori3,la_zcl_vdf.valeur()));
586
587 const double db_diffusivite = nu_mean_2_pts_(face_voisins(fac3,0),face_voisins(fac3,1));
588
589 if (ori1 == 0) // bord d'equation R = cte
590 {
591 double flux1;
592 if (ori3 == 1) // flux de tau12 a travers le bord
593 {
594 double dist3 = xv(fac3,0)-xv(fac1,0);
595 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
596 const double tau12 = (-vit_imp)/dist3;
597 flux1 = db_diffusivite*tau12*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
598 }
599 else// if (ori3 == 2) flux de tau13 a travers le bord
600 {
601 assert(ori3 == 2);
602 double dist3 = xv(fac3,0)-xv(fac1,0);
603 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
604 const double tau13 = (-vit_imp)/dist3;
605 flux1 = db_diffusivite*tau13*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
606 }
607 resu[fac3] += signe*flux1;
608 }
609 else if (ori1 == 1) // bord d'equation teta = cte
610 {
611 const double R = xv(fac3,0);
612 double d_teta = xv(fac3,1) - xv(fac1,1);
613 if (d_teta < 0) d_teta += deux_pi;
614 double dist3 = R*d_teta;
615 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
616 if (ori3 == 0) // flux de tau21 a travers le bord
617 {
618 // Terme supplementaire en axi
619 const double tau21 = (-vit_imp)/dist3;
620 const double flux2 = db_diffusivite*tau21*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
621 resu[fac3] += signe*flux2;
622
623 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
624 const double coef_laplacien_axi = 0.5*db_diffusivite*tau21;
625 resu(fac1) += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
626 resu(fac2) += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
627 }
628 else if (ori3 == 2) // flux de tau23 a travers le bord
629 {
630 const double tau23 = (-vit_imp)/dist3;
631 const double flux3 = db_diffusivite*tau23*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
632 resu[fac3] += signe*flux3;
633 }
634 }
635 else // (ori1 == 2) bord d'equation Z = cte
636 {
637 double flux4;
638 if (ori3 == 0) // flux de tau31 a travers le bord
639 {
640 double dist3 = xv(fac3,2) - xv(fac1,2);
641 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
642 const double tau31 = (-vit_imp)/dist3;
643 flux4 = db_diffusivite*tau31*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
644 }
645 else //if (ori3 == 1) flux de tau32 a travers le bord
646 {
647 assert(ori3 == 1);
648 double dist3 = xv(fac3,2)-xv(fac1,2);
649 if (n_type != TypeAreteBordVDF::PAROI_PAROI) dist3 *= 1;
650 const double tau32 = (-vit_imp)/dist3;
651 flux4 = db_diffusivite*tau32*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
652 }
653 resu[fac3] += signe*flux4;
654 }
655 break;
656 }
657 case TypeAreteBordVDF::NAVIER_NAVIER: // pas de flux diffusif calcule
658 break;
659 default :
660 {
661 Cerr << "On a rencontre un type d'arete non prevu : [ num arete : " << n_arete << " ], [ type : " << n_type << " ]" << finl;
662 Process::exit();
663 break;
664 }
665 }
666 }
667}
Champ_Face_VDF
class Champ_Face_VDF Cette classe sert a representer un champ vectoriel dont on ne calcule
Definition Champ_Face_VDF.h:42
Champ_Inc_base
Classe Champ_Inc_base.
Definition Champ_Inc_base.h:53
Domaine_Cl_VDF
class Domaine_Cl_VDF
Definition Domaine_Cl_VDF.h:63
Domaine_Cl_VDF::type_arete_bord
int type_arete_bord(int num_arete) const
Definition Domaine_Cl_VDF.h:76
Domaine_Cl_dis_base
classe Domaine_Cl_dis_base Les objets Domaine_Cl_dis_base representent les conditions aux limites
Definition Domaine_Cl_dis_base.h:37
Domaine_VDF
class Domaine_VDF
Definition Domaine_VDF.h:64
Domaine_VDF::orientation
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
Definition Domaine_VDF.h:297
Domaine_VDF::Qdm
int Qdm(int num_arete, int) const
Definition Domaine_VDF.h:231
Domaine_VF::face_surfaces
virtual const DoubleVect & face_surfaces() const
Definition Domaine_VF.h:51
Domaine_VF::volumes_entrelaces
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
Domaine_VF::xp
double xp(int num_elem, int k) const
Definition Domaine_VF.h:77
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_dis_base
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
Definition Domaine_dis_base.h:39
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Matrice_Morse
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
Definition Matrice_Morse.h:50
Matrice_Morse::get_set_tab2
auto & get_set_tab2()
Definition Matrice_Morse.h:103
Matrice_Morse::get_set_coeff
auto & get_set_coeff()
Definition Matrice_Morse.h:108
Matrice_Morse::get_set_tab1
auto & get_set_tab1()
Definition Matrice_Morse.h:98
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Op_Diff_VDF_Face_Axi_base
Definition Op_Diff_VDF_Face_Axi_base.h:25
Op_Diff_VDF_Face_Axi_base::deux_pi
static constexpr double deux_pi
Definition Op_Diff_VDF_Face_Axi_base.h:39
Op_Diff_VDF_Face_Axi_base::volumes_entrelaces
DoubleVect volumes_entrelaces
Definition Op_Diff_VDF_Face_Axi_base.h:45
Op_Diff_VDF_Face_Axi_base::xv
DoubleTab xv
Definition Op_Diff_VDF_Face_Axi_base.h:46
Op_Diff_VDF_Face_Axi_base::nu_mean_4_pts_
virtual double nu_mean_4_pts_(const int, const int) const =0
Op_Diff_VDF_Face_Axi_base::nu_
virtual double nu_(const int) const =0
Op_Diff_VDF_Face_Axi_base::face_voisins
IntTab face_voisins
Definition Op_Diff_VDF_Face_Axi_base.h:44
Op_Diff_VDF_Face_Axi_base::xp
DoubleTab xp
Definition Op_Diff_VDF_Face_Axi_base.h:46
Op_Diff_VDF_Face_Axi_base::Qdm
IntTab Qdm
Definition Op_Diff_VDF_Face_Axi_base.h:44
Op_Diff_VDF_Face_Axi_base::associer
void associer(const Domaine_dis_base &, const Domaine_Cl_dis_base &, const Champ_Inc_base &) override
Definition Op_Diff_VDF_Face_Axi_base.cpp:23
Op_Diff_VDF_Face_Axi_base::elem_faces
IntTab elem_faces
Definition Op_Diff_VDF_Face_Axi_base.h:44
Op_Diff_VDF_Face_Axi_base::calculer_dt_stab
double calculer_dt_stab() const override
Calcul dt_stab.
Definition Op_Diff_VDF_Face_Axi_base.cpp:43
Op_Diff_VDF_Face_Axi_base::orientation
IntVect orientation
Definition Op_Diff_VDF_Face_Axi_base.h:43
Op_Diff_VDF_Face_Axi_base::type_arete_bord
IntVect type_arete_bord
Definition Op_Diff_VDF_Face_Axi_base.h:43
Op_Diff_VDF_Face_Axi_base::calculer
DoubleTab & calculer(const DoubleTab &, DoubleTab &) const override
calcule la contribution de la diffusion, la range dans resu
Definition Op_Diff_VDF_Face_Axi_base.cpp:286
Op_Diff_VDF_Face_Axi_base::surface
DoubleVect surface
Definition Op_Diff_VDF_Face_Axi_base.h:45
Op_Diff_VDF_Face_Axi_base::nu_mean_2_pts_
virtual double nu_mean_2_pts_(const int, const int) const =0
Op_Diff_VDF_Face_Axi_base::porosite
DoubleVect porosite
Definition Op_Diff_VDF_Face_Axi_base.h:45
Op_Diff_VDF_Face_base
Definition Op_Diff_VDF_Face_base.h:27
Op_Diff_VDF_base::calculer_dt_stab_
double calculer_dt_stab_(const Domaine_VDF &zone_VDF) const
Definition Op_Diff_VDF_base.cpp:158
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67