TrioCFD 1.9.9_beta
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 // Compute tau11
55 const double tau11 = (inco[fx1]-inco[fx0])/(xv(fx1,0) - xv(fx0,0));
56 // Compute 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; // extra terms in axisymmetric
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 // Extra terms in the axisymmetric Laplacian: integrated as source terms
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 // Compute 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) // wall_fluid edge: we must determine which is the fluid face
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) // boundary with equation R = const
114 {
115 double flux1;
116 if (ori3 == 1) // flux of tau12 through the boundary
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 of tau13 through the boundary
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) // boundary with equation theta = const
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 of tau21 through the boundary
142 {
143 double tau21 = (inco[fac3]-vit_imp)/dist3;
144 tau21 -= 0.5*(inco[fac1]+inco[fac2])/R; // Extra term in 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 // Extra terms in the axisymmetric Laplacian
149 // They are integrated as source terms
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 of tau23 through the boundary
155 {
156 // XXX : beware if discrepancy: in the constant case it was tau23 = signe*(vit_imp-inco[fac3])/dist3 (normally the same but just in case)
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) boundary with equation Z = const
163 {
164 double flux4;
165 if (ori3 == 0) // flux of tau31 through the boundary
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 of tau32 through the boundary
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: // no diffusive flux computed
187 break;
188 default :
189 {
190 Cerr << "Encountered an unexpected edge type: [ edge num: " << 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) // (only possibility: ori3 =0) XY edge
207 {
208 double flux1;
209 // Compute 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 // Extra term in axi
216 tau21 -= 0.5*(inco[fac1]+inco[fac2])/R;
217
218 // flux of tau21 on the facet straddling faces fac1 and 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 // Extra terms in the axisymmetric Laplacian: integrated as source terms
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 // Compute tau12
229 const double tau12 = (inco(fac2)-inco(fac1))/(xv(fac2,0) - xv(fac1,0));
230 // flux of tau12 on the facet straddling faces fac3 and 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) // (only possibility: ori1 = 2) YZ edge
236 {
237 double flux2;
238 // Compute tau32
239 const double tau32 = (inco(fac4)-inco(fac3))/(xv(fac4,2) - xv(fac3,2));
240
241 // flux of tau32 on the facet straddling faces fac1 and 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 // Compute 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 of tau23 on the facet straddling faces fac3 and 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 // only possibility: ori1 = 2 and ori3 = 0: XZ edge
257 {
258 double flux3;
259 // Compute tau31
260 const double tau31 = (inco(fac4)-inco(fac3))/(xv(fac4,2) - xv(fac3,2));
261
262 // flux of tau31 on the facet straddling faces fac1 and 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 // Compute tau13
268 const double tau13 = (inco(fac2)-inco(fac1))/(xv(fac2,0) - xv(fac1,0));
269 // flux of tau13 on the facet straddling faces fac3 and 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); // Loop over elements
280 if (dimension == 3) ajouter_elem_3D(inco,resu); // Loop over extra elements if 3D
281 ajouter_aretes_bords(inco, resu); // Loop over boundary edges
282 ajouter_aretes_mixtes_internes(inco, resu); // Loop over mixed and internal edges
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 // Compute tau11
320 const double tau11 = 1/(xv(fx1,0) - xv(fx0,0));
321
322 // Compute 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 // extra terms in 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 // Extra terms in the axisymmetric Laplacian: integrated as source terms
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 // Compute 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) // boundary with equation R = const
377 {
378 double flux1;
379 if (ori3 == 1) // flux of tau12 through the boundary
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 of tau13 through the boundary
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) // boundary with equation theta = const
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 of tau21 through the boundary
408 {
409 double tau21 = 1/dist3;
410
411 // Extra term in 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 // Extra terms in the axisymmetric Laplacian: integrated as source terms
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 of tau23 through the boundary
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) boundary with equation Z = const
438 {
439 double flux4;
440 if (ori3 == 0) // flux of tau31 through the boundary
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 of tau32 through the boundary
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: // no diffusive flux computed
462 break;
463 default :
464 {
465 Cerr << "Encountered an unexpected edge type: [ edge num: " << 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) // (only possibility: ori3 =0) XY edge
485 {
486 double flux1;
487 // Compute 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 // Extra term in axi
493 tau21 -= 0.5/R;
494
495 // flux of tau21 on the facet straddling faces fac1 and 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 // Extra terms in the axisymmetric Laplacian: integrated as source terms
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 // Compute tau12
508 const double tau12 = 1/(xv(fac2,0) - xv(fac1,0));
509
510 // flux of tau12 on the facet straddling faces fac3 and 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) // (only possibility: ori1 = 2) YZ edge
515 {
516 double flux2;
517 // Compute tau32
518 const double tau32 = 1/(xv(fac4,2) - xv(fac3,2));
519
520 // flux of tau32 on the facet straddling faces fac1 and 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 // Compute 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 of tau23 on the facet straddling faces fac3 and 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 // only possibility: ori1 = 2 and ori3 = 0: XZ edge
535 {
536 double flux3;
537 // Compute tau31
538 const double tau31 = 1/(xv(fac4,2) - xv(fac3,2));
539
540 // flux of tau31 on the facet straddling faces fac1 and 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 // Compute tau13
545 const double tau13 = 1/(xv(fac2,0) - xv(fac1,0));
546
547 // flux of tau13 on the facet straddling faces fac3 and 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); // Loop over elements
557 if (dimension == 3) ajouter_contribution_elem_3D(matrice); // Loop over extra elements if 3D
558 ajouter_contribution_aretes_bords(matrice); // Loop over boundary edges
559 ajouter_contribution_aretes_mixtes_internes(matrice); // Loop over mixed and internal edges
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) // wall_fluid edge: we must determine which is the fluid face
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) // boundary with equation R = const
590 {
591 double flux1;
592 if (ori3 == 1) // flux of tau12 through the boundary
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) // boundary with equation theta = const
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 of tau21 through the boundary
617 {
618 // Extra term in 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 // Extra terms in the axisymmetric Laplacian: integrated as source terms
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 of tau23 through the boundary
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) boundary with equation Z = const
636 {
637 double flux4;
638 if (ori3 == 0) // flux of tau31 through the boundary
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: // no diffusive flux computed
658 break;
659 default :
660 {
661 Cerr << "Encountered an unexpected edge type: [ edge num: " << n_arete << " ], [ type: " << n_type << " ]" << finl;
662 Process::exit();
663 break;
664 }
665 }
666 }
667}
Champ_Face_VDF
class Champ_Face_VDF
Definition Champ_Face_VDF.h:43
Champ_Inc_base
Class 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
class Domaine_Cl_dis_base Domaine_Cl_dis_base objects represent discretized boundary conditions
Definition Domaine_Cl_dis_base.h:37
Domaine_VDF
class Domaine_VDF
Definition Domaine_VDF.h:61
Domaine_VDF::orientation
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
Definition Domaine_VDF.h:291
Domaine_VDF::Qdm
int Qdm(int num_arete, int) const
Definition Domaine_VDF.h:225
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
Returns the index of the i-th face of element num_elem; the face numbering convention is.
Definition Domaine_VF.h:542
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
Returns the neighbouring element of num_face in direction i.
Definition Domaine_VF.h:418
Domaine_dis_base
class Domaine_dis_base This class is the base of the hierarchy of discretized domains.
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
Matrice_Morse class - Represents a (sparse) matrix M, not necessarily square,.
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: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::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
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
Computes 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
Computes the diffusion contribution and stores it in 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)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Sortie
Base class for output streams.
Definition Sortie.h:52
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67