TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Force_ph.cpp
1/****************************************************************************
2* Copyright (c) 2021, 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 <Force_ph.h>
17#include <IJK_Field_vector.h>
18#include <fstream>
19#include <math.h>
20#include <IJK_Navier_Stokes_tools.h>
21#include <communications.h>
22#include <Domaine_IJK.h>
23
24// #include <fftw3.h>
25
26Implemente_instanciable( Force_ph, "Force_ph", Objet_U ) ;
27
28Sortie& Force_ph::printOn( Sortie& os ) const
29{
30 Objet_U::printOn( os );
31 return os;
32}
33
34Entree& Force_ph::readOn( Entree& is )
35{
36 Objet_U::readOn( is );
37 return is;
38}
39
40void Force_ph::initialise(int a_nproc_tot, int a_ni,int a_nj,int a_nk,int a_nl,int a_nm, int a_nn,
41 double a_Lx, double a_Ly, double a_Lz, double a_Ox,double a_Oy,double a_Oz, int a_momin, int a_momax, double a_kmin, double a_kmax,
42 std::string nom_fichier, const Domaine_IJK& splitting,
43 int a_i_offset, int a_j_offset, int a_k_offset
44 )
45
46{
47// int ind,i,j,k;
48 nproc_tot = a_nproc_tot;
49
50 ni = a_ni;
51 nj = a_nj;
52 nk = a_nk;
53 n_ijk = ni*nj*nk;
54 nl = a_nl;
55 nm = a_nm;
56 nn = a_nn;
57 n_lmn = (2*nl+1)*(2*nm+1)*(2*nn+1);
58 Lx = a_Lx;
59 Ly = a_Ly;
60 Lz = a_Lz;
61 Ox = a_Ox;
62 Oy = a_Oy;
63 Oz = a_Oz;
64 momin = a_momin;
65 momax = a_momax;
66 kmin = a_kmin;
67 kmax = a_kmax;
68 i_offset = a_i_offset;
69 j_offset = a_j_offset;
70 k_offset = a_k_offset;
71// Cout << "a_i_offset : " << a_i_offset << finl;
72
73// NB : Normalement force est juste reel, mais on verifie qd mm ...
74 force = set_dimensions(force,2,3,n_ijk);
75
76 allocate_velocity(force_, splitting, 2); //: 2 a la fin -> 2 cellules ghost a la fin
77 force_.nommer("FORCE_PH");
78///////////////////////////////////////////////////////////////
79 std::ofstream Physical_flux(nom_fichier.c_str());
80 if(Physical_flux)
81 {
82 Physical_flux << "-------- PHYSICAL_FORCE --------" << std::endl;
83 Physical_flux << std::endl << "i,j,k \t : f_x, \tf_y, \tf_z,\t";
84 }
85///////////////////////////////////////////////////////////////
86// TODO fix that path wtf
87 std::ofstream Offset_flux("/volatile/FFTW/creeping_flow/creeping_flow_bis/src/offset.txt");
88 if(Offset_flux)
89 {
90 Offset_flux << "-------- PHYSICAL_POSITION --------" << std::endl;
91 Offset_flux << std::endl << "i,j,k \t : f_x, \tf_y, \tf_z,\t";
92 }
93
94 energie = 0.;
95
96}
97
98void Force_ph::set_zero()
99{
100 /*Mise a zero de tout le champ de force physique*/
101 force_[0].data() = 0.;
102 force_[1].data() = 0.;
103 force_[2].data() = 0.;
104}
105
106void Force_ph::cheat_function()
107{
108 /*
109 Fonction de test. Ecrit directement la fonction cos(x-z) ex dans le champ
110 de force physique sans passer par le champ spectral.
111 QUESTION : Si je mets += au lieu de = a partir du premier pdt
112 j'obtiens une amplitude -6;+6 avec le time_scheme : RK3
113 */
114 double pos[3];
115 int roc(momin+0);
116 double t(2.*M_PI*roc/Lx);
117
118 for (int k=0; k<nk; k++)
119 {
120 pos[2] = Oz + k*Lz/nk;
121 for (int j=0; j<nj; j++)
122 {
123 pos[1] = Oy + j*Ly/nj;
124 for (int i=0; i<ni; i++)
125 {
126 int ind_ijk = (k*nj + j)*ni + i;
127 pos[0] = Ox + i*Lx/ni;
128
129 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
130 force[0][0][ind_ijk] = 2*cos(t*(pos[0]-pos[2]));
131 force[0][2][ind_ijk] = 2*cos(t*(pos[0]-pos[2]));
132 force[1][0][ind_ijk] = 0.;
133 force[1][2][ind_ijk] = 0.;
134 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
136 force_[0](i,j,k) = 2*cos(t*(pos[0]-pos[2]));
137 force_[2](i,j,k) = 2*cos(t*(pos[0]-pos[2]));
138 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
139 }
140 }
141 }
142}
143
144
145void Force_ph::from_spect_to_phys2(const std:: vector <double >& coeff_force)
146{
147 //////////////////////////////////////////////////////////////////////////////
148 /*
149 Version primitive de la transformee de Fourier inverse :
150 - 6 boucles imbriquees (3 d'espace physique, 3 d'espace spectral).
151 - coefficients de la force spectrale stockes dans un tableau
152 a une dimension, on reconstruit l'indice a aller chercher.
153 - coefficients de la force physique recopies dans un tableau
154 a trois dimensions.
155 */
156 //////////////////////////////////////////////////////////////////////////////
157
158 // Initialisation du champ de force
159 for (int dir = 0; dir < 3; dir++)
160 {
161 force_[dir].data() = 0.;
162 }
163
164 const int n_ind_mni((2*nm+1)*(2*nn+1)*1*ni);
165 const int n_ind_nij((2*nn+1)*1*ni*1*nj);
166 const int n_ind_lmn((2*nl+1)*(2*nm+1)*(2*nn+1));
167
168 std:: vector <double > bout_a_quart(2*3*n_ind_mni,0.);
169 std:: vector <double > bout_b_quart(2*3*n_ind_nij,0.);
170
171 double expo_a[2], expo_b[2], expo_c[2];
172 double pos[3];
173 double dz(Lz/nk), dy(Ly/nj), dx(Lx/ni);
174
175 int n_dir = 3;
176
177 //////////////////////////////////////////////////////////////
178 // TODO : remplacer mes pos[X] par un resultat de build_local_coords
179 // --> Dans IJK_Switch ou IJK_switch_FT
180 for (int k=0; k<nk; k++)
181 {
182 pos[2] = Oz + k*Lz/nk + dz/2.;
183 for (int j=0; j<nj; j++)
184 {
185 pos[1] = Oy + j*Ly/nj + dy/2.;
186 for (int i=0; i<ni; i++)
187 {
188 pos[0] = Ox + i*Lx/ni + dx/2.;
189 for (int n=0; n<2*nn; n++)
190 {
191 int ind_nij = ( (j*ni + i) * (2*nn+1) + n );
192
193 // Initialisation de la valeur de bout_b
194
195 for(int dir=0; dir<3; dir++)
196 {
197 int ind_RDInij((0*n_dir+dir) * n_ind_nij + ind_nij);
198 int ind_CDInij((1*n_dir+dir) * n_ind_nij + ind_nij);
199 bout_b_quart[ind_RDInij] = 0.0;
200 bout_b_quart[ind_CDInij] = 0.0;
201 }
202
203 double kappa_n = (- kmax + (n)*(2*kmax)/(2*nn) );
204 for (int m=0; m<2*nm; m++)
205 {
206 int ind_mni = ((i*(2*nn+1) + n) * (2*nm+1) + m );
207
208 // Initialisation de la valeur de bout_a
209
210 for(int dir=0; dir<3; dir++)
211 {
212 int ind_RDImni((0*n_dir+dir) * n_ind_mni + ind_mni);
213 int ind_CDImni((1*n_dir+dir) * n_ind_mni + ind_mni);
214 bout_a_quart[ind_RDImni] = 0.0;
215 bout_a_quart[ind_CDImni] = 0.0;
216 }
217
218 double kappa_m = (- kmax + (m)*(2*kmax)/(2*nm));
219
220 for (int l=0; l<2*nl; l++)
221 {
222 int ind_lmn = ( (n*(2*nm+1) + m) * (2*nl+1) +l);
223 double kappa_l = (- kmax + (l)*(2*kmax)/(2*nl));
224
225 expo_a[0] = cos(kappa_l*pos[0]);
226 expo_a[1] = sin(kappa_l*pos[0]);
227
228 for (int dir=0; dir<3; dir++)
229 {
230 int ind_RDImni((0*n_dir+dir)*n_ind_mni+ind_mni);
231 int ind_RDIlmn((0*n_dir+dir)*n_ind_lmn+ind_lmn);
232 int ind_CDImni((1*n_dir+dir)*n_ind_mni+ind_mni);
233 int ind_CDIlmn((1*n_dir+dir)*n_ind_lmn+ind_lmn);
234
235 bout_a_quart[ind_RDImni] += 1*(coeff_force[ind_RDIlmn]*expo_a[0] - coeff_force[ind_CDIlmn]*expo_a[1]);
236 bout_a_quart[ind_CDImni] += 1*(coeff_force[ind_RDIlmn]*expo_a[1] + coeff_force[ind_CDIlmn]*expo_a[0]);
237 }
238 }
239 expo_b[0] = cos(kappa_m*pos[1]);
240 expo_b[1] = sin(kappa_m*pos[1]);
241 for (int dir=0; dir<3; dir++)
242 {
243 int ind_RDInij((0*n_dir+dir)*n_ind_nij+ind_nij);
244 int ind_RDImni((0*n_dir+dir)*n_ind_mni+ind_mni);
245 int ind_CDInij((1*n_dir+dir)*n_ind_nij+ind_nij);
246 int ind_CDImni((1*n_dir+dir)*n_ind_mni+ind_mni);
247
248 bout_b_quart[ind_RDInij] += 1*(bout_a_quart[ind_RDImni]*expo_b[0] - bout_a_quart[ind_CDImni]*expo_b[1]);
249 bout_b_quart[ind_CDInij] += 1*(bout_a_quart[ind_RDImni]*expo_b[1] + bout_a_quart[ind_CDImni]*expo_b[0]);
250 }
251 }
252 expo_c[0] = cos(kappa_n*pos[2]);
253 expo_c[1] = sin(kappa_n*pos[2]);
254 for (int dir=0; dir<3; dir++)
255 {
256 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
257 // ~ RDI : Real,Direction,Index
258 // ~ CDI : Complex,Direction,Index
259 int ind_RDInij((0*n_dir+dir)*n_ind_nij+ind_nij);
260 int ind_CDInij((1*n_dir+dir)*n_ind_nij+ind_nij);
261
262 // force[0][dir][ind_ijk] += 1*(bout_b_quart[ind_RDInij]*expo_c[0] - bout_b_quart[ind_CDInij]*expo_c[1]);
263 // force[1][dir][ind_ijk] += 1*(bout_b_quart[ind_RDInij]*expo_c[1] + bout_b_quart[ind_CDInij]*expo_c[0]);
264
265 force_[dir](i,j,k) += 1*(bout_b_quart[ind_RDInij]*expo_c[0] - bout_b_quart[ind_CDInij]*expo_c[1]);
266 }
267 }
268 }
269 }
270 }
271}
272
273
274
275void Force_ph::from_spect_to_phys_opti2(ArrOfDouble& coeff_force )
276{
277 //////////////////////////////////////////////////////////////////////////////
278 /*
279 Version la plus avancee de la transformee de Fourier inverse :
280 - Les 6 boucles imbriquees sont eclatees en 3 blocs de complexite
281 au plus ni*nj*nk*nn.
282 - coefficients de la force spectrale stockes dans un tableau
283 a une dimension, on reconstruit l'indice a aller chercher.
284 - coefficients de la force physique recopies dans un tableau
285 a trois dimensions (a changer, facile a changer), mais aussi stockes dans un
286 IJK_Field_vector3_double
287 */
288 //////////////////////////////////////////////////////////////////////////////
289
290 const int ni_o(ni);
291 const int nj_o(nj);
292
293 // Initialisation du champ de force
294 for (int dir = 0; dir < 3; dir++)
295 {
296 force_[dir].data() = 0.;
297 }
298 int n_dir(3);
299 int n_ind_mni((2*nm+1)*(2*nn+1)*1*ni_o);
300 int n_ind_nij((2*nn+1)*1*ni_o*1*nj_o);
301 int n_ind_lmn((2*nl+1)*(2*nm+1)*(2*nn+1));
302
303 std:: vector <double > bout_a_quart(2*3*n_ind_mni,0.);
304 std:: vector <double > bout_b_quart(2*3*n_ind_nij,0.);
305
306 double expo_a[2], expo_b[2], expo_c[2];
307 double pos[3];
308
309 //////////////////////////////////////////////////////////////////////////////
310 // Construction des posisitons. Comme la force est localisee (pour le moment)
311 // aux faces, on est obliges de definir les coord pour les faces x, celles pour les faces y
312 // et puis pour les faces z.
313 // 1) Reconstruire a la main ferai economiser 6 ArrOfDouble
314 // 2) localisr aux centres (sans reconstruire a la main) ferai econmiser 6 ArrOfDouble
315 ArrOfDouble pos_x_fx, pos_y_fx, pos_z_fx;
316 build_local_coords(force_[0], pos_x_fx, pos_y_fx, pos_z_fx);
317 ArrOfDouble pos_x_fy, pos_y_fy, pos_z_fy;
318 build_local_coords(force_[1], pos_x_fy, pos_y_fy, pos_z_fy);
319 ArrOfDouble pos_x_fz, pos_y_fz, pos_z_fz;
320 build_local_coords(force_[2], pos_x_fz, pos_y_fz, pos_z_fz);
321 //////////////////////////////////////////////////////////////////////////////
322
323 for (int dir=0; dir<3; dir++)
324 {
325 ///////////////////////////////////////////////////////////////
326 for (int m=0; m<2*nm+1; m++)
327 for (int n=0; n<2*nn+1; n++)
328 for (int i=0; i<ni; i++)
329 {
330 int ind_mni = ((i*(2*nn+1) + n) * (2*nm+1) + m );
331 // pos[0] = Ox + i*Lx/(ni) - dx/2.; // Localise aux faces
332 // pos[0] = Ox + i*Lx/(ni+i_offset) - dx/2.; // Localise aux faces
333 pos[0] = pos_x_fx[i]; // Ox + i*Lx/(ni+i_offset) - dx/2.; // Localise aux faces
334 for (int l=0; l<2*nl+1; l++)
335 {
336 int ind_lmn = ( (n*(2*nm+1) + m) * (2*nl+1) +l);
337 double kappa_l = (- kmax + (l)*(2*kmax)/(2*nl));
338 expo_a[0] = cos(kappa_l*pos[0]);
339 expo_a[1] = sin(kappa_l*pos[0]);
340
341 int ind_RDImni((0*n_dir+dir) * n_ind_mni + ind_mni);
342 int ind_RDIlmn((0*n_dir+dir) * n_ind_lmn + ind_lmn);
343 int ind_CDImni((1*n_dir+dir) * n_ind_mni + ind_mni);
344 int ind_CDIlmn((1*n_dir+dir) * n_ind_lmn + ind_lmn);
345 // Complex multiplication
346 bout_a_quart[ind_RDImni] += 1*(coeff_force[ind_RDIlmn]*expo_a[0] - coeff_force[ind_CDIlmn]*expo_a[1]);
347 bout_a_quart[ind_CDImni] += 1*(coeff_force[ind_RDIlmn]*expo_a[1] + coeff_force[ind_CDIlmn]*expo_a[0]);
348 }
349 }
350 ///////////////////////////////////////////////////////////////
351
352 ///////////////////////////////////////////////////////////////
353 for (int n=0; n<2*nn+1; n++)
354 for (int i=0; i<ni; i++)
355 for (int j=0; j<nj; j++)
356 {
357 int ind_nij = ( (j*ni_o + i) * (2*nn+1) + n );
358 // pos[1] = Oy + j*Ly/(nj) - dy/2.; // Localise aux faces
359 // pos[1] = Oy + j*Ly/(nj+j_offset) - dy/2.; // Localise aux faces
360 pos[1] = pos_y_fy[j]; // Localise aux faces
361 for (int m=0; m<2*nm; m++)
362 {
363 int ind_mni = ((i*(2*nn+1) + n) * (2*nm+1) + m );
364 double kappa_m = (- kmax + (m)*(2*kmax)/(2*nm));
365 expo_b[0] = cos(kappa_m*pos[1]);
366 expo_b[1] = sin(kappa_m*pos[1]);
367
368 int ind_RDInij((0*n_dir+dir)*n_ind_nij+ind_nij);
369 int ind_RDImni((0*n_dir+dir)*n_ind_mni+ind_mni);
370 int ind_CDInij((1*n_dir+dir)*n_ind_nij+ind_nij);
371 int ind_CDImni((1*n_dir+dir)*n_ind_mni+ind_mni);
372 // Complex multiplication
373 bout_b_quart[ind_RDInij] += 1*(bout_a_quart[ind_RDImni]*expo_b[0] - bout_a_quart[ind_CDImni]*expo_b[1]);
374 bout_b_quart[ind_CDInij] += 1*(bout_a_quart[ind_RDImni]*expo_b[1] + bout_a_quart[ind_CDImni]*expo_b[0]);
375 }
376 }
377 ///////////////////////////////////////////////////////////////
378
379 ///////////////////////////////////////////////////////////////
380 // GAB, question : pour cette boucle il vaut mieux mettre la boucle en k le plus a lexterieur, on est d'acc ?
381 for (int k=0; k<nk; k++)
382 for (int j=0; j<nj; j++)
383 for (int i=0; i<ni; i++)
384 {
385 // pos[2] = Oz + k*Lz/(nk) - dz/2.; // Localise aux faces
386 // pos[2] = Oz + k*Lz/(nk+k_offset) - dz/2.; // Localise aux faces
387 pos[2] = pos_z_fz[k]; // Localise aux faces
388 for (int n=0; n<2*nn+1; n++)
389 {
390 int ind_nij = ( (j*ni_o + i) * (2*nn+1) + n );
391 double kappa_n = - kmax + (n)*(2*kmax)/(2*nn);
392 expo_c[0] = cos(kappa_n*pos[2]);
393 expo_c[1] = sin(kappa_n*pos[2]);
394
395 int ind_RDInij((0*n_dir+dir)*n_ind_nij+ind_nij);
396 int ind_CDInij((1*n_dir+dir)*n_ind_nij+ind_nij);
397 // Smetrie Hermitienne :
398 // force[0][dir][ind_ijk] += 1*(bout_b_quart[ind_RDInij]*expo_c[0] - bout_b_quart[ind_CDInij]*expo_c[1]);
399 // force[1][dir][ind_ijk] += 1*(bout_b_quart[ind_RDInij]*expo_c[1] + bout_b_quart[ind_CDInij]*expo_c[0]);
400 // Real part only
401 force_[dir](i,j,k) += 1*(bout_b_quart[ind_RDInij]*expo_c[0] - bout_b_quart[ind_CDInij]*expo_c[1]);
402 }
403 }
404 ///////////////////////////////////////////////////////////////
405 // Si on plante ici, mettre un 0 au moment du allocate
406 force_[dir].echange_espace_virtuel(force_[dir].ghost());
407 }
408}
409
410void Force_ph::from_spect_to_phys_opti2_advection(ArrOfDouble& coeff_force, const ArrOfDouble& advection_length )
411{
412 //////////////////////////////////////////////////////////////////////////////
413 /*
414 Version la plus avancee de la transformee de Fourier inverse, avec advection du champ de force :
415 - Les 6 boucles imbriquees sont eclatees en 3 blocs de complexite
416 au plus ni*nj*nk*nn.
417 - coefficients de la force spectrale stockes dans un tableau
418 a une dimension, on reconstruit l'indice a aller chercher.
419 - coefficients de la force physique recopies dans un tableau
420 a trois dimensions (a changer, facile a changer), mais aussi stockes dans un
421 IJK_Field_vector3_double
422 */
423 //////////////////////////////////////////////////////////////////////////////
424
425 const int ni_o(ni);
426 const int nj_o(nj);
427
428 // Initialisation du champ de force
429 for (int dir = 0; dir < 3; dir++)
430 {
431 force_[dir].data() = 0.;
432 }
433 int n_dir(3);
434 int n_ind_mni((2*nm+1)*(2*nn+1)*1*ni_o);
435 int n_ind_nij((2*nn+1)*1*ni_o*1*nj_o);
436 int n_ind_lmn((2*nl+1)*(2*nm+1)*(2*nn+1));
437
438 std:: vector <double > bout_a_quart(2*3*n_ind_mni,0.);
439 std:: vector <double > bout_b_quart(2*3*n_ind_nij,0.);
440
441 double expo_a[2], expo_b[2], expo_c[2];
442 double pos[3];
443
444 //////////////////////////////////////////////////////////////////////////////
445 // Construction des posisitons. Comme la force est localisee (pour le moment)
446 // aux faces, on est obliges de definir les coord pour les faces x, celles pour les faces y
447 // et puis pour les faces z.
448 // 1) Reconstruire a la main ferai economiser 6 ArrOfDouble
449 // 2) localisr aux centres (sans reconstruire a la main) ferai econmiser 6 ArrOfDouble
450 ArrOfDouble pos_x_fx, pos_y_fx, pos_z_fx;
451 build_local_coords(force_[0], pos_x_fx, pos_y_fx, pos_z_fx);
452 ArrOfDouble pos_x_fy, pos_y_fy, pos_z_fy;
453 build_local_coords(force_[1], pos_x_fy, pos_y_fy, pos_z_fy);
454 ArrOfDouble pos_x_fz, pos_y_fz, pos_z_fz;
455 build_local_coords(force_[2], pos_x_fz, pos_y_fz, pos_z_fz);
456 //////////////////////////////////////////////////////////////////////////////
457
458 for (int dir=0; dir<3; dir++)
459 {
460 ///////////////////////////////////////////////////////////////
461 for (int m=0; m<2*nm+1; m++)
462 for (int n=0; n<2*nn+1; n++)
463 for (int i=0; i<ni; i++)
464 {
465 int ind_mni = ((i*(2*nn+1) + n) * (2*nm+1) + m );
466 // pos[0] = Ox + i*Lx/(ni) - dx/2.; // Localise aux faces
467 // pos[0] = Ox + i*Lx/(ni+i_offset) - dx/2.; // Localise aux faces
468 pos[0] = pos_x_fx[i]; // Ox + i*Lx/(ni+i_offset) - dx/2.; // Localise aux faces
469 for (int l=0; l<2*nl+1; l++)
470 {
471 int ind_lmn = ( (n*(2*nm+1) + m) * (2*nl+1) +l);
472 double kappa_l = (- kmax + (l)*(2*kmax)/(2*nl));
473 expo_a[0] = cos(kappa_l*(pos[0]-advection_length[0]));
474 expo_a[1] = sin(kappa_l*(pos[0]-advection_length[0]));
475
476 int ind_RDImni((0*n_dir+dir) * n_ind_mni + ind_mni);
477 int ind_RDIlmn((0*n_dir+dir) * n_ind_lmn + ind_lmn);
478 int ind_CDImni((1*n_dir+dir) * n_ind_mni + ind_mni);
479 int ind_CDIlmn((1*n_dir+dir) * n_ind_lmn + ind_lmn);
480
481 bout_a_quart[ind_RDImni] += 1*(coeff_force[ind_RDIlmn]*expo_a[0] - coeff_force[ind_CDIlmn]*expo_a[1]);
482 bout_a_quart[ind_CDImni] += 1*(coeff_force[ind_RDIlmn]*expo_a[1] + coeff_force[ind_CDIlmn]*expo_a[0]);
483 }
484 }
485 ///////////////////////////////////////////////////////////////
486
487 ///////////////////////////////////////////////////////////////
488 for (int n=0; n<2*nn+1; n++)
489 for (int i=0; i<ni; i++)
490 for (int j=0; j<nj; j++)
491 {
492 int ind_nij = ( (j*ni_o + i) * (2*nn+1) + n );
493 // pos[1] = Oy + j*Ly/(nj) - dy/2.; // Localise aux faces
494 // pos[1] = Oy + j*Ly/(nj+j_offset) - dy/2.; // Localise aux faces
495 pos[1] = pos_y_fy[j]; // Localise aux faces
496 for (int m=0; m<2*nm; m++)
497 {
498 int ind_mni = ((i*(2*nn+1) + n) * (2*nm+1) + m );
499 double kappa_m = (- kmax + (m)*(2*kmax)/(2*nm));
500 expo_b[0] = cos(kappa_m*(pos[1]-advection_length[1]));
501 expo_b[1] = sin(kappa_m*(pos[1]-advection_length[1]));
502
503 int ind_RDInij((0*n_dir+dir)*n_ind_nij+ind_nij);
504 int ind_RDImni((0*n_dir+dir)*n_ind_mni+ind_mni);
505 int ind_CDInij((1*n_dir+dir)*n_ind_nij+ind_nij);
506 int ind_CDImni((1*n_dir+dir)*n_ind_mni+ind_mni);
507
508 bout_b_quart[ind_RDInij] += 1*(bout_a_quart[ind_RDImni]*expo_b[0] - bout_a_quart[ind_CDImni]*expo_b[1]);
509 bout_b_quart[ind_CDInij] += 1*(bout_a_quart[ind_RDImni]*expo_b[1] + bout_a_quart[ind_CDImni]*expo_b[0]);
510 }
511 }
512 ///////////////////////////////////////////////////////////////
513
514 ///////////////////////////////////////////////////////////////
515 // GAB, question : la boucle en k est exterieure. C'est le plus optimal ? ?
516 for (int k=0; k<nk; k++)
517 for (int j=0; j<nj; j++)
518 for (int i=0; i<ni; i++)
519 {
520 // pos[2] = Oz + k*Lz/(nk) - dz/2.; // Localise aux faces
521 // pos[2] = Oz + k*Lz/(nk+k_offset) - dz/2.; // Localise aux faces
522 pos[2] = pos_z_fz[k]; // Localise aux faces
523 for (int n=0; n<2*nn+1; n++)
524 {
525 int ind_nij = ( (j*ni_o + i) * (2*nn+1) + n );
526 double kappa_n = - kmax + (n)*(2*kmax)/(2*nn);
527 expo_c[0] = cos(kappa_n*(pos[2]-advection_length[2]));
528 expo_c[1] = sin(kappa_n*(pos[2]-advection_length[2]));
529
530 int ind_RDInij((0*n_dir+dir)*n_ind_nij+ind_nij);
531 int ind_CDInij((1*n_dir+dir)*n_ind_nij+ind_nij);
532 // Symetrie Hermitienne :
533 // force[0][dir][ind_ijk] += 1*(bout_b_quart[ind_RDInij]*expo_c[0] - bout_b_quart[ind_CDInij]*expo_c[1]);
534 // force[1][dir][ind_ijk] += 1*(bout_b_quart[ind_RDInij]*expo_c[1] + bout_b_quart[ind_CDInij]*expo_c[0]);
535
536 force_[dir](i,j,k) += 1*(bout_b_quart[ind_RDInij]*expo_c[0] - bout_b_quart[ind_CDInij]*expo_c[1]);
537 }
538 }
539 ///////////////////////////////////////////////////////////////
540 // Si on plante ici, mettre un 0 au moment du allocate
541 force_[dir].echange_espace_virtuel(force_[dir].ghost());
542 }
543
544}
545
546
547
548
549void Force_ph::write(std::string nom_fichier_sortie, double t)
550{
551 //////////////////////////////////////////////////////////////////////////////
552 // ATTENTION : Cette fonction ecrit force et pas force_
553 // On a supprime force maintenant. Il faut re ecrire cette fct
554 // si on souhaite l'utiliser
555 //////////////////////////////////////////////////////////////////////////////
556 std::cout << "test in force_ph/write" << std::endl;
557 std::ofstream Physical_flux(nom_fichier_sortie.c_str(), std::ios::app);
558 if (Physical_flux)
559 {
560 int i, j, k, ind;
561
562 Physical_flux << std::endl << "time : " << t << std::endl << std::endl;
563 for (k=0; k<nk; k++)
564 for (j=0; j<nj; j++)
565 for (i=0; i<ni; i++)
566 {
567 ind = (k*nj + j)*ni + i;
568 Physical_flux << i<<","<<j<<","<<k<<"\t : ";
569 Physical_flux << force[0][0][ind] << " + i" << force[1][0][ind]<<", \t";
570 Physical_flux << force[0][1][ind] << " + i" << force[1][1][ind]<<", \t";
571 Physical_flux << force[0][2][ind] << " + i" << force[1][2][ind]<<", \t";
572 Physical_flux << std::endl;
573 }
574 }
575}
576
577
578void Force_ph::write_offset_index_position( const Domaine_IJK& my_splitting)
579{
580 std::ofstream Offset_flux("/volatile/FFTW/creeping_flow/creeping_flow_bis/src/offset.txt", std::ios::app);
581 if (Offset_flux)
582 {
583
584 Offset_flux << std::endl << "of_i,of_j,of_k,i,j,k,x,y,z"<< std::endl;
585 double dz(Lz/nk), dy(Ly/nj), dx(Lx/ni);
586 for (int k=0; k<nk; k++)
587 {
588 double z = Oz + k*Lz/nk + dz/2.;
589 for (int j=0; j<nj; j++)
590 {
591 double y = Oy + j*Ly/nj + dy/2.;
592 for (int i=0; i<ni; i++)
593 {
594 double x = Ox + i*Lx/ni + dx/2.;
595 Offset_flux << my_splitting.get_offset_local(DIRECTION_I)<<",";
596 Offset_flux << my_splitting.get_offset_local(DIRECTION_J)<<",";
597 Offset_flux << my_splitting.get_offset_local(DIRECTION_K)<<";";
598 Offset_flux << i<<","<<j<<","<<k<<";";
599 Offset_flux << x<<","<<y<<","<<z<<";";
600 Offset_flux << std::endl;
601 }
602 }
603 }
604 }
605}
606
607
608
609
610void Force_ph::compute_energie()
611{
612 int i,j,k,dir,ind;
613 for (k=0; k<nk; k++)
614 for (j=0; j<nj; j++)
615 for (i=0; i<ni; i++)
616 {
617 ind = (k*nj + j)*ni + i;
618 for (dir=0; dir<3; dir++)
619 energie += force[0][dir][ind]*force[0][dir][ind] + force[1][dir][ind]*force[1][dir][ind];
620 }
621 energie *= ((1.)/(ni*nj*nk)) ;
622}
623
624double Force_ph::get_energie()
625{
626 return energie;
627}
628
629void Force_ph::write_separate(std::string nom_fichier_sortie, double t)
630{
631 //////////////////////////////////////////////////////////////////////////////
632 // ATTENTION : Cette fonction ecrit force et pas force_
633 // On a supprime force maintenant. Il faut re ecrire cette fct
634 // si on souhaite l'utiliser
635 //////////////////////////////////////////////////////////////////////////////
636 std::ofstream Physical_flux(nom_fichier_sortie.c_str());
637 if (Physical_flux)
638 {
639 int i, j, k, ind;
640 double x,y,z;
641 Physical_flux << std::endl << "i,j,k, \t x,y,z, \t rf_x, \t cf_x, \t\t rf_y, \t cf_y, \t\t rf_z, \t cf_z \t";
642 Physical_flux << std::endl;
643
644 for (k=0; k<nk; k++)
645 for (j=0; j<nj; j++)
646 for (i=0; i<ni; i++)
647 {
648 // ind = (k*nk + j)*nj +i;
649 ind = (k*nj + j)*ni + i;
650
651 x = -Lx/2 + i*Lx/ni;
652 y = -Ly/2 + j*Ly/nj;
653 z = -Lz/2 + k*Lz/nk;
654
655 Physical_flux << i<<","<<j<<","<<k<<"\t,";
656 Physical_flux << x << ",\t" << y << ",\t" << z <<", \t\t";
657 Physical_flux << force[0][0][ind] << ",\t" << force[1][0][ind]<<", \t\t";
658 Physical_flux << force[0][1][ind] << ",\t" << force[1][1][ind]<<", \t\t";
659 Physical_flux << force[0][2][ind] << ",\t" << force[1][2][ind]<<"\t";
660 Physical_flux << std::endl;
661 }
662 }
663}
664
665IJK_Field_vector3_double Force_ph::get_force_attribute()
666{
667 return force_;
668}
669
670IJK_Field_vector3_double& Force_ph::get_force_attribute2()
671{
672 return force_;
673}
674
675
676std::vector< std::vector< std:: vector <double >>> set_dimensions(std::vector< std::vector< std:: vector <double >>> the_vector,
677 int dim_one, int dim_two, int dim_three)
678{
679 int i,j;
680 the_vector.resize(dim_one);
681 for (i=0; i<dim_one; i++)
682 {
683 the_vector[i].resize(dim_two);
684 for (j=0; j<dim_two; j++)
685 {
686 the_vector[i][j].resize(dim_three,0.0);
687 }
688 }
689 return the_vector;
690}
Domaine_IJK
This class encapsulates all the information related to the eulerian mesh for TrioIJK.
Definition Domaine_IJK.h:47
Domaine_IJK::get_offset_local
int get_offset_local(int direction) const
Returns the local offset in requested direction.
Definition Domaine_IJK.h:304
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Force_ph
Definition Force_ph.h:33
Force_ph::compute_energie
void compute_energie()
Definition Force_ph.cpp:610
Force_ph::get_force_attribute2
IJK_Field_vector3_double & get_force_attribute2()
Definition Force_ph.cpp:670
Force_ph::from_spect_to_phys_opti2
void from_spect_to_phys_opti2(ArrOfDouble &coeff_force)
Definition Force_ph.cpp:275
Force_ph::write_offset_index_position
void write_offset_index_position(const Domaine_IJK &my_splitting)
Definition Force_ph.cpp:578
Force_ph::get_energie
double get_energie()
Definition Force_ph.cpp:624
Force_ph::set_zero
void set_zero()
Definition Force_ph.cpp:98
Force_ph::write
void write(std::string nom_fichier_sortie, double t)
Definition Force_ph.cpp:549
Force_ph::from_spect_to_phys2
void from_spect_to_phys2(const std::vector< double > &coeff_force)
Definition Force_ph.cpp:145
Force_ph::get_force_attribute
IJK_Field_vector3_double get_force_attribute()
Definition Force_ph.cpp:665
Force_ph::write_separate
void write_separate(std::string nom_fichier_sortie, double t)
Definition Force_ph.cpp:629
Force_ph::cheat_function
void cheat_function()
Definition Force_ph.cpp:106
Force_ph::initialise
void initialise(int nproc_tot, int ni, int nj, int nk, int nl, int nm, int nn, double Lx, double Ly, double Lz, double Ox, double Oy, double Oz, int momin, int momax, double kmin, double kmax, std::string nom_fichier, const Domaine_IJK &splitting, int a_i_offset, int a_j_offset, int a_k_offset)
Definition Force_ph.cpp:40
Force_ph::from_spect_to_phys_opti2_advection
void from_spect_to_phys_opti2_advection(ArrOfDouble &coeff_force, const ArrOfDouble &advection_length)
Definition Force_ph.cpp:410
Objet_U
classe Objet_U Cette classe est la classe de base des Objets de TRUST
Definition Objet_U.h:73
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52