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 feedback problem requires particles to be enabled.
00029 #endif
00030
00031
00032
00033
00034
00035
00036 static Vector ParticleTroj(Real t);
00037 static Vector ParticleVel(Real t);
00038
00039
00040 Real mratio,wx,wy,wz;
00041 Real x1min,x1max,x2min,x2max;
00042 Real x0[3];
00043 int Npar;
00044 long idlab;
00045 int cpuid;
00046 char name[50];
00047
00048
00049
00050
00051
00052
00053 void problem(Grid *pGrid, Domain *pDomain)
00054 {
00055 int i,j,ip,jp;
00056 long p,lab;
00057 Real x1,x2,x3,x1l,x1u,x2l,x2u,x1p,x2p,x3p,cosphi,sinphi;
00058 Real vdif_p,vdif_v,w,u,ux,uy,uz;
00059
00060 if (par_geti("grid","Nx2") == 1) {
00061 ath_error("[par_collision]: this test only works for 2D problem.\n");
00062 }
00063 if (par_geti("grid","Nx3") > 1) {
00064 ath_error("[par_collision]: this test only works for 2D problem.\n");
00065 }
00066
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
00073 vdif_p = par_getd("problem","vdif_p");
00074 vdif_v = par_getd("problem","vdif_v");
00075 cosphi = par_getd("problem","oblique");
00076 sinphi = sqrt(1.0-SQR(cosphi));
00077
00078
00079
00080
00081 if (par_geti("particle","partypes") != 1)
00082 ath_error("[par_collision]: This test only allows ONE particle species!\n");
00083
00084 Npar = (int)(sqrt(par_geti("particle","parnumcell")));
00085 pGrid->nparticle = Npar*Npar*pGrid->Nx1*pGrid->Nx2;
00086 pGrid->grproperty[0].num = pGrid->nparticle;
00087 if (pGrid->nparticle+2 > pGrid->arrsize)
00088 particle_realloc(pGrid, pGrid->nparticle+2);
00089
00090
00091 tstop0[0] = par_getd("problem","tstop");
00092 if (par_geti("particle","tsmode") != 3)
00093 ath_error("[par_collision]: This test only allow fixed stopping time!\n");
00094
00095
00096 #ifdef FEEDBACK
00097 mratio = par_getd_def("problem","mratio",0.0);
00098 if (mratio < 0.0)
00099 ath_error("[par_collision]: mratio must be positive!\n");
00100 pGrid->grproperty[0].m = mratio/SQR(Npar);
00101 #else
00102 mratio = 0.0;
00103 #endif
00104
00105 w = vdif_p/(1.0+mratio);
00106 u = -mratio*w;
00107
00108 wx = w * cosphi;
00109 wy = w * sinphi;
00110 ux = u * cosphi;
00111 uy = u * sinphi;
00112 wz = vdif_v/(1.0+mratio);
00113 uz = -mratio*wz;
00114
00115
00116
00117 for (j=pGrid->js; j<=pGrid->je; j++) {
00118 for (i=pGrid->is; i<=pGrid->ie; i++) {
00119 cc_pos(pGrid,i,j,pGrid->ks,&x1,&x2,&x3);
00120 pGrid->U[pGrid->ks][j][i].d = 1.0;
00121
00122 pGrid->U[pGrid->ks][j][i].M1 = ux;
00123 pGrid->U[pGrid->ks][j][i].M2 = uy;
00124 pGrid->U[pGrid->ks][j][i].M3 = uz;
00125 }}
00126
00127
00128 p = 0;
00129 x3p = pGrid->x3_0 + (pGrid->ks+pGrid->kdisp)*pGrid->dx3;
00130
00131 for (j=pGrid->js; j<=pGrid->je; j++)
00132 {
00133 x2l = pGrid->x2_0 + (j+pGrid->jdisp)*pGrid->dx2;
00134 x2u = pGrid->x2_0 + ((j+pGrid->jdisp)+1.0)*pGrid->dx2;
00135
00136 for (i=pGrid->is; i<=pGrid->ie; i++)
00137 {
00138 x1l = pGrid->x1_0 + (i + pGrid->idisp)*pGrid->dx1;
00139 x1u = pGrid->x1_0 + ((i + pGrid->idisp) + 1.0)*pGrid->dx1;
00140
00141 for (ip=0;ip<Npar;ip++)
00142 {
00143 x1p = x1l+(x1u-x1l)/Npar*(ip+0.5);
00144
00145 for (jp=0;jp<Npar;jp++)
00146 {
00147 x2p = x2l+(x2u-x2l)/Npar*(jp+0.5);
00148
00149 pGrid->particle[p].property = 0;
00150 pGrid->particle[p].x1 = x1p;
00151 pGrid->particle[p].x2 = x2p;
00152 pGrid->particle[p].x3 = x3p;
00153
00154 pGrid->particle[p].v1 = wx;
00155 pGrid->particle[p].v2 = wy;
00156 pGrid->particle[p].v3 = wz;
00157
00158 pGrid->particle[p].pos = 1;
00159 pGrid->particle[p].my_id = p;
00160 #ifdef MPI_PARALLEL
00161 pGrid->particle[p].init_id = pGrid->my_id;
00162 #endif
00163 p += 1;
00164 }
00165 }
00166 }
00167 }
00168
00169
00170 if (pGrid->my_id == 0) {
00171 lab = 0;
00172 x0[0] = pGrid->particle[lab].x1;
00173 x0[1] = pGrid->particle[lab].x2;
00174 x0[2] = pGrid->particle[lab].x3;
00175 idlab = pGrid->particle[lab].my_id;
00176 cpuid = 0;
00177 }
00178
00179 if (pGrid->my_id == 0) {
00180
00181 sprintf(name, "%s_partroj.dat", pGrid->outfilename);
00182 FILE *fid = fopen(name,"w");
00183 fclose(fid);
00184 #ifdef MPI_PARALLEL
00185 sprintf(name, "../%s_partroj.dat", pGrid->outfilename);
00186 #else
00187 sprintf(name, "%s_partroj.dat", pGrid->outfilename);
00188 #endif
00189 }
00190
00191 #ifdef MPI_PARALLEL
00192 MPI_Bcast(name,50,MPI_CHAR,0,MPI_COMM_WORLD);
00193 MPI_Bcast(&x0,3,MPI_DOUBLE,0,MPI_COMM_WORLD);
00194 MPI_Bcast(&idlab,1,MPI_LONG,0,MPI_COMM_WORLD);
00195 MPI_Bcast(&cpuid,1,MPI_INT,0,MPI_COMM_WORLD);
00196 #endif
00197
00198 return;
00199 }
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00212 {
00213 fwrite(&wx, sizeof(Real),1,fp);
00214 fwrite(&wy, sizeof(Real),1,fp);
00215 fwrite(&wz, sizeof(Real),1,fp);
00216 fwrite(x0, sizeof(Real),3,fp);
00217 fwrite(&idlab, sizeof(long),1,fp);
00218 fwrite(&cpuid, sizeof(int),1,fp);
00219 return;
00220 }
00221
00222 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00223 {
00224 fread(&wx, sizeof(Real),1,fp);
00225 fread(&wy, sizeof(Real),1,fp);
00226 fread(&wz, sizeof(Real),1,fp);
00227 fread(x0, sizeof(Real),3,fp);
00228 fread(&idlab, sizeof(long),1,fp);
00229 fread(&cpuid, sizeof(int),1,fp);
00230
00231 x1min = par_getd("grid","x1min");
00232 x1max = par_getd("grid","x1max");
00233 x2min = par_getd("grid","x2min");
00234 x2max = par_getd("grid","x2max");
00235 Npar = (int)(sqrt(par_geti("particle","parnumcell")));
00236 #ifdef FEEDBACK
00237 mratio = par_getd_def("problem","mratio",0.0);
00238 #else
00239 mratio = 0.0;
00240 #endif
00241
00242 sprintf(name, "%s_partroj.dat", pG->outfilename);
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 void gasvshift(const Real x1, const Real x2, const Real x3,
00262 Real *u1, Real *u2, Real *u3)
00263 {
00264 return;
00265 }
00266
00267 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00268 const Real v1, const Real v2, const Real v3)
00269 {
00270 return;
00271 }
00272 #endif
00273
00274
00275 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00276 {
00277 long p, lab;
00278 Grain *gr;
00279 Real Mg1,Mg2,Mg3,mp,t,s,ds,v,dv,Mtot;
00280 Vector pos0, vel0;
00281 FILE *fid;
00282
00283 p = 0; lab = -1;
00284 while (p<pGrid->nparticle) {
00285 #ifdef MPI_PARALLEL
00286 if ((pGrid->particle[p].my_id == idlab) &&
00287 (pGrid->particle[p].init_id == cpuid))
00288 #else
00289 if (pGrid->particle[p].my_id == idlab)
00290 #endif
00291 {
00292 lab = p;
00293 break;
00294 }
00295 p++;
00296 }
00297
00298 if (lab >= 0)
00299 {
00300
00301
00302 t = pGrid->time + pGrid->dt;
00303 pos0 = ParticleTroj(t);
00304 vel0 = ParticleVel(t);
00305 Mg1 = pGrid->U[pGrid->ks][pGrid->js][pGrid->is].M1;
00306 Mg2 = pGrid->U[pGrid->ks][pGrid->js][pGrid->is].M2;
00307 Mg3 = pGrid->U[pGrid->ks][pGrid->js][pGrid->is].M3;
00308 gr = &(pGrid->particle[lab]);
00309 #ifdef FEEDBACK
00310 mp = pGrid->grproperty[0].m*SQR(Npar);
00311 #else
00312 mp = SQR(Npar);
00313 #endif
00314
00315
00316 s = sqrt(SQR(gr->x1-x0[0])+SQR(gr->x2-x0[1])+SQR(gr->x3-x0[2]));
00317 ds = sqrt(SQR(gr->x1-pos0.x1)+SQR(gr->x2-pos0.x2)+SQR(gr->x3-pos0.x3));
00318 v = sqrt(SQR(gr->v1)+SQR(gr->v2)+SQR(gr->v3));
00319 dv = sqrt(SQR(gr->v1-vel0.x1)+SQR(gr->v2-vel0.x2)+SQR(gr->v3-vel0.x3));
00320 Mtot = sqrt(SQR(mp*gr->v1+Mg1)+SQR(mp*gr->v2+Mg2)+SQR(mp*gr->v3+Mg3));
00321
00322
00323 fid = fopen(name,"a+");
00324 fprintf(fid,"%e %e %e %e %e %e\n",
00325 t, s, ds, v, dv, Mtot);
00326 fclose(fid);
00327 }
00328
00329 return;
00330 }
00331
00332
00333
00334
00335
00336
00337
00338 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00339 {
00340 return;
00341 }
00342
00343
00344
00345
00346
00347
00348 static Vector ParticleTroj(Real t)
00349 {
00350 Vector pos;
00351 Real L1,L2,L3;
00352 Real factor=tstop0[0]/(1.0+mratio)*(1.0-exp(-(1.0+mratio)*t/tstop0[0]));
00353
00354 pos.x1 = x0[0] + wx*factor;
00355 pos.x2 = x0[1] + wy*factor;
00356 pos.x3 = x0[2];
00357
00358 L1 = x1max-x1min; L2 = x2max-x2min;
00359
00360 pos.x1 = pos.x1-floor((pos.x1-x1min)/L1)*L1;
00361 pos.x2 = pos.x2-floor((pos.x2-x2min)/L2)*L2;
00362
00363 return pos;
00364 }
00365
00366
00367
00368 static Vector ParticleVel(Real t)
00369 {
00370 Vector vel;
00371 Real factor = exp(-(1.0+mratio)*t/tstop0[0]);
00372
00373 vel.x1 = wx*factor;
00374 vel.x2 = wy*factor;
00375 vel.x3 = wz*factor;
00376
00377 return vel;
00378 }