00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include <float.h>
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include <time.h>
00021 #include "defs.h"
00022 #include "athena.h"
00023 #include "globals.h"
00024 #include "prototypes.h"
00025 #include "particles/particle.h"
00026
00027 #ifndef SHEARING_BOX
00028 #error : The epicycle problem requires shearing-box to be enabled.
00029 #endif
00030
00031 #ifndef PARTICLES
00032 #error : The epicycle problem requires particles to be enabled.
00033 #endif
00034
00035 #ifndef ISOTHERMAL
00036 #error : The epicycle problem requires isothermal equation of state.
00037 #endif
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3);
00049 static Vector ParticlePosition(const Real t);
00050 static Vector ParticleVelocity(const Vector pos, const Real t);
00051 static int ParticleLocator(const Vector pos);
00052
00053
00054 char name[50];
00055
00056 Real amp, Lx, Ly, x1min, x1max, x2min, x2max, omg;
00057
00058
00059
00060
00061
00062
00063
00064 void problem(Grid *pGrid, Domain *pDomain)
00065 {
00066 int in,i,j,k;
00067 Real x1,x2,x3;
00068 long p;
00069 Vector parpos, parvel;
00070
00071 if (par_geti("grid","Nx2") == 1) {
00072 ath_error("[par_epicycle]: par_epicycle must work in 2D or 3D.\n");
00073 }
00074
00075
00076 x1min = par_getd("grid","x1min");
00077 x1max = par_getd("grid","x1max");
00078 Lx = x1max - x1min;
00079
00080 x2min = par_getd("grid","x2min");
00081 x2max = par_getd("grid","x2max");
00082 Ly = x2max - x2min;
00083 if (par_geti("grid","Nx3") == 1) {
00084 Ly = 0.0;
00085 }
00086
00087
00088 Omega_0 = par_getd("problem","omega");
00089 qshear = par_getd_def("problem","qshear",1.5);
00090 amp = par_getd("problem","amp");
00091 omg = sqrt(2.0*(2.0-qshear))*Omega_0;
00092
00093
00094 if (par_geti("particle","partypes") != 1)
00095 ath_error("[par_epicycle]: This test only allows ONE particle species!\n");
00096
00097
00098 tstop0[0] = par_getd_def("problem","tstop",1.0e20);
00099 if (par_geti("particle","tsmode") != 3)
00100 ath_error("[par_epicycle]: This test only allows fixed stopping time!\n");
00101
00102
00103 parpos = ParticlePosition(0.0);
00104 parvel = ParticleVelocity(parpos, 0.0);
00105 in = ParticleLocator(parpos);
00106
00107 pGrid->nparticle = in;
00108 pGrid->grproperty[0].num = in;
00109
00110 if (pGrid->nparticle+2 > pGrid->arrsize)
00111 particle_realloc(pGrid, pGrid->nparticle+2);
00112
00113
00114 for (k=pGrid->ks; k<=pGrid->ke; k++) {
00115 for (j=pGrid->js; j<=pGrid->je; j++) {
00116 for (i=pGrid->is; i<=pGrid->ie; i++) {
00117 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00118 pGrid->U[k][j][i].d = 1.0;
00119 pGrid->U[k][j][i].M1 = 0.0;
00120 pGrid->U[k][j][i].M2 = 0.0;
00121 pGrid->U[k][j][i].M3 = 0.0;
00122 #ifndef FARGO
00123 if (Ly>0.0)
00124 pGrid->U[k][j][i].M2 -= qshear*Omega_0*x1;
00125 else
00126 pGrid->U[k][j][i].M3 -= qshear*Omega_0*x1;
00127 #endif
00128 }}}
00129
00130
00131 for (p=0; p<in; p++)
00132 {
00133 pGrid->particle[p].property = 0;
00134 pGrid->particle[p].x1 = parpos.x1;
00135 pGrid->particle[p].x2 = parpos.x2;
00136 pGrid->particle[p].x3 = parpos.x3;
00137 pGrid->particle[p].v1 = parvel.x1;
00138 pGrid->particle[p].v2 = parvel.x2;
00139 pGrid->particle[p].v3 = parvel.x3;
00140 pGrid->particle[p].pos = 1;
00141 pGrid->particle[p].my_id = p;
00142 #ifdef MPI_PARALLEL
00143 pGrid->particle[p].init_id = pGrid->my_id;
00144 #endif
00145 }
00146
00147
00148 StaticGravPot = ShearingBoxPot;
00149
00150 if (pGrid->my_id == 0) {
00151
00152 sprintf(name, "%s_Traj.dat", pGrid->outfilename);
00153 FILE *fid = fopen(name,"w");
00154 fclose(fid);
00155 #ifdef MPI_PARALLEL
00156 sprintf(name, "../%s_Traj.dat", pGrid->outfilename);
00157 #else
00158 sprintf(name, "%s_Traj.dat", pGrid->outfilename);
00159 #endif
00160 }
00161
00162 #ifdef MPI_PARALLEL
00163 MPI_Bcast(name,50,MPI_CHAR,0,MPI_COMM_WORLD);
00164 #endif
00165
00166 return;
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00180 {
00181 fwrite(name, sizeof(char),50,fp);
00182 return;
00183 }
00184
00185 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00186 {
00187 Omega_0 = par_getd("problem","omega");
00188 qshear = par_getd_def("problem","qshear",1.5);
00189
00190 StaticGravPot = ShearingBoxPot;
00191
00192 Ly = x2max - x2min;
00193 if (par_geti("grid","Nx3") == 1) {
00194 Ly = 0.0;
00195 }
00196
00197 amp = par_getd("problem","amp");
00198 omg = sqrt(2.0*(2.0-qshear))*Omega_0;
00199
00200 fread(name, sizeof(char),50,fp);
00201 return;
00202 }
00203
00204 Gasfun_t get_usr_expr(const char *expr)
00205 {
00206 return NULL;
00207 }
00208
00209 VGFunout_t get_usr_out_fun(const char *name){
00210 return NULL;
00211 }
00212
00213 #ifdef PARTICLES
00214 PropFun_t get_usr_par_prop(const char *name)
00215 {
00216 return NULL;
00217 }
00218
00219 void gasvshift(const Real x1, const Real x2, const Real x3,
00220 Real *u1, Real *u2, Real *u3)
00221 {
00222 return;
00223 }
00224
00225 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00226 const Real v1, const Real v2, const Real v3)
00227 {
00228 return;
00229 }
00230 #endif
00231
00232 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00233 {
00234 long p;
00235 Real t, ds, E, vshift;
00236 Vector pos0;
00237 Grain *gr;
00238 FILE *fid;
00239
00240 t = pGrid->time+pGrid->dt;
00241 pos0 = ParticlePosition(t);
00242 for (p=0; p<pGrid->nparticle; p++)
00243 {
00244 gr = &(pGrid->particle[p]);
00245 if (gr->pos == 1)
00246 {
00247
00248 ds = sqrt(SQR(gr->x1-pos0.x1)+SQR(gr->x2-pos0.x2)+SQR(gr->x3-pos0.x3));
00249
00250
00251 E = 0.5*SQR(gr->v1) - qshear*SQR(Omega_0*gr->x1);
00252 #ifdef FARGO
00253 vshift = -qshear*Omega_0*gr->x1;
00254 #else
00255 vshift = 0.0;
00256 #endif
00257 if (Ly>0.0)
00258 E += 0.5*SQR(gr->v2+vshift);
00259 else
00260 E += 0.5*SQR(gr->v3+vshift);
00261
00262
00263 fid = fopen(name,"a+");
00264 fprintf(fid,"%e %e %e %e %e %e %e %e %e\n", t, ds, E, gr->x1, gr->x2, gr->x3, pos0.x1, pos0.x2, pos0.x3);
00265 fclose(fid);
00266 }
00267 }
00268 return;
00269 }
00270
00271
00272
00273
00274
00275
00276
00277 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00278 {
00279 return;
00280 }
00281
00282
00283
00284
00285
00286 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00287 {
00288 Real phi=0.0;
00289 #ifndef FARGO
00290 phi -= qshear*Omega_0*Omega_0*x1*x1;
00291 #endif
00292 return phi;
00293 }
00294
00295
00296
00297 static Vector ParticlePosition(const Real t)
00298 {
00299 Real x,y;
00300 Vector pos;
00301 x = amp*cos(omg*t);
00302 y = -2.0*amp*Omega_0/omg*sin(omg*t);
00303
00304 x = x-floor((x-x1min)/Lx)*Lx;
00305 if (Ly > 0.0)
00306 y = y-floor((y-x2min)/Ly)*Ly;
00307 else
00308 y = 0.0;
00309
00310 pos.x1 = x; pos.x2 = y; pos.x3 = 0.0;
00311 return pos;
00312 }
00313
00314
00315
00316 static Vector ParticleVelocity(const Vector pos, const Real t)
00317 {
00318 Real vx,vy;
00319 Vector vel;
00320 vx = -amp*omg*sin(omg*t);
00321 vy = -2.0*amp*Omega_0*cos(omg*t);
00322 #ifdef FARGO
00323 vy = vy + qshear*Omega_0*pos.x1;
00324 #endif
00325 if (Ly>0.0) {
00326 vel.x1 = vx; vel.x2 = vy; vel.x3 = 0.0;
00327 } else {
00328 vel.x1 = vx; vel.x3 = vy; vel.x2 = 0.0;
00329 }
00330 return vel;
00331 }
00332
00333
00334
00335 static int ParticleLocator(const Vector pos)
00336 {
00337 if ((pos.x1<x1upar) && (pos.x1>=x1lpar) && (pos.x2<x2upar)
00338 && (pos.x2>=x2lpar) &&(pos.x3<x3upar) && (pos.x3>=x3lpar))
00339 return 1;
00340 else
00341 return 0;
00342 }