TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Tensors_Computation_VEF.cpp
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#include <Tensors_Computation_VEF.h>
17#include <TRUSTTab.h>
18#include <Domaine_VEF.h>
19#include <Domaine_Cl_VEF.h>
20#include <Champ_P1NC.h>
21
22using TCV = Tensors_Computation_VEF;
23
24// Compute the antisymetric tensor
25void TCV::compute_enstrophy(const Domaine_VEF& dom_VEF,
26 const Domaine_Cl_VEF& dom_BC_VEF,
27 const DoubleTab& velocity,
28 DoubleTab& enstrophy) const
29{
30 // Initialisation
31 const int nb_elem_tot = dom_VEF.nb_elem_tot();
32 const int dimension = Objet_U::dimension;
33 DoubleTrav gradient_elem(nb_elem_tot, dimension, dimension);
34
35 // Computation of gradient
36 Champ_P1NC::calcul_gradient(velocity, gradient_elem, dom_BC_VEF);
37
38 // Loop on edge faces
39 antisym_loop_edge_faces_enstrophy(dom_VEF, dom_BC_VEF, gradient_elem, enstrophy);
40
41 // Loop on internal faces
42 antisym_loop_internal_faces_enstrophy(dom_VEF, gradient_elem, enstrophy);
43
44}
45
47 const Domaine_Cl_VEF& dom_BC_VEF,
48 const DoubleTab& gradient_elem,
49 DoubleTab& enstrophy) const
50{
51 for (int n_edge = 0; n_edge < dom_VEF.nb_front_Cl(); ++n_edge)
52 {
53 const Cond_lim& current_BC = dom_BC_VEF.les_conditions_limites(n_edge);
54 const Front_VF& the_edge = ref_cast(Front_VF, current_BC->frontiere_dis());
55
56 if (sub_type(Periodique, current_BC.valeur()))
57 antisym_loop_edges_periodiqueBC_enstrophy(dom_VEF, the_edge, gradient_elem, enstrophy);
58 else
59 antisym_loop_edges_general_enstrophy(dom_VEF, the_edge, gradient_elem, enstrophy);
60 }
61}
62
64 const Front_VF& the_edge,
65 const DoubleTab& tab_gradient_elem,
66 DoubleTab& tab_enstrophy) const
67{
68 const IntTab& tab_neighbour_face = dom_VEF.face_voisins();
69 const DoubleVect& tab_volumes = dom_VEF.volumes();
70
71 const int nbegin = the_edge.num_premiere_face();
72 const int nend = nbegin + the_edge.nb_faces();
73 const int dim = Objet_U::dimension;
74
75 CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
76 CDoubleArrView volumes = tab_volumes.view_ro();
77 CIntTabView neighbour_face = tab_neighbour_face.view_ro();
78 DoubleArrView enstrophy = static_cast<ArrOfDouble&>(tab_enstrophy).view_wo();
79 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(nbegin, nend), KOKKOS_LAMBDA(const int nface)
80 {
81 const int f0 = neighbour_face(nface, 0);
82 const int f1 = neighbour_face(nface, 1);
83
84 const double inv_vol = 1./(volumes(f0) + volumes(f1));
85 const double vol0 = volumes(f0)*inv_vol;
86 const double vol1 = volumes(f1)*inv_vol;
87
88 // gradient_elem(elem_num, dimension, dimension)
89 const double du_dy = vol0*gradient_elem(f0, 0, 1) + vol1*gradient_elem(f1, 0, 1);
90 const double dv_dx = vol0*gradient_elem(f0, 1, 0) + vol1*gradient_elem(f1, 1, 0);
91
92 double tmp = 0.5*((du_dy - dv_dx)*(du_dy - dv_dx) +
93 (dv_dx - du_dy)*(dv_dx - du_dy));
94
95 if (dim == 3)
96 {
97 const double du_dz = vol0*gradient_elem(f0, 0, 2) + vol1*gradient_elem(f1, 0, 2);
98 const double dv_dz = vol0*gradient_elem(f0, 1, 2) + vol1*gradient_elem(f1, 1, 2);
99 const double dw_dx = vol0*gradient_elem(f0, 2, 0) + vol1*gradient_elem(f1, 2, 0);
100 const double dw_dy = vol0*gradient_elem(f0, 2, 1) + vol1*gradient_elem(f1, 2, 1);
101
102 tmp += 0.5*((du_dz - dw_dx)*(du_dz - dw_dx) +
103 (dv_dz - dw_dy)*(dv_dz - dw_dy) +
104 (dw_dx - du_dz)*(dw_dx - du_dz) +
105 (dw_dy - dv_dz)*(dw_dy - dv_dz));
106 }
107
108 enstrophy(nface) = Kokkos::sqrt(tmp);
109 });
110 end_gpu_timer(__KERNEL_NAME__);
111}
112
113
115 const Front_VF& the_edge,
116 const DoubleTab& tab_gradient_elem,
117 DoubleTab& tab_enstrophy) const
118{
119 const IntTab& tab_neighbour_face = dom_VEF.face_voisins();
120
121 const int nbegin = the_edge.num_premiere_face();
122 const int nend = nbegin + the_edge.nb_faces();
123 const int dim = Objet_U::dimension;
124
125 // Create Kokkos views
126 CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
127 CIntTabView neighbour_face = tab_neighbour_face.view_ro();
128 DoubleArrView enstrophy = static_cast<ArrOfDouble&>(tab_enstrophy).view_wo();
129
130 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(nbegin, nend), KOKKOS_LAMBDA(const int nface)
131 {
132 const int f0 = neighbour_face(nface, 0);
133
134 // gradient_elem(elem_num, dimension, dimension)
135 const double du_dy = gradient_elem(f0, 0, 1);
136 const double dv_dx = gradient_elem(f0, 1, 0);
137
138 double tmp = 0.5*((du_dy - dv_dx)*(du_dy - dv_dx) +
139 (dv_dx - du_dy)*(dv_dx - du_dy));
140
141 if (dim == 3)
142 {
143 const double du_dz = gradient_elem(f0, 0, 2);
144 const double dv_dz = gradient_elem(f0, 1, 2);
145 const double dw_dx = gradient_elem(f0, 2, 0);
146 const double dw_dy = gradient_elem(f0, 2, 1);
147
148 tmp += 0.5*((du_dz - dw_dx)*(du_dz - dw_dx) +
149 (dv_dz - dw_dy)*(dv_dz - dw_dy) +
150 (dw_dx - du_dz)*(dw_dx - du_dz) +
151 (dw_dy - dv_dz)*(dw_dy - dv_dz));
152 }
153
154 enstrophy(nface) = Kokkos::sqrt(tmp);
155 });
156 end_gpu_timer(__KERNEL_NAME__);
157}
158
159// Même fonction que périodique !
161 const DoubleTab& tab_grad_elem,
162 DoubleTab& tab_enstrophy) const
163{
164 const IntTab& tab_neighbour_face = dom_VEF.face_voisins();
165 const DoubleVect& tab_volumes = dom_VEF.volumes();
166 const int first_internal_face = dom_VEF.premiere_face_int();
167 const int nb_faces = dom_VEF.nb_faces();
168 const int dim = Objet_U::dimension;
169
170 // Create Kokkos views
171 CDoubleTabView3 grad_elem = tab_grad_elem.view_ro<3>();
172 CDoubleArrView volumes = tab_volumes.view_ro();
173 CIntTabView neighbour_face = tab_neighbour_face.view_ro();
174 DoubleArrView enstrophy = static_cast<ArrOfDouble&>(tab_enstrophy).view_wo();
175
176 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(first_internal_face, nb_faces), KOKKOS_LAMBDA(const int nface)
177 {
178 const int f0 = neighbour_face(nface, 0);
179 const int f1 = neighbour_face(nface, 1);
180
181 const double inv_vol = 1./(volumes(f0) + volumes(f1));
182 const double vol0 = volumes(f0)*inv_vol;
183 const double vol1 = volumes(f1)*inv_vol;
184
185 // format: grad_elem(elem_num, dimension, dimension)
186 const double du_dy = vol0*grad_elem(f0, 0, 1) + vol1*grad_elem(f1, 0, 1);
187 const double dv_dx = vol0*grad_elem(f0, 1, 0) + vol1*grad_elem(f1, 1, 0);
188
189 double tmp = 0.5*((du_dy - dv_dx)*(du_dy - dv_dx) +
190 (dv_dx - du_dy)*(dv_dx - du_dy));
191
192 if (dim == 3)
193 {
194 const double du_dz = vol0*grad_elem(f0, 0, 2) + vol1*grad_elem(f1, 0, 2);
195 const double dv_dz = vol0*grad_elem(f0, 1, 2) + vol1*grad_elem(f1, 1, 2);
196 const double dw_dx = vol0*grad_elem(f0, 2, 0) + vol1*grad_elem(f1, 2, 0);
197 const double dw_dy = vol0*grad_elem(f0, 2, 1) + vol1*grad_elem(f1, 2, 1);
198
199 tmp += 0.5*((du_dz - dw_dx)*(du_dz - dw_dx) +
200 (dv_dz - dw_dy)*(dv_dz - dw_dy) +
201 (dw_dx - du_dz)*(dw_dx - du_dz) +
202 (dw_dy - dv_dz)*(dw_dy - dv_dz));
203 }
204
205 enstrophy(nface) = Kokkos::sqrt(tmp);
206 });
207 end_gpu_timer(__KERNEL_NAME__);
208}
209
210// Compute the antisymetric tensor
212 const Domaine_Cl_VEF& dom_BC_VEF,
213 const DoubleTab& velocity,
214 DoubleTab& strain_invariant) const
215{
216 // Initialisation
217 const int nb_elem_tot = dom_VEF.nb_elem_tot();
218 const int dimension = Objet_U::dimension;
219 DoubleTrav gradient_elem(nb_elem_tot, dimension, dimension);
220
221 // Computation of gradient
222 Champ_P1NC::calcul_gradient(velocity, gradient_elem, dom_BC_VEF);
223
224 // Loop on edge faces
225 antisym_loop_edge_faces_strain_invariant(dom_VEF, dom_BC_VEF, gradient_elem, strain_invariant);
226
227 // Loop on internal faces
228 antisym_loop_internal_faces_strain_invariant(dom_VEF, gradient_elem, strain_invariant);
229
230}
231
233 const Domaine_Cl_VEF& dom_BC_VEF,
234 const DoubleTab& gradient_elem,
235 DoubleTab& strain_invariant) const
236{
237 for (int n_edge = 0; n_edge < dom_VEF.nb_front_Cl(); ++n_edge)
238 {
239 const Cond_lim& current_BC = dom_BC_VEF.les_conditions_limites(n_edge);
240 const Front_VF& the_edge = ref_cast(Front_VF, current_BC->frontiere_dis());
241
242 if (sub_type(Periodique, current_BC.valeur()))
243 antisym_loop_edges_periodiqueBC_strain_invariant(dom_VEF, the_edge, gradient_elem, strain_invariant);
244 else
245 antisym_loop_edges_general_strain_invariant(dom_VEF, the_edge, gradient_elem, strain_invariant);
246 }
247}
248
250 const Front_VF& the_edge,
251 const DoubleTab& tab_gradient_elem,
252 DoubleTab& tab_strain_invariant) const
253{
254 const IntTab& tab_neighbour_face = dom_VEF.face_voisins();
255 const DoubleVect& tab_volumes = dom_VEF.volumes();
256
257 const int nbegin = the_edge.num_premiere_face();
258 const int nend = nbegin + the_edge.nb_faces();
259 const int dim = Objet_U::dimension;
260
261 CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
262 CDoubleArrView volumes = tab_volumes.view_ro();
263 CIntTabView neighbour_face = tab_neighbour_face.view_ro();
264 DoubleArrView strain_invariant = static_cast<ArrOfDouble&>(tab_strain_invariant).view_wo();
265 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(nbegin, nend), KOKKOS_LAMBDA(const int nface)
266 {
267 const int f0 = neighbour_face(nface, 0);
268 const int f1 = neighbour_face(nface, 1);
269
270 const double inv_vol = 1./(volumes(f0) + volumes(f1));
271 const double vol0 = volumes(f0)*inv_vol;
272 const double vol1 = volumes(f1)*inv_vol;
273
274 // gradient_elem(elem_num, dimension, dimension)
275 const double du_dx = vol0*gradient_elem(f0, 0, 0) + vol1*gradient_elem(f1, 0, 0);
276 const double du_dy = vol0*gradient_elem(f0, 0, 1) + vol1*gradient_elem(f1, 0, 1);
277 const double dv_dx = vol0*gradient_elem(f0, 1, 0) + vol1*gradient_elem(f1, 1, 0);
278 const double dv_dy = vol0*gradient_elem(f0, 1, 1) + vol1*gradient_elem(f1, 1, 1);
279
280 double SijSij = du_dx*du_dx
281 +0.5*(du_dy+dv_dx)*(du_dy+dv_dx)
282 +dv_dy*dv_dy ;
283
284 if (dim == 3)
285 {
286 const double du_dz = vol0*gradient_elem(f0, 0, 2) + vol1*gradient_elem(f1, 0, 2);
287 const double dv_dz = vol0*gradient_elem(f0, 1, 2) + vol1*gradient_elem(f1, 1, 2);
288 const double dw_dx = vol0*gradient_elem(f0, 2, 0) + vol1*gradient_elem(f1, 2, 0);
289 const double dw_dy = vol0*gradient_elem(f0, 2, 1) + vol1*gradient_elem(f1, 2, 1);
290 const double dw_dz = vol0*gradient_elem(f0, 2, 2) + vol1*gradient_elem(f1, 2, 2);
291
292 SijSij += 0.5*(du_dz+dw_dx)*(du_dz+dw_dx)
293 + 0.5*(dv_dz+dw_dy)*(dv_dz+dw_dy)
294 + dw_dz*dw_dz ;
295 }
296
297 strain_invariant(nface) = Kokkos::sqrt(2*SijSij);
298 });
299 end_gpu_timer(__KERNEL_NAME__);
300}
301
302
304 const Front_VF& the_edge,
305 const DoubleTab& tab_gradient_elem,
306 DoubleTab& tab_strain_invariant) const
307{
308 const IntTab& tab_neighbour_face = dom_VEF.face_voisins();
309
310 const int nbegin = the_edge.num_premiere_face();
311 const int nend = nbegin + the_edge.nb_faces();
312 const int dim = Objet_U::dimension;
313
314 CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
315 CIntTabView neighbour_face = tab_neighbour_face.view_ro();
316 DoubleArrView strain_invariant = static_cast<ArrOfDouble&>(tab_strain_invariant).view_wo();
317 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(nbegin, nend), KOKKOS_LAMBDA(const int nface)
318 {
319 const int f0 = neighbour_face(nface, 0);
320
321 // gradient_elem(elem_num, dimension, dimension)
322 const double du_dx = gradient_elem(f0, 0, 0);
323 const double du_dy = gradient_elem(f0, 0, 1);
324 const double dv_dx = gradient_elem(f0, 1, 0);
325 const double dv_dy = gradient_elem(f0, 1, 1);
326
327 double SijSij = du_dx*du_dx
328 +0.5*(du_dy+dv_dx)*(du_dy+dv_dx)
329 +dv_dy*dv_dy ;
330
331 if (dim == 3)
332 {
333 const double du_dz = gradient_elem(f0, 0, 2);
334 const double dv_dz = gradient_elem(f0, 1, 2);
335 const double dw_dx = gradient_elem(f0, 2, 0);
336 const double dw_dy = gradient_elem(f0, 2, 1);
337 const double dw_dz = gradient_elem(f0, 2, 2);
338
339 SijSij += 0.5*(du_dz+dw_dx)*(du_dz+dw_dx)
340 + 0.5*(dv_dz+dw_dy)*(dv_dz+dw_dy)
341 + dw_dz*dw_dz ;
342 }
343
344 strain_invariant(nface) = Kokkos::sqrt(2*SijSij);
345 });
346 end_gpu_timer(__KERNEL_NAME__);
347}
348
349// Même fonction que périodique !
351 const DoubleTab& tab_grad_elem,
352 DoubleTab& tab_strain_invariant) const
353{
354 const IntTab& tab_neighbour_face = dom_VEF.face_voisins();
355 const DoubleVect& tab_volumes = dom_VEF.volumes();
356 const int first_internal_face = dom_VEF.premiere_face_int();
357 const int nb_faces = dom_VEF.nb_faces();
358 const int dim = Objet_U::dimension;
359
360 CDoubleTabView3 grad_elem = tab_grad_elem.view_ro<3>();
361 CDoubleArrView volumes = tab_volumes.view_ro();
362 CIntTabView neighbour_face = tab_neighbour_face.view_ro();
363 DoubleArrView strain_invariant = static_cast<ArrOfDouble&>(tab_strain_invariant).view_wo();
364 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(first_internal_face, nb_faces), KOKKOS_LAMBDA(const int nface)
365 {
366 const int f0 = neighbour_face(nface, 0);
367 const int f1 = neighbour_face(nface, 1);
368
369 const double inv_vol = 1./(volumes(f0) + volumes(f1));
370 const double vol0 = volumes(f0)*inv_vol;
371 const double vol1 = volumes(f1)*inv_vol;
372
373 // format: grad_elem(elem_num, dimension, dimension)
374 const double du_dx = vol0*grad_elem(f0, 0, 0) + vol1*grad_elem(f1, 0, 0);
375 const double du_dy = vol0*grad_elem(f0, 0, 1) + vol1*grad_elem(f1, 0, 1);
376 const double dv_dx = vol0*grad_elem(f0, 1, 0) + vol1*grad_elem(f1, 1, 0);
377 const double dv_dy = vol0*grad_elem(f0, 1, 1) + vol1*grad_elem(f1, 1, 1);
378
379 double SijSij = du_dx*du_dx
380 +0.5*(du_dy+dv_dx)*(du_dy+dv_dx)
381 +dv_dy*dv_dy ;
382
383 if (dim == 3)
384 {
385 const double du_dz = vol0*grad_elem(f0, 0, 2) + vol1*grad_elem(f1, 0, 2);
386 const double dv_dz = vol0*grad_elem(f0, 1, 2) + vol1*grad_elem(f1, 1, 2);
387 const double dw_dx = vol0*grad_elem(f0, 2, 0) + vol1*grad_elem(f1, 2, 0);
388 const double dw_dy = vol0*grad_elem(f0, 2, 1) + vol1*grad_elem(f1, 2, 1);
389 const double dw_dz = vol0*grad_elem(f0, 2, 2) + vol1*grad_elem(f1, 2, 2);
390
391 SijSij += 0.5*(du_dz+dw_dx)*(du_dz+dw_dx)
392 + 0.5*(dv_dz+dw_dy)*(dv_dz+dw_dy)
393 + dw_dz*dw_dz ;
394 }
395
396 strain_invariant(nface) = Kokkos::sqrt(2*SijSij);
397 });
398 end_gpu_timer(__KERNEL_NAME__);
399}
static DoubleTab & calcul_gradient(const DoubleTab &, DoubleTab &, const Domaine_Cl_VEF &)
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
class Domaine_VEF
Definition Domaine_VEF.h:54
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
double volumes(int i) const
Definition Domaine_VF.h:113
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
int nb_elem_tot() const
int nb_front_Cl() const
class Front_VF
Definition Front_VF.h:36
int nb_faces() const
Definition Front_VF.h:53
int num_premiere_face() const
Definition Front_VF.h:63
static int dimension
Definition Objet_U.h:99
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, ConstView< _TYPE_, _SHAPE_ > > view_ro() const
Definition TRUSTTab.h:261
void compute_enstrophy(const Domaine_VEF &, const Domaine_Cl_VEF &, const DoubleTab &velocity, DoubleTab &enstrophy) const
void antisym_loop_internal_faces_strain_invariant(const Domaine_VEF &dom_VEF, const DoubleTab &gradient_elem, DoubleTab &strain_invariant) const
void antisym_loop_edges_general_enstrophy(const Domaine_VEF &dom_VEF, const Front_VF &the_edge, const DoubleTab &gradient_elem, DoubleTab &enstrophy) const
void antisym_loop_internal_faces_enstrophy(const Domaine_VEF &dom_VEF, const DoubleTab &gradient_elem, DoubleTab &enstrophy) const
void compute_strain_invariant(const Domaine_VEF &, const Domaine_Cl_VEF &, const DoubleTab &velocity, DoubleTab &strain_invariant) const
void antisym_loop_edges_periodiqueBC_strain_invariant(const Domaine_VEF &dom_VEF, const Front_VF &the_edge, const DoubleTab &gradient_elem, DoubleTab &strain_invariant) const
void antisym_loop_edge_faces_enstrophy(const Domaine_VEF &dom_VEF, const Domaine_Cl_VEF &dom_BC_VEF, const DoubleTab &gradient_elem, DoubleTab &enstrophy) const
void antisym_loop_edges_periodiqueBC_enstrophy(const Domaine_VEF &dom_VEF, const Front_VF &the_edge, const DoubleTab &gradient_elem, DoubleTab &enstrophy) const
void antisym_loop_edge_faces_strain_invariant(const Domaine_VEF &dom_VEF, const Domaine_Cl_VEF &dom_BC_VEF, const DoubleTab &gradient_elem, DoubleTab &strain_invariant) const
void antisym_loop_edges_general_strain_invariant(const Domaine_VEF &dom_VEF, const Front_VF &the_edge, const DoubleTab &gradient_elem, DoubleTab &strain_invariant) const