PLUTO
setparticles.c
Go to the documentation of this file.
1 #include "pluto.h"
2 
3 
4 extern int nop;
5 
6 void INTERPOL_VELOCITY (real *v_interpol,struct PARTICLES *pl,
7  real ***uu[], struct GRID gxyz[]);
8 
9 
10 void INTERPOL (real *a_interpol,struct PARTICLES *pl,
11  real ***uu[], struct GRID gxyz[], int VAR);
12 
13 void DIV ( real *diver_a, int VAR, real ***uu[], struct GRID gxyz[],
14  int indici[] );
15 
16 void LOCATE_PART ( struct PARTICLES *pl, struct GRID gxyz[]);
17 
18 void SAVE_STEP ( HEAD *list, real dt, int nop );
19 
20 /* We define a linked list for PARTICLES */
21 
22 void ADD_CELL(HEAD *list,struct PARTICLES *pl){
23 
24  CELL *new_cell;
25 
26  /* we always add cells at the head: it's easy and cheap */
27  /* it means that they are sorted in the exact opposite order of input */
28 
29  new_cell =(CELL *)malloc (sizeof(CELL));
30  new_cell->part= *pl ;
31  new_cell->next=list->first;
32  list->first=new_cell;
33 
34 }
35 
36 /* At end this part should be parametrizable (in order to choose
37  the type of initialization in space) and should stand in init.c */
38 
39 
40 void SET_PARTICLES (HEAD *list, HEAD *temp_list, struct GRID gxyz[],
41  real ***uu[] )
42 {
43 
44  int i , j , k ;
45  int coords[DIMENSIONS];
46  /* variables for interpolation */
47 
48  struct PARTICLES pl, pl_temp;
49  struct PARTICLES *p_l;
50  real *v_interpol,*lx;
51  real *b_interpol,*pr_interpol,*rho_interpol; /* bocchi */
52  real *div_v;
53  real dxA,dyA; /* variables for test1 */
54  int indici[3];
55 
56 #ifdef PARALLEL
57  int TAG=100;
58  int rank_part;
59  int cart_comm;
60  MPI_Datatype type_PARTICLES;
61  MPI_Datatype *type_ptr;
62 #endif
63 
64  lx=array_1D(3);
65  v_interpol = array_1D(3);
66  b_interpol = array_1D(3); /* bocchi */
67  pr_interpol = array_1D(1);
68  rho_interpol = array_1D(1);
69  div_v = array_1D(1);
70 
71  lx[0]=1./gxyz[0].dx[0];
72  lx[1]=1./gxyz[1].dx[0];
73  lx[2]=1./gxyz[2].dx[0];
74 
75  print1(" > Particles intialization \n");
76 
77  dxA = (0.8*(gxyz[0].xf-gxyz[0].xi))/(real)NPARTICLES;
78  dyA = (gxyz[1].xf-gxyz[1].xi)/3.;
79 
80  for(j=0;j<NPARTICLES;++j){
81 
82 #ifdef PARALLEL
83 
84  /*initialize the cartesian communicator of ArrayLib */
85  AL_Get_cart_comm(SZ, &cart_comm);
86 
87  if(prank==0){
88  type_ptr=&(type_PARTICLES);
89  pl.identity=j;
90  /* the random sequence is always the SAME since the seed is unchanged */
91  for(i=0;i<DIMENSIONS;++i){
92  pl.coor[i]=gxyz[i].xi+ (gxyz[i].xf-gxyz[i].xi)*(rand()/(1.+RAND_MAX) ) ;
93  pl.cell[i]=(int)( ( (pl.coor[i]) - gxyz[i].xi)*lx[i] ) + gxyz[i].nghost;
94  coords[i]=pl.cell[i];
95  p_l=&pl;
96  }
97 
98  BUILD_PARTICLES_TYPE(&pl, type_ptr);
99  MPI_Bcast( &pl, 1, type_PARTICLES, 0, MPI_COMM_WORLD );
100 
101  /* find the processor associated to particle position */
102  MPI_Cart_rank(cart_comm, coords, &rank_part);
103 
104  /* send the particle on its processor: rank_part */
105  MPI_SSend( pl,1,type_PARTICLES , rank_part, TAG, cart_comm );
106 
107  }
108  /* receive the particle on the processor rank_part */
109  if(prank==rank_part){
110  MPI_RECV( p_l,1,type_PARTICLES , 0, TAG,cart_comm );
111  INTERPOL ( v_interpol, &pl, uu, gxyz, VX);
112  for(i=0;i<DIMENSIONS;++i){pl.speed[i] = v_interpol[i] ;}
113  ADD_CELL(list, &pl);
114  }
115 
116 
117 #else
118 
119  pl.identity=j;
120 
121  for(i=0;i<DIMENSIONS;++i){
122  /* pl.coor[i]=gxyz[i].xi+ (gxyz[i].xf-gxyz[i].xi)*(rand()/(1.+RAND_MAX) ) ;*/
123  /* pl.coor[i]=-1.+1.*(rand()/(1.+RAND_MAX) ) ;*/
124  /* pl.cell[i]=(int)( ( (pl.coor[i]) - gxyz[i].xi)*lx[i] ) + gxyz[i].nghost; */
125  }
126 
127  /* bocchi - special initial conditions */
128  /* all particles on the jet : r=1.0 */
129  /* pl.coor[0]=gxyz[0].xi + (rand()/(1.+RAND_MAX));
130  pl.coor[1]=gxyz[1].xi+ (gxyz[1].xf-gxyz[1].xi)*(rand()/(1.+RAND_MAX) ) ;*/
131 
132  /* bocchi - initial conditions for test1 */
133 
134  if(j<(NPARTICLES/2)){
135 
136  pl.coor[0]=gxyz[0].xi + 0.1*(gxyz[0].xf-gxyz[0].xi) + j*dxA;
137  pl.coor[1]=gxyz[1].xi + dyA;
138 
139  }else{
140 
141  pl.coor[0]=gxyz[0].xi + 0.1*(gxyz[0].xf-gxyz[0].xi) +
142  (j-NPARTICLES/2)*dxA;
143  pl.coor[1]=gxyz[1].xi + 2.*dyA;
144 
145  }
146 
147  /* bocchi */
148 
149  /* printf("\nix_p:%d iy_p:%d\n", pl.cell[0], pl.cell[1]);*/
150  LOCATE_PART(&pl, gxyz);
151  /* printf("\nix_d:%d iy_d:%d\n", pl.cell[0], pl.cell[1]);*/
152 
153  INTERPOL ( v_interpol, &pl, uu, gxyz, VX);
154  INTERPOL ( b_interpol, &pl, uu, gxyz, BX);
155  INTERPOL ( pr_interpol, &pl, uu, gxyz, PR);
156  INTERPOL ( rho_interpol,&pl, uu, gxyz, DN);
157 
158  /* printf("\nic_prima:%d\n", pl.cell[0]); */
159 
160  for(i=0;i<DIMENSIONS;++i){
161  indici[i] = pl.cell[i];
162  }
163 
164  DIV ( div_v, VX, uu, gxyz, indici );
165 
166  for(i=0;i<DIMENSIONS;++i){
167  pl.speed[i] = v_interpol[i];
168  pl.mag[i] = b_interpol[i];
169  }
170  pl.pression = *pr_interpol;
171  pl.density = *rho_interpol;
172  pl.divv = *div_v;
173 
174  /* bocchi - end */
175 
176  ADD_CELL(list, &pl);
177  ADD_CELL(temp_list, &pl_temp);
178 
179  }
180 
181  nop=NPARTICLES;
182 
183  SAVE_STEP(list, 0.0, nop);
184 
185 #endif
186 
187  free_array_1D(v_interpol);
188  free_array_1D(b_interpol); /* bocchi */
189  free_array_1D(pr_interpol);
190  free_array_1D(rho_interpol);
191  free_array_1D(div_v);
192  free_array_1D(lx);
193 
194 
195 }
196 
197 
198 
199 void ADVANCE_PARTICLES_predictor( HEAD *list, HEAD *temp_list, real ***uu[],
200  struct GRID gxyz[], real dt, real striction){
201 
202  int i ;
203 
204  /* variable for interpolation*/
205  real *v_interpol;
206 
207  /* bocchi - RK2 - variable for interpolation */
208  /* real *v_int_old;
209  struct PARTICLES *pl_dummy; */
210  /* bocchi - end */
211 
212  /* variables for RK2*/
213  real *a,*v1,*x1;
214  real *lx;
215  real *gravity;
216 
217  /* Variables to handle data structure*/
218  struct PARTICLES *pl,*pl_temp;
219  CELL *actual_cell, *temp_actual_cell;
220  CELL *previous, *deleted; /* bocchi */
221  CELL *temp_previous, *temp_deleted;
222 
223  gravity= array_1D(3);
224  v_interpol = array_1D(3);
225  /* v_int_old = array_1D(3); */ /* bocchi */
226  a = array_1D(3);
227  v1 = array_1D(3);
228  x1 = array_1D(3);
229  lx = array_1D(3);
230  lx[0]=1./gxyz[0].dx[0];
231  lx[1]=1./gxyz[1].dx[0];
232  lx[2]=1./gxyz[2].dx[0];
233 
234  if(list->first != NULL ){
235 
236  actual_cell=list->first;
237  temp_actual_cell=temp_list->first;
238 
239  previous=actual_cell;
240  temp_previous=temp_actual_cell;
241 
242  while(actual_cell !=NULL){
243 
244  pl=&(actual_cell->part);
245  pl_temp=&(temp_actual_cell->part);
246  INTERPOL ( v_interpol, pl, uu, gxyz, VX);
247 /*
248  GRAVITY_FORCE (pl->coor, gravity);
249 */
250 
251  /*-----------------------------------*/
252  /* */
253  /* RK 2 scheme to move particles */
254  /* */
255  /*-----------------------------------*/
256 
257  for(i=0;i<2;++i){
258  /* a[i] = - striction*( pl->speed[i] - v_interpol[i] )+ gravity[i];
259  v1[i] = pl->speed[i] + dt*a[i];
260  x1[i] = pl->coor[i] + dt*pl->speed[i];*/ /*rem by bocchi*/
261 
262  /* bocchi - euler method */
263  /* v1[i] = v_interpol[i];
264  x1[i] = pl->coor[i] + dt*pl->speed[i]; */
265  /* bocchi - end */
266 
267  /* bocchi - RK2 */
268  /* v1[i] = v_interpol[i];
269  x1[i] = pl->coor[i] + 0.5*dt*v_interpol[i]; */
270  /* bocchi - end */
271 
272  /* bocchi - RK2 - alternative */
273  v1[i] = v_interpol[i];
274  x1[i] = pl->coor[i] + dt*pl->speed[i];
275  /*bocchi -RK2 - alternative - end */
276 
277  }
278 
279 
280  /*---------------------------------------------------*/
281  /* */
282  /* bocchi - Boundary Conditions */
283  /* xleft: mirror - xright: outflow */
284  /* yleft(bottom): periodic - yright(top): periodic */
285  /* */
286  /*---------------------------------------------------*/
287 
288  /* check x outflow first */
289 
290  if(x1[0]>gxyz[0].xf){
291  /* DESTROY ACTUAL PARTICLE */
292 
293  deleted=actual_cell;
294  temp_deleted=temp_actual_cell;
295 
296  previous->next=actual_cell->next;
297  temp_previous->next=temp_actual_cell->next;
298 
299  if(actual_cell==list->first ){
300  list->first=list->first->next;
301  temp_list->first=temp_list->first->next;
302  }
303 
304  if(actual_cell->next == NULL && actual_cell==list->first ){
305  actual_cell=NULL;
306  temp_actual_cell=NULL;
307  }else{
308  actual_cell=actual_cell->next;
309  temp_actual_cell=temp_actual_cell->next;
310  }
311 
312  free(deleted);
313  free(temp_deleted);
314 
315  --nop; /* correct the number of particles */
316 
317  continue; /* jump to the end of the loop */
318  }
319 
320  else{
321 
322  if(x1[0]<gxyz[0].xi){
323  x1[0] = 2*gxyz[0].xi - x1[0];
324  v1[0] = -v1[0];
325  }
326 
327 
328  /* check y axis */
329 
330  if(x1[1]>gxyz[1].xf){ x1[1] = gxyz[1].xi + x1[1] - gxyz[1].xf; }
331 
332  if(x1[1]<gxyz[1].xi){ x1[1] = gxyz[1].xf + x1[1] - gxyz[1].xi; }
333 
334  }
335 
336  /*-----------------------*/
337  /* bocchi - Boundary end */
338  /*-----------------------*/
339 
340 
341 
342  for(i=0;i<DIMENSIONS;++i)
343  {
344  pl_temp->coor[i]=x1[i] ;
345  /* pl_temp.cell[i]=(int)( (x1[i]-gxyz[i].xi)*lx[i] ) + gxyz[i].nghost;*/
346  pl_temp->speed[i]=v1[i];
347  }
348 
349  LOCATE_PART (pl_temp, gxyz); /* bocchi */
350 
351  /* bocchi - RK2 - fluid velocity in x1 at time t */
352  /*
353  pl_dummy = &pl_temp;
354  INTERPOL ( v_int_old, pl_dummy, uu, gxyz, VX);
355  for(i=0;i<DIMENSIONS;++i)
356  {
357  pl_temp.v_old[i]=v_int_old[i];
358  }
359  */
360  /* bocchi - end */
361 
362 
363 
364 
365 
366 # if( GEOMETRY==SPHERICAL )
367 
368 # if( DIMENSIONS>1)
369  /* special conditions for theta-boundary crossing*/
370  if( x1[1]>gxyz[1].xf ){ x1[1] = x1[1]-gxyz[1].xf; }
371  if( x1[1]<gxyz[1].xi ){ x1[1] = gxyz[1].xi-x1[1]; }
372 
373  pl_temp->coor[1]=x1[1];
374  pl_temp->cell[1]=(int)( (x1[1]-gxyz[1].xi)*lx[1] ) + gxyz[1].nghost;
375  pl_temp->speed[1]=v1[1] ;
376 
377 #endif
378 # if( DIMENSIONS>2)
379  /* special conditions for phi-boundary crossing*/
380  if( x1[2]>gxyz[2].xf ){ x1[2] = x1[2]-gxyz[2].xf; }
381  if( x1[2]<gxyz[2].xi ){ x1[2] = gxyz[2].xi-x1[2]; }
382  pl_temp->coor[2]=x1[2];
383  pl_temp->cell[2]=(int)( (x1[2]-gxyz[2].xi)*lx[2] ) + gxyz[2].nghost;
384  pl_temp->speed[2]=v1[2] ;
385 #endif
386 #endif
387 
388 # if( GEOMETRY==POLAR )
389 # if( DIMENSIONS>1)
390  /* special conditions for theta-boundary crossing*/
391  if( x1[1]>gxyz[1].xf ){ x1[1] = x1[1]-gxyz[1].xf; }
392  if( x1[1]<gxyz[1].xi ){ x1[1] = gxyz[1].xi-x1[1]; }
393 
394  pl_temp->coor[1]=x1[1];
395  pl_temp->cell[1]=(int)( (x1[1]-gxyz[1].xi)*lx[1] ) + gxyz[1].nghost;
396  pl_temp->speed[1]=v1[1] ;
397 
398 
399 #endif
400 #endif
401 # if( GEOMETRY== CYLINDRICAL )
402 # if( DIMENSIONS==3)
403  /* special conditions for theta-boundary crossing*/
404  if( x1[2]>gxyz[2].xf ){ x1[2] = x1[2]-gxyz[2].xf; }
405  if( x1[2]<gxyz[2].xi ){ x1[2] = gxyz[2].xi-x1[2]; }
406 
407  pl_temp->coor[2]=x1[2] ;
408  pl_temp->cell[2]=(int)( (x1[2]-gxyz[2].xi)*lx[2] ) + gxyz[2].nghost;
409  pl_temp->speed[2]=v1[2];
410 #endif
411 #endif
412 /*
413  ADD_CELL(temp_list, &pl_temp);
414 */
415  previous=actual_cell;
416  actual_cell=actual_cell->next;
417 
418  temp_previous=temp_actual_cell;
419  temp_actual_cell=temp_actual_cell->next;
420 
421  }
422  }
423 
424  free_array_1D(gravity);
425  free_array_1D(a);
426  free_array_1D(v1);
427  free_array_1D(x1);
428  free_array_1D(lx);
429  free_array_1D(v_interpol);
430  /* free_array_1D(v_int_old); */ /* bocchi */
431 
432 }
433 
434 
435 
436 void ADVANCE_PARTICLES_corrector( HEAD *list, HEAD *temp_list, real ***uu[], struct GRID gxyz[]
437  , real dt, real striction)
438 {
439  int i ;
440 
441  /* variable for interpolation*/
442  real *v_interpol;
443  real *b_interpol,*rho_interpol,*pr_interpol; /* bocchi */
444  real *div_v;
445  int indici[3];
446 
447  /* variables for RK2*/
448  real *a,*v,*x;
449  real *lx;
450  real *gravity;
451 
452  /* Variables to handle data structure*/
453  struct PARTICLES *pl,*pl_temp;
454 
455  CELL *actual_cell, *previous, *deleted;
456  CELL *temp_actual_cell, *temp_previous, *temp_deleted;
457 
458  gravity= array_1D(3);
459  v_interpol = array_1D(3);
460  b_interpol = array_1D(3); /* bocchi */
461  rho_interpol = array_1D(1);
462  pr_interpol = array_1D(1);
463  div_v = array_1D(1);
464  a = array_1D(3);
465  v = array_1D(3);
466  x = array_1D(3);
467  lx = array_1D(3);
468 
469  lx[0]=1./gxyz[0].dx[0];
470  lx[1]=1./gxyz[1].dx[0];
471  lx[2]=1./gxyz[2].dx[0];
472 
473  nop = 0; /* bocchi */
474 
475  if(list->first != NULL ){
476 
477  actual_cell=list->first;
478  temp_actual_cell=temp_list->first;
479 
480  previous=actual_cell;
481  temp_previous=temp_actual_cell;
482 
483  while(actual_cell !=NULL){
484  pl=&(actual_cell->part);
485  pl_temp=&(temp_actual_cell->part);
486 
487  ++nop; /* bocchi */
488 
489  INTERPOL ( v_interpol, pl_temp, uu, gxyz, VX);
490 
491  /* print1("A u1:%f,u2:%f\n", v_interpol[0],v_interpol[1]);*/
492 /*
493  GRAVITY_FORCE (pl_temp->coor, gravity);
494 */
495 
496  /*-----------------------------------*/
497  /* */
498  /* RK 2 scheme to move particles */
499  /* */
500  /*-----------------------------------*/
501 
502  for(i=0;i<2;++i){
503  /* a[i] = - striction*( pl_temp->speed[i] - v_interpol[i] )+ gravity[i];
504  v[i] = 0.5*( pl_temp->speed[i] + pl->speed[i] + dt*a[i] );
505  x[i] = pl_temp->coor[i] + dt*0.5*(pl_temp->speed[i] + pl->speed[i] );*/ /*rem by bocchi*/
506 
507  /* bocchi - euler method */
508  /* v[i] = v_interpol[i];
509  x[i] = pl_temp->coor[i];*/
510  /* bocchi - end */
511 
512  /* bocchi - RK2 */
513  /* x[i] = pl->coor[i] + dt*0.5*( v_interpol[i] + pl_temp->v_old[i] );
514  v[i] = v_interpol[i]; */ /* dummy value */
515  /* bocchi - RK2 - end */
516 
517  /* bocchi - RK2 - alternative */
518  v[i] = v_interpol[i]; /* dummy value */
519  x[i] = pl->coor[i] + dt*0.5*(pl->speed[i] + v_interpol[i]);
520  /* bocchi - RK2 - alternative - end */
521 
522  }
523 
524  /*
525  if(pl->identity==10){
526  printf("\npart10:%lf,%lf",pl->speed[1],v_interpol[1]);
527  }
528  if(pl->identity==100){
529  printf("\npart100:%lf,%lf",pl->speed[1],v_interpol[1]);
530  }
531  */
532 
533 
534 
535  /*---------------------------------------------------*/
536  /* */
537  /* bocchi - Boundary Conditions */
538  /* xleft: mirror - xright: outflow */
539  /* yleft(bottom): periodic - yright(top): periodic */
540  /* */
541  /*---------------------------------------------------*/
542 
543  /* check x outflow first */
544 
545  if(x[0]>gxyz[0].xf){
546  /* DESTROY ACTUAL PARTICLE */
547 
548  deleted=actual_cell;
549  temp_deleted=temp_actual_cell;
550 
551  previous->next=actual_cell->next;
552  temp_previous->next=temp_actual_cell->next;
553 
554  if(actual_cell==list->first ){
555  list->first=list->first->next;
556  temp_list->first=temp_list->first->next;
557  }
558 
559  if(actual_cell->next == NULL && actual_cell==list->first ){
560  actual_cell=NULL;
561  temp_actual_cell=NULL;
562  }else{
563  actual_cell=actual_cell->next;
564  temp_actual_cell=temp_actual_cell->next;
565  }
566 
567  free(deleted);
568  free(temp_deleted);
569 
570  --nop; /* correct the number of particles */
571 
572  continue; /* jump to the end of the loop */
573  }
574 
575  else{
576 
577  if(x[0]<gxyz[0].xi){
578  x[0] = 2*gxyz[0].xi - x[0];
579  /* speed correction not needed: calculated after this part */
580  }
581 
582 
583  /* check y axis */
584 
585  if(x[1]>gxyz[1].xf){ x[1] = gxyz[1].xi + x[1] - gxyz[1].xf; }
586 
587  if(x[1]<gxyz[1].xi){ x[1] = gxyz[1].xf + x[1] - gxyz[1].xi; }
588 
589  }
590 
591  /*-----------------------*/
592  /* bocchi - Boundary end */
593  /*-----------------------*/
594 
595 
596 
597 
598 #if( GEOMETRY==CARTESIAN)
599 # if( DIMENSIONS==1)
600  if(x[0]>gxyz[0].xi && x[0]<gxyz[0].xf )
601 #endif
602 # if( DIMENSIONS==2)
603  if(x[0]>gxyz[0].xi && x[0]<gxyz[0].xf &&
604  x[1]>gxyz[1].xi && x[1]<gxyz[1].xf )
605 #endif
606 # if( DIMENSIONS==3)
607  if(x[0]>gxyz[0].xi && x[0]<gxyz[0].xf &&
608  x[1]>gxyz[1].xi && x[1]<gxyz[1].xf &&
609  x[2]>gxyz[2].xi && x[2]<gxyz[2].xf)
610 #endif
611  {
612  for(i=0;i<DIMENSIONS;++i)
613  {
614  pl->coor[i]=x[i] ;
615  pl->cell[i]=(int)( (x[i]-gxyz[i].xi)*lx[i] ) + gxyz[i].nghost;
616  pl->speed[i]=v[i];
617  }
618 
619 
620 
621 
622  /*print1(" > particle : %d, x1:%f, x2:%f \n",
623  actual_cell->part.identity,
624  actual_cell->part.coor[0],
625  actual_cell->part.coor[1] );*/
626 
627 
628  previous=actual_cell;
629  actual_cell=actual_cell->next;
630 
631  temp_previous=temp_actual_cell;
632  temp_actual_cell=temp_actual_cell->next;
633  }
634 
635 #endif
636 
637 
638 # if( GEOMETRY==SPHERICAL )
639 
640  if(x[0]>gxyz[0].xi && x[0]<gxyz[0].xf){
641  pl->coor[0]=x[0] ;
642  pl->cell[0]=(int)( (x[0]-gxyz[0].xi)*lx[0] ) + gxyz[0].nghost;
643  pl->speed[0]=v[0];
644 # if( DIMENSIONS>1)
645  /* special conditions for theta-boundary crossing*/
646  if( x[1]>gxyz[1].xf ){ x[1] = x[1]-gxyz[1].xf; }
647  if( x[1]<gxyz[1].xi ){ x[1] = gxyz[1].xi-x[1]; }
648  pl->coor[1]=x[1];
649  pl->cell[1]=(int)( (x[1]-gxyz[1].xi)*lx[1] ) + gxyz[1].nghost;
650  pl->speed[1]=v[1] ;
651 
652 #endif
653 # if( DIMENSIONS>2)
654  /* special conditions for phi-boundary crossing*/
655  if( x[2]>gxyz[2].xf ){ x[2] = x[2]-gxyz[2].xf; }
656  if( x[2]<gxyz[2].xi ){ x[2] = gxyz[2].xi-x[2]; }
657  pl->coor[2]=x[2];
658  pl->cell[2]=(int)( (x[2]-gxyz[2].xi)*lx[2] ) + gxyz[2].nghost;
659  pl->speed[2]=v[2] ;
660 #endif
661  previous=actual_cell;
662  actual_cell=actual_cell->next;
663  temp_previous=temp_actual_cell;
664  temp_actual_cell=temp_actual_cell->next;
665  }
666 
667 #endif
668 
669 
670 # if( GEOMETRY==POLAR )
671 
672 # if( DIMENSIONS<3)
673  if(x[0]>gxyz[0].xi && x[0]<gxyz[0].xf)
674 #endif
675 
676 # if( DIMENSIONS==3)
677  if(x[0]>gxyz[0].xi && x[0]<gxyz[0].xf &&
678  x[2]>gxyz[2].xi && x[2]<gxyz[2].xf )
679 #endif
680  {
681  pl->coor[0]=x[0] ;
682  pl->cell[0]=(int)( (x[0]-gxyz[0].xi)*lx[0] ) + gxyz[0].nghost;
683  pl->speed[0]=v[0];
684 
685 # if( DIMENSIONS==3)
686  pl->coor[2]=x[2] ;
687  pl->cell[2]=(int)( (x[2]-gxyz[2].xi)*lx[2] ) + gxyz[2].nghost;
688  pl->speed[2]=v[2];
689 #endif
690 
691 # if( DIMENSIONS>1)
692  /* special conditions for theta-boundary crossing*/
693  if( x[1]>gxyz[1].xf ){ x[1] = x[1]-gxyz[1].xf; }
694  if( x[1]<gxyz[1].xi ){ x[1] = gxyz[1].xi-x[1]; }
695 
696  pl->coor[1]=x[1];
697  pl->cell[1]=(int)( (x[1]-gxyz[1].xi)*lx[1] ) + gxyz[1].nghost;
698  pl->speed[1]=v[1] ;
699 
700  /*print1(" > particle : %d , r:%f, theta:%f \n",
701  actual_cell->part.identity,
702  actual_cell->part.coor[0],
703  57.29578*actual_cell->part.coor[1] );*/
704 
705 #endif
706  previous=actual_cell;
707  actual_cell=actual_cell->next;
708 
709  temp_previous=temp_actual_cell;
710  temp_actual_cell=temp_actual_cell->next;
711  }
712 #endif
713 
714 
715 # if( GEOMETRY==CYLINDRICAL )
716 
717  if(x[0]>gxyz[0].xi && x[0]<gxyz[0].xf)
718 
719 # if( DIMENSIONS==3)
720  if(x[0]>gxyz[0].xi && x[0]<gxyz[0].xf &&
721  x[1]>gxyz[1].xi && x[1]<gxyz[1].xf )
722 #endif
723  {
724  pl->coor[0]=x[0] ;
725  /* pl->cell[0]=(int)( (x[0]-gxyz[0].xi)*lx[0] ) + gxyz[0].nghost;*/
726  pl->speed[0]=v[0];
727 
728 # if( DIMENSIONS>1)
729  pl->coor[1]=x[1] ;
730  /* pl->cell[1]=(int)( (x[1]-gxyz[1].xi)*lx[1] ) + gxyz[1].nghost;*/
731  pl->speed[1]=v[1];
732 #endif
733 # if( DIMENSIONS>2)
734  /* special conditions for theta-boundary crossing*/
735  if( x[2]>gxyz[2].xf ){ x[2] = x[2]-gxyz[2].xf; }
736  if( x[2]<gxyz[2].xi ){ x[2] = gxyz[2].xi-x[2]; }
737  pl->coor[2]=x[2];
738  /* pl->cell[2]=(int)( (x[2]-gxyz[2].xi)*lx[2] ) + gxyz[2].nghost;*/
739  pl->speed[2]=v[2] ;
740 
741 #endif
742 
743  /* bocchi - RK2 - adjust velocity */
744 
745  LOCATE_PART ( pl, gxyz );
746 
747  INTERPOL ( v_interpol, pl, uu, gxyz, VX);
748  INTERPOL ( b_interpol, pl, uu, gxyz, BX);
749  INTERPOL ( rho_interpol,pl,uu, gxyz, DN);
750  INTERPOL ( pr_interpol, pl,uu, gxyz, PR);
751 
752  for(i=0;i<DIMENSIONS;++i){
753  indici[i] = pl->cell[i];
754  }
755 
756  DIV ( div_v, VX, uu, gxyz, indici );
757 
758  for(i=0;i<DIMENSIONS;++i)
759  {
760  pl->speed[i]=v_interpol[i];
761  pl->mag[i]=b_interpol[i];
762  }
763 
764  pl->pression=*pr_interpol;
765  pl->density=*rho_interpol;
766  pl->divv = *div_v;
767 
768 
769  /* bocchi - RK2 - end */
770 
771  /* print1("B u1:%f,u2:%f\n", v_interpol[0],v_interpol[1]);*/
772  /* print1("\n part:%d, v1:%f, v2:%f, b2:%lf, divv:%f\n",
773  actual_cell->part.identity,
774  actual_cell->part.speed[0],
775  actual_cell->part.speed[1],
776  actual_cell->part.mag[1],
777  actual_cell->part.divv); */ /* bocchi */
778 
779 
780 
781  /* print1(" > particule: %d, rayon: %f\n",
782  actual_cell->part.identity,
783  actual_cell->part.coor[0] );*/
784 
785  previous=actual_cell;
786  actual_cell=actual_cell->next;
787 
788  temp_previous=temp_actual_cell;
789  temp_actual_cell=temp_actual_cell->next;
790  }
791 
792 #endif
793  else{
794  /* enter here if only particle is out the domain*/
795  /*print1(" > Canceled particle : %d, r: %f, theta:%f\n",
796  actual_cell->part.identity,
797  x[0],
798  actual_cell->part.coor[1]);*/
799 
800  deleted=actual_cell;
801  temp_deleted=temp_actual_cell;
802 
803  previous->next=actual_cell->next;
804  temp_previous->next=temp_actual_cell->next;
805 
806  if(actual_cell==list->first ){
807  list->first=list->first->next;
808  temp_list->first=temp_list->first->next;
809  }
810 
811  if(actual_cell->next == NULL && actual_cell==list->first ){
812  actual_cell=NULL;
813  temp_actual_cell=NULL;
814  }else{
815  actual_cell=actual_cell->next;
816  temp_actual_cell=temp_actual_cell->next;
817  }
818 
819  free(deleted);
820  free(temp_deleted);
821  }
822  }
823  }
824 
825 
826  /* last thing to do: destruct all cells of the temporary list*/
827  /* if(temp_list->first != NULL ){
828 
829  temp_actual_cell=temp_list->first;
830  temp_previous=temp_actual_cell;
831 
832  while(temp_list->first !=NULL){
833  deleted=temp_list->first;
834  temp_list->first=temp_list->first->next;
835  free(deleted);
836  }
837 
838  }*/ /* now the temp_list is mem-allocked with the list so
839  no need to destroy it now */
840 
841  free_array_1D(gravity);
842  free_array_1D(a);
843  free_array_1D(v);
844  free_array_1D(x);
845  free_array_1D(lx);
846  free_array_1D(v_interpol);
847  free_array_1D(b_interpol); /* bocchi */
848  free_array_1D(rho_interpol);
849  free_array_1D(pr_interpol);
850  free_array_1D(div_v);
851 
852 
853 }
854 void SAVE_PARTICLES( HEAD *list , int nt, real t, real striction,
855  real ***uu[] ){
856 
857  FILE *pout ,*verif;
858  int j,i0,i1,count=0 ;
859  real u0, u1, v0, v1, r,theta;
860  struct PARTICLES *pl;
861  CELL *actual_cell;
862 
863  if(list->first!= NULL){
864  actual_cell=list->first;
865  if(nt==0){
866  print("WRITING Particles_0\n" );
867  pout = fopen("0.out", "w");
868  }
869  if(nt==50){
870  print("WRITING Particles_1\n" );
871  pout = fopen("100.out", "w");
872  }
873  if(nt==100){
874  print("WRITING Particles_2\n" );
875  pout = fopen("200.out", "w");
876  }
877  if(nt==150){
878  print("WRITING Particles_3\n" );
879  pout = fopen("300.out", "w");
880  }
881  while(actual_cell!= NULL) {
882 
883  count+=1;
884  pl=&(actual_cell->part);
885 
886  /*r = pl->coor[0];
887  theta = pl->coor[1];*/
888 
889  fprintf(pout," %f %f",pl->coor[0],pl->coor[1] );
890  fprintf(pout,"\n" );
891  actual_cell=actual_cell->next;
892  }
893 
894  print(">number of particles written:%d \n",count );
895 
896  fclose(pout);
897 
898  }
899 }
900 
901 void INTERPOL_VELOCITY (real *v_interpol,struct PARTICLES *pl, real ***uu[]
902  , struct GRID gxyz[]){
903 
904 
905  /* variables for interpolation*/
906  int i,il, ir, jr, jl ,kr, kl, ind0, ind1, ind2;
907  real ix, iy, iz;
908 
909  ind0=pl->cell[0];
910  ind1=pl->cell[1];
911 
912  /*-------------------------------------------------------*/
913  /* */
914  /* This Perform Lagrangian interpolation of velocity */
915  /* base of polynoms: < 1, x, y, xy > */
916  /* */
917  /*-------------------------------------------------------*/
918 
919  if(gxyz[0].x[ind0] >pl->coor[0]){
920  il = ind0-1;
921  ir = ind0;
922  }else{
923  il = ind0;
924  ir = ind0+1;
925  }
926  ix = ( pl->coor[0] - gxyz[0].x[il] )/( gxyz[0].x[ir] - gxyz[0].x[il] );
927 
928 #if ( DIMENSIONS == 2 ||DIMENSIONS ==3 )
929 
930  if(gxyz[1].x[ind1] >pl->coor[1]){
931  jl = ind1-1;
932  jr = ind1 ;
933  }else{
934  jl = ind1;
935  jr = ind1+1 ;
936  }
937  iy = ( pl->coor[1] - gxyz[1].x[jl] )/( gxyz[1].x[jr] - gxyz[1].x[jl] );
938 
939 #endif
940 
941 #if DIMENSIONS ==3
942  if(gxyz[2].x[ind2] >pl->coor[2]){
943  kl = ind2-1;
944  kr = ind2 ;
945  }else{
946  kl = ind2;
947  kr = ind2+1 ;
948  }
949  iz = ( pl->coor[2] - gxyz[2].x[kl] )/( gxyz[2].x[kr] - gxyz[2].x[kl] );
950 #endif
951 
952 #if DIMENSIONS == 1
953  v_interpol[0] = uu[VX1][0][0][il]
954  + ix*(uu[VX1][0][0][ir] - uu[VX1][0][0][il] ) ;
955 #endif
956 
957 #if DIMENSIONS == 2
958 
959  for(i=0;i<DIMENSIONS;++i){
960  v_interpol[i] =
961  (1.-ix)*( uu[VX+i][0][jl][il]+iy*(uu[VX+i][0][jr][il]-uu[VX+i][0][jl][il]))
962  + ix*(uu[VX+i][0][jl][ir]+iy*(uu[VX+i][0][jr][ir] -uu[VX+i][0][jl][ir]));
963  }
964 
965 #endif
966 
967 #if DIMENSIONS == 3
968 
969  /* We use a formulation which asks for 18 calls to variables and 8 products*/
970  for(i=0;i<DIMENSIONS;++i){
971  v_interpol[i] =
972  (1.-ix)*( uu[VX+i][kl][jl][il]
973  + iz*(uu[VX+i][kr][jl][il] - uu[VX+i][kl][jl][il] )
974  + iy *( (uu[VX+i][kl][jr][il] - uu[VX+i][kl][jl][il] )
975  +iz*( uu[VX+i][kr][jr][il] - uu[VX+i][kl][jr][il]
976  - uu[VX+i][kr][jl][il] + uu[VX+i][kl][jl][il])
977  )
978  )
979  + ix*( uu[VX+i][kl][jl][ir]
980  + iz*(uu[VX+i][kr][jl][ir] - uu[VX+i][kl][jl][ir] )
981  + iy *( (uu[VX+i][kl][jr][ir] - uu[VX+i][kl][jl][ir] )
982  +iz*( uu[VX+i][kr][jr][ir] - uu[VX+i][kl][jr][ir]
983  - uu[VX+i][kr][jl][ir] + uu[VX+i][kl][jl][ir])
984  )
985  );
986 
987  }
988 #endif
989 
990 }
991 
992 
993 
994 
995 void INTERPOL (real *a_interpol,struct PARTICLES *pl, real ***uu[]
996  , struct GRID gxyz[], int VAR){
997 
998 
999  /* variables for interpolation*/
1000  int i,il, ir, jr, jl ,kr, kl, ind0, ind1, ind2;
1001  real ix, iy, iz;
1002 
1003  ind0=pl->cell[0];
1004  ind1=pl->cell[1];
1005 
1006 
1007  /*-------------------------------------------------------*/
1008  /* */
1009  /* This Perform Lagrangian interpolation of velocity */
1010  /* base of polynoms: < 1, x, y, xy > */
1011  /* */
1012  /*-------------------------------------------------------*/
1013 
1014  if(gxyz[0].x[ind0] >pl->coor[0]){
1015  il = ind0-1;
1016  ir = ind0;
1017  }else{
1018  il = ind0;
1019  ir = ind0+1;
1020  }
1021  ix = ( pl->coor[0] - gxyz[0].x[il] )/( gxyz[0].x[ir] - gxyz[0].x[il] );
1022 
1023 #if ( DIMENSIONS == 2 ||DIMENSIONS ==3 )
1024 
1025  if(gxyz[1].x[ind1] >pl->coor[1]){
1026  jl = ind1-1;
1027  jr = ind1 ;
1028  }else{
1029  jl = ind1;
1030  jr = ind1+1 ;
1031  }
1032  iy = ( pl->coor[1] - gxyz[1].x[jl] )/( gxyz[1].x[jr] - gxyz[1].x[jl] );
1033 
1034 #endif
1035 
1036 #if DIMENSIONS ==3
1037  if(gxyz[2].x[ind2] >pl->coor[2]){
1038  kl = ind2-1;
1039  kr = ind2 ;
1040  }else{
1041  kl = ind2;
1042  kr = ind2+1 ;
1043  }
1044  iz = ( pl->coor[2] - gxyz[2].x[kl] )/( gxyz[2].x[kr] - gxyz[2].x[kl] );
1045 #endif
1046 
1047 
1048  if(VAR==VX||VAR==BX){
1049 
1050 
1051 #if DIMENSIONS == 1
1052  a_interpol[0] = uu[var][0][0][il]
1053  + ix*(uu[VAR][0][0][ir] - uu[VAR][0][0][il] ) ;
1054 #endif
1055 
1056 #if DIMENSIONS == 2
1057 
1058  for(i=0;i<DIMENSIONS;++i){
1059  a_interpol[i] =
1060  (1.-ix)*( uu[VAR+i][0][jl][il]+iy*(uu[VAR+i][0][jr][il]-uu[VAR+i][0][jl][il]))
1061  + ix*(uu[VAR+i][0][jl][ir]+iy*(uu[VAR+i][0][jr][ir] -uu[VAR+i][0][jl][ir]));
1062  }
1063 
1064 #endif
1065 
1066 #if DIMENSIONS == 3
1067 
1068  /* We use a formulation which asks for 18 calls to variables and 8 products*/
1069  for(i=0;i<DIMENSIONS;++i){
1070  a_interpol[i] =
1071  (1.-ix)*( uu[VAR+i][kl][jl][il]
1072  + iz*(uu[VAR+i][kr][jl][il] - uu[VAR+i][kl][jl][il] )
1073  + iy *( (uu[VAR+i][kl][jr][il] - uu[VAR+i][kl][jl][il] )
1074  +iz*( uu[VAR+i][kr][jr][il] - uu[VAR+i][kl][jr][il]
1075  - uu[VAR+i][kr][jl][il] + uu[VAR+i][kl][jl][il])
1076  )
1077  )
1078  + ix*( uu[VAR+i][kl][jl][ir]
1079  + iz*(uu[VAR+i][kr][jl][ir] - uu[VAR+i][kl][jl][ir] )
1080  + iy *( (uu[VAR+i][kl][jr][ir] - uu[VAR+i][kl][jl][ir] )
1081  +iz*( uu[VAR+i][kr][jr][ir] - uu[VAR+i][kl][jr][ir]
1082  - uu[VAR+i][kr][jl][ir] + uu[VAR+i][kl][jl][ir])
1083  )
1084  );
1085 
1086  }
1087 #endif
1088  }
1089  else{
1090 
1091 
1092 #if DIMENSIONS == 1
1093  *a_interpol = uu[VAR][0][0][il]
1094  + ix*(uu[VAR][0][0][ir] - uu[VAR][0][0][il] ) ;
1095 #endif
1096 
1097 #if DIMENSIONS == 2
1098 
1099 
1100  *a_interpol =
1101  (1.-ix)*( uu[VAR][0][jl][il]+iy*(uu[VAR][0][jr][il]-uu[VAR][0][jl][il]))
1102  + ix*(uu[VAR][0][jl][ir]+iy*(uu[VAR][0][jr][ir] -uu[VAR][0][jl][ir]));
1103 
1104 
1105 #endif
1106 
1107 #if DIMENSIONS == 3
1108 
1109  /* We use a formulation which asks for 18 calls to variables and 8 products*/
1110 
1111  *a_interpol =
1112  (1.-ix)*( uu[VAR][kl][jl][il]
1113  + iz*(uu[VAR][kr][jl][il] - uu[VAR][kl][jl][il] )
1114  + iy *( (uu[VAR][kl][jr][il] - uu[VAR][kl][jl][il] )
1115  +iz*( uu[VAR][kr][jr][il] - uu[VAR][kl][jr][il]
1116  - uu[VAR][kr][jl][il] + uu[VAR][kl][jl][il])
1117  )
1118  )
1119  + ix*( uu[VAR][kl][jl][ir]
1120  + iz*(uu[VAR][kr][jl][ir] - uu[VAR][kl][jl][ir] )
1121  + iy *( (uu[VAR][kl][jr][ir] - uu[VAR][kl][jl][ir] )
1122  +iz*( uu[VAR][kr][jr][ir] - uu[VAR][kl][jr][ir]
1123  - uu[VAR][kr][jl][ir] + uu[VAR][kl][jl][ir])
1124  )
1125  );
1126 #endif
1127 
1128  }
1129 
1130 }
1131 
1132 
1133 
1134 
1135 
1136 void RENAME_PARTICLES( HEAD *list ){
1137 
1138 
1139  /* This routine simply rename particles, starting from the root*/
1140 
1141  int count=0;
1142  struct PARTICLES *pl;
1143  CELL *actual_cell;
1144 
1145 
1146  if(list->first!= NULL){
1147  actual_cell=list->first;
1148 
1149  while(actual_cell!= NULL) {
1150 
1151  pl=&(actual_cell->part);
1152  pl->identity=count;
1153  count+=1;
1154  actual_cell=actual_cell->next;
1155  }
1156 
1157  print(" > particles re_named: %d \n",count );
1158  }
1159 }
1160 
1161 /*void SEND_PARTICLES(PARTICLES *pl, int proc_ini, int proc_end ){*/
1162 
1163 /* This Function sends a particle: */
1164 /* from processor proc_ini */
1165 /* to processor proc_end */
1166 
1167 
1168 
1169 
1170 /* bocchi */
1171 void DIV ( real *diver, int VAR, real ***uu[], struct GRID gxyz[],
1172  int indici[] ){
1173 
1174  /* perform divergence of the array_1D VAR: es VX */
1175 
1176  int il, ic, ir, jl, jc, jr, kl, kc, kr;
1177  real dx1, dx2, dy1, dy2, dz1, dz2;
1178  real ax, bx, cx, ay, by, cy, az, bz, cz;
1179 
1180 #if ( GEOMETRY == CYLINDRICAL )
1181 
1182 # if DIMENSIONS == 1
1183 
1184  ic = indici[0];
1185  il = ic - 1;
1186  ir = ic + 1;
1187 
1188  dx1 = gxyz[0].x[ic] - gxyz[0].x[il];
1189  dx2 = gxyz[0].x[ir] - gxyz[0].x[ic];
1190 
1191  ax = (dx2*dx2-dx1*dx1)/(dx1*dx2*(dx1+dx2));
1192  bx = -dx2/(dx1*(dx1+dx2));
1193  cx = dx1/(dx2*(dx1+dx2));
1194 
1195  *diver = ax*uu[VAR][0][0][ic] + bx*uu[VAR][0][0][il] + cx*uu[VAR][0][0][ir]
1196  + uu[VAR][0][0][ic]/gxyz[0].x[ic];
1197 
1198 #endif
1199 
1200 
1201 # if DIMENSIONS == 2
1202 
1203  ic = indici[0];
1204  il = ic - 1;
1205  ir = ic + 1;
1206 
1207  jc = indici[1];
1208  jl = jc - 1;
1209  jr = jc + 1;
1210 
1211  dx1 = gxyz[0].x[ic] - gxyz[0].x[il];
1212  dx2 = gxyz[0].x[ir] - gxyz[0].x[ic];
1213 
1214  dy1 = gxyz[1].x[jc] - gxyz[1].x[jl];
1215  dy2 = gxyz[1].x[jr] - gxyz[1].x[jc];
1216 
1217 
1218  ax = (dx2*dx2-dx1*dx1)/(dx1*dx2*(dx1+dx2));
1219  bx = -dx2/(dx1*(dx1+dx2));
1220  cx = dx1/(dx2*(dx1+dx2));
1221 
1222  ay = (dy2*dy2-dy1*dy1)/(dy1*dy2*(dy1+dy2));
1223  by = -dy2/(dy1*(dy1+dy2));
1224  cy = dy1/(dy2*(dy1+dy2));
1225 
1226  *diver = ax*uu[VAR][0][jc][ic] + bx*uu[VAR][0][jc][il]
1227  + cx*uu[VAR][0][jc][ir]
1228  + uu[VAR][0][jc][ic]/gxyz[0].x[ic]
1229  + ay*uu[VAR+1][0][jc][ic] + by*uu[VAR+1][0][jl][ic]
1230  + cy*uu[VAR+1][0][jr][ic];
1231 
1232 
1233 
1234 #endif
1235 
1236 # if DIMENSIONS == 3
1237 
1238  ic = pl->cell[0];
1239  il = ic - 1;
1240  ir = ic + 1;
1241 
1242  jc = pl->cell[1];
1243  jl = jc - 1;
1244  jr = jc + 1;
1245 
1246  kc = pl->cell[2];
1247  kl = kc - 1;
1248  kr = kc + 1;
1249 
1250  dx1 = gxyz[0].x[ic] - gxyz[0].x[il];
1251  dx2 = gxyz[0].x[ir] - gxyz[0].x[ic];
1252 
1253  dy1 = gxyz[1].x[jc] - gxyz[1].x[jl];
1254  dy2 = gxyz[1].x[jr] - gxyz[1].x[jc];
1255 
1256  dz1 = gxyz[2].x[kc] - gxyz[2].x[kl];
1257  dz2 = gxyz[2].x[kr] - gxyz[2].x[kc];
1258 
1259  ax = (dx2*dx2-dx1*dx1)/(dx1*dx2*(dx1+dx2));
1260  bx = -dx2/(dx1*(dx1+dx2));
1261  cx = dx1/(dx2*(dx1+dx2));
1262 
1263  ay = (dy2*dy2-dy1*dy1)/(dy1*dy2*(dy1+dy2));
1264  by = -dy2/(dy1*(dy1+dy2));
1265  cy = dy1/(dy2*(dy1+dy2));
1266 
1267  az = (dz2*dz2-dz1*dz1)/(dz1*dz2*(dz1+dz2));
1268  bz = -dz2/(dz1*(dz1+dz2));
1269  cz = dz1/(dz2*(dz1+dz2));
1270 
1271  *diver = ax*uu[VAR][kc][jc][ic] + bx*uu[VAR][kc][jc][il]
1272  + cx*uu[VAR][kc][jc][ir]
1273  + uu[VAR][kc][jc][ic]/gxyz[0].x[ic]
1274  + ay*uu[VAR+1][kc][jc][ic] + by*uu[VAR+1][kc][jl][ic]
1275  + cy*uu[VAR+1][kc][jr][ic]
1276  + (az*uu[VAR+2][kc][jc][ic] + bz*uu[VAR+2][kl][jc][ic]
1277  + cz*uu[VAR+2][kr][jc][ic])/gxyz[0].x[ic];
1278 
1279 #endif
1280 
1281 #endif
1282 
1283 
1284 
1285 
1286 }
1287 /* bocchi end */
1288 
1289 #ifdef PARALLEL
1290 void BUILD_PARTICLES_TYPE(struct PARTICLES* pl,MPI_Datatype* type_ptr){
1291 
1292  /* This function create a MPI type for the particles to ease */
1293  /* SEND and RECV procedures. */
1294 
1295  /* For information see "A User Guide to MPI" by Peter S. Pachero */
1296  /* (his responsible for my mistakes...) page 28-30 */
1297 
1298 
1299  int block_lengths[4];
1300  MPI_Aint displacements[4];
1301  MPI_Aint addresses[5];
1302  MPI_Datatype typelist[4];
1303 
1304  /* specify the types */
1305  typelist[0] = MPI_INT ;
1306  typelist[1] = MPI_DOUBLE ;
1307  typelist[2] = MPI_DOUBLE ;
1308  typelist[3] = MPI_INT ;
1309 
1310  /* specify the number of elements of each type */
1311 
1312  block_lengths[0] = 1;
1313  block_lengths[1] = DIMENSIONS;
1314  block_lengths[2] = DIMENSIONS;
1315  block_lengths[3] = DIMENSIONS;
1316 
1317  /* Calculate the displacements of the members */
1318 
1319  MPI_Address(pl , &addresses[0]);
1320  MPI_Address(&(pl->identity), &addresses[1]);
1321  MPI_Address(&(pl->coor) , &addresses[2]);
1322  MPI_Address(&(pl->speed) , &addresses[3]);
1323  MPI_Address(&(pl->cell) , &addresses[4]);
1324 
1325  displacements[0]=addresses[1]-addresses[0];
1326  displacements[1]=addresses[2]-addresses[0];
1327  displacements[2]=addresses[3]-addresses[0];
1328  displacements[3]=addresses[4]-addresses[0];
1329 
1330  /* Create the derived type */
1331 
1332  MPI_Type_struct(4,block_lengths,displacements,typelist,type_ptr);
1333 
1334  /* Commit it so it can be used */
1335 
1336  MPI_Type_commit(type_ptr);
1337 
1338 }
1339 
1340 #endif
1341 
1342 
1343 
1344 /* bocchi */
1345 
1346 void SAVE_STEP ( HEAD *list, real dt, int nop ){
1347 
1348  /* save structure particle, with dt and number of particles */
1349 
1350  FILE *stream;
1351  struct PARTICLES *pl;
1352  CELL *actual_cell;
1353 
1354  stream = fopen("part_data.bin.out", "a");
1355 
1356  fwrite(&nop, sizeof(int), 1, stream);
1357  fwrite(&dt, sizeof(real), 1, stream);
1358 
1359  if (list->first != NULL){
1360 
1361  actual_cell = list->first;
1362 
1363  while (actual_cell != NULL){
1364 
1365  pl = &(actual_cell->part);
1366 
1367  fwrite(pl, sizeof(particle), 1, stream);
1368 
1369  actual_cell = actual_cell->next;
1370 
1371  }
1372  }
1373 
1374  fclose(stream);
1375 
1376 }
1377 
1378 
1379 
1380 void LOCATE_PART ( struct PARTICLES *pl, struct GRID gxyz[] ){
1381 
1382 
1383  /* determine the cell of the particle, even in stretched grids */
1384 
1385 
1386  int l_ind, r_ind;
1387  int m_ind, i, j;
1388  int true = 1;
1389 
1390  for(i=0;i<DIMENSIONS;++i){
1391 
1392 
1393  l_ind = gxyz[i].nghost;
1394  r_ind = gxyz[i].np_int_glob + gxyz[i].nghost - 1;
1395 
1396 
1397  while ( true ){
1398 
1399  m_ind = l_ind + (r_ind - l_ind)/2;
1400 
1401  if (pl->coor[i] <= gxyz[i].xr[m_ind]){
1402  r_ind = m_ind;
1403  }
1404  else{
1405  l_ind = m_ind + 1;
1406  }
1407 
1408  if( l_ind == r_ind ) break;
1409 
1410  }
1411 
1412 
1413  /*now ind_left == ind_right == ind_part */
1414 
1415  pl->cell[i] = r_ind;
1416 
1417  }
1418 
1419 }
1420 
1421 /* bocchi - end */
PARTICLES::identity
int identity
Definition: particles.h:9
a
static double a
Definition: init.c:135
SET_PARTICLES
void SET_PARTICLES(HEAD *list, HEAD *temp_list, struct GRID gxyz[], real ***uu[])
Definition: setparticles.c:40
LOCATE_PART
void LOCATE_PART(struct PARTICLES *pl, struct GRID gxyz[])
Definition: setparticles.c:1380
real
double real
Definition: pluto.h:488
NPARTICLES
#define NPARTICLES
Definition: particles.h:5
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
GRID::xr
double * xr
Definition: structs.h:81
GRID::np_int_glob
int np_int_glob
Total number of points in the global domain (boundaries excluded).
Definition: structs.h:98
cell::next
struct cell * next
Definition: particles.h:26
GRID::xf
double xf
Leftmost and rightmost point in the local domain.
Definition: structs.h:79
GRID::dx
double * dx
Definition: structs.h:83
prank
int prank
Processor rank.
Definition: globals.h:33
VX1
#define VX1
Definition: mod_defs.h:28
ADVANCE_PARTICLES_predictor
void ADVANCE_PARTICLES_predictor(HEAD *list, HEAD *temp_list, real ***uu[], struct GRID gxyz[], real dt, real striction)
Definition: setparticles.c:199
SAVE_STEP
void SAVE_STEP(HEAD *list, real dt, int nop)
Definition: setparticles.c:1346
BX
#define BX
Definition: mod_defs.h:108
DIV
void DIV(real *diver_a, int VAR, real ***uu[], struct GRID gxyz[], int indici[])
Definition: setparticles.c:1171
PARTICLES::mag
real mag[DIMENSIONS]
Definition: particles.h:13
PARTICLES::density
real density
Definition: particles.h:14
AL_Get_cart_comm
int AL_Get_cart_comm(int, MPI_Comm *)
Definition: al_sz_get.c:117
ADD_CELL
void ADD_CELL(HEAD *list, struct PARTICLES *pl)
Definition: setparticles.c:22
head
Definition: particles.h:20
ADVANCE_PARTICLES_corrector
void ADVANCE_PARTICLES_corrector(HEAD *list, HEAD *temp_list, real ***uu[], struct GRID gxyz[], real dt, real striction)
Definition: setparticles.c:436
GRID
Definition: structs.h:78
cell
Definition: particles.h:24
INTERPOL
void INTERPOL(real *a_interpol, struct PARTICLES *pl, real ***uu[], struct GRID gxyz[], int VAR)
Definition: setparticles.c:995
j
int j
Definition: analysis.c:2
GRID::nghost
int nghost
Number of ghost zones.
Definition: structs.h:104
k
int k
Definition: analysis.c:2
cell::part
struct PARTICLES part
Definition: particles.h:25
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
PARTICLES::speed
real speed[DIMENSIONS]
Definition: particles.h:12
GRID::x
double * x
Definition: structs.h:80
VX
#define VX
Definition: mod_defs.h:100
INTERPOL_VELOCITY
void INTERPOL_VELOCITY(real *v_interpol, struct PARTICLES *pl, real ***uu[], struct GRID gxyz[])
Definition: setparticles.c:901
PARTICLES::cell
int cell[DIMENSIONS]
Definition: particles.h:10
PARTICLES::pression
real pression
Definition: particles.h:15
pluto.h
PLUTO main header file.
RENAME_PARTICLES
void RENAME_PARTICLES(HEAD *list)
Definition: setparticles.c:1136
i
int i
Definition: analysis.c:2
head::first
struct cell * first
Definition: particles.h:21
nop
int nop
PARTICLES::coor
real coor[DIMENSIONS]
Definition: particles.h:11
szz
Definition: al_hidden.h:38
PARTICLES::divv
real divv
Definition: particles.h:16
SAVE_PARTICLES
void SAVE_PARTICLES(HEAD *list, int nt, real t, real striction, real ***uu[])
Definition: setparticles.c:854
GRID::xi
double xi
Definition: structs.h:79
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2