00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <math.h>
00014 #include <stdio.h>
00015 #include <stdlib.h>
00016 #include "defs.h"
00017 #include "athena.h"
00018 #include "globals.h"
00019 #include "prototypes.h"
00020 #include "particles/particle.h"
00021
00022 #ifndef SHEARING_BOX
00023 #error : The shear wave problem requires shearing-box to be enabled.
00024 #endif
00025
00026 int Nx12;
00027
00028
00029
00030
00031
00032
00033
00034 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3);
00035 #ifdef PARTICLES
00036 int GetPosition(Grain *gr);
00037 static int property_mybin(const Grain *gr, const GrainAux *grsub);
00038 extern Real expr_V2par(const Grid *pG, const int i, const int j, const int k);
00039 #endif
00040
00041
00042
00043
00044
00045
00046 void problem(Grid *pGrid, Domain *pDomain)
00047 {
00048 int i=0,j=0,k=0;
00049 int is,ie,js,je,ks,ke,n,m,nwave,samp;
00050 Real x1,x2,x3,x1max,x1min,x2max,x2min;
00051 Real ky,omg,omg2,amp,v1x,v1y;
00052 #ifdef PARTICLES
00053 long p;
00054 int Npar,ip,jp,kp;
00055 Real x1p,x2p,x3p,x1l,x1u,x2l,x2u,x3l,x3u;
00056 #endif
00057
00058 if (par_geti("grid","Nx3") == 1) {
00059 ath_error("[par_shwave2d]: par_shwave2d must work in 3D grid.\n");
00060 }
00061
00062 is = pGrid->is; ie = pGrid->ie;
00063 js = pGrid->js; je = pGrid->je;
00064 ks = pGrid->ks; ke = pGrid->ke;
00065
00066
00067 amp = par_getd("problem","amp");
00068 Omega_0 = par_getd_def("problem","omega",1.0);
00069 qshear = par_getd_def("problem","qshear",1.5);
00070 nwave = par_geti("problem","nwave");
00071 samp = par_geti("problem","sample");
00072 x1min = par_getd("grid","x1min");
00073 x1max = par_getd("grid","x1max");
00074 x2min = par_getd("grid","x2min");
00075 x2max = par_getd("grid","x2max");
00076
00077 Nx12 = pGrid->Nx1*pGrid->Nx2;
00078
00079 if (nwave <= 0)
00080 ath_error("[par_shwave2d]: nwave must be positive!\n");
00081
00082
00083 ky = 2.0*(PI)*nwave/(x2max-x2min);
00084 omg2 = SQR(Iso_csound*ky)+2.0*(2.0-qshear)*SQR(Omega_0);
00085 omg = sqrt(omg2);
00086 v1y = omg*amp/ky;
00087 v1x = 2.0*Omega_0/omg*v1y;
00088
00089
00090
00091 for (k=ks; k<=ke; k++) {
00092 for (j=js; j<=je; j++) {
00093 for (i=is; i<=ie; i++) {
00094 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00095 pGrid->U[k][j][i].d = 1.0+amp*cos(ky*x2);
00096 pGrid->U[k][j][i].M1 = -pGrid->U[k][j][i].d*v1x*sin(ky*x2);
00097 pGrid->U[k][j][i].M2 = pGrid->U[k][j][i].d*v1y*cos(ky*x2);
00098 #ifndef FARGO
00099 pGrid->U[k][j][i].M2 -= pGrid->U[k][j][i].d*qshear*Omega_0*x1;
00100 #endif
00101 pGrid->U[k][j][i].M3 = 0.0;
00102 #if (NSCALARS > 0)
00103 if (samp == 1)
00104 for (n=0; n<NSCALARS; n++)
00105 pGrid->U[k][j][i].s[n] = pGrid->U[k][j][i].d;
00106 else
00107 for (n=0; n<NSCALARS; n++)
00108 pGrid->U[k][j][i].s[n] = 1.0;
00109 #endif
00110 }}}
00111
00112
00113 #ifdef PARTICLES
00114
00115
00116 if (par_geti("particle","partypes") != 1)
00117 ath_error("[par_shwave2d]: This test only allows ONE particle species!\n");
00118
00119 Npar = (int)(sqrt(par_geti("particle","parnumcell")));
00120 pGrid->nparticle = Npar*Npar*Npar*pGrid->Nx1*pGrid->Nx2*pGrid->Nx3;
00121 pGrid->grproperty[0].num = pGrid->nparticle;
00122 if (pGrid->nparticle+2 > pGrid->arrsize)
00123 particle_realloc(pGrid, pGrid->nparticle+2);
00124
00125
00126 tstop0[0] = par_getd_def("problem","tstop",0.0);
00127 if (par_geti("particle","tsmode") != 3)
00128 ath_error("[par_shwave2d]: This test only allows fixed stopping time!\n");
00129
00130
00131 p = 0;
00132
00133 for (k=ks; k<=ke; k++)
00134 {
00135 x3l = pGrid->x3_0 + (k+pGrid->kdisp)*pGrid->dx3;
00136 x3u = pGrid->x3_0 + ((k+pGrid->kdisp)+1.0)*pGrid->dx3;
00137
00138 for (j=js; j<=je; j++)
00139 {
00140 x2l = pGrid->x2_0 + (j+pGrid->jdisp)*pGrid->dx2;
00141 x2u = pGrid->x2_0 + ((j+pGrid->jdisp)+1.0)*pGrid->dx2;
00142
00143 for (i=is; i<=ie; i++)
00144 {
00145 x1l = pGrid->x1_0 + (i + pGrid->idisp)*pGrid->dx1;
00146 x1u = pGrid->x1_0 + ((i + pGrid->idisp) + 1.0)*pGrid->dx1;
00147
00148 for (ip=0;ip<Npar;ip++)
00149 {
00150 x1p = x1l+(x1u-x1l)/Npar*(ip+0.5);
00151
00152 for (jp=0;jp<Npar;jp++)
00153 {
00154 x2p = x2l+(x2u-x2l)/Npar*(jp+0.5);
00155
00156 for (kp=0;kp<Npar;kp++)
00157 {
00158 x3p = x3l+(x3u-x3l)/Npar*(kp+0.5);
00159
00160 pGrid->particle[p].property = 0;
00161
00162 pGrid->particle[p].x1 = x1p;
00163 pGrid->particle[p].x2 = x2p;
00164 if (samp == 1)
00165 pGrid->particle[p].x2 += (-amp*sin(ky*x2p));
00166
00167 pGrid->particle[p].x3 = x3p;
00168
00169 pGrid->particle[p].v1 = -v1x*sin(ky*x2p);
00170 pGrid->particle[p].v2 = v1y*cos(ky*x2p);
00171 #ifndef FARGO
00172 pGrid->particle[p].v2 -= qshear*Omega_0*x1p;
00173 #endif
00174 pGrid->particle[p].v3 = 0.0;
00175
00176 pGrid->particle[p].pos = GetPosition(&pGrid->particle[p]);
00177
00178 pGrid->particle[p].my_id = p;
00179 #ifdef MPI_PARALLEL
00180 pGrid->particle[p].init_id = pGrid->my_id;
00181 #endif
00182 p += 1;
00183 }
00184 }
00185 }
00186 }
00187 }
00188 }
00189
00190 #endif
00191
00192
00193
00194 StaticGravPot = ShearingBoxPot;
00195
00196 return;
00197 }
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00211 {
00212 return;
00213 }
00214
00215 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00216 {
00217 Omega_0 = par_getd_def("problem","omega",1.0e-3);
00218 qshear = par_getd_def("problem","qshear",1.5);
00219
00220 StaticGravPot = ShearingBoxPot;
00221 return;
00222 }
00223
00224 #if (NSCALARS > 0)
00225
00226
00227
00228 static Real ScalarDen(const Grid *pG, const int i, const int j, const int k)
00229 {
00230 return pG->U[k][j][i].s[0]-1.0;
00231 }
00232 #endif
00233
00234
00235
00236 static Real diffd(const Grid *pG, const int i, const int j, const int k)
00237 {
00238 return pG->U[k][j][i].d-1.0;
00239 }
00240
00241
00242
00243
00244 static Real expr_dV2(const Grid *pG, const int i, const int j, const int k)
00245 {
00246 Real x1,x2,x3;
00247 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00248 #ifdef FARGO
00249 return pG->U[k][j][i].M2/pG->U[k][j][i].d;
00250 #else
00251 return (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00252 #endif
00253 }
00254
00255 #ifdef PARTICLES
00256
00257
00258
00259
00260 {
00261 return pG->Coup[k][j][i].grid_d-1.0;
00262 }
00263
00264
00265
00266
00267 static Real expr_dV2par(const Grid *pG, const int i, const int j, const int k)
00268 {
00269 Real x1,x2,x3,v2par;
00270 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00271 #ifdef FARGO
00272 v2par = expr_V2par(pG, i, j, k);
00273 #else
00274 v2par = expr_V2par(pG, i, j, k) + qshear*Omega_0*x1;
00275 #endif
00276 return v2par;
00277 }
00278 #endif
00279
00280 Gasfun_t get_usr_expr(const char *expr)
00281 {
00282 #if (NSCALARS > 0)
00283 if(strcmp(expr,"scalar")==0) return ScalarDen;
00284 #endif
00285 if(strcmp(expr,"difd")==0) return diffd;
00286 if(strcmp(expr,"dVy")==0) return expr_dV2;
00287 #ifdef PARTICLES
00288 if(strcmp(expr,"dVypar")==0) return expr_dV2par;
00289 if(strcmp(expr,"difdpar")==0) return diffdpar;
00290 #endif
00291 return NULL;
00292 }
00293
00294 VGFunout_t get_usr_out_fun(const char *name){
00295 return NULL;
00296 }
00297
00298 #ifdef PARTICLES
00299 PropFun_t get_usr_par_prop(const char *name)
00300 {
00301 if (strcmp(name,"plane")==0) return property_mybin;
00302 return NULL;
00303 }
00304
00305 void gasvshift(const Real x1, const Real x2, const Real x3,
00306 Real *u1, Real *u2, Real *u3)
00307 {
00308 return;
00309 }
00310
00311 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00312 const Real v1, const Real v2, const Real v3)
00313 {
00314 return;
00315 }
00316 #endif
00317
00318 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00319 {
00320 return;
00321 }
00322
00323
00324
00325
00326
00327
00328
00329 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00330 {
00331 int i=0,j=0,k=0;
00332 int is,js,je,ks,n;
00333 Real x1,x2,x3,vshear,time;
00334 char *fname;
00335
00336 is = (pGrid->is+pGrid->ie)/2;
00337 js = pGrid->js; je = pGrid->je;
00338 ks = pGrid->ks;
00339
00340 time = pGrid->time;
00341
00342
00343 particle_to_grid(pGrid, pDomain, property_all);
00344
00345
00346 fname = ath_fname(NULL,"Par_Shwave2d-errors",0,0,NULL,"dat");
00347
00348 FILE *fid = fopen(fname,"w");
00349
00350
00351 fprintf(fid, "%f %d\n", time, je-js+1);
00352 for (j=js; j<=je; j++) {
00353 cc_pos(pGrid,is,j,ks,&x1,&x2,&x3);
00354 #ifdef FARGO
00355 vshear = 0.0;
00356 #else
00357 vshear = -qshear*Omega_0*x1;
00358 #endif
00359 fprintf(fid,"%f ",x2);
00360 fprintf(fid,"%e ",pGrid->U[ks][j][is].d-1.0);
00361 #ifdef PARTICLES
00362 fprintf(fid,"%e %e ",pGrid->Coup[ks][j][is].grid_d-1.0,
00363 pGrid->Coup[ks][j][is].grid_d-pGrid->U[ks][j][is].d);
00364 #if (NSCALARS > 0)
00365 for (n=0; n<NSCALARS; n++)
00366 fprintf(fid,"%e %e ",pGrid->U[ks][j][is].s[n]-1.0,
00367 pGrid->Coup[ks][j][is].grid_d-pGrid->U[ks][j][is].s[n]);
00368 #endif
00369 #endif
00370 fprintf(fid,"\n");
00371 }
00372
00373 fclose(fid);
00374
00375 return;
00376 }
00377
00378
00379
00380
00381
00382
00383
00384
00385 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00386 {
00387 Real phi=0.0;
00388 #ifndef FARGO
00389 phi -= qshear*SQR(Omega_0*x1);
00390 #endif
00391 return phi;
00392 }
00393
00394 #ifdef PARTICLES
00395
00396
00397 int GetPosition(Grain *gr)
00398 {
00399 if ((gr->x1>=x1upar) || (gr->x1<x1lpar) || (gr->x2>=x2upar) || (gr->x2<x2lpar) || (gr->x3>=x3upar) || (gr->x3<x3lpar))
00400 return 10;
00401 else
00402 return 1;
00403 }
00404
00405
00406
00407
00408
00409 static int property_mybin(const Grain *gr, const GrainAux *grsub)
00410 {
00411 if ((gr->my_id<Nx12) && (gr->pos == 1))
00412 return 1;
00413 else
00414 return 0;
00415 }
00416 #endif