TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Matrice_Grossiere.tpp
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#ifndef Matrice_Grossiere_H_TPP
17#define Matrice_Grossiere_H_TPP
18
19#include <Matrice_Morse_Sym.h>
20
21template <typename _TYPE_, typename _TYPE_ARRAY_>
23{
27 const Domaine_IJK& splitting = coeffs_face.get_domaine();
28
29 const int ni = splitting.get_nb_elem_local(DIRECTION_I);
30 const int nj = splitting.get_nb_elem_local(DIRECTION_J);
31 const int nk = splitting.get_nb_elem_local(DIRECTION_K);
32
33
34 renum_.resize(nk+2, nj+2, ni+2, RESIZE_OPTIONS::NOCOPY_NOINIT);
35 renum_ = -1; // init a -1
36 // plusieur vecteur renum pour le cas shear periodic ou une case peut renvoyer vers plusieurs
37 // + 4 autres vecteur contenant les ponderation associee pour interpolation 4th order
42
43 int count = 0;
44 for (int k = 0; k < nk; k++)
45 for (int j = 0; j < nj; j++)
46 for (int i = 0; i < ni; i++)
47 {
48 int index = count++;
49 renum(i,j,k) = index;
50 }
51
52
53 // initialisation du tableau d'indice
54
55 ArrOfInt pe_voisins(6);
56
57 VECT(ArrOfInt) items_to_send(6);
58 VECT(ArrOfInt) items_to_recv(6);
59 VECT(ArrOfInt) blocs_to_recv(6);
60 int npe = 0;
61
62 const int pe_imin = splitting.get_neighbour_processor(0 /* left */, DIRECTION_I);
63 const int pe_jmin = splitting.get_neighbour_processor(0 /* left */, DIRECTION_J);
64 const int pe_kmin = splitting.get_neighbour_processor(0 /* left */, DIRECTION_K);
65 const int pe_imax = splitting.get_neighbour_processor(1 /* right */, DIRECTION_I);
66 const int pe_jmax = splitting.get_neighbour_processor(1 /* right */, DIRECTION_J);
67 const int pe_kmax = splitting.get_neighbour_processor(1 /* right */, DIRECTION_K);
68
69 int pe = pe_kmin;
70 if (pe >= 0)
71 {
72 pe_voisins[npe] = pe;
73 if (pe != pe_kmax)
74 {
75 add_virt_bloc(pe, count, 0,0,-1, ni,nj,0, blocs_to_recv[npe]);
76 add_dist_bloc(pe, 0,0,0, ni,nj,1, items_to_send[npe]);
77 }
78 else
79 {
80 // un processeur voisin a gauche et a droite (par periodicite)
81 // attention a l'ordre des blocs:
82 add_virt_bloc(pe, count, 0,0,-1, ni,nj,0, blocs_to_recv[npe]);
83 add_virt_bloc(pe, count, 0,0,nk, ni,nj,nk+1, blocs_to_recv[npe]);
84 add_dist_bloc(pe, 0,0,nk-1, ni,nj,nk, items_to_send[npe]);
85 add_dist_bloc(pe, 0,0,0, ni,nj,1, items_to_send[npe]);
86 }
87 npe++;
88 }
89
90 pe = pe_jmin;
91 if (pe >= 0)
92 {
93 pe_voisins[npe] = pe;
94 if (pe != pe_jmax)
95 {
96 add_virt_bloc(pe, count, 0,-1,0, ni,0,nk, blocs_to_recv[npe]);
97 add_dist_bloc(pe, 0,0,0, ni,1,nk, items_to_send[npe]);
98 }
99 else
100 {
101 // un processeur voisin a gauche et a droite (par periodicite)
102 // attention a l'ordre des blocs:
103 add_virt_bloc(pe, count, 0,-1,0, ni,0,nk, blocs_to_recv[npe]);
104 add_virt_bloc(pe, count, 0,nj,0, ni,nj+1,nk, blocs_to_recv[npe]);
105 add_dist_bloc(pe, 0,nj-1,0, ni,nj,nk, items_to_send[npe]);
106 add_dist_bloc(pe, 0,0,0, ni,1,nk, items_to_send[npe]);
107 }
108 npe++;
109 }
110
111 pe = pe_imin;
112 if (pe >= 0)
113 {
114 pe_voisins[npe] = pe;
115 if (pe != pe_imax)
116 {
117 add_virt_bloc(pe, count, -1,0,0, 0,nj,nk, blocs_to_recv[npe]);
118 add_dist_bloc(pe, 0,0,0, 1,nj,nk, items_to_send[npe]);
119 }
120 else
121 {
122 // un processeur voisin a gauche et a droite (par periodicite)
123 // attention a l'ordre des blocs:
124 add_virt_bloc(pe, count, -1,0,0, 0,nj,nk, blocs_to_recv[npe]);
125 add_virt_bloc(pe, count, ni,0,0, ni+1,nj,nk, blocs_to_recv[npe]);
126 add_dist_bloc(pe, ni-1,0,0, ni,nj,nk, items_to_send[npe]);
127 add_dist_bloc(pe, 0,0,0, 1,nj,nk, items_to_send[npe]);
128 }
129 npe++;
130 }
131
132 pe = pe_imax;
133 if (pe >= 0 && pe != pe_imin)
134 {
135 pe_voisins[npe] = pe;
136 add_virt_bloc(pe, count, ni,0,0, ni+1,nj,nk, blocs_to_recv[npe]);
137 add_dist_bloc(pe, ni-1,0,0, ni,nj,nk, items_to_send[npe]);
138 npe++;
139 }
140
141 pe = pe_jmax;
142 if (pe >= 0 && pe != pe_jmin)
143 {
144 pe_voisins[npe] = pe;
145 add_virt_bloc(pe, count, 0,nj,0, ni,nj+1,nk, blocs_to_recv[npe]);
146 add_dist_bloc(pe, 0,nj-1,0, ni,nj,nk, items_to_send[npe]);
147 npe++;
148 }
149
150 pe = pe_kmax;
151 if (pe >= 0 && pe != pe_kmin)
152 {
153 pe_voisins[npe] = pe;
154 add_virt_bloc(pe, count, 0,0,nk, ni,nj,nk+1, blocs_to_recv[npe]);
155 add_dist_bloc(pe, 0,0,nk-1, ni,nj,nk, items_to_send[npe]);
156 npe++;
157 }
158
159 MD_Vector_std md_std(count /* nb_items_tot */,
160 ni * nj * nk /* nb_items_reels */,
161 pe_voisins,
162 items_to_send,
163 items_to_recv,
164 blocs_to_recv);
165 md_.copy(md_std);
166
167
168
169
170 const int n_reels = md_->get_nb_items_reels();
171 voisins_.dimensionner(n_reels);
172 coeffs_.dimensionner(n_reels);
173 coeff_diag_.resize(n_reels);
174 voisins_virt_.dimensionner(n_reels);
175 coeffs_virt_.dimensionner(n_reels);
176 constexpr int z_index_min = 0;
177 const int z_index = splitting.get_local_slice_index(2);
178 const int z_index_max = splitting.get_nprocessor_per_direction(2) - 1;
179
180 for (int i = 0; i < ni; i++)
181 {
182
183 const double DX = splitting.get_constant_delta(0);
184 const double shear_per_DX = shear_x_time_ / DX;
185
186 const int offset_pos = static_cast<int>(std::round(shear_per_DX));
187 interpolation_for_shear_periodicity(i , offset_pos, shear_per_DX, ni, 1.);
188
189 const int offset_neg = static_cast<int>(std::round(-shear_per_DX));
190 interpolation_for_shear_periodicity(i ,offset_neg, -shear_per_DX, ni, -1.);
191
192
193 for (int k = 0; k < nk; k++)
194 {
195 for (int j = 0; j < nj; j++)
196 {
197 if (z_index==z_index_min && defilement_==1 && k==0)
198 {
199 ajoute_coeff(i,j,k,i,j,k-1,coeffs_face(i,j,k,2), 1.);
200 }
201 else
202 {
203 ajoute_coeff(i,j,k,i,j,k-1,coeffs_face(i,j,k,2));
204 }
205 ajoute_coeff(i,j,k,i,j-1,k,coeffs_face(i,j,k,1));
206 ajoute_coeff(i,j,k,i-1,j,k,coeffs_face(i,j,k,0));
207 ajoute_coeff(i,j,k,i+1,j,k,coeffs_face(i+1,j,k,0));
208 ajoute_coeff(i,j,k,i,j+1,k,coeffs_face(i,j+1,k,1));
209 if (z_index==z_index_max && defilement_==1 && k==nk-1)
210 {
211 ajoute_coeff(i,j,k,i,j,k+1,coeffs_face(i,j,k+1,2), -1.);
212 }
213 else
214 {
215 ajoute_coeff(i,j,k,i,j,k+1,coeffs_face(i,j,k+1,2));
216 }
217 }
218 }
219 }
220
221 for (int i = 0; i < n_reels; i++)
222 {
223 if (coeff_diag_[i] == 0.)
224 coeff_diag_[i] = 1.;
225 }
226 int nnz = 0;
227 for (int i = 0; i < n_reels; i++)
228 nnz += 1 + coeffs_[i].size();
229 int nnz_virt = 0;
230 for (int i = 0; i < n_reels; i++)
231 nnz_virt += coeffs_virt_[i].size();
232
233 if (Process::me() == 0)
234 coeff_diag_[0] *= 2;
235
236 mat_.dimensionner(1, 2);
237 mat_.get_bloc(0,0).typer("Matrice_Morse_Sym");
238 mat_.get_bloc(0,1).typer("Matrice_Morse");
239 Matrice_Morse_Sym& carre = ref_cast(Matrice_Morse_Sym ,mat_.get_bloc(0,0).valeur());
240 Matrice_Morse& rect = ref_cast(Matrice_Morse , mat_.get_bloc(0,1).valeur());
241
242 carre.dimensionner(n_reels, nnz);
243 carre.remplir(voisins_, coeffs_, coeff_diag_);
244 rect.dimensionner(n_reels, md_->get_nb_items_tot() - n_reels, nnz_virt);
246
247 voisins_ = IntLists();
248 coeffs_ = DoubleLists();
249 coeff_diag_.reset();
250 voisins_virt_ = IntLists();
251 coeffs_virt_ = DoubleLists();
252
253
254}
255
256
257template <typename _TYPE_, typename _TYPE_ARRAY_>
259{
263 const Domaine_IJK& splitting = coeffs_face.get_domaine();
264
265 const int ni = splitting.get_nb_elem_local(DIRECTION_I);
266 const int nj = splitting.get_nb_elem_local(DIRECTION_J);
267 const int nk = splitting.get_nb_elem_local(DIRECTION_K);
268
269
270 renum_.resize(nk+2, nj+2, ni+2, RESIZE_OPTIONS::NOCOPY_NOINIT);
271 renum_ = -1; // init a -1
272 // plusieur vecteur renum pour le cas shear periodic ou une case peut renvoyer vers plusieurs
273 // + 4 autres vecteur contenant les ponderation associee pour interpolation 4th order
278
279 int count = 0;
280 for (int k = 0; k < nk; k++)
281 for (int j = 0; j < nj; j++)
282 for (int i = 0; i < ni; i++)
283 {
284 int index = count++;
285 renum(i,j,k) = index;
286 }
287
288
289 // initialisation du tableau d'indice
290
291 ArrOfInt pe_voisins(6);
292
293 VECT(ArrOfInt) items_to_send(6);
294 VECT(ArrOfInt) items_to_recv(6);
295 VECT(ArrOfInt) blocs_to_recv(6);
296 int npe = 0;
297
298 const int pe_imin = splitting.get_neighbour_processor(0 /* left */, DIRECTION_I);
299 const int pe_jmin = splitting.get_neighbour_processor(0 /* left */, DIRECTION_J);
300 const int pe_kmin = splitting.get_neighbour_processor(0 /* left */, DIRECTION_K);
301 const int pe_imax = splitting.get_neighbour_processor(1 /* right */, DIRECTION_I);
302 const int pe_jmax = splitting.get_neighbour_processor(1 /* right */, DIRECTION_J);
303 const int pe_kmax = splitting.get_neighbour_processor(1 /* right */, DIRECTION_K);
304
305 int pe = pe_kmin;
306 if (pe >= 0)
307 {
308 pe_voisins[npe] = pe;
309 if (pe != pe_kmax)
310 {
311 add_virt_bloc(pe, count, 0,0,-1, ni,nj,0, blocs_to_recv[npe]);
312 add_dist_bloc(pe, 0,0,0, ni,nj,1, items_to_send[npe]);
313 }
314 else
315 {
316 // un processeur voisin a gauche et a droite (par periodicite)
317 // attention a l'ordre des blocs:
318 add_virt_bloc(pe, count, 0,0,-1, ni,nj,0, blocs_to_recv[npe]);
319 add_virt_bloc(pe, count, 0,0,nk, ni,nj,nk+1, blocs_to_recv[npe]);
320 add_dist_bloc(pe, 0,0,nk-1, ni,nj,nk, items_to_send[npe]);
321 add_dist_bloc(pe, 0,0,0, ni,nj,1, items_to_send[npe]);
322 }
323 npe++;
324 }
325
326 pe = pe_jmin;
327 if (pe >= 0)
328 {
329 pe_voisins[npe] = pe;
330 if (pe != pe_jmax)
331 {
332 add_virt_bloc(pe, count, 0,-1,0, ni,0,nk, blocs_to_recv[npe]);
333 add_dist_bloc(pe, 0,0,0, ni,1,nk, items_to_send[npe]);
334 }
335 else
336 {
337 // un processeur voisin a gauche et a droite (par periodicite)
338 // attention a l'ordre des blocs:
339 add_virt_bloc(pe, count, 0,-1,0, ni,0,nk, blocs_to_recv[npe]);
340 add_virt_bloc(pe, count, 0,nj,0, ni,nj+1,nk, blocs_to_recv[npe]);
341 add_dist_bloc(pe, 0,nj-1,0, ni,nj,nk, items_to_send[npe]);
342 add_dist_bloc(pe, 0,0,0, ni,1,nk, items_to_send[npe]);
343 }
344 npe++;
345 }
346
347 pe = pe_imin;
348 if (pe >= 0)
349 {
350 pe_voisins[npe] = pe;
351 if (pe != pe_imax)
352 {
353 add_virt_bloc(pe, count, -1,0,0, 0,nj,nk, blocs_to_recv[npe]);
354 add_dist_bloc(pe, 0,0,0, 1,nj,nk, items_to_send[npe]);
355 }
356 else
357 {
358 // un processeur voisin a gauche et a droite (par periodicite)
359 // attention a l'ordre des blocs:
360 add_virt_bloc(pe, count, -1,0,0, 0,nj,nk, blocs_to_recv[npe]);
361 add_virt_bloc(pe, count, ni,0,0, ni+1,nj,nk, blocs_to_recv[npe]);
362 add_dist_bloc(pe, ni-1,0,0, ni,nj,nk, items_to_send[npe]);
363 add_dist_bloc(pe, 0,0,0, 1,nj,nk, items_to_send[npe]);
364 }
365 npe++;
366 }
367
368 pe = pe_imax;
369 if (pe >= 0 && pe != pe_imin)
370 {
371 pe_voisins[npe] = pe;
372 add_virt_bloc(pe, count, ni,0,0, ni+1,nj,nk, blocs_to_recv[npe]);
373 add_dist_bloc(pe, ni-1,0,0, ni,nj,nk, items_to_send[npe]);
374 npe++;
375 }
376
377 pe = pe_jmax;
378 if (pe >= 0 && pe != pe_jmin)
379 {
380 pe_voisins[npe] = pe;
381 add_virt_bloc(pe, count, 0,nj,0, ni,nj+1,nk, blocs_to_recv[npe]);
382 add_dist_bloc(pe, 0,nj-1,0, ni,nj,nk, items_to_send[npe]);
383 npe++;
384 }
385
386 pe = pe_kmax;
387 if (pe >= 0 && pe != pe_kmin)
388 {
389 pe_voisins[npe] = pe;
390 add_virt_bloc(pe, count, 0,0,nk, ni,nj,nk+1, blocs_to_recv[npe]);
391 add_dist_bloc(pe, 0,0,nk-1, ni,nj,nk, items_to_send[npe]);
392 npe++;
393 }
394
395 MD_Vector_std md_std(count /* nb_items_tot */,
396 ni * nj * nk /* nb_items_reels */,
397 pe_voisins,
398 items_to_send,
399 items_to_recv,
400 blocs_to_recv);
401 md_.copy(md_std);
402
403
404
405
406 const int n_reels = md_.valeur().get_nb_items_reels();
407 voisins_2_.dimensionner(n_reels);
408 coeffs_2_.dimensionner(n_reels);
409 coeff_diag_.resize(n_reels);
410 voisins_virt_2_.dimensionner(n_reels);
411 coeffs_virt_2_.dimensionner(n_reels);
412 int z_index_min = 0;
413 int z_index = splitting.get_local_slice_index(2);
414 int z_index_max = splitting.get_nprocessor_per_direction(2) - 1;
415
416 for (int i = 0; i < ni; i++)
417 {
418
419 const double DX = splitting.get_constant_delta(0);
420 const double shear_per_DX = shear_x_time_ / DX;
421
422 const int offset_pos = static_cast<int>(std::round(shear_per_DX));
423 interpolation_for_shear_periodicity(i , offset_pos, shear_per_DX, ni, 1.);
424
425 const int offset_neg = static_cast<int>(std::round(-shear_per_DX));
426 interpolation_for_shear_periodicity(i ,offset_neg, -shear_per_DX, ni, -1.);
427
428
429 for (int k = 0; k < nk; k++)
430 {
431 for (int j = 0; j < nj; j++)
432 {
433 if (z_index==z_index_min && defilement_==1 && k==0)
434 {
435 ajoute_coeff2(i,j,k,i,j,k-1,coeffs_face(i,j,k,2), 1.);
436 }
437 else
438 {
439 ajoute_coeff2(i,j,k,i,j,k-1,coeffs_face(i,j,k,2));
440 }
441 ajoute_coeff2(i,j,k,i,j-1,k,coeffs_face(i,j,k,1));
442 ajoute_coeff2(i,j,k,i-1,j,k,coeffs_face(i,j,k,0));
443 ajoute_coeff2(i,j,k,i+1,j,k,coeffs_face(i+1,j,k,0));
444 ajoute_coeff2(i,j,k,i,j+1,k,coeffs_face(i,j+1,k,1));
445 if (z_index==z_index_max && defilement_==1 && k==nk-1)
446 {
447 ajoute_coeff2(i,j,k,i,j,k+1,coeffs_face(i,j,k+1,2), -1.);
448 }
449 else
450 {
451 ajoute_coeff2(i,j,k,i,j,k+1,coeffs_face(i,j,k+1,2));
452 }
453 }
454 }
455 }
456
457 for (int i = 0; i < n_reels; i++)
458 {
459 if (coeff_diag_[i] == 0.)
460 coeff_diag_[i] = 1.;
461 }
462 int nnz = 0;
463 for (int i = 0; i < n_reels; i++)
464 nnz += 1 + coeffs_2_[i].size();
465 int nnz_virt = 0;
466 for (int i = 0; i < n_reels; i++)
467 nnz_virt += coeffs_virt_2_[i].size();
468
469 if (Process::me() == 0)
470 coeff_diag_[0] *= 2;
471
472 mat_.dimensionner(1, 2);
473 mat_.get_bloc(0,0).typer("Matrice_Morse");
474 mat_.get_bloc(0,1).typer("Matrice_Morse");
475 Matrice_Morse& carre = ref_cast(Matrice_Morse ,mat_.get_bloc(0,0).valeur());
476 Matrice_Morse& rect = ref_cast(Matrice_Morse , mat_.get_bloc(0,1).valeur());
477
478 carre.dimensionner(n_reels, nnz);
479 carre.remplir(voisins_2_, coeffs_2_, coeff_diag_);
480 rect.dimensionner(n_reels, md_.valeur().get_nb_items_tot() - n_reels, nnz_virt);
482
483 voisins_2_ = IntLists();
484 coeffs_2_ = DoubleLists();
485 coeff_diag_.reset();
486 voisins_virt_2_ = IntLists();
487 coeffs_virt_2_ = DoubleLists();
488
489
490}
491
492#endif
This class encapsulates all the information related to the eulerian mesh for TrioIJK.
Definition Domaine_IJK.h:47
int get_neighbour_processor(int previous_or_next, int direction) const
Returns the index of the requested neighbour processor (-1 if no neighbour).
int get_nb_elem_local(int direction) const
Returns the number of elements owned by this processor in the given direction.
int get_nprocessor_per_direction(int direction) const
Returns the number of slices in the given direction.
double get_constant_delta(int direction) const
Returns the size of cells in a direction.
int get_local_slice_index(int direction) const
Returns the position of the local subdomain in the requested direction.
: This class is an IJK_Field_local with parallel informations.
const Domaine_IJK & get_domaine() const
C'est le plus simple des descripteurs, utilise pour les tableaux de valeurs aux sommets,...
ArrOfDouble ponderation_shear_p_
void build_matrix_test(const IJK_Field_template< _TYPE_, _TYPE_ARRAY_ > &coeffs_face)
void build_matrix(const IJK_Field_template< _TYPE_, _TYPE_ARRAY_ > &coeffs_face)
const int & renum(int i, int j, int k) const
void interpolation_for_shear_periodicity(const int i, const int send_i, const double istmp, const int real_size_i, const double shear_perio)
void ajoute_coeff(int i, int j, int k, int i_voisin, int j_voisin, int k_voisin, const double coeff, const double shear_perio)
ajoute deux coefficients diagonal/extra-diagonal a la matrice
ArrOfDouble ponderation_shear_m_
void ajoute_coeff2(int i, int j, int k, int i_voisin, int j_voisin, int k_voisin, const double coeff, const double shear_perio)
ajoute deux coefficients diagonal/extra-diagonal a la matrice
void add_dist_bloc(int pe, int imin, int jmin, int kmin, int imax, int jmax, int kmax, ArrOfInt &items_to_send)
void add_virt_bloc(int pe, int &count, int imin, int jmin, int kmin, int imax, int jmax, int kmax, ArrOfInt &virt_blocs)
Classe Matrice_Morse_Sym Represente une matrice M (creuse) symetrique stockee au format Morse.
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
void dimensionner(int n, _SIZE_ nnz)
Size the matrix with n lines and n columns and nnz zero-values coefficients.
void remplir(const IntLists &, const DoubleLists &, const DoubleVect &)
static int me()
renvoie mon rang dans le groupe de communication courant.
Definition Process.cpp:125