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 <time.h>
00024 #include "defs.h"
00025 #include "athena.h"
00026 #include "globals.h"
00027 #include "prototypes.h"
00028 #include "particles/particle.h"
00029
00030 #ifndef PARTICLES
00031 #error : The circular motion problem requires particles to be enabled.
00032 #endif
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 static Vector ParticleTroj(Real t);
00043 static Vector ParticleVel(Vector pos);
00044 static int ParticleLocator(Real x1, Real x2, Real x3);
00045 double ran2(long int *idum);
00046 extern Vector Get_Term(Grid *pG, int type, Real x1, Real x2, Real x3,
00047 Vector cell1, Real *tstop);
00048
00049 Real x1c,x2c,x3c;
00050 Real x01,x02,x03;
00051 Real omg, omgx1, omgx2, omgx3;
00052 Real r01,r02,r03,r0dn,r0cn1,r0cn2,r0cn3,n1,n2,n3;
00053 char name[50];
00054
00055
00056
00057
00058
00059
00060 void problem(Grid *pGrid, Domain *pDomain)
00061 {
00062 int i,j,k,in;
00063 long p, seed;
00064 Real x1min,x1max,x2min,x2max,x3min,x3max;
00065 Real ranv1, ranv2, ranv3, ranvnorm, tstop;
00066 Real theta, phi, rad, vran, thetaran, phiran;
00067 Vector cell1, vterm;
00068
00069 if (par_geti("grid","Nx1") == 1 || par_geti("grid","Nx2") == 1) {
00070 ath_error("[par_circ]: this test only works with Nx1 & Nx2 > 1\n");
00071 }
00072
00073
00074 if (pGrid->Nx1 > 1) cell1.x1 = 1.0/pGrid->dx1; else cell1.x1 = 0.0;
00075 if (pGrid->Nx2 > 1) cell1.x2 = 1.0/pGrid->dx2; else cell1.x2 = 0.0;
00076 if (pGrid->Nx3 > 1) cell1.x3 = 1.0/pGrid->dx3; else cell1.x3 = 0.0;
00077
00078
00079 x1min = par_getd("grid","x1min");
00080 x1max = par_getd("grid","x1max");
00081 x2min = par_getd("grid","x2min");
00082 x2max = par_getd("grid","x2max");
00083 x3min = par_getd("grid","x3min");
00084 x3max = par_getd("grid","x3max");
00085 x1c = 0.5*(x1min+x1max);
00086 x2c = 0.5*(x2min+x2max);
00087 x3c = 0.5*(x3min+x3max);
00088
00089
00090 omg = par_getd("problem","omega");
00091 theta = par_getd("problem","theta");
00092 phi = par_getd("problem","phi");
00093 rad = par_getd("problem","rad");
00094 vran = par_getd("problem","vran");
00095 seed = par_geti_def("problem","seed",0);
00096
00097 omgx1 = omg*sin(theta)*cos(phi);
00098 omgx2 = omg*sin(theta)*sin(phi);
00099 omgx3 = omg*cos(theta);
00100 n1 = omgx1/omg;
00101 n2 = omgx2/omg;
00102 n3 = omgx3/omg;
00103
00104
00105 if (par_geti("particle","partypes") != 1)
00106 ath_error("[par_circ]: number of particle types must be 1!\n");
00107
00108
00109 tstop0[0] = par_getd_def("problem","tstop",0.0);
00110 if (par_geti("particle","tsmode") != 3)
00111 ath_error("[par_circ]: This test works only for fixed stopping time!\n");
00112
00113
00114 r01 = rad*cos(theta)*cos(phi);
00115 r02 = rad*cos(theta)*sin(phi);
00116 r03 = -rad*sin(theta);
00117
00118
00119 x01 = x1c+r01;
00120 x02 = x2c+r02;
00121 x03 = x3c+r03;
00122
00123 r0dn = r01*n1+r02*n2+r03*n3;
00124 r0cn1 = n2*r03-n3*r02;
00125 r0cn2 = n3*r01-n1*r03;
00126 r0cn3 = n1*r02-n2*r01;
00127
00128 in = ParticleLocator(x01, x02, x03);
00129
00130 pGrid->nparticle = in;
00131 pGrid->grproperty[0].num = in;
00132
00133 if (pGrid->nparticle+2 > pGrid->arrsize)
00134 particle_realloc(pGrid, pGrid->nparticle+2);
00135
00136
00137 ranv1 = ran2(&seed)-0.5;
00138 ranv2 = ran2(&seed)-0.5;
00139 ranv3 = ran2(&seed)-0.5;
00140 ranvnorm = vran/sqrt(SQR(ranv1)+SQR(ranv2)+SQR(ranv3)+1.0e-8);
00141 ranv1 = ranvnorm*ranv1;
00142 ranv2 = ranvnorm*ranv2;
00143 ranv3 = ranvnorm*ranv3;
00144
00145
00146 for (k=pGrid->ks; k<=pGrid->ke; k++) {
00147 for (j=pGrid->js; j<=pGrid->je; j++) {
00148 for (i=pGrid->is; i<=pGrid->ie; i++) {
00149 pGrid->U[k][j][i].d = 1.0;
00150 pGrid->U[k][j][i].M1 = 0.0;
00151 pGrid->U[k][j][i].M2 = 0.0;
00152 pGrid->U[k][j][i].M3 = 0.0;
00153 }}}
00154
00155
00156 for (p=0; p<in; p++)
00157 {
00158 pGrid->particle[p].property = 0;
00159 pGrid->particle[p].x1 = x01;
00160 pGrid->particle[p].x2 = x02;
00161 pGrid->particle[p].x3 = x03;
00162
00163 vterm = Get_Term(pGrid,0,x01,x02,x03,cell1,&tstop);
00164
00165
00166
00167
00168 pGrid->particle[p].v1 = ranv1+vterm.x1;
00169 pGrid->particle[p].v2 = ranv2+vterm.x2;
00170 pGrid->particle[p].v3 = ranv3+vterm.x3;
00171
00172 pGrid->particle[p].pos = 1;
00173 pGrid->particle[p].my_id = p;
00174 #ifdef MPI_PARALLEL
00175 pGrid->particle[p].init_id = pGrid->my_id;
00176 #endif
00177 }
00178
00179 if (pGrid->my_id == 0) {
00180
00181 sprintf(name, "%s_Traj.dat", pGrid->outfilename);
00182 FILE *fid = fopen(name,"w");
00183 fclose(fid);
00184 #ifdef MPI_PARALLEL
00185 sprintf(name, "../%s_Traj.dat", pGrid->outfilename);
00186 #else
00187 sprintf(name, "%s_Traj.dat", pGrid->outfilename);
00188 #endif
00189 }
00190
00191 #ifdef MPI_PARALLEL
00192 MPI_Bcast(name,50,MPI_CHAR,0,MPI_COMM_WORLD);
00193 #endif
00194
00195 return;
00196 }
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00210 {
00211 fwrite(name, sizeof(char),50,fp);
00212 return;
00213 }
00214
00215 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00216 {
00217 Real x1min,x1max,x2min,x2max,x3min,x3max;
00218 Real theta, phi, rad;
00219
00220 x1min = par_getd("grid","x1min"); x1max = par_getd("grid","x1max");
00221 x2min = par_getd("grid","x2min"); x2max = par_getd("grid","x2max");
00222 x3min = par_getd("grid","x3min"); x3max = par_getd("grid","x3max");
00223 x1c = 0.5*(x1min+x1max); x2c = 0.5*(x2min+x2max); x3c = 0.5*(x3min+x3max);
00224
00225 omg = par_getd("problem","omega");
00226 theta = par_getd("problem","theta");
00227 phi = par_getd("problem","phi");
00228 rad = par_getd("problem","rad");
00229
00230 omgx1 = omg*sin(theta)*cos(phi); r01 = rad*cos(theta)*cos(phi);
00231 omgx2 = omg*sin(theta)*sin(phi); r02 = rad*cos(theta)*sin(phi);
00232 omgx3 = omg*cos(theta); r03 = -rad*sin(theta);
00233
00234 n1 = omgx1/omg; n2 = omgx2/omg; n3 = omgx3/omg;
00235 x01 = x1c+r01; x02 = x2c+r02; x03 = x3c+r03;
00236
00237 r0dn = r01*n1+r02*n2+r03*n3;
00238 r0cn1 = n2*r03-n3*r02;
00239 r0cn2 = n3*r01-n1*r03;
00240 r0cn3 = n1*r02-n2*r01;
00241
00242 fread(name, sizeof(char),50,fp);
00243 return;
00244 }
00245
00246 Gasfun_t get_usr_expr(const char *expr)
00247 {
00248 return NULL;
00249 }
00250
00251 VGFunout_t get_usr_out_fun(const char *name){
00252 return NULL;
00253 }
00254
00255 #ifdef PARTICLES
00256 PropFun_t get_usr_par_prop(const char *name)
00257 {
00258 return NULL;
00259 }
00260
00261
00262
00263
00264 void gasvshift(const Real x1, const Real x2, const Real x3,
00265 Real *u1, Real *u2, Real *u3)
00266 {
00267 Real dx1, dx2, dx3;
00268 dx1 = x1-x1c; dx2 = x2-x2c; dx3 = x3-x3c;
00269 *u1 = omgx2*dx3-omgx3*dx2;
00270 *u2 = omgx3*dx1-omgx1*dx3;
00271 *u3 = omgx1*dx2-omgx2*dx1;
00272
00273 return;
00274 }
00275
00276
00277
00278
00279 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00280 const Real v1, const Real v2, const Real v3)
00281 {
00282 return;
00283 }
00284 #endif
00285
00286 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00287 {
00288 long p;
00289 Real t,ds,dv,r1,r0,dr;
00290 Vector pos0, vel0, pos;
00291 Grain *gr;
00292 FILE *fid;
00293
00294 t = pGrid->time+pGrid->dt;
00295 pos0 = ParticleTroj(t);
00296
00297 for (p=0; p<pGrid->nparticle; p++)
00298 {
00299 gr = &(pGrid->particle[p]);
00300 if (gr->pos == 1)
00301 {
00302 pos.x1 = gr->x1;
00303 pos.x2 = gr->x2;
00304 pos.x3 = gr->x3;
00305 vel0 = ParticleVel(pos);
00306
00307 ds = sqrt(SQR(pos0.x1-gr->x1)+SQR(pos0.x2-gr->x2)+SQR(pos0.x3-gr->x3));
00308 dv = sqrt(SQR(vel0.x1-gr->v1)+SQR(vel0.x2-gr->v2)+SQR(vel0.x3-gr->v3));
00309 r0 = sqrt(SQR(pos0.x1-x1c)+SQR(pos0.x2-x2c)+SQR(pos0.x3-x3c));
00310 r1 = sqrt(SQR(gr->x1-x1c)+SQR(gr->x2-x2c)+SQR(gr->x3-x3c));
00311 dr = r1-r0;
00312
00313 fid = fopen(name,"a+");
00314
00315 fprintf(fid,"%e %e %e\n", t, ds, dr);
00316 fclose(fid);
00317 }
00318 }
00319 return;
00320 }
00321
00322
00323
00324
00325
00326
00327
00328 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00329 {
00330 return;
00331 }
00332
00333
00334
00335
00336
00337
00338 static Vector ParticleTroj(Real t)
00339 {
00340 Vector pos;
00341 Real r1,r2,r3,omgt,sinot, cosot;
00342
00343 omgt = omg*t;
00344 sinot = sin(omgt); cosot = cos(omgt);
00345 r1 = r0dn*n1 + sinot*r0cn1 + cosot*(r01-r0dn*n1);
00346 r2 = r0dn*n2 + sinot*r0cn2 + cosot*(r02-r0dn*n2);
00347 r3 = r0dn*n3 + sinot*r0cn3 + cosot*(r03-r0dn*n3);
00348
00349 pos.x1 = x1c+r1; pos.x2 = x2c+r2; pos.x3 = x3c+r3;
00350 return pos;
00351 }
00352
00353
00354
00355 static Vector ParticleVel(Vector pos)
00356 {
00357 Vector vel;
00358 Real r1,r2,r3;
00359
00360 r1 = pos.x1-x1c; r2 = pos.x2-x2c; r3 = pos.x3-x3c;
00361 vel.x1 = omgx2*r3-omgx3*r2;
00362 vel.x2 = omgx3*r1-omgx1*r3;
00363 vel.x3 = omgx1*r2-omgx2*r1;
00364
00365 return vel;
00366 }
00367
00368
00369
00370 static int ParticleLocator(Real x1, Real x2, Real x3)
00371 {
00372 if ((x1<x1upar) && (x1>=x1lpar) && (x2<x2upar) && (x2>=x2lpar) &&(x3<x3upar) && (x3>=x3lpar))
00373 return 1;
00374 else
00375 return 0;
00376 }
00377
00378
00379
00380 #define DBL_EPSILON 1.0e-16
00381 #define IM1 2147483563
00382 #define IM2 2147483399
00383 #define AM (1.0/IM1)
00384 #define IMM1 (IM1-1)
00385 #define IA1 40014
00386 #define IA2 40692
00387 #define IQ1 53668
00388 #define IQ2 52774
00389 #define IR1 12211
00390 #define IR2 3791
00391 #define NTAB 32
00392 #define NDIV (1+IMM1/NTAB)
00393 #define RNMX (1.0-DBL_EPSILON)
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408 double ran2(long int *idum)
00409 {
00410 int j;
00411 long int k;
00412 static long int idum2=123456789;
00413 static long int iy=0;
00414 static long int iv[NTAB];
00415 double temp;
00416
00417 if (*idum <= 0) {
00418 if (-(*idum) < 1) *idum=1;
00419 else *idum = -(*idum);
00420 idum2=(*idum);
00421 for (j=NTAB+7;j>=0;j--) {
00422 k=(*idum)/IQ1;
00423 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00424 if (*idum < 0) *idum += IM1;
00425 if (j < NTAB) iv[j] = *idum;
00426 }
00427 iy=iv[0];
00428 }
00429 k=(*idum)/IQ1;
00430 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00431 if (*idum < 0) *idum += IM1;
00432 k=idum2/IQ2;
00433 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00434 if (idum2 < 0) idum2 += IM2;
00435 j=(int)(iy/NDIV);
00436 iy=iv[j]-idum2;
00437 iv[j] = *idum;
00438 if (iy < 1) iy += IMM1;
00439 if ((temp=AM*iy) > RNMX) return RNMX;
00440 else return temp;
00441 }
00442
00443 #undef IM1
00444 #undef IM2
00445 #undef AM
00446 #undef IMM1
00447 #undef IA1
00448 #undef IA2
00449 #undef IQ1
00450 #undef IQ2
00451 #undef IR1
00452 #undef IR2
00453 #undef NTAB
00454 #undef NDIV
00455 #undef RNMX
00456 #undef DBL_EPSILON