Athena: prob/par_strat3d.c Source File

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