TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Force_sp.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_sp.h>
17#include <fstream>
18#include <math.h>
19#include <Domaine_IJK.h>
20#include <communications.h>
21#include <MaillerParallel.h>
22#include <EChaine.h>
23#include <SChaine.h>
24#include <Probleme_base.h>
25#include <Domaine_VF.h>
26#include <IJK_Navier_Stokes_tools.h>
27
28Implemente_instanciable( Force_sp, "Force_sp", Objet_U ) ;
29
30Sortie& Force_sp::printOn( Sortie& os ) const
31{
32 Objet_U::printOn( os );
33 return os;
34}
35
36Entree& Force_sp::readOn( Entree& is )
37{
38 Objet_U::readOn( is );
39 return is;
40}
41
42void Force_sp::initialise(int a_nl,int a_nn,int a_nm, int a_momin, int a_momax, double a_kmin, double a_kmax, std::string nom_fichier
43 )
44{
45 nl = a_nl;
46 nm = a_nm;
47 nn = a_nn;
48 n_lmn = (2*nl+1)*(2*nm+1)*(2*nn+1);
49 momin = a_momin;
50 momax = a_momax;
51 kmin = a_kmin;
52 kmax = a_kmax;
53 amplitude = 0.;
54
55
56// force[0] -> Partie reelle;
57// force[1] -> partie imaginaire;
58 force.resize(2);
59// force[_][i] -> Composant i = x, y ou z
60 force[0].resize(3);
61 force[1].resize(3);
62 force[0][0].resize(n_lmn,0.0);
63 force[1][0].resize(n_lmn,0.0);
64 force[0][1].resize(n_lmn,0.0);
65 force[1][1].resize(n_lmn,0.0);
66 force[0][2].resize(n_lmn,0.0);
67 force[1][2].resize(n_lmn,0.0);
68
69 force_flt.resize_array(2*3*n_lmn);
70
71 std::ofstream Spectral_flux(nom_fichier.c_str());
72 if(Spectral_flux)
73 {
74 Spectral_flux << "-------- SECTRAL_FORCE --------" << std::endl;
75 Spectral_flux << std::endl << "l,m,n \t : f_x, \tf_y, \tf_z,\t";
76 }
77
78 energie = 0.;
79
80}
81
82void Force_sp::initialise(int a_nl,int a_nn,int a_nm, int a_momin, int a_momax, double a_kmin, double a_kmax, double a_amplitude, std::string nom_fichier)
83{
84
85 nl = a_nl;
86 nm = a_nm;
87 nn = a_nn;
88 n_lmn = (2*nl+1)*(2*nm+1)*(2*nn+1);
89 momin = a_momin;
90 momax = a_momax;
91 kmin = a_kmin;
92 kmax = a_kmax;
93 amplitude = a_amplitude;
94
95// force[0] -> Partie reelle;
96// force[1] -> partie imaginaire;
97 force.resize(2);
98// force[_][i] -> Composant i = x, y ou z
99 force[0].resize(3);
100 force[1].resize(3);
101
102 force[0][0].resize(n_lmn,0.0);
103 force[1][0].resize(n_lmn,0.0);
104 force[0][1].resize(n_lmn,0.0);
105 force[1][1].resize(n_lmn,0.0);
106 force[0][2].resize(n_lmn,0.0);
107 force[1][2].resize(n_lmn,0.0);
108
109 force_flt.resize_array(2*3*n_lmn);
110
111 std::ofstream Spectral_flux(nom_fichier.c_str());
112 if(Spectral_flux)
113 {
114 Spectral_flux << "-------- SECTRAL_FORCE --------" << std::endl;
115 Spectral_flux << std::endl << "l,m,n \t : f_x, \tf_y, \tf_z,\t";
116 }
117
118}
119
120void Force_sp::field_advection(const ArrOfDouble& advection_length, const double dt, const int it)
121{
122 if (it>0)
123 {
124 if (Process::je_suis_maitre())
125 {
126 /*
127 Version avancée du calcul du processus aléatoire :
128 - Les champs sont stockés dans des tableaux à une dimension,
129 l'indice est reconstruit.
130 */
131 int n_dir(3);
132 double kappa[3];
133 double coefficient_translation[2];
134 for (int n=0; n<2*nn+1; n++)
135 for (int m=0; m<2*nm+1; m++)
136 for (int l=0; l<2*nl+1; l++) //for (int l=nl; l<2*nl+1; l++) // La TF-1 de cette force est réelle. Donc cette force est à symétrie Hermitienne
137 {
138 int ind_lmn = (n*(2*nm+1) + m) * (2*nl+1) +l;
139 // Position dans l'espace spectral. k_x va de -kmin a +kmin,
140 // d'ou le l-nl
141 kappa[0] = - kmax + (l)*(2*kmax)/(2*nl);
142 kappa[1] = - kmax + (m)*(2*kmax)/(2*nm);
143 kappa[2] = - kmax + (n)*(2*kmax)/(2*nn);
144
145 if (!(std::fabs(kappa[0])<kmin && std::fabs(kappa[1])<kmin && std::fabs(kappa[2])<kmin))
146 {
147 /* Le nombre d'onde est dans le domaine force*/
148
149
150 double argument_advection(0);
151 for (int dor=0; dor<n_dir; dor++)
152 {
153 /* Produit scalaire : kappa . advection_length */
154 Cout << dor << "argument_advection" << argument_advection << finl;
155 argument_advection += kappa[dor]*advection_length[dor];
156 }
157 Cout << "argument_advection" << argument_advection << finl;
158 coefficient_translation[0] = cos(-1.*argument_advection);
159 coefficient_translation[1] = sin(-1.*argument_advection);
160
161 for (int dir=0; dir<n_dir; dir++)
162 {
163 // Translation dans le domaine physique <=> multiplication par expo cpx en spectral
164 // TF(f(x-u Dt)) = exp(-i kappa u Dt) * TF(f(x))
165
166 {
167 int ind_CDIlmn((1*n_dir+dir)*n_lmn+ind_lmn);
168 int ind_RDIlmn((0*n_dir+dir)*n_lmn+ind_lmn);
169
170 force_flt[ind_RDIlmn] = force_flt[ind_RDIlmn]*coefficient_translation[0] - force_flt[ind_CDIlmn]*coefficient_translation[1];
171 force_flt[ind_CDIlmn] = force_flt[ind_CDIlmn]*coefficient_translation[0] + force_flt[ind_RDIlmn]*coefficient_translation[1];
172
173 // Symétrie Hermitienne : Sous nos conventions, -k[ind] = k[n_lmn - ind]
174// force_flt[(0*n_dir+dir)*n_lmn+(ind_lmn_moins)] = force_flt[(0*n_dir+dir)*n_lmn+(ind_lmn)];
175// force_flt[(1*n_dir+dir)*n_lmn+(ind_lmn_moins)] = - force_flt[(1*n_dir+dir)*n_lmn+(ind_lmn)];
176 }
177
178 }
179 }
180 }
181 }
182 envoyer_broadcast(force_flt,0);
183 }
184}
185
186void Force_sp::set_zero()
187{
188 /*Mise à zero de tout le champ de force spectral*/
189 n_lmn = (2*nl+1)*(2*nm+1)*(2*nn+1);
190
191 force_flt.resize(2*3*n_lmn);//,0.0);
192}
193
194// --------------------------------------
195void Force_sp::compute_force_kappa()
196{
197 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
198 if (Process::je_suis_maitre())
199 {
200 /*
201 Fonction de test, non opérationnelle pour le moment. Peut facilement
202 ettre opérationnelle.
203 Construit une force spectrale : f_sp(k) = k
204 */
205 int cpx, dir, l, m, n, ind;
206 double kappa[3];
207 double norme_kappa;
208 for (n=0; n<2*nn+1; n++)
209 {
210 for (m=0; m<2*nm+1; m++)
211 for (l=0; l<2*nl+1; l++) // La TF-1 de cette force est réelle. Donc cette force est à symétrie Hermitienne
212 {
213 ind = (n*(2*nm+1) + m) * (2*nl+1) +l;
214
215 kappa[0] = - kmax + (l)*(2*kmax)/(2*nl);
216 kappa[1] = - kmax + (m)*(2*kmax)/(2*nm);
217 kappa[2] = - kmax + (n)*(2*kmax)/(2*nn);
218 norme_kappa = sqrt(kappa[0]*kappa[0] + kappa[1]*kappa[1] + kappa[2]*kappa[2]);
219 if (norme_kappa >= kmin)
220 {
221 for (dir=0; dir<3; dir++)
222 for (cpx=0; cpx<2; cpx++)
223 {
224 force[cpx][dir][ind] = kappa[0];
225 }
226 }
227 }
228 }
229 }
230 envoyer_broadcast(force_flt,0);
231}
232
233void Force_sp::compute_dirac_board()
234{
235 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
236 if (Process::je_suis_maitre())
237 {
238 /*
239 Fonction de test. Non opérationnelle à ce jour
240 Construit une force spectral qui vaut 1 sur un plan et 0 ailleurs
241 */
242 int ind;
243 double kappa[3];
244 double norme_kappa;
245 for (int n=0; n<2*nn+1; n++)
246 {
247
248 for (int m=0; m<2*nm+1; m++)
249 for (int l=0; l<2*nl+1; l++) // La TF-1 de cette force est réelle. Donc cette force est à symétrie Hermitienne
250 {
251 ind = (n*(2*nm+1) + m) * (2*nl+1) +l;
252 kappa[0] = - kmax + (l)*(2*kmax)/(2*nl);
253 kappa[1] = - kmax + (m)*(2*kmax)/(2*nm);
254 kappa[2] = - kmax + (n)*(2*kmax)/(2*nn);
255 norme_kappa = sqrt(kappa[0]*kappa[0] + kappa[1]*kappa[1] + kappa[2]*kappa[2]);
256
257 if (norme_kappa >= kmin)
258 {
259 for (int dir=0; dir<3; dir++)
260 {
261 // ICI ON FORCE SUR LES 2 PLAQUES QUI SONT X=Xo. Or ce n'est pas ce qu'on veut !
262 if (l==nl+nl/2 || l==nl/2)
263 {
264 force[0][dir][ind] = amplitude;
265 force[1][dir][ind] = 0;
266 }
267 }
268 }
269 }
270 }
271 }
272 envoyer_broadcast(force_flt,0);
273}
274
275void Force_sp::compute_dirac_point()
276{
277 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
278 if (Process::je_suis_maitre())
279 {
280 /*
281 Fonciton de test. Operationnelle
282 Construit la force spectrale qui vaut 1 en deux point et 0 ailleurs
283 de sorte à ce que la fonction de forçage physique soit : cos(x+y) ex
284 */
285 int dir(0);
286
287 int roc(momin+1);
288
289 int l1(nl-roc), l2(nl+roc);
290 int n1(nn-roc), n2(nn+roc);
291 int fsp_m1(nm), fsp_m2(nm);
292
293 int ind1((n1*(2*nm+1) + fsp_m1) * (2*nl+1) +l1), ind2((n2*(2*nm+1) + fsp_m2) * (2*nl+1) +l2);
294 force[0][dir][ind1]=amplitude;
295 force[0][dir][ind2]=amplitude;
296 }
297 envoyer_broadcast(force_flt,0);
298}
299
300void Force_sp::compute_dirac_point_div_nulle()
301{
302 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
303 if (Process::je_suis_maitre())
304 {
305 /*
306 Fonction de test. Operationnelle
307 Construit la force spectrale qui vaut 1 en deux point et 0 ailleurs
308 de sorte à ce que la fonction de forçage physique soit : cos(x-z) (ex + ez)
309 f_sp(k) = dirac() + dirac()
310 */
311
312 int n_dir(3);
313 int roc(momin);
314 int n_ind_lmn(n_lmn);
315
316 int l_moins(nl-roc), l_plus(nl+roc);
317 int n_moins(nn-roc), n_plus(nn+roc);
318 int m_moins(nm);
319 int m_plus(nm);
320
321 int ind_moins((n_moins*(2*nm+1) + m_moins) * (2*nl+1) +l_plus);
322 int ind_plus((n_plus*(2*nm+1) + m_plus) * (2*nl+1) +l_moins);
323
324 for (int dir=0; dir<2; dir++)
325 {
326 dir = 2*dir;
327 double coeff(0);
328 if (dir==0)
329 coeff=1;
330 else
331 coeff=-1;
332
333 int ind_RDImoins((0*n_dir+dir)*n_ind_lmn+ind_moins);
334 int ind_RDIplus((0*n_dir+dir)*n_ind_lmn+ind_plus);
335 force_flt[ind_RDImoins]=amplitude*(1./sqrt(2))*coeff;
336 force_flt[ind_RDIplus]=amplitude*(1./sqrt(2))*coeff;
337
338 ind_RDImoins = (0*n_dir+dir+2)*n_ind_lmn+ind_moins;
339 ind_RDIplus = (0*n_dir+dir+2)*n_ind_lmn+ind_plus;
340 force_flt[ind_RDImoins]=amplitude*(1./sqrt(2))*coeff;
341 force_flt[ind_RDIplus]=amplitude*(1./sqrt(2))*coeff;
342 }
343 }
344 envoyer_broadcast(force_flt,0);
345}
346
347void Force_sp::compute_dirac_point_uniZ()
348{
349
350 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
351 if (Process::je_suis_maitre())
352 {
353 // donne cos(z) (ex)
354 int dir(2);
355 int n_dir(3);
356
357 int roc(momin+0);
358 int n_ind_lmn(n_lmn);
359
360 // int l_moins(nl-roc);
361 int l_zero(nl);
362 // int l_plus(nl+roc);
363 // int m_moins(nm-roc);
364 int m_zero(nm);
365 // ints m_plus(nm+roc);
366 int n_moins(nn-roc);
367 // int n_zero(nn);
368 int n_plus(nn+roc);
369
370 int ind_moins((n_moins*(2*nm+1) + m_zero) * (2*nl+1) +l_zero);
371 int ind_plus((n_plus*(2*nm+1) + m_zero) * (2*nl+1) +l_zero);
372
373
374 int ind_RDImoins((0*n_dir+dir)*n_ind_lmn+ind_moins);
375 int ind_RDIplus((0*n_dir+dir)*n_ind_lmn+ind_plus);
376 force_flt[ind_RDImoins]=amplitude;
377 force_flt[ind_RDIplus]=amplitude;
378
379 force[0][dir][ind_moins]=amplitude;
380 force[0][dir][ind_plus]=amplitude;
381 }
382 envoyer_broadcast(force_flt,0);
383}
384
385void Force_sp::compute_dirac_point_uniY()
386{
387 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
388 if (Process::je_suis_maitre())
389 {
390 // donne cos(momin*2pi/L*(y)) (ex)
391 int dir(1);
392 int n_dir(3);
393 int roc(momin+0);
394 int n_ind_lmn(n_lmn);
395
396 // int l_moins(nl-roc);
397 int l_zero(nl);
398 // int l_plus(nl+roc);
399 int m_moins(nm-roc);
400 // int m_zero(nm);
401 int m_plus(nm+roc);
402 // int n_moins(nn-roc);
403 int n_zero(nn);
404 // int n_plus(nn+roc);
405
406 int ind_moins((n_zero*(2*nm+1) + m_moins) * (2*nl+1) +l_zero);
407 int ind_plus((n_zero*(2*nm+1) + m_plus) * (2*nl+1) +l_zero);
408
409 int ind_RDImoins((0*n_dir+dir)*n_ind_lmn+ind_moins);
410 int ind_RDIplus((0*n_dir+dir)*n_ind_lmn+ind_plus);
411 force_flt[ind_RDImoins]=amplitude;
412 force_flt[ind_RDIplus]=amplitude;
413
414
415 force[0][dir][ind_moins]=amplitude;
416 force[0][dir][ind_plus]=amplitude;
417 }
418 envoyer_broadcast(force_flt,0);
419}
420
421void Force_sp::compute_diracs_for_cos_squarred()
422{
423 // -------------------------------------------
424 // cos(ax) <-------------------------- delta(k-a) + delta(k+a)
425 // cos(ax)cos(ax) = 0.5*(1+cos(2ax)) <- 0.5*(delta(k-2a) + 2delta(k) + delta(k+2a))
426 // -------------------------------------------
427 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
428 if (Process::je_suis_maitre())
429 {
430 // donne cos(momin*2pi/L*(y)) (ex)
431 int dir(1);
432 int n_dir(3);
433 int roc(momin);
434 int n_ind_lmn(n_lmn);
435
436 // int l_moins(nl-roc);
437 int l_zero(nl);
438 // int l_plus(nl+roc);
439 int m_moins(nm-2*roc); // k-2a
440 int m_zero(nm); // k
441 int m_plus(nm+2*roc); // k+2a
442 // int n_moins(nn-roc);
443 int n_zero(nn);
444 // int n_plus(nn+roc);
445
446 // - Stands for 0.5*(delta(k-2ko) + delta(k+2a))
447 int ind_moins((n_zero*(2*nm+1) + m_moins) * (2*nl+1) +l_zero);
448 int ind_plus((n_zero*(2*nm+1) + m_plus) * (2*nl+1) +l_zero);
449
450 int ind_RDImoins((0*n_dir+dir)*n_ind_lmn+ind_moins);
451 int ind_RDIplus((0*n_dir+dir)*n_ind_lmn+ind_plus);
452
453 force_flt[ind_RDImoins]=0.5*amplitude;
454 force_flt[ind_RDIplus]=0.5*amplitude;
455
456 force[0][dir][ind_moins]=0.5*amplitude;
457 force[0][dir][ind_plus]=0.5*amplitude;
458 // -----------------------------------------------------------
459
460 // - Stands for delta(k)
461 int ind_zero((n_zero*(2*nm+1) + m_zero) * (2*nl+1) +l_zero);
462
463 int ind_RDIzero((0*n_dir+dir)*n_ind_lmn+ind_zero);
464
465 force_flt[ind_RDIzero]=amplitude;
466
467 force[0][dir][ind_zero]=amplitude;
468 // -----------------------------------------------------------
469 }
470 envoyer_broadcast(force_flt,0);
471}
472
473void Force_sp::compute_diracs_for_t_times_cos_squarred(double time)
474{
475 // -------------------------------------------
476 // cos(ax) <-------------------------- delta(k-a) + delta(k+a)
477 // cos(ax)cos(ax) = 0.5*(1+cos(2ax)) <- 0.5*(delta(k-2a) + 2delta(k) + delta(k+2a))
478 // -------------------------------------------
479 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
480 if (Process::je_suis_maitre())
481 {
482 // donne cos(momin*2pi/L*(y)) (ex)
483 int dir(1);
484 int n_dir(3);
485 int roc(momin);
486 int n_ind_lmn(n_lmn);
487
488 // int l_moins(nl-roc);
489 int l_zero(nl);
490 // int l_plus(nl+roc);
491 int m_moins(nm-2*roc); // k-2a
492 int m_zero(nm); // k
493 int m_plus(nm+2*roc); // k+2a
494 // int n_moins(nn-roc);
495 int n_zero(nn);
496 // int n_plus(nn+roc);
497
498 // - Stands for 0.5*(delta(k-2ko) + delta(k+2a))
499 int ind_moins((n_zero*(2*nm+1) + m_moins) * (2*nl+1) +l_zero);
500 int ind_plus((n_zero*(2*nm+1) + m_plus) * (2*nl+1) +l_zero);
501
502 int ind_RDImoins((0*n_dir+dir)*n_ind_lmn+ind_moins);
503 int ind_RDIplus((0*n_dir+dir)*n_ind_lmn+ind_plus);
504
505 force_flt[ind_RDImoins]=time*0.5*amplitude;
506 force_flt[ind_RDIplus]=time*0.5*amplitude;
507
508 force[0][dir][ind_moins]=time*0.5*amplitude;
509 force[0][dir][ind_plus]=time*0.5*amplitude;
510 // -----------------------------------------------------------
511
512 // - Stands for delta(k)
513 int ind_zero((n_zero*(2*nm+1) + m_zero) * (2*nl+1) +l_zero);
514
515 int ind_RDIzero((0*n_dir+dir)*n_ind_lmn+ind_zero);
516
517 force_flt[ind_RDIzero]=time*amplitude;
518
519 force[0][dir][ind_zero]=time*amplitude;
520 // -----------------------------------------------------------
521 }
522 envoyer_broadcast(force_flt,0);
523}
524
525void Force_sp::compute_dirac_point_uniX_alongX()
526{
527 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
528 if (Process::je_suis_maitre())
529 {
530 // donne cos(x) (ex)
531 int dir(0);
532 int n_dir(3);
533 int n_ind_lmn(n_lmn);
534 int roc(momin+0);
535
536 int l_moins(nl-roc);
537 // int l_zero(nl);
538 int l_plus(nl+roc);
539 // int n_moins(nn-roc);
540 int n_zero(nn);
541 // int n_plus(nn+roc);
542 // int m_moins(nm-roc);
543 int m_zero(nm);
544 // ints m_plus(nm+roc);
545
546 int ind_moins((n_zero*(2*nm+1) + m_zero) * (2*nl+1) +l_moins);
547 int ind_plus((n_zero*(2*nm+1) + m_zero) * (2*nl+1) +l_plus);
548 std::cout << "ind_moins : " << ind_moins << std::endl;
549 std::cout << "ind_plus : " << ind_plus << std::endl;
550 std::cout << "in force_sp : " << Process::me() << std::endl;
551 int ind_RDImoins((0*n_dir+dir)*n_ind_lmn+ind_moins);
552 int ind_RDIplus((0*n_dir+dir)*n_ind_lmn+ind_plus);
553 force_flt[ind_RDImoins]=amplitude;
554 force_flt[ind_RDIplus]=amplitude;
555
556 force[0][dir][ind_moins]=amplitude;
557 force[0][dir][ind_plus]=amplitude;
558 }
559 envoyer_broadcast(force_flt,0);
560}
561
562void Force_sp::compute_dirac_point_uniX_alongY()
563{
564 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
565 if (Process::je_suis_maitre())
566 {
567 // donne cos(x) (ex)
568 int dir(1);
569 int n_dir(3);
570 int n_ind_lmn(n_lmn);
571 int roc(momin+0);
572
573 int l_moins(nl-roc);
574 // int l_zero(nl);
575 int l_plus(nl+roc);
576 // int n_moins(nn-roc);
577 int n_zero(nn);
578 // int n_plus(nn+roc);
579 // int m_moins(nm-roc);
580 int m_zero(nm);
581 // ints m_plus(nm+roc);
582
583 int ind_moins((n_zero*(2*nm+1) + m_zero) * (2*nl+1) +l_moins);
584 int ind_plus((n_zero*(2*nm+1) + m_zero) * (2*nl+1) +l_plus);
585 std::cout << "ind_moins : " << ind_moins << std::endl;
586 std::cout << "ind_plus : " << ind_plus << std::endl;
587 std::cout << "in force_sp : " << Process::me() << std::endl;
588 int ind_RDImoins((0*n_dir+dir)*n_ind_lmn+ind_moins);
589 int ind_RDIplus((0*n_dir+dir)*n_ind_lmn+ind_plus);
590 force_flt[ind_RDImoins]=amplitude;
591 force_flt[ind_RDIplus]=amplitude;
592
593 force[0][dir][ind_moins]=amplitude;
594 force[0][dir][ind_plus]=amplitude;
595 }
596 envoyer_broadcast(force_flt,0);
597}
598
599void Force_sp::compute_door_rope()
600{
601 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
602 if (Process::je_suis_maitre())
603 {
604 /*
605 Fonction de test. Met le champ de force spectral à 1 sur un segment défini.
606 Met le champ de force spectral à 0 ailleurs
607 */
608 int dir, l, m, n, ind;
609 int roc(momin+0);
610
611 // int l1(nl-roc), l2(nl+roc);
612 int fsp_m1(nm-roc), fsp_m2(nm+roc);
613 int n1(nn-roc), n2(nn+roc);
614
615 int l1(nl), l2(nl);
616 // int fsp_m1(nm), fsp_m2(nm);
617 // int n1(nn), n2(nn);
618
619 for (n=n1; n<n2+1; n++)
620 for (m=fsp_m1; m<fsp_m2+1; m++)
621 for (l=l1; l<l2+1; l++) // La TF-1 de cette force est réelle. Donc cette force est à symétrie Hermitienne
622 {
623 ind = (n*(2*nm+1) + m) * (2*nl+1) +l;
624 for (dir=0; dir<3; dir++)
625 {
626 force[0][dir][ind] = amplitude;
627 force[1][dir][ind] = 0.0;
628 }
629 }
630 }
631 envoyer_broadcast(force_flt,0);
632}
633
634void Force_sp::compute_door_cube()
635{
636 // Seul le proc 0 fait le calcul, comme c'est le mm clc pour tous les procs
637 if (Process::je_suis_maitre())
638 {
639 int cpx, dir, l, m, n, ind;
640
641 int roc(momin);
642
643 int l1(nl-roc), l2(nl+roc);
644 // int fsp_m1(nm-roc), fsp_m2(nm+roc);
645 int n1(nn-roc), n2(nn+roc);
646
647 // int l1(nl), l2(nl);
648 int fsp_m1(nm), fsp_m2(nm);
649 // int n1(nn), n2(nn);
650
651 for (n=n1; n<n2+1; n++)
652 for (m=fsp_m1; m<fsp_m2+1; m++)
653 for (l=l1; l<l2+1; l++) // La TF-1 de cette force est réelle. Donc cette force est à symétrie Hermitienne
654 {
655 ind = (n*(2*nm+1) + m) * (2*nl+1) +l;
656 for (dir=0; dir<3; dir++)
657 {
658 for (cpx=0; cpx<2; cpx++)
659 {
660 force[0][dir][ind] = amplitude;
661 force[1][dir][ind] = 0;
662 }
663 }
664 }
665 }
666 envoyer_broadcast(force_flt,0);
667}
668
669void Force_sp::compute_step2(ArrOfDouble& process_flt)
670{
671 if (Process::je_suis_maitre())
672 {
673 /*
674 Version avancée du calcul du processus aléatoire :
675 - Les champs sont stockés dans des tableaux à une dimension,
676 l'indice est reconstruit.
677 */
678 int n_dir(3);
679 double kappa[3];
680 double terme[2];
681 for (int n=0; n<2*nn+1; n++)
682 for (int m=0; m<2*nm+1; m++)
683 for (int l=nl; l<2*nl+1; l++) // La TF-1 de cette force est réelle. Donc cette force est à symétrie Hermitienne
684 {
685 int ind_lmn = (n*(2*nm+1) + m) * (2*nl+1) +l;
686 int ind_lmn_moins = (n*(2*nm+1) + m) * (2*nl+1) + (2*nl-l);
687 // Position dans l'espace spectral. k_x va de -kmin a +kmin,
688 // d'ou le l-nl
689 kappa[0] = - kmax + (l)*(2*kmax)/(2*nl);
690 kappa[1] = - kmax + (m)*(2*kmax)/(2*nm);
691 kappa[2] = - kmax + (n)*(2*kmax)/(2*nn);
692
693 if (std::fabs(kappa[0])<kmin && std::fabs(kappa[1])<kmin && std::fabs(kappa[2])<kmin)
694 {
695 /* Le nombre d'onde est trop petit pour etre dans le domaine force*/
696 for (int cpx=0; cpx<2; cpx++)
697 {
698 for (int dir=0; dir<3; dir++)
699 {
700 int ind_CDI((cpx*n_dir+dir)*n_lmn+ind_lmn);
701 force_flt[ind_CDI] = 0;
702 force_flt[(0*n_dir+dir)*n_lmn+(ind_lmn_moins)] = force_flt[(0*n_dir+dir)*n_lmn+(ind_lmn)];
703 force_flt[(1*n_dir+dir)*n_lmn+(ind_lmn_moins)] = - force_flt[(1*n_dir+dir)*n_lmn+(ind_lmn)];
704 }
705 }
706 }
707 else
708 {
709 /* On est dans le domaine force */
710 double norme_kappa = sqrt(kappa[0]*kappa[0] + kappa[1]*kappa[1] + kappa[2]*kappa[2]);
711 for (int cpx=0; cpx<2; cpx++)
712 {
713 terme[cpx] = kappa[0]*process_flt[(cpx*n_dir+0)*n_lmn+ind_lmn];
714 terme[cpx] += kappa[1]*process_flt[(cpx*n_dir+1)*n_lmn+ind_lmn];
715 terme[cpx] += kappa[2]*process_flt[(cpx*n_dir+2)*n_lmn+ind_lmn];
716 terme[cpx] /= pow(norme_kappa,2);
717 for (int dir=0; dir<3; dir++)
718 {
719 int ind_CDI((cpx*n_dir+dir)*n_lmn+ind_lmn);
720 force_flt[ind_CDI] = process_flt[ind_CDI] - kappa[dir]*terme[cpx];
721 // Symétrie Hermitienne : Sous nos conventions, -k[ind] = k[n_lmn - ind]
722 force_flt[(0*n_dir+dir)*n_lmn+(ind_lmn_moins)] = force_flt[(0*n_dir+dir)*n_lmn+(ind_lmn)];
723 force_flt[(1*n_dir+dir)*n_lmn+(ind_lmn_moins)] = - force_flt[(1*n_dir+dir)*n_lmn+(ind_lmn)];
724 }
725 }
726 }
727 }
728 }
729 envoyer_broadcast(force_flt,0);
730}
731
732// --------------------------------------
733void Force_sp::write(std::string nom_fichier_sortie, double t)
734{
735 // TODO : il faut certainement mettre du envoyer_broadcast ou quoi
736 std::ofstream Spectral_flux(nom_fichier_sortie.c_str(), std::ios::app);
737 if (Spectral_flux)
738 {
739 int l, m, n, ind;
740 double norme_kappa;
741 double kappa[3];
742 Spectral_flux << std::endl << "time : " << t << std::endl << std::endl;
743 for (n=0; n<2*nn+1; n++)
744 for (m=0; m<2*nm+1; m++)
745 for (l=0; l<2*nl+1; l++)
746 {
747 ind = (n*(2*nm+1) + m) * (2*nl+1) +l;
748 kappa[0] = - kmax + (l)*(2*kmax)/(2*nl);
749 kappa[1] = - kmax + (m)*(2*kmax)/(2*nm);
750 kappa[2] = - kmax + (n)*(2*kmax)/(2*nn);
751 norme_kappa = sqrt(kappa[0]*kappa[0] + kappa[1]*kappa[1] + kappa[2]*kappa[2]);
752 if (norme_kappa >= kmin)
753 {
754 Spectral_flux << l<<","<<m<<","<<n<<"\t : ";
755 Spectral_flux << force[0][0][ind] << " + i" << force[1][0][ind]<<", \t";
756 Spectral_flux << force[0][1][ind] << " + i" << force[1][1][ind]<<", \t";
757 Spectral_flux << force[0][2][ind] << " + i" << force[1][2][ind]<<", \t";
758 Spectral_flux << std::endl;
759 }
760 }
761 }
762}
763
764void Force_sp::write_separate(std::string nom_fichier_sortie, double t)
765{
766 std::ofstream Spectral_flux(nom_fichier_sortie.c_str());
767 if (Spectral_flux)
768 {
769 double kappa[3];
770 Spectral_flux << std::endl << "l,m,n, \t kx,ky,kz, \t rf_x, \t cf_x, \t\t rf_y, \t cf_y, \t\t rf_z, \t cf_z \t";
771 Spectral_flux << std::endl;
772
773 for (int n=0; n<2*nn+1; n++)
774 for (int m=0; m<2*nm+1; m++)
775 {
776 for (int l=0; l<2*nl+1; l++)
777 {
778 int ind = (n*(2*nm+1) + m) * (2*nl+1) +l;
779
780 kappa[0] = - kmax + (l)*(2*kmax)/(2*nl);
781 kappa[1] = - kmax + (m)*(2*kmax)/(2*nm);
782 kappa[2] = - kmax + (n)*(2*kmax)/(2*nn);
783 double norme_kappa = sqrt(kappa[0]*kappa[0] + kappa[1]*kappa[1] + kappa[2]*kappa[2]);
784 if (norme_kappa >= kmin)
785 {
786 Spectral_flux << l<<","<<m<<","<<n<<"\t,";
787 Spectral_flux << kappa[0] << ",\t" << kappa[1] << ",\t" << kappa[2] <<", \t\t";
788 Spectral_flux << force[0][0][ind] << ",\t" << force[1][0][ind]<<", \t\t";
789 Spectral_flux << force[0][1][ind] << ",\t" << force[1][1][ind]<<", \t\t";
790 Spectral_flux << force[0][2][ind] << ",\t" << force[1][2][ind]<<"\t";
791 Spectral_flux << std::endl;
792 }
793 }
794
795 }
796 }
797}
798
799// --------------------------------------
800void Force_sp::compute_energie()
801{
802 // Pas vérifié pour parallèle
803 int l,m,n,dir,ind;
804 energie = 0;
805 for (n=0; n<2*nn+1; n++)
806 for (m=0; m<2*nm+1; m++)
807 for (l=0; l<2*nl+1; l++)
808 {
809 ind = (n*(2*nm+1) + m) * (2*nl+1) +l;
810 for (dir=0; dir<3; dir++)
811 energie += force[0][dir][ind]*force[0][dir][ind] + force[1][dir][ind]*force[1][dir][ind];
812 }
813}
814
815// --------------------------------------
816double Force_sp::get_energie()
817{
818 return energie;
819}
820
821double Force_sp::get_force(int cpx, int dir, int ind)
822{
823 return force[cpx][dir][ind];
824}
825
826std::vector< std::vector< std:: vector <double >>> Force_sp::get_coeff()
827{
828 return force;
829}
830
831
832ArrOfDouble& Force_sp::get_coeff_flt()
833{
834 return force_flt;
835}
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Force_sp
Definition Force_sp.h:36
Force_sp::compute_dirac_point_div_nulle
void compute_dirac_point_div_nulle()
Definition Force_sp.cpp:300
Force_sp::initialise
void initialise(int nl, int nn, int nm, int momin, int momax, double kmin, double kmax, std::string nom_fichier)
Definition Force_sp.cpp:42
Force_sp::compute_force_kappa
void compute_force_kappa()
Definition Force_sp.cpp:195
Force_sp::set_zero
void set_zero()
Definition Force_sp.cpp:186
Force_sp::compute_dirac_point_uniX_alongY
void compute_dirac_point_uniX_alongY()
Definition Force_sp.cpp:562
Force_sp::get_coeff_flt
ArrOfDouble & get_coeff_flt()
Definition Force_sp.cpp:832
Force_sp::get_coeff
std::vector< std::vector< std::vector< double > > > get_coeff()
Definition Force_sp.cpp:826
Force_sp::get_force
double get_force(int cpx, int dir, int ind)
Definition Force_sp.cpp:821
Force_sp::compute_diracs_for_t_times_cos_squarred
void compute_diracs_for_t_times_cos_squarred(double time)
Definition Force_sp.cpp:473
Force_sp::compute_diracs_for_cos_squarred
void compute_diracs_for_cos_squarred()
Definition Force_sp.cpp:421
Force_sp::compute_dirac_board
void compute_dirac_board()
Definition Force_sp.cpp:233
Force_sp::compute_door_cube
void compute_door_cube()
Definition Force_sp.cpp:634
Force_sp::compute_dirac_point_uniX_alongX
void compute_dirac_point_uniX_alongX()
Definition Force_sp.cpp:525
Force_sp::compute_step2
void compute_step2(ArrOfDouble &process_flt)
Definition Force_sp.cpp:669
Force_sp::get_energie
double get_energie()
Definition Force_sp.cpp:816
Force_sp::write
void write(std::string nom_fichier_sortie, double time)
Definition Force_sp.cpp:733
Force_sp::compute_dirac_point_uniY
void compute_dirac_point_uniY()
Definition Force_sp.cpp:385
Force_sp::field_advection
void field_advection(const ArrOfDouble &advection_length, const double dt, const int it)
Definition Force_sp.cpp:120
Force_sp::compute_energie
void compute_energie()
Definition Force_sp.cpp:800
Force_sp::write_separate
void write_separate(std::string nom_fichier_sortie, double t)
Definition Force_sp.cpp:764
Force_sp::compute_door_rope
void compute_door_rope()
Definition Force_sp.cpp:599
Force_sp::compute_dirac_point
void compute_dirac_point()
Definition Force_sp.cpp:275
Force_sp::compute_dirac_point_uniZ
void compute_dirac_point_uniZ()
Definition Force_sp.cpp:347
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
Process::me
static int me()
renvoie mon rang dans le groupe de communication courant.
Definition Process.cpp:125
Process::je_suis_maitre
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Definition Process.cpp:86
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52