Athena: prob/par_strat2d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file par_strat2d.c
00004  *  \brief Problem generator for non-linear streaming instability in
00005  *   stratified disks.
00006  *
00007  * PURPOSE: Problem generator for non-linear streaming instability in
00008  *   stratified disks. This code works in 2D ONLY. Isothermal eos is assumed,
00009  *   and the value etavk/iso_sound is fixed. MPI domain decomposition in x is
00010  *   allowed, but not in z.
00011  *
00012  * Perturbation modes:
00013  *  - ipert = 0: multi-nsh equilibtium
00014  *  - ipert = 1: white noise within the entire grid
00015  *  - ipert = 2: non-nsh velocity
00016  *
00017  *  Must be configured using --enable-shearing-box and --with-eos=isothermal.
00018  *  FARGO is recommended.
00019  *
00020  * Reference:
00021  *   Johansen & Youdin, 2007, ApJ, 662, 627
00022  *   Bai & Stone, 2009, in preparation                                        */
00023 /*============================================================================*/
00024 
00025 #include <float.h>
00026 #include <math.h>
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include <string.h>
00030 #include <time.h>
00031 #include "defs.h"
00032 #include "athena.h"
00033 #include "globals.h"
00034 #include "prototypes.h"
00035 #include "particles/particle.h"
00036 
00037 #ifndef SHEARING_BOX
00038 #error : The streaming2d problem requires shearing-box to be enabled.
00039 #endif /* SHEARING_BOX */
00040 
00041 #ifndef PARTICLES
00042 #error : The streaming2d problem requires particles to be enabled.
00043 #endif /* PARTICLES */
00044 
00045 #ifndef ISOTHERMAL
00046 #error : The streaming2d problem requires isothermal equation of state.
00047 #endif /* ISOTHERMAL */
00048 
00049 /*------------------------ filewide global variables -------------------------*/
00050 /* flow properties */
00051 Real vsc1,vsc2;
00052 int ipert;
00053 /* domain size variables */
00054 Real x1min,x1max,x2min,x2max,Lx,Lz,Lg;
00055 long Npar;
00056 /* output variables */
00057 long ntrack;   /* number of tracer particles */
00058 long nlis;     /* number of output particles for list output */
00059 int mytype;    /* specific particle type for output */
00060 Real dpar_thresh; /* threshold particle density */
00061 
00062 /*==============================================================================
00063  * PRIVATE FUNCTION PROTOTYPES:
00064  * ran2()            - random number generator
00065  * Normal()          - normal distribution generator
00066  * Erf()             - error function
00067  * MultiNSH()        - multiple component NSH equilibrium solver
00068  * ShearingBoxPot()  - shearing box tidal gravitational potential
00069  * hst_rho_Vx_dVy()  - total Reynolds stress for history dump
00070  * property_???()    - particle property selection function
00071  *============================================================================*/
00072 double ran2(long int *idum);
00073 double Normal(long int *idum);
00074 Real Erf(Real z);
00075 
00076 void MultiNSH(int n, Real *tstop, Real *mratio, Real etavk,
00077                      Real *uxNSH, Real *uyNSH,  Real *wxNSH, Real *wyNSH);
00078 static Real hst_rho_Vx_dVy(const Grid *pG,const int i,const int j,const int k);
00079 static void close_ix2(Grid *pGrid);
00080 static void close_ox2(Grid *pGrid);
00081 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3);
00082 
00083 static int property_limit(const Grain *gr, const GrainAux *grsub);
00084 static int property_trace(const Grain *gr, const GrainAux *grsub);
00085 static int property_type(const Grain *gr, const GrainAux *grsub);
00086 static int property_dense(const Grain *gr, const GrainAux *grsub);
00087 
00088 extern Real expr_dpar(const Grid *pG, const int i, const int j, const int k);
00089 extern Real expr_V1par(const Grid *pG, const int i, const int j, const int k);
00090 extern Real expr_V2par(const Grid *pG, const int i, const int j, const int k);
00091 extern Real expr_V3par(const Grid *pG, const int i, const int j, const int k);
00092 extern Real expr_V2(const Grid *pG, const int i, const int j, const int k);
00093 
00094 /*=========================== PUBLIC FUNCTIONS =================================
00095  *============================================================================*/
00096 /*----------------------------------------------------------------------------*/
00097 /* problem:   */
00098 
00099 void problem(Grid *pGrid, Domain *pDomain)
00100 {
00101   int i,j,js,pt,tsmode,vertical_bc;
00102   long p,q;
00103   Real ScaleHg,tsmin,tsmax,tscrit,amin,amax,Hparmin,Hparmax;
00104   Real *ep,*ScaleHpar,epsum,mratio,pwind,rhoaconv,etavk;
00105   Real *epsilon,*uxNSH,*uyNSH,**wxNSH,**wyNSH;
00106   Real rhog,h,x1,x2,x3,t,x1p,x2p,x3p,zmin,zmax,dx2_1,b;
00107   long int iseed = -1; /* Initialize on the first call to ran2 */
00108 
00109   if (par_geti("grid","Nx2") == 1) {
00110     ath_error("[par_strat2d]: par_strat2d only works for 2D problem.\n");
00111   }
00112   if (par_geti("grid","Nx3") > 1) {
00113     ath_error("[par_strat2d]: par_strat2d only works for 2D problem.\n");
00114   }
00115 #ifdef MPI_PARALLEL
00116   if (par_geti("parallel","NGrid_x2") > 2) {
00117     ath_error("[par_strat2d]: Domain decomposition in x2 (z) must be no more\
00118  than 2!\n");
00119   }
00120 #endif
00121 
00122 /* Initialize boxsize */
00123   x1min = pGrid->x1_0 + (pGrid->is + pGrid->idisp)*pGrid->dx1;
00124   x1max = pGrid->x1_0 + (pGrid->ie + pGrid->idisp + 1.0)*pGrid->dx1;
00125   Lx = x1max - x1min;
00126 
00127   x2min = par_getd("grid","x2min");
00128   x2max = par_getd("grid","x2max");
00129   Lz = x2max - x2min;
00130 
00131   Lg = nghost*pGrid->dx2; /* size of the ghost zone */
00132 
00133   js = pGrid->js;
00134 
00135 /* Read initial conditions */
00136   Omega_0 = par_getd("problem","omega");
00137   qshear = par_getd_def("problem","qshear",1.5);
00138   ipert = par_geti_def("problem","ipert",1);
00139   vsc1 = par_getd_def("problem","vsc1",0.05); /* in unit of iso_sound (N.B.!) */
00140   vsc2 = par_getd_def("problem","vsc2",0.0);
00141 
00142   vsc1 = vsc1 * Iso_csound;
00143   vsc2 = vsc2 * Iso_csound;
00144 
00145   ScaleHg = Iso_csound/Omega_0;
00146 
00147   /* particle number */
00148   Npar  = (long)(par_geti("particle","parnumgrid"));
00149 
00150   pGrid->nparticle = Npar*pGrid->partypes;
00151   for (i=0; i<pGrid->partypes; i++)
00152     pGrid->grproperty[i].num = Npar;
00153 
00154   if (pGrid->nparticle+2 > pGrid->arrsize)
00155     particle_realloc(pGrid, pGrid->nparticle+2);
00156 
00157   ep = (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00158   ScaleHpar = (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00159 
00160   epsilon = (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00161   wxNSH   = (Real**)calloc_2d_array(pGrid->Nx2+1, pGrid->partypes,sizeof(Real));
00162   wyNSH   = (Real**)calloc_2d_array(pGrid->Nx2+1, pGrid->partypes,sizeof(Real));
00163   uxNSH   = (Real*)calloc_1d_array(pGrid->Nx2+1, sizeof(Real));
00164   uyNSH   = (Real*)calloc_1d_array(pGrid->Nx2+1, sizeof(Real));
00165 
00166   /* particle stopping time */
00167   tsmode = par_geti("particle","tsmode");
00168   if (tsmode == 3) {/* fixed stopping time */
00169     tsmin = par_getd("problem","tsmin"); /* in code unit */
00170     tsmax = par_getd("problem","tsmax");
00171     tscrit= par_getd("problem","tscrit");
00172 
00173     for (i=0; i<pGrid->partypes; i++) {
00174       tstop0[i] = tsmin*exp(i*log(tsmax/tsmin)/MAX(pGrid->partypes-1,1.0));
00175       pGrid->grproperty[i].rad = tstop0[i];
00176       /* use fully implicit integrator for well coupled particles */
00177       if (tstop0[i] < tscrit) pGrid->grproperty[i].integrator = 3;
00178     }
00179   }
00180   else { 
00181     amin = par_getd("problem","amin");
00182     amax = par_getd("problem","amax");
00183 
00184     for (i=0; i<pGrid->partypes; i++)
00185       pGrid->grproperty[i].rad = amin*
00186                            exp(i*log(amax/amin)/MAX(pGrid->partypes-1,1.0));
00187 
00188     if (tsmode >= 2) {/* Epstein/General regime */
00189       /* conversion factor for rhoa */
00190       rhoaconv = par_getd_def("problem","rhoaconv",1.0);
00191       for (i=0; i<pGrid->partypes; i++)
00192         grrhoa[i]=pGrid->grproperty[i].rad*rhoaconv;
00193     }
00194 
00195     if (tsmode == 3)  /* General drag formula */
00196       alamcoeff = par_getd("problem","alamcoeff");
00197   }
00198 
00199   /* particle scale height */
00200   Hparmax = par_getd("problem","hparmax"); /* in unit of gas scale height */
00201   Hparmin = par_getd("problem","hparmin");
00202   for (i=0; i<pGrid->partypes; i++) 
00203     ScaleHpar[i] = Hparmax*
00204                    exp(-i*log(Hparmax/Hparmin)/MAX(pGrid->partypes-1,1.0));
00205 
00206 #ifdef FEEDBACK
00207   mratio = par_getd_def("problem","mratio",0.0); /* total mass fraction */
00208   pwind = par_getd_def("problem","pwind",0.0);   /* power law index */
00209   if (mratio < 0.0)
00210     ath_error("[par_strat2d]: mratio must be positive!\n");
00211 
00212   epsum = 0.0;
00213   for (i=0; i<pGrid->partypes; i++)
00214   {
00215     ep[i] = pow(pGrid->grproperty[i].rad,pwind);
00216     epsum += ep[i];
00217   }
00218 
00219   for (i=0; i<pGrid->partypes; i++)
00220   {
00221     ep[i] = mratio*ep[i]/epsum;
00222     pGrid->grproperty[i].m = sqrt(2.0*PI)*ScaleHg/Lz*ep[i]*
00223                                                   pGrid->Nx1*pGrid->Nx2/Npar;
00224   }
00225 #else
00226   mratio = 0.0;
00227   for (i=0; i<pGrid->partypes; i++)
00228     ep[i] = 0.0;
00229 #endif
00230 
00231   /* NSH equilibrium */
00232   for (j=pGrid->js; j<=pGrid->je+1; j++) {
00233 
00234     h = pGrid->x2_0 + (j+pGrid->jdisp)*pGrid->dx2;
00235     q = j - js;
00236     etavk = fabs(vsc1+vsc2*SQR(h));
00237 
00238     for (i=0; i<pGrid->partypes; i++) {
00239       epsilon[i] = ep[i]/ScaleHpar[i]/erf(Lz/(sqrt(8.0)*ScaleHpar[i]*ScaleHg))*
00240                    exp(-0.5*SQR(h/ScaleHg)*(SQR(1.0/ScaleHpar[i])-1.0));
00241 
00242       if (tsmode != 3)
00243         tstop0[i] = get_ts(pGrid,i,exp(-0.5*SQR(h/ScaleHg)),Iso_csound,etavk);
00244     }
00245 
00246     MultiNSH(pGrid->partypes, tstop0, epsilon, etavk,
00247                               &uxNSH[q], &uyNSH[q], wxNSH[q], wyNSH[q]);
00248   }
00249 
00250 /* Now set initial conditions for the gas */
00251   for (j=pGrid->js; j<=pGrid->je; j++) {
00252   for (i=pGrid->is; i<=pGrid->ie; i++) {
00253     cc_pos(pGrid,i,j,pGrid->ks,&x1,&x2,&x3);
00254 
00255     rhog = exp(-0.5*SQR(x2/ScaleHg));
00256     pGrid->U[pGrid->ks][j][i].d = rhog;
00257 
00258     if (ipert != 1) {/* NSH velocity */
00259       pGrid->U[pGrid->ks][j][i].M1 = 0.5*rhog*(uxNSH[j-js]+uxNSH[j-js+1]);
00260       pGrid->U[pGrid->ks][j][i].M3 = 0.5*rhog*(uyNSH[j-js]+uyNSH[j-js+1]);
00261     } else {
00262       pGrid->U[pGrid->ks][j][i].M1 = 0.0;
00263       pGrid->U[pGrid->ks][j][i].M3 = 0.0;
00264     }
00265 
00266     pGrid->U[pGrid->ks][j][i].M2 = 0.0;
00267 #ifndef FARGO
00268     pGrid->U[pGrid->ks][j][i].M3 -= qshear*rhog*Omega*x1;
00269 #endif
00270 
00271   }}
00272 
00273 /* Now set initial conditions for the particles */
00274   p = 0;
00275   dx2_1 = 1.0/pGrid->dx2;
00276   x3p = pGrid->x3_0 + (pGrid->ks+pGrid->kdisp)*pGrid->dx3;
00277   zmin = pGrid->x2_0 + (pGrid->js+pGrid->jdisp)*pGrid->dx2;
00278   zmax = pGrid->x2_0 + (pGrid->je+pGrid->jdisp+1.0)*pGrid->dx2;
00279 
00280   for (q=0; q<Npar; q++) {
00281 
00282     for (pt=0; pt<pGrid->partypes; pt++) {
00283 
00284       x1p = x1min + Lx*ran2(&iseed);
00285       x2p = ScaleHpar[pt]*ScaleHg*Normal(&iseed);
00286       while ((x2p >= zmax) || (x2p < zmin))
00287         x2p = ScaleHpar[pt]*ScaleHg*Normal(&iseed);
00288 
00289       pGrid->particle[p].property = pt;
00290       pGrid->particle[p].x1 = x1p;
00291       pGrid->particle[p].x2 = x2p;
00292       pGrid->particle[p].x3 = x3p;
00293 
00294       if (ipert != 1) {/* NSH velocity */
00295 
00296         cellj(pGrid, x2p, dx2_1, &j, &b);
00297         j = j-pGrid->js;  b = b - pGrid->js;
00298 
00299         pGrid->particle[p].v1 = (j+1-b)*wxNSH[j][pt]+(b-j)*wxNSH[j+1][pt];
00300         pGrid->particle[p].v3 = (j+1-b)*wyNSH[j][pt]+(b-j)*wyNSH[j+1][pt];
00301 
00302       } else {
00303 
00304         pGrid->particle[p].v1 = 0.0;
00305         pGrid->particle[p].v3 = vsc1+vsc2*SQR(x2p);
00306 
00307       }
00308 
00309       pGrid->particle[p].v2 = 0.0;
00310 #ifndef FARGO
00311       pGrid->particle[p].v3 -= qshear*Omega*x1p;
00312 #endif
00313 
00314       pGrid->particle[p].pos = 1; /* grid particle */
00315       pGrid->particle[p].my_id = p;
00316 #ifdef MPI_PARALLEL
00317       pGrid->particle[p].init_id = pGrid->my_id;
00318 #endif
00319       p++;
00320   }}
00321 
00322 /* enroll gravitational potential function, shearing sheet BC functions */
00323   StaticGravPot = ShearingBoxPot;
00324 
00325   dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00326 
00327   /* set the # of the particles in list output
00328    * (by default, output 1 particle per cell)
00329    */
00330   nlis = par_geti_def("problem","nlis",pGrid->Nx1*pGrid->Nx2);
00331 
00332   /* set the number of particles to keep track of */
00333   ntrack = par_geti_def("problem","ntrack",2000);
00334 
00335   /* set the threshold particle density */
00336   dpar_thresh = par_geti_def("problem","dpar_thresh",10.0);
00337 
00338   /* set vertical boundary conditions (by default, periodic) */
00339   vertical_bc = par_geti_def("problem","vertical_bc",0);
00340 
00341   if (vertical_bc == 1) /* closed BC */
00342   {
00343     set_bvals_mhd_fun(left_x2,  close_ix2);
00344     set_bvals_mhd_fun(right_x2, close_ox2);
00345   }
00346 
00347   /* finalize */
00348   free(ep);  free(ScaleHpar);
00349   free(epsilon);
00350   free_2d_array(wxNSH);  free_2d_array(wyNSH);
00351   free(uxNSH);           free(uyNSH);
00352 
00353   return;
00354 }
00355 
00356 /*==============================================================================
00357  * PROBLEM USER FUNCTIONS:
00358  * problem_write_restart() - writes problem-specific user data to restart files
00359  * problem_read_restart()  - reads problem-specific user data from restart files
00360  * get_usr_expr()          - sets pointer to expression for special output data
00361  * get_usr_out_fun()       - returns a user defined output function pointer
00362  * Userwork_in_loop        - problem specific work IN     main loop
00363  * Userwork_after_loop     - problem specific work AFTER  main loop
00364  *----------------------------------------------------------------------------*/
00365 
00366 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00367 {
00368   return;
00369 }
00370 
00371 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00372 {
00373 
00374   StaticGravPot = ShearingBoxPot;
00375 
00376   Omega_0 = par_getd("problem","omega");
00377   qshear = par_getd_def("problem","qshear",1.5);
00378   ipert = par_geti_def("problem","ipert",1);
00379 
00380   x1min = pG->x1_0 + (pG->is + pG->idisp)*pG->dx1;
00381   x1max = pG->x1_0 + (pG->ie + pG->idisp + 1.0)*pG->dx1;
00382   Lx = x1max - x1min;
00383 
00384   x2min = par_getd("grid","x2min");
00385   x2max = par_getd("grid","x2max");
00386   Lz = x2max - x2min;
00387 
00388   Lg = nghost*pG->dx2; /* size of the ghost zone */
00389 
00390   vsc1 = par_getd_def("problem","vsc1",0.05); /* in unit of iso_sound (N.B.!) */
00391   vsc2 = par_getd_def("problem","vsc2",0.0);
00392 
00393   vsc1 = vsc1 * Iso_csound;
00394   vsc2 = vsc2 * Iso_csound;
00395 
00396   Npar  = (int)(sqrt(par_geti("particle","parnumgrid")));
00397   nlis = par_geti_def("problem","nlis",pG->Nx1*pG->Nx2);
00398   ntrack = par_geti_def("problem","ntrack",2000);
00399 
00400   dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00401 
00402   return;
00403 }
00404 
00405 static Real hst_rho_Vx_dVy(const Grid *pG,
00406                            const int i, const int j, const int k)
00407 {
00408   Real x1,x2,x3;
00409   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00410 #ifdef FARGO
00411   return pG->U[k][j][i].M1*pG->U[k][j][i].M3/pG->U[k][j][i].d;
00412 #else
00413   return pG->U[k][j][i].M1*(pG->U[k][j][i].M3/pG->U[k][j][i].d
00414                              + qshear*Omega_0*x1);
00415 #endif
00416 }
00417 
00418 static void close_ix2(Grid *pGrid)
00419 {
00420   int js = pGrid->js;
00421   int ks = pGrid->ks, ke = pGrid->ke;
00422   int i,j,k,il,iu; /* i-lower/upper */
00423 
00424   iu = pGrid->ie + nghost;
00425   il = pGrid->is - nghost;
00426 
00427   for (k=ks; k<=ke; k++)
00428     for (j=1; j<=nghost; j++)
00429       for (i=il; i<=iu; i++) {
00430         pGrid->U[k][js-j][i] = pGrid->U[k][js][i];
00431         pGrid->U[k][js-j][i].M2 = 0.0;
00432       }
00433 
00434   return;
00435 }
00436 
00437 static void close_ox2(Grid *pGrid)
00438 {
00439   int je = pGrid->je;
00440   int ks = pGrid->ks, ke = pGrid->ke;
00441   int i,j,k,il,iu; /* i-lower/upper */
00442 
00443   iu = pGrid->ie + nghost;
00444   il = pGrid->is - nghost;
00445 
00446   for (k=ks; k<=ke; k++)
00447     for (j=1; j<=nghost; j++)
00448       for (i=il; i<=iu; i++) {
00449         pGrid->U[k][je+j][i] = pGrid->U[k][je][i];
00450         pGrid->U[k][je+j][i].M2 = 0.0;
00451       }
00452 
00453   return;
00454 }
00455 
00456 Gasfun_t get_usr_expr(const char *expr)
00457 {
00458   return NULL;
00459 }
00460 
00461 VGFunout_t get_usr_out_fun(const char *name){
00462   return NULL;
00463 }
00464 
00465 #ifdef PARTICLES
00466 PropFun_t get_usr_par_prop(const char *name)
00467 {
00468   if (strcmp(name,"limit")==0)    return property_limit;
00469   if (strcmp(name,"trace")==0)    return property_trace;
00470   if (strcmp(name,"type")==0)    return property_type;
00471 
00472   return NULL;
00473 }
00474 
00475 void gasvshift(const Real x1, const Real x2, const Real x3,
00476                                     Real *u1, Real *u2, Real *u3)
00477 {
00478   return;
00479 }
00480 
00481 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00482                                     const Real v1, const Real v2, const Real v3)
00483 {
00484   Real z,fac;
00485 
00486   ft->x1 -= 2.0*(vsc1 + vsc2*SQR(x2))*Omega_0;
00487 
00488   if(x2 > x2max)
00489     z = x2-Lz;
00490   else if (x2 < x2min)
00491     z = x2+Lz;
00492   else
00493     z = x2;
00494 
00495   fac = Lg/(0.5*Lz+Lg-fabs(z));
00496   ft->x2 -= SQR(Omega_0)*z*(1.0-SQR(fac)*fac); /* 3rd order sharp */
00497 //  *ft.x2 -= SQR(Omega_0)*z*(1.0-SQR(fac));  /* 2nd order sharp */
00498 
00499   return;
00500 }
00501 #endif
00502 
00503 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00504 {
00505   return;
00506 }
00507 
00508 /*---------------------------------------------------------------------------
00509  * Userwork_after_loop: computes L1-error in linear waves,
00510  * ASSUMING WAVE HAS PROPAGATED AN INTEGER NUMBER OF PERIODS
00511  * Must set parameters in input file appropriately so that this is true
00512  */
00513 
00514 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00515 {
00516   return;
00517 }
00518  
00519 /*=========================== PRIVATE FUNCTIONS ==============================*/
00520 /*--------------------------------------------------------------------------- */
00521 /*! \fn static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00522  *  \brief Shearing box tidal gravitational potential */
00523 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00524 {
00525   Real z,z0,phi=0.0;
00526 
00527 #ifndef FARGO
00528   phi -= qshear*SQR(Omega_0*x1);
00529 #endif
00530   /* Ensure vertical periodicity in ghost zones */
00531   if(x2 > x2max)
00532     z=fabs(x2-Lz);
00533   else if (x2 < x2min)
00534     z=fabs(x2+Lz);
00535   else
00536     z=fabs(x2);
00537 
00538   phi += 0.5*SQR(Omega_0*z);
00539 
00540   /* smooth the potential at domain edges */
00541   z0 = 0.5*Lz+Lg;
00542   phi -= SQR(Omega_0*Lg)*Lg*(2.0*z-z0)/(2.0*SQR(z0-z)); /* 3rd order sharp */
00543 
00544 //  phi -= SQR(Omega_0*Lg)*(z/(z0-z)+log(z0/(z0-z))); /* 2nd order sharp */
00545 
00546   return phi;
00547 }
00548 
00549 /*! \fn static int property_limit(const Grain *gr, const GrainAux *grsub)
00550  *  \brief User defined particle selection function (1: true; 0: false) */
00551 static int property_limit(const Grain *gr, const GrainAux *grsub)
00552 {
00553   if ((gr->pos == 1) && (gr->my_id<nlis))
00554     return 1;
00555   else
00556     return 0;
00557 }
00558 
00559 /*! \fn static int property_trace(const Grain *gr, const GrainAux *grsub)
00560  *  \brief User defined particle selection function (1: true; 0: false) */
00561 static int property_trace(const Grain *gr, const GrainAux *grsub)
00562 {
00563   if ((gr->pos == 1) && (gr->my_id<ntrack))
00564     return 1;
00565   else
00566     return 0;
00567 }
00568 
00569 /*! \fn static int property_type(const Grain *gr, const GrainAux *grsub)
00570  *  \brief User defined particle selection function (1: true; 0: false) */
00571 static int property_type(const Grain *gr, const GrainAux *grsub)
00572 {
00573   if ((gr->pos == 1) && (gr->property == mytype))
00574     return 1;
00575   else
00576     return 0;
00577 }
00578 
00579 /*! \fn static int property_dense(const Grain *gr, const GrainAux *grsub)
00580  *  \brief User defined particle selection function (1: true; 0: false) */
00581 static int property_dense(const Grain *gr, const GrainAux *grsub)
00582 {
00583   if ((gr->pos == 1) && (grsub->dpar > dpar_thresh))
00584     return 1;
00585   else
00586     return 0;
00587 }
00588 
00589 /*----------------------------------------------------------------------------*/
00590 /*! \fn void MultiNSH(int n, Real *tstop, Real *mratio, Real etavk,
00591  *                   Real *uxNSH, Real *uyNSH, Real *wxNSH, Real *wyNSH)
00592  *  \brief Multi-species NSH equilibrium
00593  *
00594  * Input: # of particle types (n), dust stopping time and mass ratio array, and
00595  *        drift speed etavk.
00596  * Output: gas NSH equlibrium velocity (u), and dust NSH equilibrium velocity
00597  *         array (w).
00598  */
00599 void MultiNSH(int n, Real *tstop, Real *mratio, Real etavk,
00600                      Real *uxNSH, Real *uyNSH, Real *wxNSH, Real *wyNSH)
00601 {
00602   int i,j;
00603   Real *Lambda1,**Lam1GamP1, **A, **B, **Tmp;
00604 
00605   Lambda1 = (Real*)calloc_1d_array(n, sizeof(Real));     /* Lambda^{-1} */
00606   Lam1GamP1=(Real**)calloc_2d_array(n, n, sizeof(Real)); /* Lambda1*(1+Gamma) */
00607   A       = (Real**)calloc_2d_array(n, n, sizeof(Real));
00608   B       = (Real**)calloc_2d_array(n, n, sizeof(Real));
00609   Tmp     = (Real**)calloc_2d_array(n, n, sizeof(Real));
00610 
00611   /* definitions */
00612   for (i=0; i<n; i++){
00613     for (j=0; j<n; j++)
00614       Lam1GamP1[i][j] = mratio[j];
00615     Lam1GamP1[i][i] += 1.0;
00616     Lambda1[i] = 1.0/(tstop[i]+1.0e-16);
00617     for (j=0; j<n; j++)
00618       Lam1GamP1[i][j] *= Lambda1[i];
00619   }
00620 
00621   /* Calculate A and B */
00622   MatrixMult(Lam1GamP1, Lam1GamP1, n,n,n, Tmp);
00623   for (i=0; i<n; i++) Tmp[i][i] += 1.0;
00624   InverseMatrix(Tmp, n, B);
00625   for (i=0; i<n; i++)
00626   for (j=0; j<n; j++)
00627     B[i][j] *= Lambda1[j];
00628   MatrixMult(Lam1GamP1, B, n,n,n, A);
00629 
00630   /* Obtain NSH velocities */
00631   *uxNSH = 0.0;  *uyNSH = 0.0;
00632   for (i=0; i<n; i++){
00633     wxNSH[i] = 0.0;
00634     wyNSH[i] = 0.0;
00635     for (j=0; j<n; j++){
00636       wxNSH[i] -= B[i][j];
00637       wyNSH[i] -= A[i][j];
00638     }
00639     wxNSH[i] *= 2.0*etavk;
00640     wyNSH[i] *= etavk;
00641     *uxNSH -= mratio[i]*wxNSH[i];
00642     *uyNSH -= mratio[i]*wyNSH[i];
00643     wyNSH[i] += etavk;
00644   }
00645 
00646   free(Lambda1);
00647   free_2d_array(A);         free_2d_array(B);
00648   free_2d_array(Lam1GamP1); free_2d_array(Tmp);
00649 
00650   return;
00651 }
00652 
00653 /*------------------------------------------------------------------------------
00654  */
00655 
00656 #define IM1 2147483563
00657 #define IM2 2147483399
00658 #define AM (1.0/IM1)
00659 #define IMM1 (IM1-1)
00660 #define IA1 40014
00661 #define IA2 40692
00662 #define IQ1 53668
00663 #define IQ2 52774
00664 #define IR1 12211
00665 #define IR2 3791
00666 #define NTAB 32
00667 #define NDIV (1+IMM1/NTAB)
00668 #define RNMX (1.0-DBL_EPSILON)
00669 
00670 /*! \fn double ran2(long int *idum)
00671  *  \brief Extracted from the Numerical Recipes in C (version 2) code.  Modified
00672  *   to use doubles instead of floats. -- T. A. Gardiner -- Aug. 12, 2003
00673  *
00674  * Long period (> 2 x 10^{18}) random number generator of L'Ecuyer
00675  * with Bays-Durham shuffle and added safeguards.  Returns a uniform
00676  * random deviate between 0.0 and 1.0 (exclusive of the endpoint
00677  * values).  Call with idum = a negative integer to initialize;
00678  * thereafter, do not alter idum between successive deviates in a
00679  * sequence.  RNMX should appriximate the largest floating point value
00680  * that is less than 1.
00681  */
00682 
00683 double ran2(long int *idum)
00684 {
00685   int j;
00686   long int k;
00687   static long int idum2=123456789;
00688   static long int iy=0;
00689   static long int iv[NTAB];
00690   double temp;
00691 
00692   if (*idum <= 0) { /* Initialize */
00693     if (-(*idum) < 1) *idum=1; /* Be sure to prevent idum = 0 */
00694     else *idum = -(*idum);
00695     idum2=(*idum);
00696     for (j=NTAB+7;j>=0;j--) { /* Load the shuffle table (after 8 warm-ups) */
00697       k=(*idum)/IQ1;
00698       *idum=IA1*(*idum-k*IQ1)-k*IR1;
00699       if (*idum < 0) *idum += IM1;
00700       if (j < NTAB) iv[j] = *idum;
00701     }
00702     iy=iv[0];
00703   }
00704   k=(*idum)/IQ1;                 /* Start here when not initializing */
00705   *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without */
00706   if (*idum < 0) *idum += IM1;   /* overflows by Schrage's method */
00707   k=idum2/IQ2;
00708   idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise */
00709   if (idum2 < 0) idum2 += IM2;
00710   j=(int)(iy/NDIV);              /* Will be in the range 0...NTAB-1 */
00711   iy=iv[j]-idum2;                /* Here idum is shuffled, idum and idum2 */
00712   iv[j] = *idum;                 /* are combined to generate output */
00713   if (iy < 1) iy += IMM1;
00714   if ((temp=AM*iy) > RNMX) return RNMX; /* No endpoint values */
00715   else return temp;
00716 }
00717 
00718 /*--------------------------------------------------------------------------- */
00719 /*! \fn double Normal(long int *idum)
00720  *  \brief Normal distribution random number generator */
00721 double Normal(long int *idum)
00722 {
00723   double Y,X1,X2;
00724 
00725   X1 = ran2(idum);
00726   X2 = ran2(idum);
00727 
00728   Y = sqrt(-2.0*log(X1+TINY_NUMBER))*cos(2*PI*X2);
00729 
00730   return Y;
00731 }
00732 
00733 /*--------------------------------------------------------------------------- */
00734 /*! \fn Real Erf(Real z)
00735  *  \brief Error function  */
00736 Real Erf(Real z)
00737 {
00738   /* arrays of the error function y=erf(x) */
00739   static double x[101]={
00740         0.000000e+000,  3.783387e-003,  7.709914e-003,  1.178500e-002,  1.601425e-002,  2.040352e-002,
00741         2.495885e-002,  2.968653e-002,  3.459307e-002,  3.968525e-002,  4.497008e-002,  5.045486e-002,
00742         5.614715e-002,  6.205480e-002,  6.818596e-002,  7.454909e-002,  8.115295e-002,  8.800667e-002,
00743         9.511969e-002,  1.025018e-001,  1.101632e-001,  1.181145e-001,  1.263667e-001,  1.349310e-001,
00744         1.438193e-001,  1.530440e-001,  1.626176e-001,  1.725534e-001,  1.828652e-001,  1.935671e-001,
00745         2.046738e-001,  2.162008e-001,  2.281639e-001,  2.405796e-001,  2.534651e-001,  2.668380e-001,
00746         2.807169e-001,  2.951209e-001,  3.100699e-001,  3.255844e-001,  3.416859e-001,  3.583966e-001,
00747         3.757395e-001,  3.937386e-001,  4.124186e-001,  4.318054e-001,  4.519256e-001,  4.728071e-001,
00748         4.944786e-001,  5.169701e-001,  5.403124e-001,  5.645379e-001,  5.896800e-001,  6.157732e-001,
00749         6.428537e-001,  6.709587e-001,  7.001271e-001,  7.303990e-001,  7.618162e-001,  7.944220e-001,
00750         8.282614e-001,  8.633812e-001,  8.998296e-001,  9.376570e-001,  9.769156e-001,  1.017659e+000,
00751         1.059945e+000,  1.103830e+000,  1.149376e+000,  1.196644e+000,  1.245701e+000,  1.296614e+000,
00752         1.349454e+000,  1.404292e+000,  1.461205e+000,  1.520272e+000,  1.581573e+000,  1.645193e+000,
00753         1.711221e+000,  1.779746e+000,  1.850864e+000,  1.924673e+000,  2.001274e+000,  2.080774e+000,
00754         2.163281e+000,  2.248909e+000,  2.337778e+000,  2.430008e+000,  2.525728e+000,  2.625070e+000,
00755         2.728170e+000,  2.835170e+000,  2.946219e+000,  3.061469e+000,  3.181080e+000,  3.305216e+000,
00756         3.434048e+000,  3.567755e+000,  3.706521e+000,  3.850536e+000,  4.000000e+000
00757   };
00758   static double y[101]={
00759         0.00000000e+000,        4.26907434e-003,        8.69953340e-003,        1.32973284e-002,        1.80686067e-002,        2.30197153e-002,
00760         2.81572033e-002,        3.34878242e-002,        3.90185379e-002,        4.47565113e-002,        5.07091186e-002,        5.68839404e-002,
00761         6.32887618e-002,        6.99315688e-002,        7.68205444e-002,        8.39640613e-002,        9.13706742e-002,        9.90491090e-002,
00762         1.07008250e-001,        1.15257124e-001,        1.23804883e-001,        1.32660778e-001,        1.41834139e-001,        1.51334337e-001,
00763         1.61170754e-001,        1.71352743e-001,        1.81889576e-001,        1.92790394e-001,        2.04064148e-001,        2.15719527e-001,
00764         2.27764884e-001,        2.40208149e-001,        2.53056730e-001,        2.66317410e-001,        2.79996226e-001,        2.94098338e-001,
00765         3.08627885e-001,        3.23587825e-001,        3.38979770e-001,        3.54803790e-001,        3.71058224e-001,        3.87739454e-001,
00766         4.04841688e-001,        4.22356710e-001,        4.40273635e-001,        4.58578645e-001,        4.77254725e-001,        4.96281391e-001,
00767         5.15634428e-001,        5.35285634e-001,        5.55202571e-001,        5.75348359e-001,        5.95681482e-001,        6.16155658e-001,
00768         6.36719759e-001,        6.57317799e-001,        6.77889021e-001,        6.98368078e-001,        7.18685336e-001,        7.38767318e-001,
00769         7.58537287e-001,        7.77916009e-001,        7.96822665e-001,        8.15175962e-001,        8.32895397e-001,        8.49902691e-001,
00770         8.66123358e-001,        8.81488386e-001,        8.95935967e-001,        9.09413237e-001,        9.21877939e-001,        9.33299942e-001,
00771         9.43662512e-001,        9.52963249e-001,        9.61214608e-001,        9.68443923e-001,        9.74692870e-001,        9.80016358e-001,
00772         9.84480847e-001,        9.88162149e-001,        9.91142807e-001,        9.93509200e-001,        9.95348535e-001,        9.96745927e-001,
00773         9.97781755e-001,        9.98529475e-001,        9.99054014e-001,        9.99410828e-001,        9.99645625e-001,        9.99794704e-001,
00774         9.99885782e-001,        9.99939162e-001,        9.99969080e-001,        9.99985060e-001,        9.99993164e-001,        9.99997050e-001,
00775         9.99998805e-001,        9.99999548e-001,        9.99999841e-001,        9.99999948e-001,        9.99999985e-001
00776   };
00777   double z1, g;
00778   int il, iu, im;
00779   z1 = fabs(z);
00780   /* look for the location of z1 in the x array */
00781   il=0; iu=100;
00782   while(iu-il>1) {
00783     im = (iu+il)/2;
00784     if (x[im] > z1) iu = im;
00785     else il = im;
00786   }
00787   /* linear interpolation */
00788   g = (y[iu]*(z1-x[il])+y[il]*(x[iu]-z1))/(x[iu]-x[il]);
00789 
00790   g = MIN(g,1.0);
00791   if (z<0.0) g = -g;
00792   return g;
00793 }
00794 
00795 #undef IM1
00796 #undef IM2
00797 #undef AM
00798 #undef IMM1
00799 #undef IA1
00800 #undef IA2
00801 #undef IQ1
00802 #undef IQ2
00803 #undef IR1
00804 #undef IR2
00805 #undef NTAB
00806 #undef NDIV
00807 #undef RNMX