00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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
00031
00032
00033
00034
00035
00036
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
00044
00045
00046
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
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
00085
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
00096 tstop0[0] = par_getd_def("problem","tstop",0.0);
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
00103
00104
00105 if (ipert == 1)
00106 {
00107 if (nwave <= 0)
00108 ath_error("[par_shwave1d]: nwave must be positive!\n");
00109
00110
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
00118
00119
00120 #endif
00121
00122
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
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
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
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
00238
00239
00240
00241 StaticGravPot = ShearingBoxPot;
00242
00243 return;
00244 }
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
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
00273
00274
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
00283
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
00295
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
00309
00310
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
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
00366
00367
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
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
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
00404 particle_to_grid(pGrid, pDomain, property_all);
00405 #endif
00406
00407
00408 fname = ath_fname(NULL,"Par_Shwave1d-errors",0,0,NULL,"dat");
00409
00410 FILE *fid = fopen(fname,"w");
00411
00412
00413 fprintf(fid, "%f %d\n", time, ie-is+1);
00414 for (i=is; i<=ie; i++) {
00415
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
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
00434 #else
00435
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
00452
00453
00454 #ifdef PARTICLES
00455
00456
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;
00461 else
00462 return 1;
00463 }
00464 #endif
00465
00466
00467
00468
00469
00470
00471
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 }