PLUTO
Classes | Macros | Typedefs | Functions
particles.h File Reference

Go to the source code of this file.

Classes

struct  PARTICLES
 
struct  head
 
struct  cell
 

Macros

#define NPARTICLES   256
 

Typedefs

typedef struct PARTICLES particle
 
typedef struct head HEAD
 
typedef struct cell CELL
 

Functions

void SET_PARTICLES (HEAD *list, HEAD *temp_list, struct GRID gxyz[], real ***uu[])
 
void ADVANCE_PARTICLES_predictor (HEAD *list, HEAD *temp_list, real ***uu[], struct GRID gxyz[], real dt, real striction)
 
void ADVANCE_PARTICLES_corrector (HEAD *list, HEAD *temp_list, real ***uu[], struct GRID gxyz[], real dt, real striction)
 
void SAVE_PARTICLES (HEAD *list, int nt, real t, real striction, real ***uu[])
 
void SAVE_STEP (HEAD *list, real dt, int nop)
 
void RENAME_PARTICLES (HEAD *list)
 
void ADD_CELL (HEAD *list, struct PARTICLES *pl)
 
void DELETE_CELL (HEAD *list, int i)
 
void DISPLAY_LIST (HEAD *list)
 

Macro Definition Documentation

#define NPARTICLES   256

Definition at line 5 of file particles.h.

Typedef Documentation

typedef struct cell CELL
typedef struct head HEAD
typedef struct PARTICLES particle

Function Documentation

void ADD_CELL ( HEAD *  list,
struct PARTICLES *  pl 
)

Definition at line 22 of file setparticles.c.

22  {
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 }
cell::next
struct cell * next
Definition: particles.h:26
cell
Definition: particles.h:24
cell::part
struct PARTICLES part
Definition: particles.h:25
head::first
struct cell * first
Definition: particles.h:21

Here is the caller graph for this function:

void ADVANCE_PARTICLES_corrector ( HEAD *  list,
HEAD *  temp_list,
real ***  uu[],
struct GRID  gxyz[],
real  dt,
real  striction 
)

Definition at line 436 of file setparticles.c.

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 }
a
static double a
Definition: init.c:135
LOCATE_PART
void LOCATE_PART(struct PARTICLES *pl, struct GRID gxyz[])
Definition: setparticles.c:1380
real
double real
Definition: pluto.h:488
cell::next
struct cell * next
Definition: particles.h:26
func_param::v
double v[NVAR]
Definition: eos.h:106
GRID::xf
double xf
Leftmost and rightmost point in the local domain.
Definition: structs.h:79
GRID::dx
double * dx
Definition: structs.h:83
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
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
cell::part
struct PARTICLES part
Definition: particles.h:25
PARTICLES::speed
real speed[DIMENSIONS]
Definition: particles.h:12
VX
#define VX
Definition: mod_defs.h:100
PARTICLES::cell
int cell[DIMENSIONS]
Definition: particles.h:10
PARTICLES::pression
real pression
Definition: particles.h:15
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
PARTICLES::divv
real divv
Definition: particles.h:16
GRID::xi
double xi
Definition: structs.h:79
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2

Here is the call graph for this function:

void ADVANCE_PARTICLES_predictor ( HEAD *  list,
HEAD *  temp_list,
real ***  uu[],
struct GRID  gxyz[],
real  dt,
real  striction 
)

Definition at line 199 of file setparticles.c.

200  {
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 }
a
static double a
Definition: init.c:135
LOCATE_PART
void LOCATE_PART(struct PARTICLES *pl, struct GRID gxyz[])
Definition: setparticles.c:1380
real
double real
Definition: pluto.h:488
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
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
cell::part
struct PARTICLES part
Definition: particles.h:25
PARTICLES::speed
real speed[DIMENSIONS]
Definition: particles.h:12
VX
#define VX
Definition: mod_defs.h:100
PARTICLES::cell
int cell[DIMENSIONS]
Definition: particles.h:10
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
GRID::xi
double xi
Definition: structs.h:79
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2

Here is the call graph for this function:

void DELETE_CELL ( HEAD *  list,
int  i 
)
void DISPLAY_LIST ( HEAD *  list)
void RENAME_PARTICLES ( HEAD *  list)

Definition at line 1136 of file setparticles.c.

1136  {
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 }
PARTICLES::identity
int identity
Definition: particles.h:9
cell::next
struct cell * next
Definition: particles.h:26
cell
Definition: particles.h:24
cell::part
struct PARTICLES part
Definition: particles.h:25
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
head::first
struct cell * first
Definition: particles.h:21

Here is the call graph for this function:

void SAVE_PARTICLES ( HEAD *  list,
int  nt,
real  t,
real  striction,
real ***  uu[] 
)

Definition at line 854 of file setparticles.c.

855  {
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 }
real
double real
Definition: pluto.h:488
cell::next
struct cell * next
Definition: particles.h:26
cell
Definition: particles.h:24
j
int j
Definition: analysis.c:2
cell::part
struct PARTICLES part
Definition: particles.h:25
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
head::first
struct cell * first
Definition: particles.h:21
PARTICLES::coor
real coor[DIMENSIONS]
Definition: particles.h:11

Here is the call graph for this function:

void SAVE_STEP ( HEAD *  list,
real  dt,
int  nop 
)

Definition at line 1346 of file setparticles.c.

1346  {
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 }
real
double real
Definition: pluto.h:488
cell::next
struct cell * next
Definition: particles.h:26
cell
Definition: particles.h:24
cell::part
struct PARTICLES part
Definition: particles.h:25
head::first
struct cell * first
Definition: particles.h:21
nop
int nop

Here is the caller graph for this function:

void SET_PARTICLES ( HEAD *  list,
HEAD *  temp_list,
struct GRID  gxyz[],
real ***  uu[] 
)

Definition at line 40 of file setparticles.c.

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 }
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::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
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
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
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
VX
#define VX
Definition: mod_defs.h:100
i
int i
Definition: analysis.c:2
nop
int nop
szz
Definition: al_hidden.h:38
GRID::xi
double xi
Definition: structs.h:79
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2

Here is the call graph for this function: