Athena: prob/par_shwave1d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file par_shwave1d.c
00004  *  \brief Problem generator radial epicyclic wave test with Lagrangian
00005  *   particles.
00006  *
00007  * PURPOSE: Problem generator radial epicyclic wave test with Lagrangian
00008  *   particles.
00009  *
00010  *   Only works for 2D, and the wavevector is in x1 direction. Particles can be
00011  *   initialized either with uniform density or with density profile same as the
00012  *   gas. SHEARING_BOX must be turned on, FARGO is optional.
00013  *   ipert: 1 for linear wave; 2 for non-linear wave from Fromang & Papaloizou.
00014  *
00015  * Reference: Fromang & Papaloizou, A&A, 468, 1-18 (2007).
00016  */
00017 /*============================================================================*/
00018 
00019 #include <math.h>
00020 #include <stdio.h>
00021 #include <stdlib.h>
00022 #include "defs.h"
00023 #include "athena.h"
00024 #include "globals.h"
00025 #include "prototypes.h"
00026 #include "particles/particle.h"
00027 
00028 #ifndef SHEARING_BOX
00029 #error : The shear wave problem requires shearing-box to be enabled.
00030 #endif /* SHEARING_BOX */
00031 
00032 
00033 /*==============================================================================
00034  * PRIVATE FUNCTION PROTOTYPES:
00035  * GetPosition()    - get particle status (grid/crossing)
00036  * ShearingBoxPot() - shearing box tidal potential
00037  *============================================================================*/
00038 #ifdef PARTICLES
00039 int GetPosition(Grain *gr);
00040 #endif
00041 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3);
00042 
00043 /*=========================== PUBLIC FUNCTIONS =================================
00044  *============================================================================*/
00045 /*----------------------------------------------------------------------------*/
00046 /* problem:   */
00047 
00048 void problem(Grid *pGrid, Domain *pDomain)
00049 {
00050   FILE *fp;
00051   Real xFP[160],dFP[160],vxFP[160],vyFP[160];
00052   int i=0,j=0,k=0,ipert;
00053   int is,ie,js,je,ks,ke,n,m,nwave,samp;
00054   Real x1,x2,x3,x1max,x1min,x2max,x2min;
00055   Real kx,omg,omg2,amp,v1x,v1y;
00056 #ifdef PARTICLES
00057   long p;
00058   int Npar,ip,jp;
00059   Real x1p,x2p,x3p,x1l,x1u,x2l,x2u,par_amp,factor2;
00060 #endif
00061 
00062   if ((par_geti("grid","Nx2") == 1) || (par_geti("grid","Nx3") > 1)) {
00063     ath_error("[par_shwave1d]: par_shwave1d must work in 2D grid.\n");
00064   }
00065 
00066   is = pGrid->is; ie = pGrid->ie;
00067   js = pGrid->js; je = pGrid->je;
00068   ks = pGrid->ks; ke = pGrid->ke;
00069 
00070 /* Read initial conditions  */
00071   amp = par_getd("problem","amp");
00072   Omega_0 = par_getd_def("problem","omega",1.0);
00073   qshear  = par_getd_def("problem","qshear",1.5);
00074   ipert = par_geti_def("problem","ipert", 1);
00075   nwave = par_geti("problem","nwave");
00076   samp = par_geti("problem","sample");
00077   x1min = par_getd("grid","x1min");
00078   x1max = par_getd("grid","x1max");
00079   x2min = par_getd("grid","x2min");
00080   x2max = par_getd("grid","x2max");
00081   if (ipert != 1) samp = 2;
00082 
00083 #ifdef PARTICLES
00084 /* Read initial conditions for the particles */
00085   /* basic parameters */
00086   if (par_geti("particle","partypes") != 1)
00087     ath_error("[par_shwave1d]: This test only allows ONE particle species!\n");
00088 
00089   Npar = (int)(sqrt(par_geti("particle","parnumcell")));
00090   pGrid->nparticle = Npar*Npar*pGrid->Nx1*pGrid->Nx2;
00091   pGrid->grproperty[0].num = pGrid->nparticle;
00092   if (pGrid->nparticle+2 > pGrid->arrsize)
00093     particle_realloc(pGrid, pGrid->nparticle+2);
00094 
00095   /* particle stopping time */
00096   tstop0[0] = par_getd_def("problem","tstop",0.0); /* in code unit */
00097   if (par_geti("particle","tsmode") != 3)
00098     ath_error("[par_shwave1d]: This test only allows fixed stopping time!\n");
00099 
00100   par_amp = amp;
00101   factor2 = 0.5;
00102 #endif /* PARTICLES */
00103 
00104 /* Initialize gas and particles */
00105   if (ipert == 1)
00106   {
00107     if (nwave <= 0)
00108       ath_error("[par_shwave1d]: nwave must be positive!\n");
00109 
00110   /* calculate dispersion relation and find eigen vectors */
00111     kx = 2.0*(PI)*nwave/(x1max-x1min);
00112     omg2 = SQR(Iso_csound*kx)+2.0*(2.0-qshear)*SQR(Omega_0);
00113     omg = sqrt(omg2);
00114     v1x = omg*amp/kx;
00115     v1y = (2.0-qshear)*Omega_0/omg*v1x;
00116 #ifdef PARTICLES
00117     /* particle perturbation amplitude */
00118 //    par_amp = amp*kx*pGrid->dx1/sin(kx*pGrid->dx1);
00119 //    factor2 = 0.5*tan(kx*pGrid->dx1)/(kx*pGrid->dx1);
00120 #endif
00121 
00122   /* Now set initial conditions to wave solution */ 
00123     for (k=ks; k<=ke; k++) {
00124     for (j=js; j<=je; j++) {
00125     for (i=is; i<=ie; i++) {
00126       cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00127       pGrid->U[k][j][i].d = 1.0+amp*cos(kx*x1);
00128       pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d*v1x*cos(kx*x1);
00129       pGrid->U[k][j][i].M3 = pGrid->U[k][j][i].d*v1y*sin(kx*x1);
00130 #ifndef FARGO
00131       pGrid->U[k][j][i].M3 -= pGrid->U[k][j][i].d*qshear*Omega_0*x1;
00132 #endif
00133       pGrid->U[k][j][i].M2 = 0.0;
00134 #if (NSCALARS > 0)
00135       if (samp == 1)
00136         for (n=0; n<NSCALARS; n++)
00137           pGrid->U[k][j][i].s[n] = pGrid->U[k][j][i].d;
00138       else
00139         for (n=0; n<NSCALARS; n++)
00140           pGrid->U[k][j][i].s[n] = 1.0;
00141 #endif
00142     }}}
00143   }
00144   else /* ipert != 1 (read FP solution */
00145   {
00146     if (pGrid->Nx1 == 160) {
00147       if((fp = fopen("Data-160-FPwave.dat","r")) == NULL)
00148          ath_error("Error opening Data-160-FPwave.dat\n");
00149       for (i=0; i<160; i++) {
00150         fscanf(fp,"%lf %lf %lf %lf",&xFP[i],&dFP[i],&vxFP[i],&vyFP[i]);
00151       }
00152     }
00153 
00154     if (pGrid->Nx1 == 40) {
00155       if((fp = fopen("Data-40-FPwave.dat","r")) == NULL)
00156          ath_error("Error opening Data-40-FPwave.dat\n");
00157       for (i=0; i<40; i++) {
00158         fscanf(fp,"%lf %lf %lf %lf",&xFP[i],&dFP[i],&vxFP[i],&vyFP[i]);
00159       }
00160     }
00161 
00162     if (x1min != -4.7965)
00163       ath_error("[par_shwave1d]: ipert=2 requires x1min=-4.7965\n");
00164     if (x1max != 4.7965)
00165       ath_error("[par_shwave1d]: ipert=2 requires x1max=4.7965\n");
00166 
00167   /* Now set initial conditions to wave solution */
00168 
00169     for (k=ks; k<=ke; k++) {
00170     for (j=js; j<=je; j++) {
00171     for (i=is; i<=ie; i++) {
00172       cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00173       pGrid->U[k][j][i].d = dFP[i+pGrid->idisp];
00174       pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d*vxFP[i+pGrid->idisp];
00175       pGrid->U[k][j][i].M3 = pGrid->U[k][j][i].d*vyFP[i+pGrid->idisp];
00176 #ifdef FARGO
00177       pGrid->U[k][j][i].M3 += pGrid->U[k][j][i].d*qshear*Omega_0*x1;
00178 #endif
00179       pGrid->U[k][j][i].M2 = 0.0;
00180 #if (NSCALARS > 0)
00181       for (n=0; n<NSCALARS; n++)
00182         pGrid->U[k][j][i].s[n] = 1.0;
00183 #endif
00184     }}}
00185   }
00186 
00187 #ifdef PARTICLES
00188   /* Now set initial conditions for the particles */
00189   p = 0;
00190   x3p = pGrid->x3_0 + (pGrid->ks+pGrid->kdisp)*pGrid->dx3;
00191 
00192   for (j=pGrid->js; j<=pGrid->je; j++)
00193   {
00194     x2l = pGrid->x2_0 + (j+pGrid->jdisp)*pGrid->dx2;
00195     x2u = pGrid->x2_0 + ((j+pGrid->jdisp)+1.0)*pGrid->dx2;
00196 
00197     for (i=pGrid->is; i<=pGrid->ie; i++)
00198     {
00199       x1l = pGrid->x1_0 + (i + pGrid->idisp)*pGrid->dx1;
00200       x1u = pGrid->x1_0 + ((i + pGrid->idisp) + 1.0)*pGrid->dx1;
00201 
00202       for (ip=0;ip<Npar;ip++)
00203       {
00204         x1p = x1l+(x1u-x1l)/Npar*(ip+0.5);
00205 
00206         for (jp=0;jp<Npar;jp++)
00207         {
00208           x2p = x2l+(x2u-x2l)/Npar*(jp+0.5);
00209 
00210           pGrid->particle[p].property = 0;
00211 
00212           pGrid->particle[p].x1 = x1p;
00213           pGrid->particle[p].x2 = x2p;
00214           if (samp == 1)
00215             pGrid->particle[p].x1 += -par_amp*sin(kx*x1p)/kx
00216                                      +factor2*SQR(amp)*sin(2.0*kx*x1p)/kx;
00217           pGrid->particle[p].x3 = x3p;
00218 
00219           pGrid->particle[p].v1 = v1x*cos(kx*x1p);
00220           pGrid->particle[p].v3 = v1y*sin(kx*x1p);
00221 #ifndef FARGO
00222           pGrid->particle[p].v3 -= qshear*Omega_0*x1;
00223 #endif
00224           pGrid->particle[p].v2 = 0.0;
00225 
00226           pGrid->particle[p].pos = GetPosition(&pGrid->particle[p]);
00227 
00228           pGrid->particle[p].my_id = p;
00229 #ifdef MPI_PARALLEL
00230           pGrid->particle[p].init_id = pGrid->my_id;
00231 #endif
00232           p += 1;
00233         }
00234       }
00235     }
00236   }
00237 #endif /* PARTICLES */
00238 
00239 /* enroll gravitational potential function, shearing sheet BC functions */
00240 
00241   StaticGravPot = ShearingBoxPot;
00242 
00243   return;
00244 }
00245 
00246 /*==============================================================================
00247  * PROBLEM USER FUNCTIONS:
00248  * problem_write_restart() - writes problem-specific user data to restart files
00249  * problem_read_restart()  - reads problem-specific user data from restart files
00250  * get_usr_expr()          - sets pointer to expression for special output data
00251  * get_usr_out_fun()       - returns a user defined output function pointer
00252  * get_usr_par_prop()      - returns a user defined particle selection function
00253  * Userwork_in_loop        - problem specific work IN     main loop
00254  * Userwork_after_loop     - problem specific work AFTER  main loop
00255  *----------------------------------------------------------------------------*/
00256 
00257 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00258 {
00259   return;
00260 }
00261 
00262 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00263 {
00264   Omega_0 = par_getd_def("problem","omega",1.0);
00265   qshear  = par_getd_def("problem","qshear",1.5);
00266 
00267   StaticGravPot = ShearingBoxPot;
00268   return;
00269 }
00270 
00271 #if (NSCALARS > 0)
00272 /*! \fn static Real ScalarDen(const Grid *pG,const int i,
00273  *                            const int j,const int k)
00274  *  \brief Scalar density */
00275 static Real ScalarDen(const Grid *pG, const int i, const int j, const int k)
00276 {
00277   return pG->U[k][j][i].s[0];
00278 }
00279 #endif
00280 
00281 #ifdef PARTICLES
00282 /*! \fn static Real dratio(const Grid *pG, const int i, const int j,const int k)
00283  *  \brief Density ratio */
00284 static Real dratio(const Grid *pG, const int i, const int j, const int k)
00285 {
00286 #if (NSCALARS > 0)
00287   return pG->Coup[k][j][i].grid_d/pG->U[k][j][i].s[0];
00288 #else
00289   return pG->Coup[k][j][i].grid_d/pG->U[k][j][i].d;
00290 #endif
00291 }
00292 #endif
00293 
00294 /*! \fn static Real expr_dV3(const Grid *pG,const int i,const int j,const int k)
00295  *  \brief 3-component of velocity */
00296 static Real expr_dV3(const Grid *pG, const int i, const int j, const int k)
00297 {
00298   Real x1,x2,x3;
00299   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00300 #ifdef FARGO
00301   return pG->U[k][j][i].M3/pG->U[k][j][i].d;
00302 #else
00303   return (pG->U[k][j][i].M3/pG->U[k][j][i].d + qshear*Omega_0*x1);
00304 #endif
00305 }
00306 
00307 #ifdef PARTICLES
00308 /*! \fn static Real expr_dV3par(const Grid *pG, const int i, const int j, 
00309  *                              const int k)
00310  *  \brief 3-component of particle velocity */
00311 static Real expr_dV3par(const Grid *pG, const int i, const int j, const int k)
00312 {
00313   Real x1,x2,x3;
00314   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00315 #ifdef FARGO
00316   return expr_V3par(pG, i, j, k);
00317 #else
00318   return expr_V3par(pG, i, j, k) + qshear*Omega_0*x1;
00319 #endif
00320 }
00321 #endif /* PARTICLES */
00322 
00323 Gasfun_t get_usr_expr(const char *expr)
00324 {
00325 #if (NSCALARS > 0)
00326   if(strcmp(expr,"scalar")==0) return ScalarDen;
00327 #endif
00328 #ifdef PARTICLES
00329   if(strcmp(expr,"dratio")==0) return dratio;
00330   if(strcmp(expr,"dVypar")==0) return expr_dV3par;
00331 #endif
00332   if(strcmp(expr,"dVy")==0) return expr_dV3;
00333   return NULL;
00334 }
00335 
00336 VGFunout_t get_usr_out_fun(const char *name){
00337   return NULL;
00338 }
00339 
00340 #ifdef PARTICLES
00341 PropFun_t get_usr_par_prop(const char *name)
00342 {
00343   return NULL;
00344 }
00345 
00346 void gasvshift(const Real x1, const Real x2, const Real x3,
00347                                     Real *u1, Real *u2, Real *u3)
00348 {
00349   return;
00350 }
00351 
00352 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00353                                     const Real v1, const Real v2, const Real v3)
00354 {
00355   return;
00356 }
00357 #endif
00358 
00359 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00360 {
00361   return;
00362 }
00363 
00364 /*---------------------------------------------------------------------------
00365  * Userwork_after_loop: computes L1-error in linear waves,
00366  * ASSUMING WAVE HAS PROPAGATED AN INTEGER NUMBER OF PERIODS
00367  * Must set parameters in input file appropriately so that this is true
00368  */
00369 
00370 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00371 {
00372   int i=0,j=0,k=0;
00373   int is,ie,js,je,ks,ke,n,m,nwave,samp;
00374   Real x1,x2,x3,x1max,x1min,x2max,x2min;
00375   Real time,kx,omg,omg2,v1x,v1y,amp,SolGasd,SolLagd,SolGasM1,SolGasM3;
00376   char *fname;
00377 
00378   is = pGrid->is; ie = pGrid->ie;
00379   js = pGrid->js; je = pGrid->je;
00380   ks = pGrid->ks; ke = pGrid->ke;
00381 
00382   /* Read initial conditions  */
00383   amp = par_getd("problem","amp");
00384   Omega_0 = par_getd_def("problem","omega",1.0);
00385   qshear  = par_getd_def("problem","qshear",1.5);
00386   nwave = par_geti("problem","nwave");
00387   samp = par_geti("problem","sample");
00388   x1min = par_getd("grid","x1min");
00389   x1max = par_getd("grid","x1max");
00390   x2min = par_getd("grid","x2min");
00391   x2max = par_getd("grid","x2max");
00392 
00393 /* calculate dispersion relation and find eigen vectors */
00394   kx = 2.0*(PI)*nwave/(x1max-x1min);
00395   omg2 = SQR(Iso_csound*kx)+2.0*(2.0-qshear)*SQR(Omega_0);
00396   omg = sqrt(omg2);
00397   v1x = omg*amp/kx;
00398   v1y = (2.0-qshear)*Omega_0/omg*v1x;
00399 
00400   time = pGrid->time;
00401 
00402 #ifdef PARTICLES
00403   /* Bin particles to grid */
00404   particle_to_grid(pGrid, pDomain, property_all);
00405 #endif
00406 
00407   /* Print error to file "Par_LinWave-errors.#.dat", where #=wavedir  */
00408   fname = ath_fname(NULL,"Par_Shwave1d-errors",0,0,NULL,"dat");
00409   /* Open output file in write mode */
00410   FILE *fid = fopen(fname,"w");
00411 
00412   /* Calculate the error and output */
00413   fprintf(fid, "%f      %d\n", time, ie-is+1);
00414   for (i=is; i<=ie; i++) {
00415     /* calculate analytic solution */
00416     cc_pos(pGrid,i,js,ks,&x1,&x2,&x3);
00417     SolGasd = 1.0+amp*cos(kx*x1-omg*time);
00418     SolGasM1 = SolGasd*v1x*cos(kx*x1);
00419     SolGasM3 = SolGasd*v1y*sin(kx*x1);
00420     if (samp == 1)
00421       SolLagd = SolGasd;
00422     else
00423       SolLagd = SolGasd-amp*cos(kx*x1);
00424     /* output */
00425     fprintf(fid,"%f     ",x1);
00426     fprintf(fid,"%e     %e      %e      ",pGrid->U[ks][js][i].d-1.0,
00427                                 SolGasd-1.0,pGrid->U[ks][js][i].d-SolGasd);
00428 #ifdef PARTICLES
00429     fprintf(fid,"%e     %e      %e      ", pG->Coup[ks][js][i].grid_d-1.0,
00430                         SolLagd-1.0, pG->Coup[ks][js][i].grid_d-SolLagd);
00431 #endif
00432 #ifdef FARGO
00433 //    fprintf(fid,"%e     %e      %e      ",pGrid->U[ks][js][i].d, pGrid->U[ks][js][i].M1, pGrid->U[ks][js][i].M3);
00434 #else
00435 //    fprintf(fid,"%e     %e      %e      ",pGrid->U[ks][js][i].d, pGrid->U[ks][js][i].M1, pGrid->U[ks][js][i].M3+pGrid->U[ks][js][i].d*qshear*Omega_0*x1);
00436 #endif
00437 #if (NSCALARS > 0)
00438     for (n=0; n<NSCALARS; n++)
00439       fprintf(fid,"%e   %e      ",pGrid->U[ks][js][i].s[n]-1.0,
00440                                   pGrid->U[ks][js][i].s[n]-SolLagd);
00441 #endif
00442     fprintf(fid,"\n");
00443   }
00444 
00445   fclose(fid);
00446 
00447   return;
00448 }
00449 
00450 
00451 /*=========================== PRIVATE FUNCTIONS ==============================*/
00452 /*--------------------------------------------------------------------------- */
00453 
00454 #ifdef PARTICLES
00455 /*! \fn int GetPosition(Grain *gr)
00456  *  \brief Get particle status (grid/crossing) */
00457 int GetPosition(Grain *gr)
00458 {
00459   if ((gr->x1>=x1upar) || (gr->x1<x1lpar) || (gr->x2>=x2upar) || (gr->x2<x2lpar) || (gr->x3>=x3upar) || (gr->x3<x3lpar))
00460     return 10; /* crossing particle */
00461   else
00462     return 1;  /* grid particle */
00463 }
00464 #endif
00465 
00466 /*------------------------------------------------------------------------------
00467  * ShearingBoxPot: 
00468  */
00469 
00470 /*! \fn static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00471  *  \brief shearing box tidal potential */
00472 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00473 {
00474   Real phi=0.0;
00475 #ifndef FARGO
00476   phi -= qshear*SQR(Omega_0*x1);
00477 #endif
00478   return phi;
00479 }