00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
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
00029
00030
00031
00032
00033
00034 int GetPosition(Grain *gr);
00035
00036
00037
00038
00039
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
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
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
00115 #ifdef PARTICLES
00116
00117
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
00128 tstop0[0] = par_getd_def("problem","tstop",0.0);
00129 if (par_geti("particle","tsmode") != 3)
00130 ath_error("[par_linwave1d]: This test only allows fixed stopping time!\n");
00131
00132
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
00147
00148
00149
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
00216
00217 return;
00218 }
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
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
00243
00244
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
00253
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
00305
00306
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
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
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
00342 particle_to_grid(pGrid, pDomain, property_all);
00343
00344
00345 fname = ath_fname(NULL,"Par_LinWave1d-errors",1,wavedir,NULL,"dat");
00346
00347
00348 FILE *fid = fopen(fname,"w");
00349
00350
00351 switch(wavedir){
00352 case 1:
00353 fprintf(fid, "%f %d\n", time, ie-is+1);
00354 for (i=is; i<=ie; i++) {
00355
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
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
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
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
00413
00414
00415
00416
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;
00421 else
00422 return 1;
00423 }