TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
Op_Diff_VEF_Face.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_VEF_Face.h>
17#include <Champ_P1NC.h>
18#include <Champ_Uniforme.h>
19#include <Periodique.h>
20#include <Symetrie.h>
21#include <Neumann_homogene.h>
22#include <Neumann_paroi.h>
23#include <Echange_externe_impose.h>
24#include <Echange_externe_radiatif.h>
25#include <Neumann_sortie_libre.h>
26#include <Milieu_base.h>
27#include <TRUSTTrav.h>
28#include <Probleme_base.h>
29#include <Navier_Stokes_std.h>
30#include <Porosites_champ.h>
31#include <Device.h>
32#include <Echange_couplage_thermique.h>
33#include <Champ_front_calc_interne.h>
34#include <Robin_VEF.h>
35
36Implemente_instanciable_sans_constructeur(Op_Diff_VEF_Face,"Op_Diff_VEF_P1NC",Op_Diff_VEF_base);
37
42
43Sortie& Op_Diff_VEF_Face::printOn(Sortie& s ) const
44{
45 return s << que_suis_je() ;
46}
47
48Entree& Op_Diff_VEF_Face::readOn(Entree& s )
49{
50 return s ;
51}
52
53/*! @brief Associate the diffusivity field.
54 *
55 */
56void Op_Diff_VEF_Face::associer_diffusivite(const Champ_base& diffu)
57{
58 diffusivite_ = diffu;
59}
60
65
66const Champ_base& Op_Diff_VEF_Face::diffusivite() const
67{
68 return diffusivite_.valeur();
69}
70
71int ma_func_qui_renvoie_int()
72{
73 return 1;
74}
75
76void Op_Diff_VEF_Face::ajouter_cas_scalaire(const DoubleTab& tab_inconnue,
77 DoubleTab& tab_resu, DoubleTab& tab_flux_bords,
78 DoubleTab& tab_nu,
79 const Domaine_Cl_VEF& domaine_Cl_VEF,
80 const Domaine_VEF& domaine_VEF ) const
81{
82 int nb_faces = domaine_VEF.nb_faces();
83 int nb_faces_elem = domaine_VEF.domaine().nb_faces_elem();
84 int nb_bords=domaine_VEF.nb_front_Cl();
85 const int premiere_face_int=domaine_VEF.premiere_face_int();
86
87 {
88 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
89 CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
90 CDoubleTabView face_normale = domaine_VEF.face_normales().view_ro();
91 CDoubleArrView inverse_volumes = domaine_VEF.inverse_volumes().view_ro();
92 CDoubleTabView nu = tab_nu.view_ro();
93 CDoubleTabView inconnue = tab_inconnue.view_ro();
94 DoubleArrView flux_bords = static_cast<ArrOfDouble&>(tab_flux_bords).view_rw();
95 DoubleArrView resu = static_cast<ArrOfDouble&>(tab_resu).view_rw();
96 // Process boundary faces
97 for (int n_bord = 0; n_bord < nb_bords; n_bord++)
98 {
99 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
100 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
101 int num1 = 0;
102 int num2 = le_bord.nb_faces_tot();
103 int nb_faces_bord_reel = le_bord.nb_faces();
104 CIntArrView le_bord_num_face = le_bord.num_face().view_ro();
105 if (sub_type(Periodique, la_cl.valeur()))
106 {
107 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
108 CIntArrView face_associee = la_cl_perio.face_associee().view_ro();
109 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
110 Kokkos::RangePolicy<>(num1, nb_faces_bord_reel), KOKKOS_LAMBDA(
111 const int ind_face)
112 {
113 int num_face = le_bord_num_face(ind_face);
114 int fac_asso = face_associee(ind_face);
115 fac_asso = le_bord_num_face(fac_asso);
116 for (int kk = 0; kk < 2; kk++)
117 {
118 int elem = face_voisins(num_face, kk);
119 for (int i = 0; i < nb_faces_elem; i++)
120 {
121 int j = elem_faces(elem, i);
122 if (j > num_face && j != fac_asso)
123 {
124 double valA = viscA(num_face, j, elem, nu(elem, 0), face_voisins, face_normale,
125 inverse_volumes);
126 double flux = valA * (inconnue(j, 0) - inconnue(num_face, 0));
127 Kokkos::atomic_add(&resu(num_face), +flux);
128 if (j < nb_faces) // face reelle
129 Kokkos::atomic_add(&resu(j), -0.5 * flux);
130 }
131 }
132 }
133 });
134 end_gpu_timer(__KERNEL_NAME__);
135 }
136 else // There is only one component, so we handle
137 // a scalar equation (not velocity): no need to use
138 // the tangential tau (wall thermal laws do not compute
139 // turbulent exchange at the wall for now)
140 {
141 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
142 Kokkos::RangePolicy<>(num1, num2), KOKKOS_LAMBDA(
143 const int ind_face)
144 {
145 int num_face = le_bord_num_face(ind_face);
146 int elem = face_voisins(num_face, 0);
147 for (int i = 0; i < nb_faces_elem; i++)
148 {
149 int j = elem_faces(elem, i);
150 if (j > num_face || num_face >= nb_faces)
151 {
152 double valA = viscA(num_face, j, elem, nu(elem, 0), face_voisins, face_normale,
153 inverse_volumes);
154 double flux = valA * (inconnue(j, 0) - inconnue(num_face, 0));
155 if (num_face < nb_faces) // face reelle
156 {
157 Kokkos::atomic_add(&resu(num_face), +flux);
158 Kokkos::atomic_add(&flux_bords(num_face), -flux);
159 }
160 if (j < nb_faces) // face reelle
161 {
162 Kokkos::atomic_add(&resu(j), -flux);
163 if (j < premiere_face_int)
164 Kokkos::atomic_add(&flux_bords(j), +flux);
165 }
166 }
167 }
168 });
169 end_gpu_timer(__KERNEL_NAME__);
170 }
171 }
172
173 // Internal faces:
174 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
175 Kokkos::MDRangePolicy<Kokkos::Rank<2>>({premiere_face_int, 0}, {nb_faces, 2}),
176 KOKKOS_LAMBDA(const int num_face, const int k)
177 {
178 int elem = face_voisins(num_face, k);
179 for (int i = 0; i < nb_faces_elem; i++)
180 {
181 int j = elem_faces(elem, i);
182 if (j > num_face)
183 {
184 int contrib = 1;
185
186 if (j >= nb_faces) // This is a virtual face
187 {
188 int el1 = face_voisins(j, 0);
189 int el2 = face_voisins(j, 1);
190 if ((el1 == -1) || (el2 == -1))
191 contrib = 0;
192 }
193
194 if (contrib)
195 {
196 double valA = viscA(num_face, j, elem, nu(elem, 0), face_voisins,
197 face_normale, inverse_volumes);
198 double flux = valA * (inconnue(j, 0) - inconnue(num_face, 0));
199 Kokkos::atomic_add(&resu(num_face), flux);
200 if (j < nb_faces) // Process real faces
201 Kokkos::atomic_add(&resu(j), -flux);
202 }
203 }
204 }
205 });
206 end_gpu_timer(__KERNEL_NAME__);
207 }
208
209 // Neumann :
210 for (int n_bord=0; n_bord<nb_bords; n_bord++)
211 {
212 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
213 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
214 int ndeb = le_bord.num_premiere_face();
215 int nfin = ndeb + le_bord.nb_faces();
216 if (sub_type(Neumann_paroi,la_cl.valeur()))
217 {
218 const Neumann_paroi& la_cl_paroi = ref_cast(Neumann_paroi, la_cl.valeur());
219 CDoubleArrView surface = domaine_VEF.face_surfaces().view_ro();
220 CDoubleTabView flux_impose = la_cl_paroi.flux_impose().view_ro();
221 DoubleArrView flux_bords = static_cast<ArrOfDouble&>(tab_flux_bords).view_rw();
222 DoubleArrView resu = static_cast<ArrOfDouble&>(tab_resu).view_rw();
223 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int face)
224 {
225 double flux = flux_impose(face-ndeb, 0) * surface(face);
226 resu(face) += flux;
227 flux_bords(face) = flux;
228 });
229 end_gpu_timer(__KERNEL_NAME__);
230 }
231 else if (sub_type(Echange_externe_impose,la_cl.valeur()))
232 {
233 const Echange_externe_impose& la_cl_paroi = ref_cast(Echange_externe_impose, la_cl.valeur());
234 const double coeff = COEFF_STEFAN_BOLTZMANN;
235 const bool has_emissivity = la_cl_paroi.has_emissivite();
236 CDoubleArrView surface = domaine_VEF.face_surfaces().view_ro();
237 CDoubleArrView text = static_cast<const ArrOfDouble&>(la_cl_paroi.tab_T_ext()).view_ro();
238 CDoubleArrView himp = static_cast<const ArrOfDouble&>(la_cl_paroi.tab_h_imp()).view_ro();
239 CDoubleArrView eps;
240 if (has_emissivity) eps = static_cast<const ArrOfDouble&>(la_cl_paroi.tab_emissivite()).view_ro();
241 CDoubleArrView inconnue = static_cast<const ArrOfDouble&>(tab_inconnue).view_ro();
242 DoubleArrView resu = static_cast<ArrOfDouble&>(tab_resu).view_rw();
243 DoubleArrView flux_bords = static_cast<ArrOfDouble&>(tab_flux_bords).view_wo();
244 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA (const int face)
245 {
246 int ind_face = face - ndeb;
247 double flux = himp(ind_face)*(text(ind_face)-inconnue(face))*surface(face);
248 resu[face] += flux;
249 flux_bords(face) = flux;
250
251 if (has_emissivity)
252 {
253 double T = inconnue(face);
254 double t_ext = text(ind_face);
255 flux = coeff * eps(ind_face) * (t_ext * t_ext * t_ext * t_ext - T * T * T * T) * surface(face);
256 resu[face] += flux;
257 flux_bords(face) += flux;
258 }
259 });
260 end_gpu_timer(__KERNEL_NAME__);
261 }
262 else if (sub_type(Echange_couplage_thermique, la_cl.valeur()))
263 {
264 ToDo_Kokkos("critical");
265 const Echange_couplage_thermique& la_cl_paroi = ref_cast(Echange_couplage_thermique, la_cl.valeur());
266 const DoubleVect& surface = domaine_VEF.face_surfaces();
267 for (int face=ndeb; face<nfin; face++)
268 {
269 double h=la_cl_paroi.h_imp(face-ndeb);
270 double Text=la_cl_paroi.T_ext(face-ndeb);
271 double phiext=la_cl_paroi.flux_exterieur_impose(face-ndeb);
272 double flux=(phiext+h*(Text-tab_inconnue(face)))*surface(face);
273 tab_resu[face] += flux;
274 tab_flux_bords(face) = flux;
275 }
276 }
277 else if (sub_type(Neumann_homogene,la_cl.valeur())
278 || sub_type(Symetrie,la_cl.valeur())
279 || sub_type(Neumann_sortie_libre,la_cl.valeur()))
280 {
281 DoubleArrView flux_bords = static_cast<ArrOfDouble&>(tab_flux_bords).view_wo();
282 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int face)
283 {
284 flux_bords(face) = 0.;
285 });
286 end_gpu_timer(__KERNEL_NAME__);
287 }
288 }
289}
290
291void Op_Diff_VEF_Face::ajouter_cas_vectoriel(const DoubleTab& inconnue,
292 DoubleTab& resu, DoubleTab& tab_flux_bords,
293 DoubleTab& nu,
294 const Domaine_Cl_VEF& domaine_Cl_VEF,
295 const Domaine_VEF& domaine_VEF,
296 int nb_comp) const
297{
298 assert(nb_comp==dimension);
299
300 // Build grad_ array if necessary
301 if(!grad_.get_md_vector())
302 {
303 grad_.resize(0, Objet_U::dimension, Objet_U::dimension);
304 domaine_VEF.domaine().creer_tableau_elements(grad_);
305 }
306 Champ_P1NC::calcul_gradient(inconnue,grad_,domaine_Cl_VEF);
307
308 /* ToDo OpenMP: refactor with Op_Dift_VEF_Face.cpp into a template class
309 if (le_modele_turbulence->utiliser_loi_paroi())
310 {
311 Champ_P1NC::calcul_duidxj_paroi(grad_,nu,nu_turb,tau_tan_,domaine_Cl_VEF);
312 grad_.echange_espace_virtuel(); // gradient_elem up to date on virtual elements
313 }
314 DoubleTab Re;
315 Re.resize(0, Objet_U::dimension, Objet_U::dimension);
316 domaine_VEF.domaine().creer_tableau_elements(Re);
317 Re = 0.;
318 if (le_modele_turbulence->calcul_tenseur_Re(nu_turb, grad_, Re))
319 {
320 Cerr << "Using nonlinear turbulent diffusion in NS" << finl;
321 for (int elem=0; elem<nb_elem; elem++)
322 for (int i=0; i<nbr_comp; i++)
323 for (int j=0; j<nbr_comp; j++)
324 Re(elem,i,j) *= nu_turb[elem];
325 }
326 else
327 {
328 for (int elem=0; elem<nb_elem; elem++)
329 for (int i=0; i<nbr_comp; i++)
330 for (int j=0; j<nbr_comp; j++)
331 Re(elem,i,j) = nu_turb[elem]*(grad_(elem,i,j) + grad_(elem,j,i));
332 }
333 Re.echange_espace_virtuel();
334 */
335
336 int nb_faces = domaine_VEF.nb_faces();
337 int nb_faces_bord = domaine_VEF.premiere_face_int();
338 CIntTabView face_voisins_v = domaine_VEF.face_voisins().view_ro();
339 CDoubleTabView face_normales_v = domaine_VEF.face_normales().view_ro();
340 CDoubleTabView nu_v = nu.view_ro();
341 CDoubleTabView3 grad_v = grad_.view_ro<3>();
342 DoubleTabView resu_v = resu.view_rw();
343 DoubleTabView tab_flux_bords_v = tab_flux_bords.view_rw();
344
345 auto kern_ajouter = KOKKOS_LAMBDA(int
346 num_face, int k)
347 {
348 int elem = face_voisins_v(num_face, k);
349 if (elem >= 0)
350 {
351 int ori = 1 - 2 * k;
352 double nu_elem = nu_v(elem, 0);
353 for (int i = 0; i < nb_comp; i++)
354 for (int j = 0; j < nb_comp; j++)
355 {
356 double grad_ij = grad_v(elem, i, j);
357 double grad_ji = grad_v(elem, j, i);
358 double fn = face_normales_v(num_face, j);
359 double flux = ori * fn * (nu_elem * grad_ij /* + Re(elem, i, j) */ );
360 Kokkos::atomic_sub(&resu_v(num_face, i), flux);
361
362 if (num_face < nb_faces_bord)
363 {
364 double flux_bord = ori * fn * (nu_elem * (grad_ij + grad_ji));
365 Kokkos::atomic_sub(&tab_flux_bords_v(num_face, i), flux_bord);
366 }
367 }
368 }
369 };
370 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0,0}, {nb_faces,2}) , kern_ajouter);
371 end_gpu_timer(__KERNEL_NAME__);
372
373
374 const int nb_bords=domaine_VEF.nb_front_Cl();
375 for (int n_bord=0; n_bord<nb_bords; n_bord++)
376 {
377 // Update flux_bords on symmetry:
378 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
379 if (sub_type(Symetrie,la_cl.valeur()))
380 {
381 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
382 int ndeb = le_bord.num_premiere_face();
383 int nfin = ndeb + le_bord.nb_faces();
384 DoubleTabView flux_bords = tab_flux_bords.view_wo();
385 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int face)
386 {
387 flux_bords(face, 0) = 0.;
388 });
389 end_gpu_timer(__KERNEL_NAME__);
390 }
391
392 else if (sub_type(Robin_VEF, la_cl.valeur()))
393 {
394#ifdef TRUST_USE_GPU
395 Cerr << "Warning not tested on GPU" << finl;
396#endif
397 ToDo_Kokkos("critical");
398
399 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
400 const Robin_VEF& la_cl_robin = ref_cast(Robin_VEF,la_cl.valeur());
401 int marq = phi_psi_diffuse(equation());
402 const DoubleVect& porosite_face = equation().milieu().porosite_face();
403 double scale_factor_is_one = 1.;
404 double inv_alpha = 1./la_cl_robin.get_alpha_cl() ;
405 double inv_beta = 1./la_cl_robin.get_beta_cl();
406 double inv_alpha_minus_inv_beta = 1./la_cl_robin.get_alpha_cl() - 1./la_cl_robin.get_beta_cl();
407 DoubleTab normal_vector;
408 int ndeb = le_bord.num_premiere_face();
409 int nfin = ndeb +le_bord.nb_faces();
410 for (int face=ndeb; face<nfin; face++)
411 {
412 int id_face_bord = face -ndeb;
413 double face_surface = domaine_VEF.face_surfaces(face);
414 normal_vector = domaine_VEF.normalized_boundaries_outward_vector(face, scale_factor_is_one);
415
416 for (int nc1=0; nc1<nb_comp; nc1++)
417 {
418 double flux_robin_uu = 0. ;
419 double flux_robin_rhs = 0.;
420 double flux_tot ;
421
422 // forward term for velocity
423 for (int nc2 = 0; nc2<nb_comp; nc2++)
424 {
425 const double normal2 = normal_vector(nc1)*normal_vector(nc2);
426 flux_robin_uu += (inv_beta*(nc1==nc2) + inv_alpha_minus_inv_beta*normal2)* (face_surface);
427 }
428 flux_robin_uu *= inconnue(face,nc1);
429
430 // rhs for robin bc
431 double val;
432 if (dimension == 2)
433 {
434 // add normal component rhs
435 val = inv_alpha * normal_vector(nc1) * la_cl_robin.flux_normal_imp(id_face_bord);
436
437 // add tangential component rhs
438 double tgte = (2*nc1-1)*normal_vector(1-nc1);
439 val += inv_beta * la_cl_robin.flux_tangentiel_imp(id_face_bord, 0)*tgte;
440 }
441 else
442 {
443 // add normal component rhs
444 val = inv_alpha * normal_vector(nc1) * la_cl_robin.flux_normal_imp(id_face_bord);
445
446 // add tangential component rhs
447 val += inv_beta * la_cl_robin.flux_tangentiel_imp(id_face_bord, nc1) ;
448 }
449 flux_robin_rhs = val*face_surface;
450 flux_tot = (flux_robin_rhs - flux_robin_uu)* (marq ? porosite_face(face) : 1);
451 resu(face,nc1) += flux_tot;
452 tab_flux_bords(face,nc1) += flux_tot;
453
454
455 }
456 }
457 }
458 }
459}
460
461void Op_Diff_VEF_Face::ajouter_cas_multi_scalaire(const DoubleTab& inconnue,
462 DoubleTab& resu, DoubleTab& tab_flux_bords,
463 DoubleTab& nu,
464 const Domaine_Cl_VEF& domaine_Cl_VEF,
465 const Domaine_VEF& domaine_VEF,
466 int nb_comp) const
467{
468 ToDo_Kokkos("critical");
469 const IntTab& elemfaces = domaine_VEF.elem_faces();
470 const IntTab& face_voisins = domaine_VEF.face_voisins();
471 int i0,j,num_face;
472 int nb_faces = domaine_VEF.nb_faces();
473 int nb_faces_elem = domaine_VEF.domaine().nb_faces_elem();
474 int n_bord;
475 double flux0;
476 //DoubleVect n(Objet_U::dimension);
477 //DoubleTrav Tgrad(Objet_U::dimension,Objet_U::dimension);
478
479 assert(nb_comp>1);
480 int nb_bords=domaine_VEF.nb_front_Cl();
481 int ind_face;
482
483 for (n_bord=0; n_bord<nb_bords; n_bord++)
484 {
485 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
486 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
487 // const IntTab& elemfaces = domaine_VEF.elem_faces();
488 int num1 = 0;
489 int num2 = le_bord.nb_faces_tot();
490 int nb_faces_bord_reel = le_bord.nb_faces();
491
492 if (sub_type(Periodique,la_cl.valeur()))
493 {
494 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
495 int fac_asso;
496 for (ind_face=num1; ind_face<nb_faces_bord_reel; ind_face++)
497 {
498 fac_asso = la_cl_perio.face_associee(ind_face);
499 fac_asso = le_bord.num_face(fac_asso);
500 num_face = le_bord.num_face(ind_face);
501 for (int kk=0; kk<2; kk++)
502 {
503 int elem = face_voisins(num_face, kk);
504 for (i0=0; i0<nb_faces_elem; i0++)
505 {
506 if ( ( (j= elemfaces(elem,i0)) > num_face ) && (j != fac_asso ) )
507 {
508 for (int nc=0; nc<nb_comp; nc++)
509 {
510 double valA = viscA(num_face,j,elem,nu(elem,nc));
511 resu(num_face,nc)+=valA*inconnue(j,nc);
512 resu(num_face,nc)-=valA*inconnue(num_face,nc);
513 if(j<nb_faces) // face reelle
514 {
515 ////WARNING: NUM_face differs from the reference version
516 resu(j,nc)+=0.5*valA*inconnue(num_face,nc);
517 resu(j,nc)-=0.5*valA*inconnue(j,nc);
518 }
519 }
520 }
521 }
522 }
523 }
524 }// end if periodic
525 else
526 {
527 for (ind_face=num1; ind_face<num2; ind_face++)
528 {
529 num_face = le_bord.num_face(ind_face);
530 int elem=face_voisins(num_face,0);
531
532 // Loop over faces:
533 for (int i=0; i<nb_faces_elem; i++)
534 if (( (j= elemfaces(elem,i)) > num_face ) || (ind_face>=nb_faces_bord_reel))
535 {
536 for (int nc=0; nc<nb_comp; nc++)
537 {
538 double valA = viscA(num_face,j,elem,nu(elem,nc));
539 if (ind_face<nb_faces_bord_reel)
540 {
541 double flux=valA*(inconnue(j,nc)-inconnue(num_face,nc));
542 resu(num_face,nc)+=flux;
543 tab_flux_bords(num_face,nc)-=flux;
544 }
545
546 if(j<nb_faces) // face reelle
547 {
548 resu(j,nc)+=valA*inconnue(num_face,nc);
549 resu(j,nc)-=valA*inconnue(j,nc);
550 }
551 }
552 }
553 }
554 }
555 }// End for n_bord
556
557 // Process internal faces
558
559 for (num_face=domaine_VEF.premiere_face_int(); num_face<nb_faces; num_face++)
560 {
561 for (int k=0; k<2; k++)
562 {
563 int elem = face_voisins(num_face,k);
564 for (i0=0; i0<nb_faces_elem; i0++)
565 {
566 if ( (j= elemfaces(elem,i0)) > num_face )
567 {
568 int el1,el2;
569 int contrib=1;
570 if(j>=nb_faces) // This is a virtual face
571 {
572 el1 = face_voisins(j,0);
573 el2 = face_voisins(j,1);
574 if((el1==-1)||(el2==-1))
575 contrib=0;
576 }
577 if(contrib)
578 {
579 for (int nc=0; nc<nb_comp; nc++)
580 {
581 double valA = viscA(num_face,j,elem,nu(elem,nc));
582 resu(num_face,nc)+=valA*inconnue(j,nc);
583 resu(num_face,nc)-=valA*inconnue(num_face,nc);
584 if(j<nb_faces) // Process real faces only
585 {
586 resu(j,nc)+=valA*inconnue(num_face,nc);
587 resu(j,nc)-=valA*inconnue(j,nc);
588 }
589 else
590 {
591 // Face j is virtual
592 }
593 }
594 }
595 }
596 }
597 }
598 }// End internal faces
599
600
601 //Based on what is done for the scalar case
602 for (n_bord=0; n_bord<nb_bords; n_bord++)
603 {
604 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
605
606 if (sub_type(Neumann_paroi,la_cl.valeur()))
607 {
608 const Neumann_paroi& la_cl_paroi = ref_cast(Neumann_paroi, la_cl.valeur());
609 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
610 int ndeb = le_bord.num_premiere_face();
611 int nfin = ndeb + le_bord.nb_faces();
612 for (int face=ndeb; face<nfin; face++)
613 {
614 for (int nc=0; nc<nb_comp; nc++)
615 {
616 flux0=la_cl_paroi.flux_impose(face-ndeb,nc)*domaine_VEF.surface(face);
617 resu(face,nc) += flux0;
618 tab_flux_bords(face,nc) = flux0;
619 }
620 }
621 }
622 else if (sub_type(Echange_externe_impose,la_cl.valeur()))
623 {
624 throw;
625 const Echange_externe_impose& la_cl_paroi = ref_cast(Echange_externe_impose, la_cl.valeur());
626 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
627 int ndeb = le_bord.num_premiere_face();
628 int nfin = ndeb + le_bord.nb_faces();
629 for (int face=ndeb; face<nfin; face++)
630 {
631 for (int nc=0; nc<nb_comp; nc++)
632 {
633 flux0=la_cl_paroi.h_imp(face-ndeb,nc)*(la_cl_paroi.T_ext(face-ndeb,nc)-inconnue(face,nc))*domaine_VEF.surface(face);
634 resu(face,nc) += flux0;
635 tab_flux_bords(face,nc) = flux0;
636
637 if (la_cl_paroi.has_emissivite())
638 {
639 const double text = la_cl_paroi.T_ext(face - ndeb, nc), T = inconnue(face, nc);
640 flux0 = COEFF_STEFAN_BOLTZMANN * la_cl_paroi.emissivite(face - ndeb, nc) * (text * text * text * text - T * T * T * T) * domaine_VEF.face_surfaces(face);
641 resu(face, nc) += flux0;
642 tab_flux_bords(face, nc) += flux0;
643 }
644 }
645 }
646 }
647 else if (sub_type(Echange_couplage_thermique,la_cl.valeur()))
648 {
649 const Echange_couplage_thermique& la_cl_paroi = ref_cast(Echange_couplage_thermique, la_cl.valeur());
650 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
651 int ndeb = le_bord.num_premiere_face();
652 int nfin = ndeb + le_bord.nb_faces();
653 for (int face=ndeb; face<nfin; face++)
654 {
655 for (int nc=0; nc<nb_comp; nc++)
656 {
657 double phiext = la_cl_paroi.flux_exterieur_impose(face-ndeb,nc);
658 flux0 = (phiext + la_cl_paroi.h_imp(face-ndeb,nc)*(la_cl_paroi.T_ext(face-ndeb,nc)-inconnue(face,nc)))*domaine_VEF.surface(face);
659 resu(face,nc) += flux0;
660 tab_flux_bords(face,nc) = flux0;
661 }
662 }
663 }
664 else if (sub_type(Neumann_homogene,la_cl.valeur())
665 || sub_type(Symetrie,la_cl.valeur())
666 || sub_type(Neumann_sortie_libre,la_cl.valeur()))
667 {
668 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
669 int ndeb = le_bord.num_premiere_face();
670 int nfin = ndeb + le_bord.nb_faces();
671 for (int face=ndeb; face<nfin; face++)
672 for (int nc=0; nc<nb_comp; nc++)
673 tab_flux_bords(face,nc) = 0.;
674 }
675 }
676}
677
678
679DoubleTab& Op_Diff_VEF_Face::ajouter(const DoubleTab& inconnue_org, DoubleTab& resu) const
680{
681 remplir_nu(nu_);
682 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
683 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
684
685 int nb_comp = 1;
686 int nb_dim = resu.nb_dim();
687 if(nb_dim==2)
688 nb_comp=resu.dimension(1);
689 DoubleTab nu;
690 DoubleTab tab_inconnue;
691 int marq=phi_psi_diffuse(equation());
692 const DoubleVect& porosite_face = equation().milieu().porosite_face();
693 const DoubleVect& porosite_elem = equation().milieu().porosite_elem();
694 // either div(phi nu grad inco)
695 // or div(nu grad phi inco)
696 // depending on whether phi_psi or psi is diffused
697 modif_par_porosite_si_flag(nu_,nu,!marq,porosite_elem);
698 const DoubleTab& inconnue=modif_par_porosite_si_flag(inconnue_org,tab_inconnue,marq,porosite_face);
699
700 const Champ_base& inco = equation().inconnue();
701 const Nature_du_champ nature_champ = inco.nature_du_champ();
702
703 // Size and initialize the flux balance array:
704 if (flux_bords_.size_array()!=domaine_VEF.nb_faces_bord()) flux_bords_.resize(domaine_VEF.nb_faces_bord(),nature_champ==scalaire ? 1 : nb_comp);
705 flux_bords_=0.;
706
707 if(nature_champ==scalaire)
708 ajouter_cas_scalaire(inconnue, resu, flux_bords_, nu, domaine_Cl_VEF, domaine_VEF);
709 else if (nature_champ==vectoriel)
710 ajouter_cas_vectoriel(inconnue, resu, flux_bords_, nu, domaine_Cl_VEF, domaine_VEF,nb_comp);
711 else if (nature_champ==multi_scalaire)
712 ajouter_cas_multi_scalaire(inconnue, resu, flux_bords_, nu, domaine_Cl_VEF, domaine_VEF,nb_comp);
713 modifier_flux(*this);
714
715 return resu;
716}
717
718DoubleTab& Op_Diff_VEF_Face::calculer(const DoubleTab& inconnue, DoubleTab& resu) const
719{
720 resu = 0;
721 return ajouter(inconnue,resu);
722}
723
724void Op_Diff_VEF_Face::ajouter_contribution(const DoubleTab& tab_transporte, Matrice_Morse& tab_matrice) const
725{
726 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
727 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
728
729 modifier_matrice_pour_periodique_avant_contribuer(tab_matrice,equation());
730
731 // Fill the nu array because matrix assembly with ajouter_contribution
732 // may be performed before the first time step
733 remplir_nu(nu_);
734 DoubleTrav tab_nu;
735
736 // either div(phi nu grad inco)
737 // or div(nu grad phi inco)
738 // depending on whether phi_psi or psi is diffused
739 int marq = phi_psi_diffuse(equation());
740 modif_par_porosite_si_flag(nu_,tab_nu,!marq,equation().milieu().porosite_elem());
741
742 int nb_dim = tab_transporte.nb_dim();
743 int nb_comp = (nb_dim==2 ? tab_transporte.dimension(1) : 1);
744 int nb_faces_elem = domaine_VEF.domaine().nb_faces_elem();
745 int nb_bords = domaine_VEF.nb_front_Cl();
746
747
748 IntTrav tab_face_associee(domaine_VEF.premiere_face_int());
749 IntTrav tab_fac2b_idx(domaine_VEF.nb_faces());
750
751
752 // Retrieve the indices of periodic boundary faces and the
753 // associated faces into arrays beforehand to
754 // allow structured access in the kernel afterwards
755 for (int n_bord = 0; n_bord < nb_bords; n_bord++)
756 {
757 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
758 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
759 int num1 = le_bord.num_premiere_face();
760
761 if (sub_type(Periodique, la_cl.valeur()))
762 {
763 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
764
765 int nb_faces = le_bord.nb_faces();
766 int num2b = num1 + nb_faces / 2;
767
768 // only iterate over half the periodic faces
769 // the result will be copied to the associated face at the end...
770 ToDo_Kokkos("critical");
771 for (int fac = num1; fac < num2b; fac++)
772 {
773 int fac_asso = la_cl_perio.face_associee(fac - num1) + num1;
774
775 tab_face_associee(fac) = fac_asso;
776 tab_fac2b_idx(fac) = fac;
777 tab_fac2b_idx(fac+nb_faces/2) = fac;
778 }
779 }
780 }
781
782 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
783 CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
784 CDoubleArrView porosite_face = equation().milieu().porosite_face().view_ro();
785 CDoubleArrView inverse_volumes = domaine_VEF.inverse_volumes().view_ro();
786 CDoubleTabView face_normale = domaine_VEF.face_normales().view_ro();
787 CDoubleTabView nu = tab_nu.view_ro();
788 CIntArrView est_face_bord = domaine_VEF.est_face_bord().view_ro();
789 CIntArrView face_associee = static_cast<ArrOfInt&>(tab_face_associee).view_ro();
790 CIntArrView fac2b_idx = static_cast<ArrOfInt&>(tab_fac2b_idx).view_ro();
791 Matrice_Morse_View matrice;
792 matrice.set(tab_matrice);
793
794 auto ajouter_contrib = KOKKOS_LAMBDA(const int fac)
795 {
796 int type_face = est_face_bord(fac);
797 int fac2b = fac2b_idx(fac);
798
799 if (type_face == 2 && fac == fac2b) // faces perio
800 {
801 int elem1 = face_voisins(fac,0);
802 int elem2 = face_voisins(fac,1);
803
804 int fac_asso = face_associee(fac);
805 for (int i = 0; i < nb_faces_elem; i++)
806 {
807 int j = elem_faces(elem1,i);
808 if (j > fac)
809 {
810 double val = viscA(fac,j,elem1,nu(elem1,0), face_voisins, face_normale, inverse_volumes);
811 double coeff_face1 = val * (marq ? porosite_face(fac) : 1);
812 double coeff_face2 = val * (marq ? porosite_face(j) : 1);
813
814 for (int nc = 0; nc < nb_comp; nc++)
815 {
816 int n0 = fac*nb_comp + nc;
817 int j0 = j*nb_comp + nc;
818
819 matrice.atomic_add(n0, n0, +coeff_face1);
820 matrice.atomic_add(n0, j0, -coeff_face2);
821 matrice.atomic_add(j0, n0, -coeff_face1);
822 matrice.atomic_add(j0, j0, +coeff_face2);
823 }
824 }
825 if (elem2 != -1)
826 {
827 j = elem_faces(elem2,i);
828 if (j > fac)
829 {
830 double val = viscA(fac,j,elem2,nu(elem2,0), face_voisins, face_normale, inverse_volumes);
831 double coeff_face1 = val * (marq ? porosite_face(fac) : 1);
832 double coeff_face2 = val * (marq ? porosite_face(j) : 1);
833
834 for (int nc = 0; nc < nb_comp; nc++)
835 {
836 int n0 = fac*nb_comp + nc;
837 int j0 = j*nb_comp + nc;
838 int n1 = fac_asso*nb_comp+nc;
839
840 matrice.atomic_add(n0, n0, +coeff_face1);
841 matrice.atomic_add(n0, j0, -coeff_face2);
842 matrice.atomic_add(j0, n1, -coeff_face1);
843 matrice.atomic_add(j0, j0, +coeff_face2);
844 }
845 }
846 }
847 }
848 }
849 else if (type_face == 1) // non-periodic boundary faces
850 {
851 int elem1 = face_voisins(fac,0);
852 for (int i = 0; i < nb_faces_elem; i++)
853 {
854 int j = elem_faces(elem1,i);
855 if (j > fac)
856 {
857 double val = viscA(fac,j,elem1,nu(elem1,0), face_voisins, face_normale, inverse_volumes);
858 double coeff_face1 = val * (marq ? porosite_face(fac) : 1);
859 double coeff_face2 = val * (marq ? porosite_face(j) : 1);
860
861 for (int nc = 0; nc < nb_comp; nc++)
862 {
863 int n0 = fac*nb_comp + nc;
864 int j0 = j*nb_comp + nc;
865
866 matrice.atomic_add(n0, n0, +coeff_face1);
867 matrice.atomic_add(n0, j0, -coeff_face2);
868 matrice.atomic_add(j0, n0, -coeff_face1);
869 matrice.atomic_add(j0, j0, +coeff_face2);
870 }
871 }
872 }
873 }
874 else if (type_face == 0) // faces internes
875 {
876 for (int k=0; k<2; k++)
877 {
878 int elem = face_voisins(fac,k);
879 if (elem!=-1)
880 {
881 for (int i = 0; i < nb_faces_elem; i++)
882 {
883 int j = elem_faces(elem, i);
884 if (j > fac)
885 {
886 double val = viscA(fac, j, elem, nu(elem, 0), face_voisins, face_normale, inverse_volumes);
887 double coeff_face1 = val * (marq ? porosite_face(fac) : 1);
888 double coeff_face2 = val * (marq ? porosite_face(j) : 1);
889
890 for (int nc = 0; nc < nb_comp; nc++)
891 {
892 int n0 = fac * nb_comp + nc;
893 int j0 = j * nb_comp + nc;
894
895 matrice.atomic_add(n0, n0, +coeff_face1);
896 matrice.atomic_add(n0, j0, -coeff_face2);
897 matrice.atomic_add(j0, n0, -coeff_face1);
898 matrice.atomic_add(j0, j0, +coeff_face2);
899 }
900 }
901 }
902 }
903 }
904 }
905 };
906 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),domaine_VEF.nb_faces(), ajouter_contrib);
907 end_gpu_timer(__KERNEL_NAME__);
908
909 int premiere_face_int = domaine_VEF.premiere_face_int();
910 DoubleTrav tab_h_impose(premiere_face_int);
911 DoubleTrav tab_derivee_flux_exterieur_imposee(premiere_face_int);
912
913 // Neumann: fill arrays with boundary conditions for the Kokkos kernel
914 for (int n_bord = 0; n_bord < nb_bords; n_bord++)
915 {
916 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
917
918 if (sub_type(Echange_externe_impose,la_cl.valeur()))
919 {
920 const Echange_externe_impose& la_cl_paroi = ref_cast(Echange_externe_impose, la_cl.valeur());
921 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
922 int ndeb = le_bord.num_premiere_face();
923 int nfin = ndeb + le_bord.nb_faces();
924 const double coeff = COEFF_STEFAN_BOLTZMANN;
925 const bool has_emissivity = la_cl_paroi.has_emissivite();
926 CDoubleArrView inconnue = static_cast<const ArrOfDouble&>(equation().inconnue().valeurs()).view_ro();
927 CDoubleArrView himp = static_cast<const ArrOfDouble&>(la_cl_paroi.tab_h_imp()).view_ro();
928 CDoubleArrView eps;
929 if (has_emissivity) eps = static_cast<const ArrOfDouble&>(la_cl_paroi.tab_emissivite()).view_ro();
930 DoubleArrView h_impose = static_cast<ArrOfDouble&>(tab_h_impose).view_wo();
931 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int face)
932 {
933 int ind_face = face - ndeb;
934 h_impose(face) = himp(ind_face);
935 if (has_emissivity)
936 {
937 const double T = inconnue(face);
938 h_impose(face) = 4 * coeff * eps(ind_face) * T * T * T;
939 }
940 });
941 end_gpu_timer(__KERNEL_NAME__);
942 }
943 else if (sub_type(Echange_couplage_thermique, la_cl.valeur()))
944 {
945 const Echange_couplage_thermique& la_cl_paroi = ref_cast(Echange_couplage_thermique, la_cl.valeur());
946 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
947
948 int ndeb = le_bord.num_premiere_face();
949 int nfin = ndeb + le_bord.nb_faces();
950 ToDo_Kokkos("critical");
951 for (int face = ndeb; face < nfin; face++)
952 {
953 tab_h_impose(face) = la_cl_paroi.h_imp(face-ndeb);
954 tab_derivee_flux_exterieur_imposee(face) = la_cl_paroi.derivee_flux_exterieur_imposee(face-ndeb);
955 }
956 }
957 else if (sub_type(Robin_VEF, la_cl.valeur()))
958 {
959 ToDo_Kokkos("critical");
960 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
961 const Robin_VEF& la_cl_robin = ref_cast(Robin_VEF,la_cl.valeur());
962 double inv_alpha_minus_inv_beta = 1./la_cl_robin.get_alpha_cl() - 1./la_cl_robin.get_beta_cl();
963 double inv_beta = 1./la_cl_robin.get_beta_cl();
964 double scale_factor_is_one = 1. ;
965 DoubleTab normal_vector ;
966 int ndeb = le_bord.num_premiere_face();
967 int nfin = ndeb + le_bord.nb_faces();
968 const DoubleVect& tab_porosite_face = equation().milieu().porosite_face();
969 for (int face = ndeb; face < nfin; face++)
970 {
971 double face_surface = domaine_VEF.face_surfaces(face);
972 normal_vector = domaine_VEF.normalized_boundaries_outward_vector(face, scale_factor_is_one);
973 //int elem = face_voisins(face, 0) ;
974 for (int nc1 = 0; nc1 < nb_comp; nc1++)
975 {
976 const int i = face * nb_comp + nc1;
977
978 // diagonal term
979 double val = (inv_beta + inv_alpha_minus_inv_beta*normal_vector(nc1)*normal_vector(nc1))*face_surface;
980 double robin_contribution= val * (marq ? tab_porosite_face(face) : 1) ;
981 tab_matrice(i,i) += robin_contribution ;
982
983 // extradiagonal term
984 for (int nc2 = 0; nc2<nc1; nc2++)
985 {
986 const int j = face * nb_comp + nc2;
987 const double normal2 = normal_vector(nc1)*normal_vector(nc2);
988 val = inv_alpha_minus_inv_beta*normal2* (face_surface);
989 robin_contribution= val * (marq ? tab_porosite_face(face) : 1) ;
990 tab_matrice(i, j) += robin_contribution;
991 tab_matrice(j ,i) += robin_contribution;
992 }
993 }
994 }
995 }
996 }
997
998 CDoubleArrView h_impose = static_cast<const ArrOfDouble&>(tab_h_impose).view_ro();
999 CDoubleArrView derivee_flux_exterieur_imposee = static_cast<const ArrOfDouble&>(tab_derivee_flux_exterieur_imposee).view_ro();
1000 CDoubleArrView face_surfaces = domaine_VEF.face_surfaces().view_ro();
1001
1002 // Neumann: compute contributions on boundaries
1003 auto neumann = KOKKOS_LAMBDA (const int face)
1004 {
1005 double h = h_impose(face);
1006 double dphi_dT = derivee_flux_exterieur_imposee(face);
1007 matrice.add(face,face, + (h + dphi_dT) * face_surfaces(face));
1008 };
1009
1010 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), premiere_face_int, neumann);
1011 end_gpu_timer(__KERNEL_NAME__);
1012
1013 modifier_matrice_pour_periodique_apres_contribuer(tab_matrice,equation());
1014}
1015
1016void Op_Diff_VEF_Face::ajouter_contribution_multi_scalaire(const DoubleTab& tab_transporte, Matrice_Morse& tab_matrice) const
1017{
1018 modifier_matrice_pour_periodique_avant_contribuer(tab_matrice, equation());
1019
1020 // Fill the nu array because matrix assembly with ajouter_contribution
1021 // may be performed before the first time step
1022 remplir_nu(nu_);
1023
1024 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1025 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1026 const IntTab& tab_elem_faces = domaine_VEF.elem_faces();
1027 const IntTab& tab_face_voisins = domaine_VEF.face_voisins();
1028 const ArrOfInt& tab_est_face_bord = domaine_VEF.est_face_bord();
1029
1030 int nb_dim = tab_transporte.nb_dim();
1031 int nb_comp = (nb_dim == 2 ? tab_transporte.dimension(1) : 1);
1032
1033 DoubleTab tab_nu;
1034 int marq = phi_psi_diffuse(equation());
1035 const DoubleVect& porosite_elem = equation().milieu().porosite_elem();
1036
1037 // either div(phi nu grad inco)
1038 // or div(nu grad phi inco)
1039 // depending on whether phi_psi or psi is diffused
1040 modif_par_porosite_si_flag(nu_, tab_nu, !marq, porosite_elem);
1041 DoubleVect tab_porosite_eventuelle(equation().milieu().porosite_face());
1042 if (!marq)
1043 tab_porosite_eventuelle = 1;
1044
1045 int nb_faces_elem = domaine_VEF.domaine().nb_faces_elem();
1046 int nb_bords = domaine_VEF.nb_front_Cl();
1047 int nb_faces_tot = domaine_VEF.nb_faces_tot();
1048
1049 IntVect tab_face_associee(domaine_VEF.premiere_face_int());
1050 IntVect tab_fac2b_idx(domaine_VEF.nb_faces_tot());
1051
1052 // Retrieve the indices of periodic boundary faces and the
1053 // associated faces into arrays beforehand for structured
1054 // access in the kernel afterwards
1055 for (int n_bord = 0; n_bord < nb_bords; n_bord++)
1056 {
1057 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1058 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1059 int num1 = le_bord.num_premiere_face();
1060
1061 if (sub_type(Periodique, la_cl.valeur()))
1062 {
1063 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
1064
1065 int nb_faces = le_bord.nb_faces();
1066 int num2b = num1 + nb_faces / 2;
1067
1068 // only iterate over half the periodic faces
1069 // the result will be copied to the associated face at the end...
1070 ToDo_Kokkos("critical");
1071 for (int fac = num1; fac < num2b; fac++)
1072 {
1073 int fac_asso = la_cl_perio.face_associee(fac - num1) + num1;
1074
1075 tab_face_associee(fac) = fac_asso;
1076 tab_fac2b_idx(fac) = fac;
1077 tab_fac2b_idx(fac+nb_faces/2) = fac;
1078 }
1079 }
1080 }
1081
1082 CIntArrView est_face_bord = tab_est_face_bord.view_ro();
1083 CIntTabView face_voisins = tab_face_voisins.view_ro();
1084 CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
1085 CDoubleArrView inverse_volumes = domaine_VEF.inverse_volumes().view_ro();
1086 CIntTabView elem_faces = tab_elem_faces.view_ro();
1087
1088 CDoubleTabView nu = tab_nu.view_ro();
1089 CDoubleArrView porosite_eventuelle = tab_porosite_eventuelle.view_ro();
1090
1091 CIntArrView fac2b_idx = tab_fac2b_idx.view_ro();
1092 CIntArrView face_associee = tab_face_associee.view_ro();
1093
1094 Matrice_Morse_View matrice;
1095 matrice.set(tab_matrice);
1096
1097 auto kern_elem_faces = KOKKOS_LAMBDA (const int fac)
1098 {
1099 int type_face = est_face_bord(fac);
1100 int fac2b = fac2b_idx(fac);
1101
1102 // Periodic boundary faces
1103 if (type_face == 2 && fac == fac2b)
1104 {
1105 int fac_asso = face_associee(fac);
1106
1107 int elem1 = face_voisins(fac, 0);
1108 int elem2 = face_voisins(fac, 1);
1109
1110 for (int i = 0; i < nb_faces_elem; i++)
1111 {
1112 int j = elem_faces(elem1, i);
1113 if (j > fac)
1114 {
1115 for (int nc = 0; nc < nb_comp; nc++)
1116 {
1117 double val = viscA(fac, j, elem1, nu(elem1, nc), face_voisins, face_normales, inverse_volumes);
1118
1119 int n0 = fac * nb_comp + nc;
1120 int j0 = j * nb_comp + nc;
1121
1122 matrice.atomic_add(n0, n0, + val * porosite_eventuelle(fac));
1123 matrice.atomic_add(n0, j0, - val * porosite_eventuelle(j));
1124 matrice.atomic_add(j0, n0, - val * porosite_eventuelle(fac));
1125 matrice.atomic_add(j0, j0, + val * porosite_eventuelle(j));
1126 }
1127 }
1128 if (elem2 != -1)
1129 {
1130 j = elem_faces(elem2, i);
1131 if (j > fac)
1132 {
1133 for (int nc = 0; nc < nb_comp; nc++)
1134 {
1135 double val = viscA(fac, j, elem2, nu(elem1, nc), face_voisins, face_normales, inverse_volumes);
1136
1137 int n0 = fac * nb_comp + nc;
1138 int j0 = j * nb_comp + nc;
1139 int n0perio = fac_asso * nb_comp + nc;
1140
1141 matrice.atomic_add(n0, n0, + val * porosite_eventuelle(fac));
1142 matrice.atomic_add(n0, j0, - val * porosite_eventuelle(j));
1143 matrice.atomic_add(j0, n0perio, - val * porosite_eventuelle(fac));
1144 matrice.atomic_add(j0, j0, + val * porosite_eventuelle(j));
1145 }
1146 }
1147 }
1148 }
1149 }
1150 // Non-periodic boundary faces
1151 else if (type_face == 1)
1152 {
1153 int elem1 = face_voisins(fac, 0);
1154
1155 for (int i = 0; i < nb_faces_elem; i++)
1156 {
1157 int j = elem_faces(elem1, i);
1158 if (j > fac)
1159 {
1160 for (int nc = 0; nc < nb_comp; nc++)
1161 {
1162 double val = viscA(fac, j, elem1, nu(elem1, nc), face_voisins, face_normales, inverse_volumes);
1163
1164 int n0 = fac * nb_comp + nc;
1165 int j0 = j * nb_comp + nc;
1166
1167 matrice.atomic_add(n0, n0, + val * porosite_eventuelle(fac));
1168 matrice.atomic_add(n0, j0, - val * porosite_eventuelle(j));
1169 matrice.atomic_add(j0, n0, - val * porosite_eventuelle(fac));
1170 matrice.atomic_add(j0, j0, + val * porosite_eventuelle(j));
1171 }
1172 }
1173 }
1174 }
1175 // Faces internes
1176 else if (type_face == 0)
1177 {
1178 int elem1 = face_voisins(fac, 0);
1179 int elem2 = face_voisins(fac, 1);
1180
1181 for (int i = 0; i < nb_faces_elem; i++)
1182 {
1183 int j = elem_faces(elem1, i);
1184 if (j > fac)
1185 {
1186 for (int nc = 0; nc < nb_comp; nc++)
1187 {
1188 double val = viscA(fac, j, elem1, nu(elem1, nc), face_voisins, face_normales, inverse_volumes);
1189
1190 int n0 = fac * nb_comp + nc;
1191 int j0 = j * nb_comp + nc;
1192
1193 matrice.atomic_add(n0, n0, + val * porosite_eventuelle(fac));
1194 matrice.atomic_add(n0, j0, - val * porosite_eventuelle(j));
1195 matrice.atomic_add(j0, n0, - val * porosite_eventuelle(fac));
1196 matrice.atomic_add(j0, j0, + val * porosite_eventuelle(j));
1197 }
1198 }
1199
1200 if (elem2 != -1)
1201 {
1202 j = elem_faces(elem2, i);
1203 if (j > fac)
1204 {
1205 for (int nc = 0; nc < nb_comp; nc++)
1206 {
1207 double val = viscA(fac, j, elem2, nu(elem2, nc), face_voisins, face_normales, inverse_volumes);
1208 int n0 = fac * nb_comp + nc;
1209 int j0 = j * nb_comp + nc;
1210
1211 matrice.atomic_add(n0, n0, + val * porosite_eventuelle(fac));
1212 matrice.atomic_add(n0, j0, - val * porosite_eventuelle(j));
1213 matrice.atomic_add(j0, n0, - val * porosite_eventuelle(fac));
1214 matrice.atomic_add(j0, j0, + val * porosite_eventuelle(j));
1215 }
1216 }
1217 }
1218 }
1219 }
1220 };
1221
1222 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_faces_tot, kern_elem_faces);
1223 end_gpu_timer(__KERNEL_NAME__);
1224
1225 for (int n_bord = 0; n_bord < nb_bords; n_bord++)
1226 {
1227 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1228 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1229
1230 if (sub_type(Echange_externe_impose, la_cl.valeur()))
1231 {
1232 const Echange_externe_impose& la_cl_paroi = ref_cast(Echange_externe_impose, la_cl.valeur());
1233 int ndeb = le_bord.num_premiere_face();
1234 int nfin = ndeb + le_bord.nb_faces();
1235 ToDo_Kokkos("critical");
1236 for (int face = ndeb; face < nfin; face++)
1237 for (int nc = 0; nc < nb_comp; nc++)
1238 {
1239 const int i = face * nb_comp + nc;
1240 tab_matrice(i, i) += la_cl_paroi.h_imp(face - ndeb, nc) * domaine_VEF.surface(face);
1241 }
1242 }
1243 if (sub_type(Echange_externe_radiatif, la_cl.valeur()))
1244 {
1245 throw;
1246 }
1247 }
1248
1249 modifier_matrice_pour_periodique_apres_contribuer(tab_matrice, equation());
1250}
1251
1252void Op_Diff_VEF_Face::contribue_au_second_membre(DoubleTab& resu ) const
1253{
1254 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1255 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1256 int nb_comp = 1;
1257 int nb_dim = resu.nb_dim();
1258
1259 int nb_bords=domaine_VEF.nb_front_Cl();
1260
1261 if(nb_dim==2)
1262 nb_comp=resu.dimension(1);
1263
1264 // Partie imposee :
1265
1266 if (nb_dim == 1)
1267 {
1268 for (int n_bord=0; n_bord<nb_bords; n_bord++)
1269 {
1270 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1271
1272 if (sub_type(Neumann_paroi,la_cl.valeur()))
1273 {
1274 const Neumann_paroi& la_cl_paroi = ref_cast(Neumann_paroi, la_cl.valeur());
1275 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1276 int ndeb = le_bord.num_premiere_face();
1277 int nfin = ndeb + le_bord.nb_faces();
1278 ToDo_Kokkos("critical");
1279 for (int face=ndeb; face<nfin; face++)
1280 resu[face] += la_cl_paroi.flux_impose(face-ndeb)*domaine_VEF.surface(face);
1281 }
1282 else if (sub_type(Echange_externe_impose,la_cl.valeur()))
1283 {
1284 Cerr << "Non code pour Echange_externe_impose" << finl;
1285 assert(0);
1286 }
1287 else if (sub_type(Echange_externe_radiatif,la_cl.valeur()))
1288 {
1289 Cerr << "Non code pour Echange_externe_radiatif" << finl;
1290 assert(0);
1291 }
1292
1293 }
1294 }
1295 else
1296 {
1297 for (int n_bord=0; n_bord<nb_bords; n_bord++)
1298 {
1299 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1300
1301 if (sub_type(Neumann_paroi,la_cl.valeur()))
1302 {
1303 const Neumann_paroi& la_cl_paroi = ref_cast(Neumann_paroi, la_cl.valeur());
1304 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1305 int ndeb = le_bord.num_premiere_face();
1306 int nfin = ndeb + le_bord.nb_faces();
1307 ToDo_Kokkos("critical");
1308 for (int face=ndeb; face<nfin; face++)
1309 for (int comp=0; comp<nb_comp; comp++)
1310 resu(face,comp) += la_cl_paroi.flux_impose(face-ndeb,comp)*domaine_VEF.surface(face);
1311 }
1312 }
1313 }
1314}
1315
1316void Op_Diff_VEF_Face::verifier() const
1317{
1318 static int testee=0;
1319 if(testee)
1320 return;
1321 testee=1;
1322 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1323 // const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1324 // const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
1325 const DoubleVect& volumes_entrelaces = domaine_VEF.volumes_entrelaces();
1326
1327 const DoubleTab& xv=domaine_VEF.xv();
1328 DoubleTab vit(equation().inconnue().valeurs());
1329 DoubleTab resu(vit);
1330 int i, comp;
1331 if(dimension==2)
1332 {
1333 const int nbf = vit.dimension(0);
1334 Cerr << " Verification of delta(x,0) " << finl;
1335 for(i=0; i<nbf; i++)
1336 {
1337 vit(i,0)=xv(i,0);
1338 vit(i,1)=0;
1339 }
1340 calculer(vit, resu);
1341 for(i=0; i<nbf; i++)
1342 for(comp=0; comp<dimension; comp++)
1343 resu(i,comp)/=(volumes_entrelaces(i));
1344 for(i=0; i<nbf; i++)
1345 {
1346 if(std::fabs(resu(i,0))>1.e-10)
1347 {
1348 Cerr << " delta(x,0) ("<<i<<") = "
1349 << resu(i,0);
1350 Cerr << finl;
1351 }
1352 }
1353 Cerr << " Verification of delta(y(1-y),0) " << finl;
1354 for(i=0; i<nbf; i++)
1355 {
1356 vit(i,0)=xv(i,1)*(1-xv(i,1));
1357 vit(i,1)=0;
1358 }
1359 calculer(vit, resu);
1360 for(i=0; i<nbf; i++)
1361 for(comp=0; comp<dimension; comp++)
1362 resu(i,comp)/=(volumes_entrelaces(i));
1363 for(i=0; i<nbf; i++)
1364 {
1365 if(std::fabs(2-resu(i,0))>1.e-10)
1366 {
1367 Cerr << " delta(y(1-y),0) ("<<i<<") = "
1368 << resu(i,0);
1369 Cerr << finl;
1370 }
1371 }
1372 }
1373}
1374
1375
1376
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Returns the array of field values at the current time.
Definition Champ_Inc_base.h:77
Champ_P1NC::calcul_gradient
static DoubleTab & calcul_gradient(const DoubleTab &, DoubleTab &, const Domaine_Cl_VEF &)
Definition Champ_P1NC.cpp:912
Champ_base
class Champ_base This class is the base of the fields hierarchy.
Definition Champ_base.h:43
Cond_lim
class Cond_lim Generic class used to represent any class
Definition Cond_lim.h:31
Domaine_32_64::creer_tableau_elements
virtual void creer_tableau_elements(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
Creates a parallel array of values at elements.
Definition Domaine.cpp:850
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Returns the number of faces of type i of the geometric elements that make up the domain.
Definition Domaine.h:484
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:38
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Returns the i-th boundary condition.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_VF::face_surfaces
virtual const DoubleVect & face_surfaces() const
Definition Domaine_VF.h:51
Domaine_VF::nb_faces
int nb_faces() const
Returns the total number of faces.
Definition Domaine_VF.h:471
Domaine_VF::volumes_entrelaces
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
Domaine_VF::nb_faces_tot
int nb_faces_tot() const
Returns the total number of faces.
Definition Domaine_VF.h:481
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::surface
virtual double surface(int i) const
Definition Domaine_VF.h:53
Domaine_VF::normalized_boundaries_outward_vector
DoubleTab normalized_boundaries_outward_vector(int global_face_number, double scale_factor) const
Compute the normalized boundary outward vector associated to the face global_face_number and eventual...
Definition Domaine_VF.cpp:715
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::est_face_bord
ArrOfInt & est_face_bord()
Definition Domaine_VF.h:85
Domaine_VF::premiere_face_int
int premiere_face_int() const
A face is internal if and only if it separates two elements.
Definition Domaine_VF.h:463
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_VF::nb_faces_bord
int nb_faces_bord() const
Returns the number of faces on which boundary conditions are applied:
Definition Domaine_VF.h:512
Domaine_VF::inverse_volumes
double inverse_volumes(int i) const
Definition Domaine_VF.h:114
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Echange_couplage_thermique
: class Echange_couplage_thermique
Definition Echange_couplage_thermique.h:26
Echange_couplage_thermique::flux_exterieur_impose
double flux_exterieur_impose(int i) const override
Definition Echange_couplage_thermique.h:40
Echange_externe_impose
Classe Echange_externe_impose: This class represents the special case of the class.
Definition Echange_externe_impose.h:46
Echange_externe_radiatif
Echange_externe_radiatif: Combines radiative (sigma * eps * (T^4 - T_ext^4)) and convective (h * (T -...
Definition Echange_externe_radiatif.h:31
Echange_global_impose::derivee_flux_exterieur_imposee
virtual double derivee_flux_exterieur_imposee(int i) const
Definition Echange_global_impose.cpp:244
Echange_impose_base::h_imp
virtual double h_imp(int num) const
Returns the value of the imposed heat exchange coefficient on the i-th component.
Definition Echange_impose_base.cpp:112
Echange_impose_base::has_emissivite
bool has_emissivite() const
Definition Echange_impose_base.h:45
Echange_impose_base::T_ext
virtual double T_ext(int num) const
Returns the value of the imposed temperature on the i-th component of the boundary field.
Definition Echange_impose_base.cpp:76
Echange_impose_base::tab_h_imp
const DoubleTab & tab_h_imp(double temps=DMAXFLOAT) const
Definition Echange_impose_base.cpp:284
Echange_impose_base::emissivite
double emissivite(int num) const
Returns the value of the imposed emissivity on the i-th component.
Definition Echange_impose_base.cpp:150
Echange_impose_base::tab_T_ext
const DoubleTab & tab_T_ext(double temps=DMAXFLOAT) const
Definition Echange_impose_base.cpp:258
Echange_impose_base::tab_emissivite
const DoubleTab & tab_emissivite(double temps=DMAXFLOAT) const
Definition Echange_impose_base.cpp:310
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Equation_base::inconnue
virtual const Champ_Inc_base & inconnue() const =0
Field_base::nature_du_champ
virtual Nature_du_champ nature_du_champ() const
Definition Field_base.h:77
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
Front_VF::nb_faces_tot
int nb_faces_tot() const
Definition Front_VF.h:58
Front_VF::num_face
int num_face(const int) const
Definition Front_VF.h:68
Matrice_Morse
Matrice_Morse class - Represents a (sparse) matrix M, not necessarily square,.
Definition Matrice_Morse.h:50
Milieu_base::porosite_elem
DoubleVect & porosite_elem()
Definition Milieu_base.h:58
Milieu_base::porosite_face
DoubleVect & porosite_face()
Definition Milieu_base.h:62
MorEqn::equation
const Equation_base & equation() const
Returns the reference to the equation pointed to by MorEqn::mon_equation.
Definition MorEqn.h:62
Neumann_homogene
Classe Neumann_homogene This class is the base class of the hierarchy of homogeneous Neumann-type bou...
Definition Neumann_homogene.h:31
Neumann_paroi
Classe Neumann_paroi This boundary condition corresponds to an imposed flux for the.
Definition Neumann_paroi.h:29
Neumann_sortie_libre
Neumann_sortie_libre This class represents an open boundary without imposed velocity.
Definition Neumann_sortie_libre.h:34
Neumann::flux_impose
virtual double flux_impose(int i) const
Returns the value of the imposed flux on the i-th component of the field representing the flux at the...
Definition Neumann.cpp:35
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_VEF_Face
class Op_Diff_VEF_Face
Definition Op_Diff_VEF_Face.h:32
Op_Diff_VEF_Face::diffusivite
const Champ_base & diffusivite() const override
Definition Op_Diff_VEF_Face.cpp:66
Op_Diff_VEF_Face::Op_Diff_VEF_Face
Op_Diff_VEF_Face()
Definition Op_Diff_VEF_Face.cpp:38
Op_Diff_VEF_Face::ajouter
DoubleTab & ajouter(const DoubleTab &, DoubleTab &) const override
Definition Op_Diff_VEF_Face.cpp:679
Op_Diff_VEF_Face::contribue_au_second_membre
void contribue_au_second_membre(DoubleTab &) const
Definition Op_Diff_VEF_Face.cpp:1252
Op_Diff_VEF_Face::ajouter_cas_scalaire
void ajouter_cas_scalaire(const DoubleTab &inconnue, DoubleTab &resu, DoubleTab &flux_bords, DoubleTab &nu, const Domaine_Cl_VEF &domaine_Cl_VEF, const Domaine_VEF &domaine_VEF) const
Definition Op_Diff_VEF_Face.cpp:76
Op_Diff_VEF_Face::calculer
DoubleTab & calculer(const DoubleTab &, DoubleTab &) const override
Definition Op_Diff_VEF_Face.cpp:718
Op_Diff_VEF_Face::completer
void completer() override
Associates the operator with the domaine_dis, the domaine_Cl_dis, and the unknown of its equation.
Definition Op_Diff_VEF_Face.cpp:61
Op_Diff_VEF_Face::ajouter_cas_vectoriel
void ajouter_cas_vectoriel(const DoubleTab &inconnue, DoubleTab &resu, DoubleTab &flux_bords, DoubleTab &nu, const Domaine_Cl_VEF &domaine_Cl_VEF, const Domaine_VEF &domaine_VEF, int nb_comp) const
Definition Op_Diff_VEF_Face.cpp:291
Op_Diff_VEF_Face::grad_
DoubleTab grad_
Definition Op_Diff_VEF_Face.h:78
Op_Diff_VEF_Face::verifier
void verifier() const
Definition Op_Diff_VEF_Face.cpp:1316
Op_Diff_VEF_Face::ajouter_cas_multi_scalaire
void ajouter_cas_multi_scalaire(const DoubleTab &inconnue, DoubleTab &resu, DoubleTab &flux_bords, DoubleTab &nu, const Domaine_Cl_VEF &domaine_Cl_VEF, const Domaine_VEF &domaine_VEF, int nb_comp) const
Definition Op_Diff_VEF_Face.cpp:461
Op_Diff_VEF_Face::ajouter_contribution_multi_scalaire
void ajouter_contribution_multi_scalaire(const DoubleTab &, Matrice_Morse &) const
Definition Op_Diff_VEF_Face.cpp:1016
Op_Diff_VEF_Face::associer_diffusivite
void associer_diffusivite(const Champ_base &) override
Associate the diffusivity field.
Definition Op_Diff_VEF_Face.cpp:56
Op_Diff_VEF_Face::ajouter_contribution
void ajouter_contribution(const DoubleTab &, Matrice_Morse &) const
Definition Op_Diff_VEF_Face.cpp:724
Op_Diff_VEF_base
class Op_Diff_VEF_base
Definition Op_Diff_VEF_base.h:39
Op_Diff_VEF_base::phi_psi_diffuse
int phi_psi_diffuse(const Equation_base &eq) const
Determine whether to compute div(phi nu grad Psi) or div(nu grad Phi psi).
Definition Op_Diff_VEF_base.cpp:44
Op_Diff_VEF_base::nu_
DoubleTab nu_
Definition Op_Diff_VEF_base.h:66
Op_Diff_VEF_base::remplir_nu
virtual void remplir_nu(DoubleTab &) const
Definition Op_Diff_VEF_base.cpp:251
Op_Diff_VEF_base::viscA
double viscA(int face_i, int face_j, int num_elem, const _TYPE_ &diffu) const
Definition Op_Diff_VEF_base.h:92
Op_VEF_Face::modifier_matrice_pour_periodique_apres_contribuer
void modifier_matrice_pour_periodique_apres_contribuer(Matrice_Morse &matrice, const Equation_base &) const
Sums the 2 rows of the associated periodic faces, allowing computations in the code without needing t...
Definition Op_VEF_Face.cpp:738
Op_VEF_Face::modifier_matrice_pour_periodique_avant_contribuer
void modifier_matrice_pour_periodique_avant_contribuer(Matrice_Morse &matrice, const Equation_base &) const
Divides the coefficients on the periodic face rows by 2 in preparation for applying modifier_matrice_...
Definition Op_VEF_Face.cpp:725
Op_VEF_Face::modifier_flux
void modifier_flux(const Operateur_base &) const
Definition Op_VEF_Face.cpp:338
Operateur_base::flux_bords_
DoubleTab flux_bords_
Definition Operateur_base.h:140
Operateur_base::completer
virtual void completer()
Associates the operator with the domaine_dis, the domaine_Cl_dis, and the unknown of its equation.
Definition Operateur_base.cpp:89
Operateur_base::flux_bords
DoubleTab & flux_bords()
Definition Operateur_base.h:107
Periodique
class Periodique This class represents a periodic boundary condition.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Robin_VEF
Class Robin_VEF for Robin boundary conditions.
Definition Robin_VEF.h:44
Robin_VEF::get_beta_cl
double get_beta_cl() const
Definition Robin_VEF.h:53
Robin_VEF::flux_tangentiel_imp
double flux_tangentiel_imp(int i, int j) const
Returns the value of the imposed tangential flux on the (i,j)-th component of the field representing ...
Definition Robin_VEF.cpp:111
Robin_VEF::flux_normal_imp
double flux_normal_imp(int i) const
Returns the value of the imposed normal flux on the i-th component of the field representing the flux...
Definition Robin_VEF.cpp:103
Robin_VEF::get_alpha_cl
double get_alpha_cl() const
Definition Robin_VEF.h:52
Sortie
Base class for output streams.
Definition Sortie.h:52
Support_Champ_Masse_Volumique::declare_support_masse_volumique
virtual void declare_support_masse_volumique(int ok)
The constructor of a derived class that uses the density field must call this function with the value...
Definition Support_Champ_Masse_Volumique.cpp:46
Symetrie
Symetrie On symmetry faces, the following properties hold:
Definition Symetrie.h:37
TRUSTTab::view_wo
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, View< _TYPE_, _SHAPE_ > > view_wo()
Definition TRUSTTab.h:276
TRUSTTab::nb_dim
int nb_dim() const
Definition TRUSTTab.h:199
TRUSTTab::view_ro
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, ConstView< _TYPE_, _SHAPE_ > > view_ro() const
Definition TRUSTTab.h:261
TRUSTTab::view_rw
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, View< _TYPE_, _SHAPE_ > > view_rw()
Definition TRUSTTab.h:291
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133