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 <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 PARTICLES
00028 #error : The friction problem requires particles to be enabled.
00029 #endif
00030
00031
00032
00033
00034
00035
00036
00037
00038 static Vector ParticleTroj(Real t);
00039 static Vector ParticleVel(Real t);
00040 static int ParticleLocator(Real x1, Real x2, Real x3);
00041
00042
00043 Real x1c,x2c,x3c,v01,v02,v03;
00044 Real x1min,x1max,x2min,x2max,x3min,x3max;
00045 char name[50];
00046
00047
00048
00049
00050
00051
00052 void problem(Grid *pGrid, Domain *pDomain)
00053 {
00054 int i,j,k;
00055 long p,in;
00056
00057
00058 if (par_geti("grid","Nx1") == 1 || par_geti("grid","Nx2") == 1) {
00059 ath_error("[par_fric]: this test only works with Nx1 & Nx2 > 1\n");
00060 }
00061
00062
00063 x1min = par_getd("grid","x1min");
00064 x1max = par_getd("grid","x1max");
00065 x2min = par_getd("grid","x2min");
00066 x2max = par_getd("grid","x2max");
00067 x3min = par_getd("grid","x3min");
00068 x3max = par_getd("grid","x3max");
00069 x1c = 0.5*(x1min+x1max);
00070 x2c = 0.5*(x2min+x2max);
00071 x3c = 0.5*(x3min+x3max);
00072
00073
00074 v01 = par_getd("problem","v1");
00075 v02 = par_getd("problem","v2");
00076 v03 = par_getd("problem","v3");
00077
00078
00079 if (par_geti("particle","partypes") != 1)
00080 ath_error("[par_fric]: number of particle types must be 1!\n");
00081
00082
00083 tstop0 = par_getd("problem","tstop");
00084 if (par_geti("particle","tsmode") != 3)
00085 ath_error("[par_fric]: This test works only for fixed stopping time!\n");
00086
00087
00088 in = ParticleLocator(x1c, x2c, x3c);
00089
00090 pGrid->nparticle = in;
00091 pGrid->grproperty[0].num = in;
00092
00093 if (pGrid->nparticle+2 > pGrid->arrsize)
00094 particle_realloc(pGrid, pGrid->nparticle+2);
00095
00096
00097
00098 for (k=pGrid->ks; k<=pGrid->ke; k++) {
00099 for (j=pGrid->js; j<=pGrid->je; j++) {
00100 for (i=pGrid->is; i<=pGrid->ie; i++) {
00101 pGrid->U[k][j][i].d = 1.0;
00102 pGrid->U[k][j][i].M1 = 0.0;
00103 pGrid->U[k][j][i].M2 = 0.0;
00104 pGrid->U[k][j][i].M3 = 0.0;
00105 }}}
00106
00107
00108 for (p=0; p<in; p++)
00109 {
00110 pGrid->particle[p].property = 0;
00111 pGrid->particle[p].x1 = x1c;
00112 pGrid->particle[p].x2 = x2c;
00113 pGrid->particle[p].x3 = x3c;
00114 pGrid->particle[p].v1 = v01;
00115 pGrid->particle[p].v2 = v02;
00116 pGrid->particle[p].v3 = v03;
00117 pGrid->particle[p].pos = 1;
00118 pGrid->particle[p].my_id = p;
00119 #ifdef MPI_PARALLEL
00120 pGrid->particle[p].init_id = pGrid->my_id;
00121 #endif
00122 }
00123
00124 if (pGrid->my_id == 0) {
00125
00126 sprintf(name, "%s_Err.dat", pGrid->outfilename);
00127 FILE *fid = fopen(name,"w");
00128 fclose(fid);
00129 #ifdef MPI_PARALLEL
00130 sprintf(name, "../%s_Err.dat", pGrid->outfilename);
00131 #else
00132 sprintf(name, "%s_Err.dat", pGrid->outfilename);
00133 #endif
00134 }
00135
00136 #ifdef MPI_PARALLEL
00137 MPI_Bcast(name,50,MPI_CHAR,0,MPI_COMM_WORLD);
00138 #endif
00139
00140 return;
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00154 {
00155 fwrite(name, sizeof(char),50,fp);
00156 return;
00157 }
00158
00159 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00160 {
00161
00162 x1min = par_getd("grid","x1min");
00163 x1max = par_getd("grid","x1max");
00164 x2min = par_getd("grid","x2min");
00165 x2max = par_getd("grid","x2max");
00166 x3min = par_getd("grid","x3min");
00167 x3max = par_getd("grid","x3max");
00168 x1c = 0.5*(x1min+x1max);
00169 x2c = 0.5*(x2min+x2max);
00170 x3c = 0.5*(x3min+x3max);
00171
00172
00173 v01 = par_getd("problem","v1");
00174 v02 = par_getd("problem","v2");
00175 v03 = par_getd("problem","v3");
00176
00177 fread(name, sizeof(char),50,fp);
00178
00179 return;
00180 }
00181
00182 Gasfun_t get_usr_expr(const char *expr)
00183 {
00184 return NULL;
00185 }
00186
00187 VGFunout_t get_usr_out_fun(const char *name){
00188 return NULL;
00189 }
00190
00191 #ifdef PARTICLES
00192 PropFun_t get_usr_par_prop(const char *name)
00193 {
00194 return NULL;
00195 }
00196
00197 void gasvshift(const Real x1, const Real x2, const Real x3,
00198 Real *u1, Real *u2, Real *u3)
00199 {
00200 return;
00201 }
00202
00203 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00204 const Real v1, const Real v2, const Real v3)
00205 {
00206 return;
00207 }
00208 #endif
00209
00210 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00211 {
00212 long p;
00213 Real t,ds,dv;
00214 Vector pos0,vel0;
00215 Grain *gr;
00216 FILE *fid;
00217
00218 t = pGrid->time + pGrid->dt;
00219
00220 pos0 = ParticleTroj(t);
00221 vel0 = ParticleVel(t);
00222
00223 for (p=0; p<pGrid->nparticle; p++)
00224 {
00225 if (pGrid->particle[p].pos == 1)
00226 {
00227 gr = &(pGrid->particle[p]);
00228 ds = sqrt(SQR(pos0.x1-gr->x1)+SQR(pos0.x2-gr->x2)+SQR(pos0.x3-gr->x3));
00229 dv = sqrt(SQR(vel0.x1-gr->v1)+SQR(vel0.x2-gr->v2)+SQR(vel0.x3-gr->v3));
00230
00231 fid = fopen(name,"a+");
00232 fprintf(fid,"%e %e %e\n", t, ds, dv);
00233 fclose(fid);
00234 }
00235 }
00236
00237 return;
00238 }
00239
00240
00241
00242
00243
00244
00245
00246 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00247 {
00248 return;
00249 }
00250
00251
00252
00253
00254
00255
00256 static Vector ParticleTroj(Real t)
00257 {
00258 Vector pos;
00259 Real L1,L2,L3;
00260
00261 pos.x1 = x1c+v01*tstop0[0]*(1.0-exp(-t/tstop0[0]));
00262 pos.x2 = x2c+v02*tstop0[0]*(1.0-exp(-t/tstop0[0]));
00263 pos.x3 = x3c+v03*tstop0[0]*(1.0-exp(-t/tstop0[0]));
00264
00265 L1 = x1max-x1min; L2 = x2max-x2min; L3 = x3max-x3min;
00266
00267 pos.x1 = pos.x1-floor((pos.x1-x1min)/L1)*L1;
00268 pos.x2 = pos.x2-floor((pos.x2-x2min)/L2)*L2;
00269 pos.x3 = pos.x3-floor((pos.x3-x3min)/L3)*L3;
00270
00271 return pos;
00272 }
00273
00274
00275
00276 static Vector ParticleVel(Real t)
00277 {
00278 Vector vel;
00279
00280 vel.x1 = v01*exp(-t/tstop0[0]);
00281 vel.x2 = v02*exp(-t/tstop0[0]);
00282 vel.x3 = v03*exp(-t/tstop0[0]);
00283
00284 return vel;
00285 }
00286
00287
00288
00289 static int ParticleLocator(Real x1, Real x2, Real x3)
00290 {
00291 if ((x1<x1upar) && (x1>=x1lpar) && (x2<x2upar) && (x2>=x2lpar) &&(x3<x3upar) && (x3>=x3lpar))
00292 return 1;
00293 else
00294 return 0;
00295 }