Athena: particles/init_particle.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*===========================================================================*/
00003 /*! \file init_particle.c
00004  *  \brief Initialize particle related structures and functions.
00005  *
00006  * PURPOSE: Initialize particle related structures and functions. Particle
00007  *   integrator is enrolled by calling integrate_particle_init(int type). In
00008  *   init_particle(Grid *pG, Domain *pD), all the particle related arrays are
00009  *   allocated, interpolation and stopping time calculation functions are
00010  *   enrolled, most global variables for particles are evaluated. Also contains
00011  *   functions for particle array reallocation and destruction.
00012  *
00013  * CONTAINS PUBLIC FUNCTIONS:
00014  * - integrate_particle_init();
00015  * - init_particle();
00016  * - particle_destruct();
00017  * - particle_realloc();
00018  *
00019  * History:
00020  * - Written by  Xuening Bai      Apr. 2009                                   */
00021 /*============================================================================*/
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <math.h>
00025 #include "../defs.h"
00026 #include "../athena.h"
00027 #include "../prototypes.h"
00028 #include "prototypes.h"
00029 #include "particle.h"
00030 #include "../globals.h"
00031 
00032 #ifdef PARTICLES         /* endif at the end of the file */
00033 
00034 /*==============================================================================
00035  * PRIVATE FUNCTION PROTOTYPES:
00036  *   grid_limit()   - get the limit coordinates of the grid
00037  *============================================================================*/
00038 void grid_limit(Grid *pG, Domain *pD);
00039 
00040 
00041 /*=========================== PUBLIC FUNCTIONS ===============================*/
00042 /*----------------------------------------------------------------------------*/
00043 
00044 /*----------------------------------------------------------------------------*/
00045 /*! \fn void init_particle(Grid *pG, Domain *pD)
00046  *  \brief Initialization for particles
00047  *
00048  * Allocate memory for the gas velocity/sound speed array, feedback array.
00049  * Note we enforce that each type has equal number of particles to ensure equal
00050  * resolution.
00051  */
00052 void init_particle(Grid *pG, Domain *pD)
00053 {
00054   int i, N1T, N2T, N3T, integratortype, interp, tsmode;
00055   Grain *GrArray;
00056   long size = 1000, size1 = 1, size2 = 1;
00057 
00058 /* get coordinate limit */
00059   grid_limit(pG, pD);
00060   N1T = iup-ilp+1;
00061   N2T = jup-jlp+1;
00062   N3T = kup-klp+1;
00063 
00064 /* check particle types */
00065   pG->partypes = par_geti_def("particle","partypes",1);
00066 
00067   if (pG->partypes < 0)
00068     ath_error("[init_particle]: Particle types must not be negative!\n");
00069 
00070 /* initialize the particle array */
00071   if(par_exist("particle","parnumcell"))
00072   {
00073     /* if we consider number of particles per cell */
00074     size1 = N1T*N2T*N3T*(long)(pG->partypes*par_geti("particle","parnumcell"));
00075     if (size1 < 0)
00076       ath_error("[init_particle]: Particle number must not be negative!\n");
00077   }
00078 
00079   if(par_exist("particle","parnumgrid"))
00080   {
00081     /* if we consider number of particles per grid */
00082     size2 = (long)(pG->partypes*par_geti("particle","parnumgrid"));
00083     if (size2 < 0)
00084       ath_error("[init_particle]: Particle number must not be negative!\n");
00085     /* account for the ghost cells */
00086     size2 = (long)(size2/((double)(pG->Nx1*pG->Nx2*pG->Nx3))*N1T*N2T*N3T);
00087   }
00088 
00089   size = MAX(size, MAX(size1, size2));
00090   pG->arrsize = (long)(1.2*size);       /* account for number fluctuations */
00091 
00092   pG->particle = (Grain*)calloc_1d_array(pG->arrsize, sizeof(Grain));
00093   if (pG->particle == NULL) goto on_error;
00094 
00095   pG->parsub   = (GrainAux*)calloc_1d_array(pG->arrsize, sizeof(GrainAux));
00096   if (pG->parsub == NULL) goto on_error;
00097 
00098 /* allocate memory for particle properties */
00099   pG->grproperty = (Grain_Property*)calloc_1d_array(pG->partypes,
00100                                                     sizeof(Grain_Property));
00101   if (pG->grproperty == NULL) goto on_error;
00102 
00103   grrhoa = (Real*)calloc_1d_array(pG->partypes, sizeof(Real));
00104   if (grrhoa == NULL) goto on_error;
00105 
00106   tstop0 = (Real*)calloc_1d_array(pG->partypes, sizeof(Real));
00107   if (tstop0 == NULL) goto on_error;
00108 
00109 /* by default these global values are zero */
00110   for (i=0; i<pG->partypes; i++) {
00111     tstop0[i] = 0.0;
00112     grrhoa[i] = 0.0;
00113   }
00114   alamcoeff = 0.0;
00115 
00116 /* Set particle integrator to according to particle types */
00117   for (i=0; i<pG->partypes; i++)
00118     pG->grproperty[i].integrator = par_geti_def("particle","integrator",2);
00119 
00120 /* set the interpolation function pointer */
00121   interp = par_geti_def("particle","interp",2);
00122   if (interp == 1)
00123   { /* linear interpolation */
00124     getweight = getwei_linear;
00125     ncell = 2;
00126   }
00127   else if (interp == 2)
00128   { /* TSC interpolation */
00129     getweight = getwei_TSC;
00130     ncell = 3;
00131   }
00132   else if (interp == 3)
00133   { /* Quadratic polynomial interpolation */
00134     getweight = getwei_QP;
00135     ncell = 3;
00136   }
00137   else
00138     ath_error("[init_particle]: Invalid interp value (should be 1, 2 or 3)!\n");
00139 
00140 /* set the stopping time function pointer */
00141   tsmode = par_geti("particle","tsmode");
00142   if (tsmode == 1)
00143     get_ts = get_ts_general;
00144   else if (tsmode == 2)
00145     get_ts = get_ts_epstein;
00146   else if (tsmode == 3)
00147     get_ts = get_ts_fixed;
00148   else
00149     ath_error("[init_particle]: tsmode must be 1, 2 or 3!\n");
00150 
00151 /* allocate the memory for gas-particle coupling array */
00152   pG->Coup = (GPCouple***)calloc_3d_array(N3T, N2T, N1T, sizeof(GPCouple));
00153   if (pG->Coup == NULL) goto on_error;
00154 
00155   return;
00156 
00157   on_error:
00158     ath_error("[init_particle]: Error allocating memory.\n");
00159 }
00160 
00161 /*----------------------------------------------------------------------------*/
00162 /*! \fn void particle_destruct(Grid *pG)
00163  *  \brief Finalization for particles
00164  */
00165 void particle_destruct(Grid *pG)
00166 {
00167   free_1d_array(pG->particle);
00168   free_1d_array(pG->parsub);
00169 
00170   free_1d_array(pG->grproperty);
00171   free_1d_array(grrhoa);
00172 
00173   /* free memory for gas and feedback arrays */
00174   if (pG->Coup != NULL) free_3d_array(pG->Coup);
00175 
00176   return;
00177 }
00178 
00179 /*----------------------------------------------------------------------------*/
00180 /*! \fn void particle_realloc(Grid *pG, long n)
00181  *  \brief Enlarge the particle array
00182  */
00183 void particle_realloc(Grid *pG, long n)
00184 {
00185   pG->arrsize = MAX((long)(1.2*pG->arrsize), n);
00186 
00187   /* for the main particle array */
00188   if ((pG->particle = (Grain*)realloc(pG->particle,
00189                                       pG->arrsize*sizeof(Grain))) == NULL)
00190   {
00191     ath_error("[init_particle]: Error re-allocating memory with array size\
00192  %ld.\n", n);
00193   }
00194 
00195   /* for the auxilary array */
00196   if ((pG->parsub = (GrainAux*)realloc(pG->parsub,
00197                                       pG->arrsize*sizeof(GrainAux))) == NULL)
00198   {
00199     ath_error("[init_particle]: Error re-allocating memory with array size\
00200  %ld.\n", n);
00201   }
00202 
00203   return;
00204 }
00205 
00206 /*============================================================================*/
00207 /*----------------------------- Private Functions ----------------------------*/
00208 
00209 /*----------------------------------------------------------------------------*/
00210 /*! \fn void grid_limit(Grid *pG, Domain *pD)
00211  *  \brief Calculate the left and right grid limit
00212  *
00213  * Input: pG: grid;
00214  * Output: ilp,iup,jlp,jup,klp,kup: grid limit indices;
00215  *         x1lpar,x1upar,x2lpar,x2upar,x3lpar,x3upar: grid boundary coordinates
00216  */
00217 void grid_limit(Grid *pG, Domain *pD)
00218 {
00219   int m1, m2, m3;       /* dimension flags */
00220   int my_iproc, my_jproc, my_kproc;
00221 
00222   if (pG->Nx1 > 1) m1 = 1;
00223   else m1 = 0;
00224 
00225   if (pG->Nx2 > 1) m2 = 1;
00226   else m2 = 0;
00227 
00228   if (pG->Nx3 > 1) m3 = 1;
00229   else m3 = 0;
00230 
00231 /* set left and right grid indices */
00232   ilp = pG->is - m1*nghost;
00233   iup = pG->ie + m1*nghost;
00234 
00235   jlp = pG->js - m2*nghost;
00236   jup = pG->je + m2*nghost;
00237 
00238   klp = pG->ks - m3*nghost;
00239   kup = pG->ke + m3*nghost;
00240 
00241 /* set left and right boundary for removing particles
00242  * Note: for outflow B.C. (ibc=2), we only remove the particles in
00243  * the outermost layer of the ghost cells
00244  */
00245   get_myGridIndex(pD, pG->my_id, &my_iproc, &my_jproc, &my_kproc);
00246 
00247   if ((par_geti_def("grid","ibc_x1",4) == 2) && (my_iproc == 0))
00248     x1lpar = pG->x1_0 + (ilp+m1 + pG->idisp)*pG->dx1;
00249   else
00250     x1lpar = pG->x1_0 + (pG->is + pG->idisp)*pG->dx1;
00251 
00252   if ((par_geti_def("grid","obc_x1",4) == 2) && (my_iproc == pD->NGrid_x1-1))
00253     x1upar = pG->x1_0 + (iup + pG->idisp)*pG->dx1;
00254   else
00255     x1upar = pG->x1_0 + (pG->ie + 1 + pG->idisp)*pG->dx1;
00256 
00257   if ((par_geti_def("grid","ibc_x2",4) == 2) && (my_jproc == 0))
00258     x2lpar = pG->x2_0 + (jlp+m2 + pG->jdisp)*pG->dx2;
00259   else
00260     x2lpar = pG->x2_0 + (pG->js + pG->jdisp)*pG->dx2;
00261 
00262   if ((par_geti_def("grid","obc_x2",4) == 2) && (my_jproc == pD->NGrid_x2-1))
00263     x2upar = pG->x2_0 + (jup + pG->jdisp)*pG->dx2;
00264   else
00265     x2upar = pG->x2_0 + (pG->je + 1 + pG->jdisp)*pG->dx2;
00266 
00267   if ((par_geti_def("grid","ibc_x3",4) == 2) && (my_kproc == 0))
00268     x3lpar = pG->x3_0 + (klp+m3 + pG->kdisp)*pG->dx3;
00269   else
00270     x3lpar = pG->x3_0 + (pG->ks + pG->kdisp)*pG->dx3;
00271 
00272   if ((par_geti_def("grid","obc_x3",4) == 2) && (my_kproc == pD->NGrid_x3-1))
00273     x3upar = pG->x3_0 + (kup + pG->kdisp)*pG->dx3;
00274   else
00275     x3upar = pG->x3_0 + (pG->ke + 1 + pG->kdisp)*pG->dx3;
00276 
00277   if (pG->Nx1 == 1) x1upar += MAX(0.1*fabs(pG->x1_0), 1.);
00278   if (pG->Nx2 == 1) x2upar += MAX(0.1*fabs(pG->x2_0), 1.);
00279   if (pG->Nx3 == 1) x3upar += MAX(0.1*fabs(pG->x3_0), 1.);
00280 
00281   return;
00282 }
00283 
00284 #endif /*PARTICLES*/