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