TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Eval_Diff_VDF_Face_Gen.tpp
1/****************************************************************************
2* Copyright (c) 2023, 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#ifndef Eval_Diff_VDF_Face_Gen_TPP_included
17#define Eval_Diff_VDF_Face_Gen_TPP_included
18
19#include <Schema_Temps_base.h>
20#include <Probleme_base.h>
21#include <Equation_base.h>
22#include <Domaine_Cl_VDF.h>
23
24/* ************************************** *
25 * ********* POUR L'EXPLICITE ********** *
26 * ************************************** */
27
28template <typename DERIVED_T> template<Type_Flux_Fa7 Fa7_Type, typename Type_Double> inline std::enable_if_t< Fa7_Type == Type_Flux_Fa7::ELEM, void>
29Eval_Diff_VDF_Face_Gen<DERIVED_T>::flux_fa7(const DoubleTab& inco, const DoubleTab*, int elem, int fac1, int fac2, Type_Double& flux) const
30{
31 const int ori = orientation(fac1), ncomp = flux.size_array();
32 const double dist = dist_face(fac1,fac2,ori), surf = 0.5*(surface(fac1)*porosite(fac1)+surface(fac2)*porosite(fac2));
33 for (int k = 0; k < ncomp; k++)
34 {
35 const double tau = (inco(fac2,k)-inco(fac1,k))/dist, tau_tr = ACTIVATE_TAU_TR ? tau : 0.0;
36 const double visc_lam = nu_lam(elem, k), visc_turb = DERIVED_T::IS_TURB ? nu_turb(elem, k) : 0.;
37 flux[k] = ((tau + tau_tr) * (visc_lam + visc_turb)) * surf;
38 }
39}
40
41template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::INTERNE, void>
42Eval_Diff_VDF_Face_Gen<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab*, int fac1, int fac2, int fac3, int fac4, Type_Double& flux) const
43{
44 const int ori1 = orientation(fac1), ori3 = orientation(fac3), elem1 = elem_(fac3,0), elem2 = elem_(fac3,1), elem3 = elem_(fac4,0), elem4 = elem_(fac4,1), ncomp = flux.size_array();
45 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2);
46
47 for (int k = 0; k < ncomp; k++)
48 {
49 const double tau = (inco(fac4,k)-inco(fac3,k))/dist_face(fac3,fac4,ori1), tau_tr = ACTIVATE_TAU_TR ? (inco(fac2,k)-inco(fac1,k))/dist_face(fac1,fac2,ori3) : 0.0;
50 const int ind = DERIVED_T::IS_ANISO ? ori3 : k;
51 const double visc_lam = nu_lam_mean_4pts(elem1,elem2,elem3,elem4,ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_4pts(elem1,elem2,elem3,elem4,ind) : 0.0;
52 flux[k] = ((tau + tau_tr) * (visc_lam + visc_turb)) * surf * poros;
53 }
54}
55
56template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::MIXTE, void>
57Eval_Diff_VDF_Face_Gen<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab*, int fac1, int fac2, int fac3, int fac4, Type_Double& flux) const
58{
59 const int N = flux.size_array();
60 int elem[4], ori1 = orientation(fac1), ori3 = orientation(fac3);
61 elem[0] = elem_(fac3,0), elem[1] = elem_(fac3,1), elem[2] = elem_(fac4,0), elem[3] = elem_(fac4,1);
62 std::vector<double> visc_lam_temp(N), visc_turb_temp(N);
63 for (int k = 0; k < N; k++)
64 for (int i = 0; i < 4; i++)
65 if (elem[i] != -1)
66 {
67 visc_lam_temp[k] += nu_lam(elem[i], k);
68 visc_turb_temp[k] += nu_turb(elem[i], k);
69 }
70 for (int k = 0; k < N; k++)
71 {
72 visc_lam_temp[k] /= 3.0;
73 visc_turb_temp[k] /= 3.0;
74 }
75 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2);
76
77 for (int k = 0; k < N; k++)
78 if (inco(fac4,k)*inco(fac3,k) != 0)
79 {
80 const double visc_lam = visc_lam_temp[k], visc_turb = DERIVED_T::IS_TURB ? visc_turb_temp[k] : 0.0;
81 const double tau = (inco(fac4,k)-inco(fac3,k))/dist_face(fac3,fac4,ori1), tau_tr = ACTIVATE_TAU_TR ? (inco(fac2,k)-inco(fac1,k))/dist_face(fac1,fac2,ori3) : 0.0;
82 flux[k] = ((tau + tau_tr) * (visc_lam + visc_turb)) * surf * poros;
83 }
84}
85
86template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
87inline std::enable_if_t<Arete_Type == Type_Flux_Arete::PAROI || Arete_Type == Type_Flux_Arete::NAVIER_PAROI, void>
88Eval_Diff_VDF_Face_Gen<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab*, int fac1, int fac2, int fac3, int signe, Type_Double& flux) const
89{
90 constexpr bool is_PAROI = (Arete_Type == Type_Flux_Arete::PAROI);
91 const int rang1 = (fac1-premiere_face_bord), rang2 = (fac2-premiere_face_bord), ori = orientation(fac3), ncomp = flux.size_array();
92 if ( !uses_wall_law() )
93 {
94 int elem1 = elem_(fac3,0), elem2 = elem_(fac3,1);
95 if (is_PAROI)
96 {
97 if (elem1 == -1) elem1 = elem2;
98 else if (elem2 == -1) elem2 = elem1;
99 }
100
101 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2), dist = dist_norm_bord(fac1), tps = inconnue->temps();
102
103 const double vit_imp = is_PAROI ? 0.5*(Champ_Face_get_val_imp_face_bord(tps,rang1,ori,la_zcl.valeur())+Champ_Face_get_val_imp_face_bord(tps,rang2,ori,la_zcl.valeur())) :
104 0.5*(Champ_Face_get_val_imp_face_bord_sym(inco,tps,rang1,ori,la_zcl.valeur())+Champ_Face_get_val_imp_face_bord_sym(inco,tps,rang2,ori,la_zcl.valeur()));
105
106 double coeff = 0.0;
107 for (int k = 0; k < ncomp; k++)
108 {
109 if (!is_PAROI) // NAVIER_PAROI
110 coeff = 0.5 * (Champ_Face_coeff_frottement_face_bord(fac1, k, la_zcl.valeur()) + Champ_Face_coeff_frottement_face_bord(fac2, k, la_zcl.valeur()));
111
112 const int ind = DERIVED_T::IS_ANISO ? ori : k;
113 const double visc_lam = nu_lam_mean_2pts(elem1,elem2,ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_2pts(elem1,elem2,ind) : 0.0;
114 const double tau = (signe*(vit_imp-inco(fac3,k))/dist) - (signe * coeff * inco(fac3, k)), tau_tr = 0.;
115 flux[k] = ((tau + tau_tr) * (visc_lam + visc_turb)) * surf * poros;
116 }
117 }
118 else
119 {
120 double tau1 = tau_tan(rang1,ori)*0.5*surface(fac1), tau2 = tau_tan(rang2,ori)*0.5*surface(fac2);
121 for (int k = 0; k < ncomp; k++) flux[k] = tau1 + tau2;
122 }
123}
124
125template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
126inline std::enable_if_t< (Arete_Type == Type_Flux_Arete::NAVIER), void>
127Eval_Diff_VDF_Face_Gen<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab*, int fac1, int fac2, int fac3, int signe, Type_Double& flux) const
128{
129 // |
130 // fac 3 | fac 2
131 // --------
132 // | fac 1
133 // |
134
135 // fac3 est la face interne et fac1 et fac2 sont au bord Navier
136 // XXX : WARNING : nu/nu_turb deja dans coeff
137 const int ncomp = flux.size_array();
138 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2);
139 for (int k = 0; k < ncomp; k++)
140 {
141 const double coeff = 0.5 * (Champ_Face_coeff_frottement_face_bord(fac1, k, la_zcl.valeur()) + Champ_Face_coeff_frottement_face_bord(fac2, k, la_zcl.valeur()));
142 const double tau = - signe * coeff * inco(fac3, k), tau_tr = 0.;
143 flux[k] = (tau + tau_tr) * surf * poros;
144 }
145}
146
147template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
148inline std::enable_if_t< Arete_Type == Type_Flux_Arete::FLUIDE || Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE || Arete_Type == Type_Flux_Arete::PAROI_FLUIDE, void>
149Eval_Diff_VDF_Face_Gen<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab*, int fac1, int fac2, int fac3, int signe, Type_Double& flux3, Type_Double& flux1_2) const
150{
151 assert (flux3.size_array() == flux1_2.size_array());
152 constexpr bool is_NAV_FL = (Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE), is_PAR_FL = (Arete_Type == Type_Flux_Arete::PAROI_FLUIDE);
153 const int rang1 = (fac1-premiere_face_bord), rang2 = (fac2-premiere_face_bord), elem1 = elem_(fac3,0), elem2 = elem_(fac3,1), ori= orientation(fac3), ncomp = flux3.size_array();
154 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2), surfporos = surface(fac3)*porosite(fac3), tps = inconnue->temps(),
155 dist1 = dist_norm_bord(fac1), dist2 = dist_face(fac1,fac2,ori);
156
157 double vit_imp, coeff = 0.0;
158
159 for (int k = 0; k < ncomp; k++)
160 {
161 const int ind = DERIVED_T::IS_ANISO ? ori : k;
162 const double visc_lam = nu_lam_mean_2pts(elem1,elem2,ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_2pts(elem1,elem2,ind) : 0.0;
163 if (is_NAV_FL)
164 {
165 vit_imp = 0.5*(Champ_Face_get_val_imp_face_bord_sym(inco,tps,rang1,ori,la_zcl.valeur())+ Champ_Face_get_val_imp_face_bord_sym(inco,tps,rang2,ori,la_zcl.valeur()));
166 coeff = 0.5 * (Champ_Face_coeff_frottement_face_bord(fac1, k, la_zcl.valeur()) + Champ_Face_coeff_frottement_face_bord(fac2, k, la_zcl.valeur()));
167 }
168 else if (is_PAR_FL) // On ne sait pas qui de fac1 ou de fac2 est la face de paroi
169 {
170 if (est_egal(inco(fac1,k),0)) vit_imp = Champ_Face_get_val_imp_face_bord(tps,rang2,ori,la_zcl.valeur()); // fac1 est la face de paroi
171 else vit_imp = Champ_Face_get_val_imp_face_bord(tps,rang1,ori,la_zcl.valeur()); // fac2 est la face de paroi
172 }
173 else vit_imp = 0.5*(Champ_Face_get_val_imp_face_bord(tps,rang1,ori,la_zcl.valeur())+Champ_Face_get_val_imp_face_bord(tps,rang2,ori,la_zcl.valeur()));
174
175 // |
176 // fac 3 | fac 2
177 // --------
178 // | fac 1
179 // |
180
181 // fac3 est la face interne et fac1 et fac2 sont au bord
182 const double tau_3 = (signe*(vit_imp-inco(fac3,k))/dist1) -(signe * coeff * inco(fac3, k)),
183 tau_12 = (inco(fac2,k)-inco(fac1,k))/dist2, tau_tr_3 = ACTIVATE_TAU_TR ? tau_12 : 0.0, tau_tr_12 = ACTIVATE_TAU_TR ? tau_3 : 0.0;
184
185 flux3[k] = ((tau_3 + tau_tr_3) * (visc_lam + visc_turb)) * surf * poros;
186 flux1_2[k] = ((tau_12 + tau_tr_12) * (visc_lam + visc_turb)) * surfporos;
187 }
188}
189
190template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::PERIODICITE, void>
191Eval_Diff_VDF_Face_Gen<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab*, int fac1, int fac2, int fac3, int fac4, Type_Double& flux3_4, Type_Double& flux1_2) const
192{
193 assert (flux3_4.size_array() == flux1_2.size_array());
194 const int ori1 = orientation(fac1), ori3 = orientation(fac3), elem1 = elem_(fac3,0), elem2 = elem_(fac3,1), elem3 = elem_(fac4,0), elem4 = elem_(fac4,1), ncomp = flux3_4.size_array();
195 const double dist3_4 = dist_face_period(fac3,fac4,ori1), dist1_2 = dist_face_period(fac1,fac2,ori3),
196 surf1_2 = surface_(fac1,fac2), poros1_2 = porosity_(fac1, fac2), surf3_4 = surface_(fac3,fac4), poros3_4 = porosity_(fac3, fac4);
197
198 for (int k = 0; k < ncomp; k++)
199 {
200 const int ind = DERIVED_T::IS_ANISO ? ori3 : k;
201 const double visc_lam = nu_lam_mean_4pts(elem1,elem2,elem3,elem4,ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_4pts(elem1,elem2,elem3,elem4,ind) : 0.0;
202 const double tau_34 = (inco(fac4,k)-inco(fac3,k))/dist3_4, tau_12 = (inco(fac2,k)-inco(fac1,k))/dist1_2, tau_tr_34 = ACTIVATE_TAU_TR ? tau_12 : 0.0, tau_tr_12 = ACTIVATE_TAU_TR ? tau_34 : 0.0;
203 flux3_4[k] = ((tau_34 + tau_tr_34) * (visc_lam + visc_turb)) * surf1_2 * poros1_2;
204 flux1_2[k] = ((tau_12 + tau_tr_12) * (visc_lam + visc_turb)) * surf3_4 * poros3_4;
205 }
206}
207
208/* ************************************** *
209 * ********* POUR L'IMPLICITE ********** *
210 * ************************************** */
211
212template <typename DERIVED_T> template<Type_Flux_Fa7 Fa7_Type, typename Type_Double> inline std::enable_if_t< Fa7_Type == Type_Flux_Fa7::ELEM, void>
213Eval_Diff_VDF_Face_Gen<DERIVED_T>::coeffs_fa7(const DoubleTab*, int elem,int fac1, int fac2, Type_Double& f1, Type_Double& f2) const
214{
215 assert(f1.size_array() == f2.size_array());
216 const int ori = orientation(fac1), ncomp = f1.size_array();
217 const double dist = dist_face(fac1,fac2,ori), surf = 0.5*(surface(fac1)*porosite(fac1)+surface(fac2)*porosite(fac2)),
218 d_tau = 1. / dist, d_tau_tr = ACTIVATE_TAU_TR ? d_tau : 0.0; // On derive par rapport a fac1 et fac2
219
220 for (int k = 0; k < ncomp; k++)
221 {
222 const double visc_lam = nu_lam(elem, k), visc_turb = DERIVED_T::IS_TURB ? nu_turb(elem, k) : 0.;
223 f1[k] = f2[k] = ((d_tau + d_tau_tr) * (visc_lam + visc_turb)) * surf;
224 }
225
226 if (TEST_COEFFS) test_coeffs_fa7<Fa7_Type,Type_Double>(elem,fac1,fac2,f1);
227}
228
229template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::INTERNE, void>
230Eval_Diff_VDF_Face_Gen<DERIVED_T>::coeffs_arete(const DoubleTab*, int fac1, int fac2, int fac3, int fac4, Type_Double& aii, Type_Double& ajj) const
231{
232 assert(aii.size_array() == ajj.size_array());
233 const int ori1 = orientation(fac1), ori3 = orientation(fac3), elem1 = elem_(fac3,0), elem2 = elem_(fac3,1), elem3 = elem_(fac4,0), elem4 = elem_(fac4,1), ncomp = aii.size_array();
234 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2),
235 d_tau = 1.0 / dist_face(fac3,fac4,ori1), d_tau_tr = 0.; // On derive par rapport a fac3 et fac4
236
237 for (int k = 0; k < ncomp; k++)
238 {
239 const int ind = DERIVED_T::IS_ANISO ? ori3 : k;
240 const double visc_lam = nu_lam_mean_4pts(elem1,elem2,elem3,elem4,ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_4pts(elem1,elem2,elem3,elem4,ind) : 0.0;
241 aii[k] = ajj[k] = ((d_tau + d_tau_tr) * (visc_lam + visc_turb)) * surf * poros;
242 }
243
244 if (TEST_COEFFS) test_coeffs_arete<Arete_Type,Type_Double>(fac1,fac2,fac3,fac4,aii);
245}
246
247template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::MIXTE, void>
248Eval_Diff_VDF_Face_Gen<DERIVED_T>::coeffs_arete(const DoubleTab*,int fac1, int fac2, int fac3, int fac4, Type_Double& aii, Type_Double& ajj) const
249{
250 assert(aii.size_array() == ajj.size_array());
251 const int ori1 = orientation(fac1), ncomp = aii.size_array();
252 int elem[4];
253 elem[0] = elem_(fac3,0), elem[1] = elem_(fac3,1), elem[2] = elem_(fac4,0), elem[3] = elem_(fac4,1);
254
255 std::vector<double> visc_lam_temp(ncomp), visc_turb_temp(ncomp);
256 for (int k = 0; k < ncomp; k++)
257 for (int i = 0; i < 4; i++)
258 if (elem[i] != -1)
259 {
260 visc_lam_temp[k] += nu_lam(elem[i], k);
261 visc_turb_temp[k] += nu_turb(elem[i], k);
262 }
263 for (int k = 0; k < ncomp; k++)
264 {
265 visc_lam_temp[k] /= 3.0;
266 visc_turb_temp[k] /= 3.0;
267 }
268
269 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2),
270 d_tau = 1. / dist_face(fac3,fac4,ori1), d_tau_tr = 0.; // On derive par rapport a fac3 et fac4
271
272 const DoubleTab& inco = inconnue->valeurs();
273 for (int k = 0; k < ncomp; k++)
274 if (inco(fac4,k) * inco(fac3,k) != 0)
275 {
276 const double visc_lam = visc_lam_temp[k], visc_turb = DERIVED_T::IS_TURB ? visc_turb_temp[k] : 0.;
277 aii[k] = ajj[k] = ((d_tau + d_tau_tr) * (visc_lam + visc_turb)) * surf * poros;
278 }
279
280 if (TEST_COEFFS) test_coeffs_arete<Arete_Type,Type_Double>(fac1,fac2,fac3,fac4,aii);
281}
282
283template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
284inline std::enable_if_t< Arete_Type == Type_Flux_Arete::PAROI || Arete_Type == Type_Flux_Arete::NAVIER_PAROI, void>
285Eval_Diff_VDF_Face_Gen<DERIVED_T>::coeffs_arete(const DoubleTab*, int fac1, int fac2, int fac3, int signe, Type_Double& aii1_2, Type_Double& aii3_4, Type_Double& ajj1_2) const
286{
287 assert(aii1_2.size_array() == aii3_4.size_array() && aii1_2.size_array() == ajj1_2.size_array());
288 constexpr bool is_PAROI = (Arete_Type == Type_Flux_Arete::PAROI);
289 const int ncomp = aii1_2.size_array();
290 if ( !uses_wall_law())
291 {
292 int elem1 = elem_(fac3,0), elem2 = elem_(fac3,1), ori = orientation(fac3);
293 if(is_PAROI)
294 {
295 if (elem1 == -1) elem1 = elem2;
296 else if (elem2 == -1) elem2 = elem1;
297 }
298
299 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2), dist = dist_norm_bord(fac1);
300
301 double coeff = 0.0;
302 for (int k = 0; k < ncomp; k++)
303 {
304 if (!is_PAROI) // NAVIER_PAROI
305 coeff = 0.5 * (Champ_Face_coeff_frottement_face_bord(fac1, k, la_zcl.valeur()) + Champ_Face_coeff_frottement_face_bord(fac2, k, la_zcl.valeur()));
306
307 const double d_tau = signe / dist - (signe * coeff), d_tau_tr = 0.; // On a pas derive ... deja nul dans le flux !
308 const int ind = DERIVED_T::IS_ANISO ? ori : k;
309 const double visc_lam = nu_lam_mean_2pts(elem1,elem2,ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_2pts(elem1,elem2,ind) : 0.;
310 aii1_2[k] = ajj1_2[k] = 0.;
311 aii3_4[k] = ((d_tau + d_tau_tr) * (visc_lam + visc_turb)) * surf * poros;
312 }
313 }
314 else for (int k = 0; k < ncomp; k++) aii3_4[k] = aii1_2[k] = ajj1_2[k] = 0.;
315
316 if (TEST_COEFFS) test_coeffs_arete<Arete_Type,Type_Double>(fac1,fac2,fac3,signe,aii3_4);
317}
318
319template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
320inline std::enable_if_t< Arete_Type == Type_Flux_Arete::FLUIDE || Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE || Arete_Type == Type_Flux_Arete::PAROI_FLUIDE, void>
321Eval_Diff_VDF_Face_Gen<DERIVED_T>::coeffs_arete(const DoubleTab*, int fac1, int fac2, int fac3, int signe, Type_Double& aii1_2, Type_Double& aii3, Type_Double& ajj1_2) const
322{
323 assert(aii1_2.size_array() == aii3.size_array() && aii1_2.size_array() == ajj1_2.size_array());
324 constexpr bool is_NAV_FL = (Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE);
325 const int elem1 = elem_(fac3,0), elem2 = elem_(fac3,1), ori= orientation(fac3), ncomp = aii1_2.size_array();
326 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2), surfporos = surface(fac3)*porosite(fac3),
327 dist1 = dist_norm_bord(fac1), dist2 = dist_face(fac1,fac2,ori);
328
329 double coeff = 0.0;
330 for (int k = 0; k < ncomp; k++)
331 {
332 if (is_NAV_FL)
333 coeff = 0.5 * (Champ_Face_coeff_frottement_face_bord(fac1, k, la_zcl.valeur()) + Champ_Face_coeff_frottement_face_bord(fac2, k, la_zcl.valeur()));
334
335 const double d_tau_3 = (signe / dist1) - (signe * coeff), d_tau_tr_3 = 0., // On derive par rapport a fac3
336 d_tau_12 = 1. / dist2, d_tau_tr_12 = 0.; // On derive par rapport a fac1 et fac2
337
338 const int ind = DERIVED_T::IS_ANISO ? ori : k;
339 const double visc_lam = nu_lam_mean_2pts(elem1,elem2,ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_2pts(elem1,elem2,ind) : 0.;
340 aii3[k] = ((d_tau_3 + d_tau_tr_3) * (visc_lam + visc_turb)) * surf * poros;
341 aii1_2[k] = ajj1_2[k] = ((d_tau_12 + d_tau_tr_12) * (visc_lam + visc_turb)) * surfporos;
342 }
343
344 if (TEST_COEFFS) test_coeffs_arete<Arete_Type,Type_Double>(fac1,fac2,fac3,signe,aii1_2,aii3);
345}
346
347template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
348inline std::enable_if_t< (Arete_Type == Type_Flux_Arete::NAVIER), void>
349Eval_Diff_VDF_Face_Gen<DERIVED_T>::coeffs_arete(const DoubleTab*, int fac1, int fac2, int fac3, int signe, Type_Double& aii1_2, Type_Double& aii3, Type_Double& ajj1_2) const
350{
351 assert(aii1_2.size_array() == aii3.size_array() && aii1_2.size_array() == ajj1_2.size_array());
352 const int elem1 = elem_(fac3,0), elem2 = elem_(fac3,1), ncomp = aii1_2.size_array(), ori = orientation(fac3);
353 const double surf = surface_(fac1,fac2), poros = porosity_(fac1,fac2), dist2 = dist_face(fac1,fac2,ori), surfporos = surface(fac3)*porosite(fac3);
354
355 for (int k = 0; k < ncomp; k++)
356 {
357 const int ind = DERIVED_T::IS_ANISO ? ori : k;
358 const double visc_lam = nu_lam_mean_2pts(elem1, elem2, ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_2pts(elem1, elem2, ind) : 0.0;
359 const double coeff = 0.5 * (Champ_Face_coeff_frottement_face_bord(fac1, k, la_zcl.valeur()) + Champ_Face_coeff_frottement_face_bord(fac2, k, la_zcl.valeur()));
360 const double d_tau_3 = - (signe * coeff), d_tau_tr_3 = 0., // On derive par rapport a fac3
361 d_tau_12 = 1. / dist2, d_tau_tr_12 = 0.; // On derive par rapport a fac1 et fac2
362
363 aii3[k] = ((d_tau_3 + d_tau_tr_3) * (visc_lam + visc_turb)) * surf * poros;
364 aii1_2[k] = ajj1_2[k] = ((d_tau_12 + d_tau_tr_12) * (visc_lam + visc_turb)) * surfporos;
365 }
366}
367
368template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::PERIODICITE, void>
369Eval_Diff_VDF_Face_Gen<DERIVED_T>::coeffs_arete(const DoubleTab*,int fac1, int fac2, int fac3, int fac4, Type_Double& aii, Type_Double& ajj) const
370{
371 assert(aii.size_array() == ajj.size_array());
372 const int ori1 = orientation(fac1), ori3 = orientation(fac3), elem1 = elem_(fac3,0), elem2 = elem_(fac3,1), elem3 = elem_(fac4,0), elem4 = elem_(fac4,1), ncomp = aii.size_array();
373 const double dist3_4 = dist_face_period(fac3,fac4,ori1), surf = surface_(fac1,fac2), poros = porosity_(fac1, fac2),
374 d_tau = 1. / dist3_4, d_tau_tr = 0.; // On derive par rapport a fac3 et fac4
375
376 for (int k = 0; k < ncomp; k++)
377 {
378 const int ind = DERIVED_T::IS_ANISO ? ori3 : k;
379 const double visc_lam = nu_lam_mean_4pts(elem1,elem2,elem3,elem4,ind), visc_turb = DERIVED_T::IS_TURB ? nu_mean_4pts(elem1,elem2,elem3,elem4,ind) : 0.;
380 aii[k] = ajj[k] = ((d_tau + d_tau_tr) * (visc_lam + visc_turb)) * surf * poros;
381 }
382
383 if (TEST_COEFFS) test_coeffs_arete<Arete_Type,Type_Double>(fac1,fac2,fac3,fac4,aii);
384}
385
386/* ************************************** *
387 * ********* For checking *********** *
388 * ************************************** */
389
390template <typename DERIVED_T> template<typename Type_Double>
391void Eval_Diff_VDF_Face_Gen<DERIVED_T>::check_error(const char * nom_funct, const int Type_Flux, const int ncomp, const Type_Double& f1, const Type_Double& flux_p, const Type_Double& flux_m) const
392{
393 for (int k = 0; k < ncomp; k++)
394 {
395 const double err = std::fabs(((flux_p[0] - flux_m[0]) / (2.0 * EPS)) - f1[0]);
396 if (err > EPS)
397 {
398 std::cerr << "Error in function : " << nom_funct << "_Type_Flux_" << Type_Flux << " : " << f1[0] << " " << (flux_p[0] - flux_m[0]) / (2.0 * EPS) << std::endl;
399 Process::exit();
400 }
401 }
402
403 if (la_zcl->equation().probleme().schema_temps().nb_pas_dt() > 10)
404 {
405 Cerr << "All OK ! Coeffs and fluxes implementations are coherent !" << finl;
406 Process::exit();
407 }
408}
409
410template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
411inline std::enable_if_t< Arete_Type == Type_Flux_Arete::INTERNE || Arete_Type == Type_Flux_Arete::MIXTE, void>
412Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_common(const int fac1, const int fac2, const int fac3, const int fac4, Type_Double& flux_p, Type_Double& flux_m) const
413{
414 DoubleTab inco_pert = inconnue->valeurs();
415
416 const int ncomp = flux_p.size_array();
417 for (int k = 0; k < ncomp; k++) inco_pert(fac4,k) += EPS; // XXX : ATTENTION SIGNE
418 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,fac4,flux_p);
419
420 for (int k = 0; k < ncomp; k++) inco_pert(fac4,k) -= 2.0 * EPS; // XXX : ATTENTION SIGNE
421 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,fac4,flux_m);
422}
423
424template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
425inline std::enable_if_t<Arete_Type == Type_Flux_Arete::PAROI || Arete_Type == Type_Flux_Arete::NAVIER_PAROI, void>
426Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_common(const int fac1, const int fac2, const int fac3, const int signe, Type_Double& flux_p, Type_Double& flux_m) const
427{
428 DoubleTab inco_pert = inconnue->valeurs();
429
430 const int ncomp = flux_p.size_array();
431 for (int k = 0; k < ncomp; k++) inco_pert(fac3,k) -= EPS; // XXX : ATTENTION SIGNE
432
433 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,signe,flux_p);
434
435 for (int k = 0; k < ncomp; k++) inco_pert(fac3,k) += 2.0 * EPS; // XXX : ATTENTION SIGNE
436 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,signe,flux_m);
437}
438
439template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
440inline std::enable_if_t<Arete_Type == Type_Flux_Arete::FLUIDE || Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE || Arete_Type == Type_Flux_Arete::PAROI_FLUIDE, void>
441Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_common(const int fac1, const int fac2, const int fac3, const int signe, Type_Double& flux_p3, Type_Double& flux_m3, Type_Double& flux_p12, Type_Double& flux_m12) const
442{
443 DoubleTab inco_pert = inconnue->valeurs();
444 Type_Double poubelle(flux_p3.size_array());
445
446 const int ncomp = flux_p3.size_array();
447
448 // pour flux3
449 for (int k = 0; k < ncomp; k++) inco_pert(fac3,k) -= EPS; // XXX : ATTENTION SIGNE
450 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,signe,flux_p3,poubelle);
451
452 for (int k = 0; k < ncomp; k++) inco_pert(fac3,k) += 2.0 * EPS; // XXX : ATTENTION SIGNE
453 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,signe,flux_m3,poubelle);
454
455 // pour flux1_2
456 inco_pert = inconnue->valeurs(); // back to real values
457 for (int k = 0; k < ncomp; k++) inco_pert(fac2,k) += EPS; // XXX : ATTENTION SIGNE
458 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,signe,poubelle,flux_p12);
459
460 for (int k = 0; k < ncomp; k++) inco_pert(fac2,k) -= 2.0 * EPS; // XXX : ATTENTION SIGNE
461 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,signe,poubelle,flux_m12);
462}
463
464template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
465inline std::enable_if_t<Arete_Type == Type_Flux_Arete::PERIODICITE, void>
466Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_common(const int fac1, const int fac2, const int fac3, const int fac4, Type_Double& flux_p, Type_Double& flux_m) const
467{
468 DoubleTab inco_pert = inconnue->valeurs();
469 Type_Double poubelle(flux_p.size_array());
470 const int ncomp = flux_p.size_array();
471
472 for (int k = 0; k < ncomp; k++) inco_pert(fac4,k) += EPS; // XXX : ATTENTION SIGNE
473 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,fac4,flux_p,poubelle);
474
475 for (int k = 0; k < ncomp; k++) inco_pert(fac4,k) -= 2.0 * EPS; // XXX : ATTENTION SIGNE
476 flux_arete<Arete_Type,Type_Double>(inco_pert,nullptr,fac1,fac2,fac3,fac4,flux_m,poubelle);
477}
478
479template <typename DERIVED_T> template<Type_Flux_Fa7 Fa7_Type, typename Type_Double> inline std::enable_if_t< Fa7_Type == Type_Flux_Fa7::ELEM, void>
480Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_fa7(const int elem, const int fac1, const int fac2, const Type_Double& f1) const
481{
482 const int ncomp = f1.size_array();
483 Type_Double flux_p(ncomp), flux_m(ncomp);
484 DoubleTab inco_pert = inconnue->valeurs();
485
486 for (int k = 0; k < ncomp; k++) inco_pert(fac2,k) += EPS; // XXX : ATTENTION SIGNE
487 flux_fa7<Fa7_Type,Type_Double>(inco_pert,nullptr,elem,fac1,fac2,flux_p);
488
489 for (int k = 0; k < ncomp; k++) inco_pert(fac2,k) -= 2.0 * EPS; // XXX : ATTENTION SIGNE
490 flux_fa7<Fa7_Type,Type_Double>(inco_pert,nullptr,elem,fac1,fac2,flux_m);
491
492 check_error<Type_Double>(__func__,(int)Fa7_Type,ncomp,f1,flux_p,flux_m);
493}
494
495template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
496inline std::enable_if_t< Arete_Type == Type_Flux_Arete::INTERNE, void>
497Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_arete(const int fac1, const int fac2, const int fac3, const int fac4, const Type_Double& aii) const
498{
499 const int ncomp = aii.size_array();
500 Type_Double flux_p(ncomp), flux_m(ncomp);
501 test_coeffs_common<Arete_Type,Type_Double>(fac1,fac2,fac3,fac4,flux_p,flux_m);
502 check_error<Type_Double>(__func__,(int)Arete_Type,ncomp,aii,flux_p,flux_m);
503}
504
505template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
506inline std::enable_if_t<Arete_Type == Type_Flux_Arete::MIXTE, void>
507Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_arete(const int fac1, const int fac2, const int fac3, const int fac4, const Type_Double& aii) const
508{
509 const int ncomp = aii.size_array();
510 Type_Double flux_p(ncomp), flux_m(ncomp);
511 test_coeffs_common<Arete_Type,Type_Double>(fac1,fac2,fac3,fac4,flux_p,flux_m);
512 if (inconnue->valeurs()(fac4,0) * inconnue->valeurs()(fac3,0) != 0) check_error<Type_Double>(__func__,(int)Arete_Type,ncomp,aii,flux_p,flux_m);
513}
514
515template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
516inline std::enable_if_t<Arete_Type == Type_Flux_Arete::PAROI || Arete_Type == Type_Flux_Arete::NAVIER_PAROI, void>
517Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_arete(const int fac1, const int fac2, const int fac3, const int signe, const Type_Double& aii3_4) const
518{
519 const int ncomp = aii3_4.size_array();
520 Type_Double flux_p(ncomp), flux_m(ncomp);
521 test_coeffs_common<Arete_Type,Type_Double>(fac1,fac2,fac3,signe,flux_p,flux_m);
522 if ( !uses_wall_law() ) check_error<Type_Double>(__func__,(int)Arete_Type,ncomp,aii3_4,flux_p,flux_m);
523}
524
525template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
526inline std::enable_if_t<Arete_Type == Type_Flux_Arete::FLUIDE || Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE || Arete_Type == Type_Flux_Arete::PAROI_FLUIDE, void>
527Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_arete(const int fac1, const int fac2, const int fac3, const int signe, const Type_Double& aii1_2, const Type_Double& aii3) const
528{
529 const int ncomp = aii1_2.size_array();
530 Type_Double flux_p3(ncomp), flux_m3(ncomp), flux_p12(ncomp), flux_m12(ncomp);
531 test_coeffs_common<Arete_Type,Type_Double>(fac1,fac2,fac3,signe,flux_p3,flux_m3,flux_p12,flux_m12);
532 check_error<Type_Double>(__func__,(int)Arete_Type,ncomp,aii3,flux_p3,flux_m3);
533 check_error<Type_Double>(__func__,(int)Arete_Type,ncomp,aii1_2,flux_p12,flux_m12);
534}
535
536template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
537inline std::enable_if_t<Arete_Type == Type_Flux_Arete::PERIODICITE, void>
538Eval_Diff_VDF_Face_Gen<DERIVED_T>::test_coeffs_arete(const int fac1, const int fac2, const int fac3, const int fac4, const Type_Double& aii) const
539{
540 const int ncomp = aii.size_array();
541 Type_Double flux_p(ncomp), flux_m(ncomp);
542 test_coeffs_common<Arete_Type,Type_Double>(fac1,fac2,fac3,fac4,flux_p,flux_m);
543 check_error<Type_Double>(__func__,(int)Arete_Type,ncomp,aii,flux_p,flux_m);
544}
545
546#endif /* Eval_Diff_VDF_Face_Gen_TPP_included */
Eval_Diff_VDF_Face_Gen::flux_fa7
std::enable_if_t< Fa7_Type==Type_Flux_Fa7::SORTIE_LIBRE, void > flux_fa7(const DoubleTab &, const DoubleTab *, int, const Neumann_sortie_libre &, int, Type_Double &) const
Definition Eval_Diff_VDF_Face_Gen.h:46
Eval_Diff_VDF_Face_Gen::coeffs_fa7
std::enable_if_t< Fa7_Type==Type_Flux_Fa7::SORTIE_LIBRE, void > coeffs_fa7(const DoubleTab *, int, const Neumann_sortie_libre &, Type_Double &, Type_Double &) const
Definition Eval_Diff_VDF_Face_Gen.h:82
Eval_Diff_VDF_Face_Gen::TEST_COEFFS
static constexpr bool TEST_COEFFS
Definition Eval_Diff_VDF_Face_Gen.h:35
Eval_Diff_VDF_Face_Gen::ACTIVATE_TAU_TR
static constexpr bool ACTIVATE_TAU_TR
Definition Eval_Diff_VDF_Face_Gen.h:35
Eval_Diff_VDF_Face_Gen::coeffs_arete
std::enable_if_t< Arete_Type==Type_Flux_Arete::INTERNE, void > coeffs_arete(const DoubleTab *, int, int, int, int, Type_Double &, Type_Double &) const
Definition Eval_Diff_VDF_Face_Gen.tpp:230
Eval_Diff_VDF_Face_Gen::flux_arete
std::enable_if_t< Arete_Type==Type_Flux_Arete::INTERNE, void > flux_arete(const DoubleTab &, const DoubleTab *, int, int, int, int, Type_Double &) const
Definition Eval_Diff_VDF_Face_Gen.tpp:42
Evaluateur_VDF::elem_
IntTab elem_
Definition Evaluateur_VDF.h:53
Evaluateur_VDF::porosite
DoubleVect porosite
Definition Evaluateur_VDF.h:56
Evaluateur_VDF::orientation
IntVect orientation
Definition Evaluateur_VDF.h:55
Evaluateur_VDF::dist_norm_bord
double dist_norm_bord(int) const
Definition Evaluateur_VDF.cpp:51
Evaluateur_VDF::dist_face
double dist_face(int fac1, int fac2, int k) const
Definition Evaluateur_VDF.h:43
Evaluateur_VDF::premiere_face_bord
int premiere_face_bord
Definition Evaluateur_VDF.h:52
Evaluateur_VDF::surface
DoubleVect surface
Definition Evaluateur_VDF.h:54
Evaluateur_VDF::dist_face_period
double dist_face_period(int fac1, int fac2, int k) const
Definition Evaluateur_VDF.h:42
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
TRUSTArray::size_array
_SIZE_ size_array() const
Definition TRUSTArray.tpp:187