Athena: prob/par_linearwave1d.c Source File

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