TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Grid_Level_Data_template.tpp
1/****************************************************************************
2* Copyright (c) 2025, 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 Grid_Level_Data_TPP_H
17#define Grid_Level_Data_TPP_H
18
19#ifdef DUMP_FACES_COEFFICIENTS
20#include <IJK_lata_dump.h>
21#endif
22
23template<typename _TYPE_>
24Grid_Level_Data_template<_TYPE_>::Grid_Level_Data_template()
25{
26 local_delta_xyz_.dimensionner(3);
27 ghost_size_ = 0;
28}
29
30
31// Initialize the data structures (allocate memory and setup temporary storage structures)
32template<typename _TYPE_>
33void Grid_Level_Data_template<_TYPE_>::initialize(const Domaine_IJK& dom, int ghost, int additional_k_layers)
34{
35 domaine_ijk_ = dom;
36 perio_k_= dom.get_periodic_flag(DIRECTION_K);
37 ghost_size_ = ghost;
38 if (IJK_Shear_Periodic_helpler::defilement_==1)
39 {
40 ijk_rho_.allocate(domaine_ijk_, Domaine_IJK::ELEM, ghost, 0 , 1, "IJK_RHO", false, 2, IJK_Shear_Periodic_helpler::rho_vap_ref_for_poisson_, IJK_Shear_Periodic_helpler::rho_liq_ref_for_poisson_);
41 }
42 else
43 {
44 ijk_rho_.allocate(domaine_ijk_, Domaine_IJK::ELEM, ghost);
45 }
46 ijk_rho_.data() = 1.;
47 // Allocate the array of coefficients at faces with size "elements".
48 // Therefore, if the domain is not periodic, at the right end of the domain,
49 // the wall coefficient is stored in a "ghost" cell.
50 // This trick allows to have the same stride in j and k for all arrays.
51 ijk_faces_coefficients_.allocate(domaine_ijk_, Domaine_IJK::ELEM, ghost, 0 /* add.k layers */, 4 /* components */);
52 ijk_faces_coefficients_.data() = 1.;
53 ijk_x_.allocate(domaine_ijk_, Domaine_IJK::ELEM, ghost, additional_k_layers);
54 ijk_x_.data() = 0.;
55 ijk_rhs_.allocate(domaine_ijk_, Domaine_IJK::ELEM, ghost);
56 ijk_rhs_.data() = 0.;
57 ijk_residue_.allocate(domaine_ijk_, Domaine_IJK::ELEM, ghost);
58 ijk_residue_.data() = 0.;
59
60 for (int dir = 0; dir < 3; dir++)
61 domaine_ijk_.get_local_mesh_delta(dir, ghost, local_delta_xyz_[dir]);
62 for (int i = 0; i < 3; i++)
63 uniform_[i] = domaine_ijk_.is_uniform(i);
64
65}
66
67// Fill the coeffs_faces field (see ijk_faces_coefficients_) for the requested grid_level
68// nb_ghost_layers is the number of valid ghost layers provided in ijk_rho_levels_ and as input to the jacobi
69// smoother (at least 1).
70template<typename _TYPE_>
71void Grid_Level_Data_template<_TYPE_>::compute_faces_coefficients_from_rho()
72{
73 int flag = uniform_[0] ? 1 : 0;
74 flag += uniform_[1] ? 2 : 0;
75 flag += uniform_[2] ? 4 : 0;
76 switch(flag)
77 {
78 case 3: // variable in k
79 compute_faces_coefficients_from_rho_cst_i_cst_j_var_k();
80 break;
81 case 7: // completely uniform
82 compute_faces_coefficients_from_rho_cst_i_cst_j_cst_k();
83 break;
84 default:
85 Cerr << "Grid_Level_Data::compute_faces_coefficients_from_rho() flag=" << flag << " not implemented" << finl;
86 Process::exit();
87 }
88#ifdef DUMP_FACES_COEFFICIENTS
89 {
90 IJK_Field_template<_TYPE_,TRUSTArray<_TYPE_>> c[3];
91 for (int dir = 0; dir < 3; dir++)
92 {
93 Domaine_IJK::Localisation loc = (dir==0)?Domaine_IJK::FACES_I:((dir==1)?Domaine_IJK::FACES_J:Domaine_IJK::FACES_K);
94 IJK_Field_template<_TYPE_,TRUSTArray<_TYPE_>>& f = c[dir];
95 f.allocate(domaine_ijk_, loc, 1);
96 for (int k = 0; k < f.nk(); k++)
97 for (int j = 0; j < f.nj(); j++)
98 for (int i = 0; i < f.ni(); i++)
99 f(i,j,k)=ijk_faces_coefficients_(i,j,k,dir);
100 f.echange_espace_virtuel(1); // to update periodic faces
101 }
102 IJK_Field_template<_TYPE_,TRUSTArray<_TYPE_>> e;
103 e.allocate(domaine_ijk_, Domaine_IJK::ELEM, 0);
104 for (int k = 0; k < e.nk(); k++)
105 for (int j = 0; j < e.nj(); j++)
106 for (int i = 0; i < e.ni(); i++)
107 e(i,j,k)=ijk_faces_coefficients_(i,j,k,3);
108 static int step = 0;
109 Nom prefix = Nom("Grid_coefficients_")+Nom(typeid(_TYPE_).name())
110 +Nom(domaine_ijk_.get_nb_elem_tot(DIRECTION_I))+Nom("_")
111 + Nom(domaine_ijk_.get_nb_elem_tot(DIRECTION_J))+Nom("_")
112 + Nom(domaine_ijk_.get_nb_elem_tot(DIRECTION_K));
113
114 dumplata_vector(prefix+Nom("_faces.lata"), "val", c[0],c[1],c[2],step);
115 dumplata_scalar(prefix+Nom("_elem.lata"), "val", e,step);
116 step++;
117 }
118#endif
119}
120
121template<typename _TYPE_>
122void Grid_Level_Data_template<_TYPE_>::compute_faces_coefficients_from_inv_rho()
123{
124 int flag = uniform_[0] ? 1 : 0;
125 flag += uniform_[1] ? 2 : 0;
126 flag += uniform_[2] ? 4 : 0;
127 switch(flag)
128 {
129 case 3: // variable in k
130 compute_faces_coefficients_from_inv_rho_cst_i_cst_j_var_k();
131 break;
132 case 7: // completely uniform
133 compute_faces_coefficients_from_inv_rho_cst_i_cst_j_cst_k();
134 break;
135 default:
136 Cerr << "Grid_Level_Data::compute_faces_coefficients_from_inv_rho() flag=" << flag << " not implemented" << finl;
137 Process::exit();
138 }
139#ifdef DUMP_FACES_COEFFICIENTS
140 {
141 IJK_Field_template<_TYPE_,TRUSTArray<_TYPE_>> c[3];
142 for (int dir = 0; dir < 3; dir++)
143 {
144 Domaine_IJK::Localisation loc = (dir==0)?Domaine_IJK::FACES_I:((dir==1)?Domaine_IJK::FACES_J:Domaine_IJK::FACES_K);
145 IJK_Field_template<_TYPE_,TRUSTArray<_TYPE_>>& f = c[dir];
146 f.allocate(domaine_ijk_, loc, 1);
147 for (int k = 0; k < f.nk(); k++)
148 for (int j = 0; j < f.nj(); j++)
149 for (int i = 0; i < f.ni(); i++)
150 f(i,j,k)=ijk_faces_coefficients_(i,j,k,dir);
151 f.echange_espace_virtuel(1); // to update periodic faces
152 }
153 IJK_Field_template<_TYPE_,TRUSTArray<_TYPE_>> e;
154 e.allocate(domaine_ijk_, Domaine_IJK::ELEM, 0);
155 for (int k = 0; k < e.nk(); k++)
156 for (int j = 0; j < e.nj(); j++)
157 for (int i = 0; i < e.ni(); i++)
158 e(i,j,k)=ijk_faces_coefficients_(i,j,k,3);
159 static int step = 0;
160 Nom prefix = Nom("Grid_coefficients_")++Nom(typeid(_TYPE_).name())
161 +Nom(domaine_ijk_.get_nb_elem_tot(DIRECTION_I))+Nom("_")
162 + Nom(domaine_ijk_.get_nb_elem_tot(DIRECTION_J))+Nom("_")
163 + Nom(domaine_ijk_.get_nb_elem_tot(DIRECTION_K));
164
165 dumplata_vector(prefix+Nom("_faces.lata"), "val", c[0],c[1],c[2],step);
166 dumplata_scalar(prefix+Nom("_elem.lata"), "val", e,step);
167 step++;
168 }
169#endif
170}
171
172template<typename _TYPE_>
173void Grid_Level_Data_template<_TYPE_>::compute_faces_coefficients_from_rho_cst_i_cst_j_cst_k()
174{
175 // mesh is uniform in i direction: define a constant value
176 const _TYPE_ delta_i = (_TYPE_)local_delta_xyz_[0][0];
177 const _TYPE_ i_delta_i = (_TYPE_)(1. / local_delta_xyz_[0][0]);
178#define DELTA_i delta_i
179#define invDELTA_i i_delta_i
180 // mesh is uniform in j direction: define a constant value
181 const _TYPE_ delta_j = (_TYPE_)local_delta_xyz_[1][0];
182 const _TYPE_ i_delta_j = (_TYPE_)(1. / local_delta_xyz_[1][0]);
183#define DELTA_j delta_j
184#define invDELTA_j i_delta_j
185 // mesh is uniform in k direction: define a constant value
186 const _TYPE_ delta_k = (_TYPE_)local_delta_xyz_[2][0];
187 const _TYPE_ i_delta_k = (_TYPE_)(1. / local_delta_xyz_[2][0]);
188#define DELTA_k delta_k
189#define invDELTA_k i_delta_k
190
191 IJ_layout layout(ijk_rho_);
192
193 const int nb_ghost_layers_rho = ijk_rho_.ghost();
194 const int ghost = nb_ghost_layers_rho - 1;
195 const int imin = - ghost;
196 const int imax = ijk_rho_.ni() + ghost;
197 const int jmin = - ghost;
198 const int jmax = ijk_rho_.nj() + ghost;
199 const int kmin = - ghost;
200 const int kmax = ijk_rho_.nk() + ghost;
201
202 for (int k = kmin; k < kmax + 1; k++)
203 {
204
205 const bool compute_k_direction_only = (k == kmax);
206 if (!compute_k_direction_only)
207 {
208
209 // Compute layer of i faces
210 const _TYPE_ *src_rho = ijk_rho_.k_layer(k);
211 {
212 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 0);
213 for (int j = jmin; j < jmax; j++)
214 {
215 for (int i = imin; i < imax+1; i++)
216 {
217 const _TYPE_ rho = (_TYPE_)(2. / (layout(src_rho, i-1, j) + layout(src_rho, i, j)));
218 const _TYPE_ f = invDELTA_i * DELTA_j * DELTA_k;
219 layout(coeff, i, j) = rho * f;
220 }
221 }
222 }
223 // Compute layer of j faces
224 {
225 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 1);
226 for (int j = jmin; j < jmax+1; j++)
227 {
228 for (int i = imin; i < imax; i++)
229 {
230 const _TYPE_ rho = _TYPE_(2. / (layout(src_rho, i, j-1) + layout(src_rho, i, j)));
231 const _TYPE_ f = DELTA_i * invDELTA_j * DELTA_k;
232 layout(coeff, i, j) = rho * f;
233 }
234 }
235 }
236 }
237
238 {
239 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 2);
240 // Wall boundary condition:
241 const int global_k_pos = k + domaine_ijk_.get_offset_local(DIRECTION_K);
242 // Number of elements in the z direction (equal to index in k direction of the faces on the wall)
243 const int global_k_size = domaine_ijk_.get_nb_elem_tot(DIRECTION_K);
244 // the face that are located on the walls in k direction have position 0 and global_k_size:
245 if (!perio_k_ && (global_k_pos <= 0 || global_k_pos >= global_k_size))
246 {
247 // layer below or equal to first layer in the mesh
248 for (int i = imin; i < imax; i++)
249 for (int j = jmin; j < jmax; j++)
250 layout(coeff, i, j) = 0.;
251 }
252 else
253 {
254 const _TYPE_ *src_rho_left = ijk_rho_.k_layer(k-1);
255 const _TYPE_ *src_rho_right = ijk_rho_.k_layer(k);
256 for (int j = jmin; j < jmax; j++)
257 {
258 for (int i = imin; i < imax; i++)
259 {
260 const _TYPE_ rho = (_TYPE_)(2. / (layout(src_rho_left, i, j) + layout(src_rho_right, i, j)));
261 const _TYPE_ f = DELTA_i * DELTA_j * invDELTA_k;
262 layout(coeff, i, j) = rho * f;
263 }
264 }
265 }
266 }
267 if (k > kmin)
268 {
269 const _TYPE_ *coeff_i = ijk_faces_coefficients_.k_layer(k-1, 0);
270 const _TYPE_ *coeff_j = ijk_faces_coefficients_.k_layer(k-1, 1);
271 const _TYPE_ *coeff_k = ijk_faces_coefficients_.k_layer(k-1, 2);
272 const _TYPE_ *coeff_kplus = ijk_faces_coefficients_.k_layer(k, 2);
273 _TYPE_ *coeff_sum = ijk_faces_coefficients_.k_layer(k-1, 3);
274
275 // Sum of coefficients for 6 faces of each element
276 for (int i = imin; i < imax; i++)
277 for (int j = jmin; j < jmax; j++)
278 {
279 _TYPE_ s =
280 layout(coeff_i, i, j) + layout(coeff_i, i+1, j)
281 + layout(coeff_j, i, j) + layout(coeff_j, i, j+1)
282 + layout(coeff_k, i, j) + layout(coeff_kplus, i, j);
283 // ghost cells beyond wall boundaries will have zero sum, crashes the division in jacobi,
284 // so replace zero with dummy non zero value:
285 if (s == 0.)
286 s = 1.;
287 layout(coeff_sum, i, j) = s;
288 }
289 }
290 }
291#undef DELTA_i
292#undef DELTA_j
293#undef DELTA_k
294#undef invDELTA_i
295#undef invDELTA_j
296#undef invDELTA_k
297}
298
299template<typename _TYPE_>
300void Grid_Level_Data_template<_TYPE_>::compute_faces_coefficients_from_inv_rho_cst_i_cst_j_cst_k()
301{
302 // mesh is uniform in i direction: define a constant value
303 const _TYPE_ delta_i = (_TYPE_)local_delta_xyz_[0][0];
304 const _TYPE_ i_delta_i = (_TYPE_)(1. / local_delta_xyz_[0][0]);
305#define DELTA_i delta_i
306#define invDELTA_i i_delta_i
307 // mesh is uniform in j direction: define a constant value
308 const _TYPE_ delta_j = (_TYPE_)local_delta_xyz_[1][0];
309 const _TYPE_ i_delta_j = (_TYPE_)(1. / local_delta_xyz_[1][0]);
310#define DELTA_j delta_j
311#define invDELTA_j i_delta_j
312 // mesh is uniform in k direction: define a constant value
313 const _TYPE_ delta_k = (_TYPE_)local_delta_xyz_[2][0];
314 const _TYPE_ i_delta_k = (_TYPE_)(1. / local_delta_xyz_[2][0]);
315#define DELTA_k delta_k
316#define invDELTA_k i_delta_k
317
318 IJ_layout layout(ijk_rho_);
319
320 const int nb_ghost_layers_rho = ijk_rho_.ghost();
321 const int ghost = nb_ghost_layers_rho - 1;
322 const int imin = - ghost;
323 const int imax = ijk_rho_.ni() + ghost;
324 const int jmin = - ghost;
325 const int jmax = ijk_rho_.nj() + ghost;
326 const int kmin = - ghost;
327 const int kmax = ijk_rho_.nk() + ghost;
328
329 for (int k = kmin; k < kmax + 1; k++)
330 {
331
332 const bool compute_k_direction_only = (k == kmax);
333 if (!compute_k_direction_only)
334 {
335
336 // Compute layer of i faces
337 const _TYPE_ *src_rho = ijk_rho_.k_layer(k);
338 {
339 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 0);
340 for (int j = jmin; j < jmax; j++)
341 {
342 for (int i = imin; i < imax+1; i++)
343 {
344 const _TYPE_ rho = (_TYPE_)(0.5 * (layout(src_rho, i-1, j) + layout(src_rho, i, j)));
345 const _TYPE_ f = invDELTA_i * DELTA_j * DELTA_k;
346 layout(coeff, i, j) = rho * f;
347 }
348 }
349 }
350 // Compute layer of j faces
351 {
352 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 1);
353 for (int j = jmin; j < jmax+1; j++)
354 {
355 for (int i = imin; i < imax; i++)
356 {
357 const _TYPE_ rho = (_TYPE_)(0.5 * (layout(src_rho, i, j-1) + layout(src_rho, i, j)));
358 const _TYPE_ f = DELTA_i * invDELTA_j * DELTA_k;
359 layout(coeff, i, j) = rho * f;
360 }
361 }
362 }
363 }
364
365 {
366 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 2);
367 // Wall boundary condition:
368 const int global_k_pos = k + domaine_ijk_.get_offset_local(DIRECTION_K);
369 // Number of elements in the z direction (equal to index in k direction of the faces on the wall)
370 const int global_k_size = domaine_ijk_.get_nb_elem_tot(DIRECTION_K);
371 // the face that are located on the walls in k direction have position 0 and global_k_size:
372 if (!perio_k_ && (global_k_pos <= 0 || global_k_pos >= global_k_size))
373 {
374 // layer below or equal to first layer in the mesh
375 for (int i = imin; i < imax; i++)
376 for (int j = jmin; j < jmax; j++)
377 layout(coeff, i, j) = 0.;
378 }
379 else
380 {
381 const _TYPE_ *src_rho_left = ijk_rho_.k_layer(k-1);
382 const _TYPE_ *src_rho_right = ijk_rho_.k_layer(k);
383 for (int j = jmin; j < jmax; j++)
384 {
385 for (int i = imin; i < imax; i++)
386 {
387 const _TYPE_ rho = (_TYPE_)(0.5 * (layout(src_rho_left, i, j) + layout(src_rho_right, i, j)));
388 const _TYPE_ f = DELTA_i * DELTA_j * invDELTA_k;
389 layout(coeff, i, j) = rho * f;
390 }
391 }
392 }
393 }
394 if (k > kmin)
395 {
396 const _TYPE_ *coeff_i = ijk_faces_coefficients_.k_layer(k-1, 0);
397 const _TYPE_ *coeff_j = ijk_faces_coefficients_.k_layer(k-1, 1);
398 const _TYPE_ *coeff_k = ijk_faces_coefficients_.k_layer(k-1, 2);
399 const _TYPE_ *coeff_kplus = ijk_faces_coefficients_.k_layer(k, 2);
400 _TYPE_ *coeff_sum = ijk_faces_coefficients_.k_layer(k-1, 3);
401
402 // Sum of coefficients for 6 faces of each element
403 for (int i = imin; i < imax; i++)
404 for (int j = jmin; j < jmax; j++)
405 {
406 _TYPE_ s =
407 layout(coeff_i, i, j) + layout(coeff_i, i+1, j)
408 + layout(coeff_j, i, j) + layout(coeff_j, i, j+1)
409 + layout(coeff_k, i, j) + layout(coeff_kplus, i, j);
410 // ghost cells beyond wall boundaries will have zero sum, crashes the division in jacobi,
411 // so replace zero with dummy non zero value:
412 if (s == 0.)
413 s = 1.;
414 layout(coeff_sum, i, j) = s;
415 }
416 }
417 }
418#undef DELTA_i
419#undef DELTA_j
420#undef DELTA_k
421#undef invDELTA_i
422#undef invDELTA_j
423#undef invDELTA_k
424}
425
426template<typename _TYPE_>
427void Grid_Level_Data_template<_TYPE_>::compute_faces_coefficients_from_rho_cst_i_cst_j_var_k()
428{
429 // mesh is uniform in i direction: define a constant value
430 const _TYPE_ delta_i = (_TYPE_)local_delta_xyz_[0][0];
431 const _TYPE_ i_delta_i = (_TYPE_)(1. / local_delta_xyz_[0][0]);
432#define DELTA_i delta_i
433#define invDELTA_i i_delta_i
434 // mesh is uniform in j direction: define a constant value
435 const _TYPE_ delta_j = (_TYPE_)local_delta_xyz_[1][0];
436 const _TYPE_ i_delta_j = (_TYPE_)(1. / local_delta_xyz_[1][0]);
437#define DELTA_j delta_j
438#define invDELTA_j i_delta_j
439 // mesh has variable size in k direction: compute the inverse of
440 // space between cell centers and create macro to point to appropriate
441 // value in the arrays
442 IJKArray_with_ghost<_TYPE_,TRUSTArray<_TYPE_>> delta_k_array;
443 IJKArray_with_ghost<_TYPE_,TRUSTArray<_TYPE_>> i_delta_k_array;
444 {
445 const int n = ijk_rho_.nk();
446 const int ghost = ijk_rho_.ghost();
447 delta_k_array.resize(n, ghost);
448 i_delta_k_array.resize(n, ghost);
449 const int imin = -ghost + 1;
450 const int imax = n + ghost;
451 delta_k_array[-ghost] = (_TYPE_)local_delta_xyz_[2][-ghost];
452 for (int i = imin; i < imax; i++)
453 {
454 delta_k_array[i] = (_TYPE_)local_delta_xyz_[2][i];
455 _TYPE_ x = delta_k_array[i-1];
456 _TYPE_ y = delta_k_array[i];
457 i_delta_k_array[i] = (_TYPE_)(2. / (x+y));
458 }
459 }
460#define DELTA_k delta_k_array[k]
461#define invDELTA_k i_delta_k_array[k]
462
463 IJ_layout layout(ijk_rho_);
464
465 const int nb_ghost_layers_rho = ijk_rho_.ghost();
466 const int ghost = nb_ghost_layers_rho - 1;
467 const int imin = - ghost;
468 const int imax = ijk_rho_.ni() + ghost;
469 const int jmin = - ghost;
470 const int jmax = ijk_rho_.nj() + ghost;
471 const int kmin = - ghost;
472 const int kmax = ijk_rho_.nk() + ghost;
473
474 for (int k = kmin; k < kmax + 1; k++)
475 {
476
477 const bool compute_k_direction_only = (k == kmax);
478 if (!compute_k_direction_only)
479 {
480
481 // Compute layer of i faces
482 const _TYPE_ *src_rho = ijk_rho_.k_layer(k);
483 {
484 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 0);
485 for (int j = jmin; j < jmax; j++)
486 {
487 for (int i = imin; i < imax+1; i++)
488 {
489 const _TYPE_ rho = (_TYPE_)(2. / (layout(src_rho, i-1, j) + layout(src_rho, i, j)));
490 const _TYPE_ f = invDELTA_i * DELTA_j * DELTA_k;
491 layout(coeff, i, j) = rho * f;
492 }
493 }
494 }
495 // Compute layer of j faces
496 {
497 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 1);
498 for (int j = jmin; j < jmax+1; j++)
499 {
500 for (int i = imin; i < imax; i++)
501 {
502 const _TYPE_ rho = _TYPE_(2. / (layout(src_rho, i, j-1) + layout(src_rho, i, j)));
503 const _TYPE_ f = DELTA_i * invDELTA_j * DELTA_k;
504 layout(coeff, i, j) = rho * f;
505 }
506 }
507 }
508 }
509
510 {
511 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 2);
512 // Wall boundary condition:
513 const int global_k_pos = k + domaine_ijk_.get_offset_local(DIRECTION_K);
514 // Number of elements in the z direction (equal to index in k direction of the faces on the wall)
515 const int global_k_size = domaine_ijk_.get_nb_elem_tot(DIRECTION_K);
516 // the face that are located on the walls in k direction have position 0 and global_k_size:
517 if (!perio_k_ && (global_k_pos <= 0 || global_k_pos >= global_k_size))
518 {
519 // layer below or equal to first layer in the mesh
520 for (int i = imin; i < imax; i++)
521 for (int j = jmin; j < jmax; j++)
522 layout(coeff, i, j) = 0.;
523 }
524 else
525 {
526 const _TYPE_ *src_rho_left = ijk_rho_.k_layer(k-1);
527 const _TYPE_ *src_rho_right = ijk_rho_.k_layer(k);
528 for (int j = jmin; j < jmax; j++)
529 {
530 for (int i = imin; i < imax; i++)
531 {
532 const _TYPE_ rho = (_TYPE_)(2. / (layout(src_rho_left, i, j) + layout(src_rho_right, i, j)));
533 const _TYPE_ f = DELTA_i * DELTA_j * invDELTA_k;
534 layout(coeff, i, j) = rho * f;
535 }
536 }
537 }
538 }
539 if (k > kmin)
540 {
541 const _TYPE_ *coeff_i = ijk_faces_coefficients_.k_layer(k-1, 0);
542 const _TYPE_ *coeff_j = ijk_faces_coefficients_.k_layer(k-1, 1);
543 const _TYPE_ *coeff_k = ijk_faces_coefficients_.k_layer(k-1, 2);
544 const _TYPE_ *coeff_kplus = ijk_faces_coefficients_.k_layer(k, 2);
545 _TYPE_ *coeff_sum = ijk_faces_coefficients_.k_layer(k-1, 3);
546
547 // Sum of coefficients for 6 faces of each element
548 for (int i = imin; i < imax; i++)
549 for (int j = jmin; j < jmax; j++)
550 {
551 _TYPE_ s =
552 layout(coeff_i, i, j) + layout(coeff_i, i+1, j)
553 + layout(coeff_j, i, j) + layout(coeff_j, i, j+1)
554 + layout(coeff_k, i, j) + layout(coeff_kplus, i, j);
555 // ghost cells beyond wall boundaries will have zero sum, crashes the division in jacobi,
556 // so replace zero with dummy non zero value:
557 if (s == 0.)
558 s = 1.;
559 layout(coeff_sum, i, j) = s;
560 }
561 }
562 }
563#undef DELTA_i
564#undef DELTA_j
565#undef DELTA_k
566#undef invDELTA_i
567#undef invDELTA_j
568#undef invDELTA_k
569}
570
571template<typename _TYPE_>
572void Grid_Level_Data_template<_TYPE_>::compute_faces_coefficients_from_inv_rho_cst_i_cst_j_var_k()
573{
574 // mesh is uniform in i direction: define a constant value
575 const _TYPE_ delta_i = (_TYPE_)local_delta_xyz_[0][0];
576 const _TYPE_ i_delta_i = (_TYPE_)(1. / local_delta_xyz_[0][0]);
577#define DELTA_i delta_i
578#define invDELTA_i i_delta_i
579 // mesh is uniform in j direction: define a constant value
580 const _TYPE_ delta_j = (_TYPE_)local_delta_xyz_[1][0];
581 const _TYPE_ i_delta_j = (_TYPE_)(1. / local_delta_xyz_[1][0]);
582#define DELTA_j delta_j
583#define invDELTA_j i_delta_j
584 // mesh has variable size in k direction: compute the inverse of
585 // space between cell centers and create macro to point to appropriate
586 // value in the arrays
587 IJKArray_with_ghost<_TYPE_,TRUSTArray<_TYPE_>> delta_k_array;
588 IJKArray_with_ghost<_TYPE_,TRUSTArray<_TYPE_>> i_delta_k_array;
589 {
590 const int n = ijk_rho_.nk();
591 const int ghost = ijk_rho_.ghost();
592 delta_k_array.resize(n, ghost);
593 i_delta_k_array.resize(n, ghost);
594 const int imin = -ghost + 1;
595 const int imax = n + ghost;
596 delta_k_array[-ghost] = (_TYPE_)local_delta_xyz_[2][-ghost];
597 for (int i = imin; i < imax; i++)
598 {
599 delta_k_array[i] = (_TYPE_)local_delta_xyz_[2][i];
600 _TYPE_ x = delta_k_array[i-1];
601 _TYPE_ y = delta_k_array[i];
602 i_delta_k_array[i] = (_TYPE_)(2. / (x+y));
603 }
604 }
605#define DELTA_k delta_k_array[k]
606#define invDELTA_k i_delta_k_array[k]
607
608 IJ_layout layout(ijk_rho_);
609
610 const int nb_ghost_layers_rho = ijk_rho_.ghost();
611 const int ghost = nb_ghost_layers_rho - 1;
612 const int imin = - ghost;
613 const int imax = ijk_rho_.ni() + ghost;
614 const int jmin = - ghost;
615 const int jmax = ijk_rho_.nj() + ghost;
616 const int kmin = - ghost;
617 const int kmax = ijk_rho_.nk() + ghost;
618
619 for (int k = kmin; k < kmax + 1; k++)
620 {
621
622 const bool compute_k_direction_only = (k == kmax);
623 if (!compute_k_direction_only)
624 {
625
626 // Compute layer of i faces
627 const _TYPE_ *src_rho = ijk_rho_.k_layer(k);
628 {
629 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 0);
630 for (int j = jmin; j < jmax; j++)
631 {
632 for (int i = imin; i < imax+1; i++)
633 {
634 const _TYPE_ rho = (_TYPE_)(0.5 * (layout(src_rho, i-1, j) + layout(src_rho, i, j)));
635 const _TYPE_ f = invDELTA_i * DELTA_j * DELTA_k;
636 layout(coeff, i, j) = rho * f;
637 }
638 }
639 }
640 // Compute layer of j faces
641 {
642 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 1);
643 for (int j = jmin; j < jmax+1; j++)
644 {
645 for (int i = imin; i < imax; i++)
646 {
647 const _TYPE_ rho = (_TYPE_)(0.5 * (layout(src_rho, i, j-1) + layout(src_rho, i, j)));
648 const _TYPE_ f = DELTA_i * invDELTA_j * DELTA_k;
649 layout(coeff, i, j) = rho * f;
650 }
651 }
652 }
653 }
654
655 {
656 _TYPE_ *coeff = ijk_faces_coefficients_.k_layer(k, 2);
657 // Wall boundary condition:
658 const int global_k_pos = k + domaine_ijk_.get_offset_local(DIRECTION_K);
659 // Number of elements in the z direction (equal to index in k direction of the faces on the wall)
660 const int global_k_size = domaine_ijk_.get_nb_elem_tot(DIRECTION_K);
661 // the face that are located on the walls in k direction have position 0 and global_k_size:
662 if (!perio_k_ && (global_k_pos <= 0 || global_k_pos >= global_k_size))
663 {
664 // layer below or equal to first layer in the mesh
665 for (int i = imin; i < imax; i++)
666 for (int j = jmin; j < jmax; j++)
667 layout(coeff, i, j) = 0.;
668 }
669 else
670 {
671 const _TYPE_ *src_rho_left = ijk_rho_.k_layer(k-1);
672 const _TYPE_ *src_rho_right = ijk_rho_.k_layer(k);
673 for (int j = jmin; j < jmax; j++)
674 {
675 for (int i = imin; i < imax; i++)
676 {
677 const _TYPE_ rho = (_TYPE_)(0.5 * (layout(src_rho_left, i, j) + layout(src_rho_right, i, j)));
678 const _TYPE_ f = DELTA_i * DELTA_j * invDELTA_k;
679 layout(coeff, i, j) = rho * f;
680 }
681 }
682 }
683 }
684 if (k > kmin)
685 {
686 const _TYPE_ *coeff_i = ijk_faces_coefficients_.k_layer(k-1, 0);
687 const _TYPE_ *coeff_j = ijk_faces_coefficients_.k_layer(k-1, 1);
688 const _TYPE_ *coeff_k = ijk_faces_coefficients_.k_layer(k-1, 2);
689 const _TYPE_ *coeff_kplus = ijk_faces_coefficients_.k_layer(k, 2);
690 _TYPE_ *coeff_sum = ijk_faces_coefficients_.k_layer(k-1, 3);
691
692 // Sum of coefficients for 6 faces of each element
693 for (int i = imin; i < imax; i++)
694 for (int j = jmin; j < jmax; j++)
695 {
696 _TYPE_ s =
697 layout(coeff_i, i, j) + layout(coeff_i, i+1, j)
698 + layout(coeff_j, i, j) + layout(coeff_j, i, j+1)
699 + layout(coeff_k, i, j) + layout(coeff_kplus, i, j);
700 // ghost cells beyond wall boundaries will have zero sum, crashes the division in jacobi,
701 // so replace zero with dummy non zero value:
702 if (s == 0.)
703 s = 1.;
704 layout(coeff_sum, i, j) = s;
705 }
706 }
707 }
708#undef DELTA_i
709#undef DELTA_j
710#undef DELTA_k
711#undef invDELTA_i
712#undef invDELTA_j
713#undef invDELTA_k
714}
715
716
717// Copy double precision coefficients to single precision:
718template<class _TYPE_> template<class MY_TYPE> std::enable_if_t<std::is_same<MY_TYPE, float>::value, void>
719Grid_Level_Data_template<_TYPE_>::compute_faces_coefficients_from_double_coeffs(const Grid_Level_Data_double& dbl_coeffs)
720{
721 // Warn: strides might differ due to padding for different simd vector sizes.
722 // Therefore, linear data indices might not match, must explicitely loop on i, j, k indices
723 const IJK_Field_double& src = dbl_coeffs.get_faces_coefficients();
724
725 const int nb_ghost_layers_rho = ijk_rho_.ghost();
726 const int ghost = nb_ghost_layers_rho;
727 const int ni = ijk_rho_.ni();
728 const int imin = - ghost + 1;
729 const int jmin = - ghost + 1;
730 const int nj = ijk_rho_.nj();
731 const int kmin = - ghost + 1;
732 const int kmax = ijk_rho_.nk() + ghost;
733 for (int k = kmin; k < kmax; k++)
734 for (int compo = 0; compo < 4; compo++)
735 {
736 int imax = ni + ghost - 1;
737 int jmax = nj + ghost - 1;
738 if (k == (kmax-1) && compo != 2 )
739 {
740 continue; // this data is unused
741 }
742 if (compo == 0)
743 imax++;
744 else if (compo == 1)
745 jmax++;
746
747 for (int j = jmin; j < jmax; j++)
748 for (int i = imin; i < imax; i++)
749 ijk_faces_coefficients_(i, j, k, compo) = (float)src(i,j,k,compo);
750 }
751}
752
753
754#endif
Domaine_IJK
This class encapsulates all the information related to the eulerian mesh for TrioIJK.
Definition Domaine_IJK.h:47
Domaine_IJK::get_periodic_flag
bool get_periodic_flag(int direction) const
Method returns true if periodic in this direction.
Definition Domaine_IJK.h:575
Domaine_IJK::Localisation
Localisation
Localisation sub class.
Definition Domaine_IJK.h:53
Domaine_IJK::FACES_J
@ FACES_J
Definition Domaine_IJK.h:53
Domaine_IJK::FACES_K
@ FACES_K
Definition Domaine_IJK.h:53
Domaine_IJK::FACES_I
@ FACES_I
Definition Domaine_IJK.h:53
Domaine_IJK::ELEM
@ ELEM
Definition Domaine_IJK.h:53
Grid_Level_Data_template::ijk_residue_
IJK_Field_template< _TYPE_, TRUSTArray< _TYPE_ > > ijk_residue_
Definition Grid_Level_Data_template.h:91
Grid_Level_Data_template::compute_faces_coefficients_from_rho_cst_i_cst_j_cst_k
void compute_faces_coefficients_from_rho_cst_i_cst_j_cst_k()
Definition Grid_Level_Data_template.tpp:173
Grid_Level_Data_template::Grid_Level_Data_template
Grid_Level_Data_template()
Definition Grid_Level_Data_template.tpp:24
Grid_Level_Data_template::perio_k_
bool perio_k_
Definition Grid_Level_Data_template.h:69
Grid_Level_Data_template::compute_faces_coefficients_from_rho
void compute_faces_coefficients_from_rho()
Definition Grid_Level_Data_template.tpp:71
Grid_Level_Data_template::ijk_faces_coefficients_
IJK_Field_template< _TYPE_, TRUSTArray< _TYPE_ > > ijk_faces_coefficients_
Definition Grid_Level_Data_template.h:83
Grid_Level_Data_template::compute_faces_coefficients_from_inv_rho_cst_i_cst_j_var_k
void compute_faces_coefficients_from_inv_rho_cst_i_cst_j_var_k()
Definition Grid_Level_Data_template.tpp:572
Grid_Level_Data_template::ijk_rhs_
IJK_Field_template< _TYPE_, TRUSTArray< _TYPE_ > > ijk_rhs_
Definition Grid_Level_Data_template.h:89
Grid_Level_Data_template::uniform_
bool uniform_[3]
Definition Grid_Level_Data_template.h:72
Grid_Level_Data_template::compute_faces_coefficients_from_rho_cst_i_cst_j_var_k
void compute_faces_coefficients_from_rho_cst_i_cst_j_var_k()
Definition Grid_Level_Data_template.tpp:427
Grid_Level_Data_template::compute_faces_coefficients_from_inv_rho
void compute_faces_coefficients_from_inv_rho()
Definition Grid_Level_Data_template.tpp:122
Grid_Level_Data_template::domaine_ijk_
Domaine_IJK domaine_ijk_
Definition Grid_Level_Data_template.h:67
Grid_Level_Data_template::ijk_rho_
IJK_Field_template< _TYPE_, TRUSTArray< _TYPE_ > > ijk_rho_
Definition Grid_Level_Data_template.h:75
Grid_Level_Data_template::ijk_x_
IJK_Field_template< _TYPE_, TRUSTArray< _TYPE_ > > ijk_x_
Definition Grid_Level_Data_template.h:87
Grid_Level_Data_template::compute_faces_coefficients_from_inv_rho_cst_i_cst_j_cst_k
void compute_faces_coefficients_from_inv_rho_cst_i_cst_j_cst_k()
Definition Grid_Level_Data_template.tpp:300
Grid_Level_Data_template::get_faces_coefficients
const IJK_Field_template< _TYPE_, TRUSTArray< _TYPE_ > > & get_faces_coefficients() const
Definition Grid_Level_Data_template.h:58
Grid_Level_Data_template::ghost_size_
int ghost_size_
Definition Grid_Level_Data_template.h:68
Grid_Level_Data_template::compute_faces_coefficients_from_double_coeffs
std::enable_if_t< std::is_same< MY_TYPE, float >::value, void > compute_faces_coefficients_from_double_coeffs(const Grid_Level_Data_template< double > &)
Definition Grid_Level_Data_template.tpp:719
Grid_Level_Data_template::initialize
void initialize(const Domaine_IJK &, int ghost, int additional_k_layers)
Definition Grid_Level_Data_template.tpp:33
IJKArray_with_ghost
This is an array with [] operator.
Definition IJKArray_with_ghost.h:34
IJKArray_with_ghost::resize
void resize(int n, int new_ghost)
Definition IJKArray_with_ghost.h:51
IJK_Field_local_template::ni
int ni() const
Definition IJK_Field_local_template.h:167
IJK_Field_local_template::nj
int nj() const
Definition IJK_Field_local_template.h:168
IJK_Field_local_template::nk
int nk() const
Definition IJK_Field_local_template.h:169
IJK_Field_template
: This class is an IJK_Field_local with parallel informations.
Definition IJK_Field_template.h:38
IJK_Field_template::allocate
void allocate(const Domaine_IJK &d, Domaine_IJK::Localisation l, int ghost_size, int additional_k_layers=0, int nb_compo=1, const Nom &name=Nom(), bool external_storage=false, int monofluide=0, double rov=0., double rol=0., int use_inv_rho_in_pressure_solver=0)
Definition IJK_Field_template.tpp:562
IJK_Field_template::echange_espace_virtuel
void echange_espace_virtuel(int ghost)
Exchange data over "ghost" number of cells.
Definition IJK_Field_template.tpp:287
IJK_Shear_Periodic_helpler::rho_vap_ref_for_poisson_
static double rho_vap_ref_for_poisson_
Definition IJK_Shear_Periodic_helpler.h:30
IJK_Shear_Periodic_helpler::rho_liq_ref_for_poisson_
static double rho_liq_ref_for_poisson_
Definition IJK_Shear_Periodic_helpler.h:29
IJK_Shear_Periodic_helpler::defilement_
static int defilement_
Definition IJK_Shear_Periodic_helpler.h:27
IJ_layout
: This class describes an IJ plane of an IJK_Field_local.
Definition IJ_layout.h:27
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455