00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include <float.h>
00027 #include <math.h>
00028 #include <stdio.h>
00029 #include <stdlib.h>
00030 #include <string.h>
00031 #include <time.h>
00032 #include "defs.h"
00033 #include "athena.h"
00034 #include "globals.h"
00035 #include "prototypes.h"
00036 #include "particles/particle.h"
00037
00038 #ifndef SHEARING_BOX
00039 #error : The streaming2d problem requires shearing-box to be enabled.
00040 #endif
00041
00042 #ifndef PARTICLES
00043 #error : The streaming2d problem requires particles to be enabled.
00044 #endif
00045
00046 #ifndef ISOTHERMAL
00047 #error : The streaming2d problem requires isothermal equation of state.
00048 #endif
00049
00050
00051
00052 Real rho0, mratio, etavk, *epsilon, uxNSH, uyNSH, *wxNSH, *wyNSH;
00053 int ipert;
00054
00055 int Npar,Npar2;
00056
00057 long ntrack;
00058 long nlis;
00059 int mytype;
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069 double ran2(long int *idum);
00070
00071 void MultiNSH(int n, Real *tstop, Real *epsilon, Real etavk,
00072 Real *uxNSH, Real *uyNSH, Real *wxNSH, Real *wyNSH);
00073
00074 static Real hst_rho_Vx_dVy(const Grid *pG,const int i,const int j,const int k);
00075
00076 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3);
00077
00078 static int property_limit(const Grain *gr, const GrainAux *grsub);
00079 static int property_trace(const Grain *gr, const GrainAux *grsub);
00080 static int property_type(const Grain *gr, const GrainAux *grsub);
00081
00082 extern Real expr_dpar(const Grid *pG, const int i, const int j, const int k);
00083 extern Real expr_V1par(const Grid *pG, const int i, const int j, const int k);
00084 extern Real expr_V2par(const Grid *pG, const int i, const int j, const int k);
00085 extern Real expr_V3par(const Grid *pG, const int i, const int j, const int k);
00086 extern Real expr_V2(const Grid *pG, const int i, const int j, const int k);
00087
00088
00089
00090
00091
00092
00093 void problem(Grid *pGrid, Domain *pDomain)
00094 {
00095 int i,j,ip,jp,pt;
00096 long p;
00097 Real x1,x2,x3,t,x1l,x1u,x2l,x2u,x1p,x2p,x3p;
00098 Real x1min,x1max,x2min,x2max,Lx,Lz;
00099 Real tsmin,tsmax,tscrit,pwind,epsum;
00100 long int iseed = pGrid->my_id;
00101
00102 if (par_geti("grid","Nx2") == 1) {
00103 ath_error("[streaming2d]: streaming2D only works for 2D problem.\n");
00104 }
00105 if (par_geti("grid","Nx3") > 1) {
00106 ath_error("[streaming2d]: streaming2D only works for 2D problem.\n");
00107 }
00108
00109
00110
00111 x1min = pGrid->x1_0 + (pGrid->is + pGrid->idisp)*pGrid->dx1;
00112 x1max = pGrid->x1_0 + (pGrid->ie + pGrid->idisp + 1.0)*pGrid->dx1;
00113 Lx = x1max - x1min;
00114
00115 x2min = pGrid->x2_0 + (pGrid->js + pGrid->jdisp)*pGrid->dx2;
00116 x2max = pGrid->x2_0 + (pGrid->je + pGrid->jdisp + 1.0)*pGrid->dx2;
00117 Lz = x2max - x2min;
00118
00119
00120 rho0 = 1.0;
00121 Omega_0 = par_getd("problem","omega");
00122 qshear = par_getd_def("problem","qshear",1.5);
00123 ipert = par_geti_def("problem","ipert", 2);
00124 etavk = par_getd_def("problem","etavk",0.05);
00125
00126
00127 Npar = (int)(sqrt(par_geti("particle","parnumcell")));
00128 Npar2 = SQR(Npar);
00129
00130 pGrid->nparticle = Npar2*pGrid->Nx1*pGrid->Nx2;
00131 for (i=0; i<pGrid->partypes; i++)
00132 pGrid->grproperty[i].num = pGrid->nparticle;
00133 pGrid->nparticle = pGrid->partypes*pGrid->nparticle;
00134
00135 if (pGrid->nparticle+2 > pGrid->arrsize)
00136 particle_realloc(pGrid, pGrid->nparticle+2);
00137
00138
00139 tsmin = par_getd("problem","tsmin");
00140 tsmax = par_getd("problem","tsmax");
00141 tscrit= par_getd("problem","tscrit");
00142 if (par_geti("particle","tsmode") != 3)
00143 ath_error("[streaming2d]: This test works only for fixed stopping time!\n");
00144
00145 for (i=0; i<pGrid->partypes; i++) {
00146 tstop0[i] = tsmin*exp(i*log(tsmax/tsmin)/MAX(pGrid->partypes-1,1.0));
00147
00148
00149 if (tstop0[i] < tscrit) pGrid->grproperty[i].integrator = 3;
00150 }
00151
00152
00153 epsilon= (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00154 wxNSH = (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00155 wyNSH = (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00156
00157 #ifdef FEEDBACK
00158 mratio = par_getd_def("problem","mratio",0.0);
00159 pwind = par_getd_def("problem","pwind",0.0);
00160 if (mratio < 0.0)
00161 ath_error("[streaming2d]: mratio must be positive!\n");
00162
00163 epsum = 0.0;
00164 for (i=0; i<pGrid->partypes; i++)
00165 {
00166 epsilon[i] = pow(tstop0[i],pwind); epsum += epsilon[i];
00167 }
00168
00169 for (i=0; i<pGrid->partypes; i++)
00170 {
00171 epsilon[i] = mratio*epsilon[i]/epsum;
00172 pGrid->grproperty[i].m = rho0*epsilon[i]/Npar2;
00173 }
00174 #else
00175 mratio = 0.0;
00176 for (i=0; i<pGrid->partypes; i++)
00177 epsilon[i] = 0.0;
00178 #endif
00179
00180 etavk = etavk * Iso_csound;
00181
00182
00183 MultiNSH(pGrid->partypes, tstop0, epsilon, etavk,
00184 &uxNSH, &uyNSH, wxNSH, wyNSH);
00185
00186
00187 for (j=pGrid->js; j<=pGrid->je; j++) {
00188 for (i=pGrid->is; i<=pGrid->ie; i++) {
00189 cc_pos(pGrid,i,j,pGrid->ks,&x1,&x2,&x3);
00190
00191 pGrid->U[pGrid->ks][j][i].d = rho0;
00192
00193 if (ipert != 3) {
00194 pGrid->U[pGrid->ks][j][i].M1 = rho0 * uxNSH;
00195 pGrid->U[pGrid->ks][j][i].M3 = rho0 * uyNSH;
00196 } else {
00197 pGrid->U[pGrid->ks][j][i].M1 = 0.0;
00198 pGrid->U[pGrid->ks][j][i].M3 = 0.0;
00199 }
00200
00201 pGrid->U[pGrid->ks][j][i].M2 = 0.0;
00202 #ifndef FARGO
00203 pGrid->U[pGrid->ks][j][i].M3 -= qshear*rho0*Omega*x1;
00204 #endif
00205
00206 }}
00207
00208
00209 p = 0;
00210 x3p = pGrid->x3_0 + (pGrid->ks+pGrid->kdisp)*pGrid->dx3;
00211
00212 for (j=pGrid->js; j<=pGrid->je; j++)
00213 {
00214 x2l = pGrid->x2_0 + (j+pGrid->jdisp)*pGrid->dx2;
00215 x2u = pGrid->x2_0 + ((j+pGrid->jdisp)+1.0)*pGrid->dx2;
00216
00217 for (i=pGrid->is; i<=pGrid->ie; i++)
00218 {
00219 x1l = pGrid->x1_0 + (i + pGrid->idisp)*pGrid->dx1;
00220 x1u = pGrid->x1_0 + ((i + pGrid->idisp) + 1.0)*pGrid->dx1;
00221
00222 for (pt=0; pt<pGrid->partypes; pt++)
00223 {
00224 for (ip=0;ip<Npar;ip++)
00225 {
00226
00227 if (ipert == 0)
00228 x1p = x1l+pGrid->dx1/Npar*(ip+0.5);
00229
00230 for (jp=0;jp<Npar;jp++)
00231 {
00232 if (ipert == 0)
00233 x2p = x2l+pGrid->dx2/Npar*(jp+0.5);
00234 else if (ipert == 1)
00235 {
00236 x1p = x1l + pGrid->dx1*ran2(&iseed);
00237 x2p = x2l + pGrid->dx2*ran2(&iseed);
00238 }
00239 else
00240 {
00241 x1p = x1min + Lx*ran2(&iseed);
00242 x2p = x2min + Lz*ran2(&iseed);
00243 }
00244
00245 pGrid->particle[p].property = pt;
00246 pGrid->particle[p].x1 = x1p;
00247 pGrid->particle[p].x2 = x2p;
00248 pGrid->particle[p].x3 = x3p;
00249
00250 if (ipert != 3) {
00251 pGrid->particle[p].v1 = wxNSH[pt];
00252 pGrid->particle[p].v3 = wyNSH[pt];
00253 } else {
00254 pGrid->particle[p].v1 = 0.0;
00255 pGrid->particle[p].v3 = etavk;
00256 }
00257
00258 pGrid->particle[p].v2 = 0.0;
00259 #ifndef FARGO
00260 pGrid->particle[p].v3 -= qshear*Omega*x1p;
00261 #endif
00262 pGrid->particle[p].pos = 1;
00263 pGrid->particle[p].my_id = p;
00264 #ifdef MPI_PARALLEL
00265 pGrid->particle[p].init_id = pGrid->my_id;
00266 #endif
00267 p += 1;
00268 }
00269 }
00270 }
00271 }
00272 }
00273
00274
00275 StaticGravPot = ShearingBoxPot;
00276
00277 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00278
00279
00280
00281
00282 nlis = par_geti_def("problem","nlis",pGrid->Nx1*pGrid->Nx2);
00283
00284
00285 ntrack = par_geti_def("problem","ntrack",2000);
00286
00287 return;
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00301 {
00302 fwrite(&rho0, sizeof(Real),1,fp);
00303 fwrite(&mratio, sizeof(Real),1,fp); fwrite(epsilon, sizeof(Real),pG->partypes,fp);
00304 fwrite(&etavk, sizeof(Real),1,fp);
00305 fwrite(&uxNSH, sizeof(Real),1,fp);
00306 fwrite(&uyNSH, sizeof(Real),1,fp);
00307 fwrite(wxNSH, sizeof(Real),pG->partypes,fp);
00308 fwrite(wyNSH, sizeof(Real),pG->partypes,fp);
00309
00310 return;
00311 }
00312
00313 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00314 {
00315 Real singleintvl,tracetime,traceintvl;
00316
00317 StaticGravPot = ShearingBoxPot;
00318
00319 Omega_0 = par_getd("problem","omega");
00320 qshear = par_getd_def("problem","qshear",1.5);
00321 ipert = par_geti_def("problem","ipert", 1);
00322
00323 Npar = (int)(sqrt(par_geti("particle","parnumcell")));
00324 Npar2 = SQR(Npar);
00325 nlis = par_geti_def("problem","nlis",pG->Nx1*pG->Nx2);
00326 ntrack = par_geti_def("problem","ntrack",2000);
00327
00328
00329 epsilon= (Real*)calloc_1d_array(pG->partypes, sizeof(Real));
00330 wxNSH = (Real*)calloc_1d_array(pG->partypes, sizeof(Real));
00331 wyNSH = (Real*)calloc_1d_array(pG->partypes, sizeof(Real));
00332
00333 fread(&rho0, sizeof(Real),1,fp);
00334 fread(&mratio, sizeof(Real),1,fp);
00335 fread(epsilon,sizeof(Real),pG->partypes,fp);
00336 fread(&etavk, sizeof(Real),1,fp);
00337 fread(&uxNSH, sizeof(Real),1,fp);
00338 fread(&uyNSH, sizeof(Real),1,fp);
00339 fread(wxNSH, sizeof(Real),pG->partypes,fp);
00340 fread(wyNSH, sizeof(Real),pG->partypes,fp);
00341
00342 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00343
00344 return;
00345 }
00346
00347
00348
00349
00350 static Real expr_rhopardif(const Grid *pG,
00351 const int i, const int j, const int k)
00352 {
00353 Real x1,x2,x3;
00354 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00355 return pG->Coup[k][j][i].grid_d - rho0*mratio;
00356 }
00357
00358
00359
00360
00361
00362
00363
00364 static Real hst_rho_Vx_dVy(const Grid *pG,
00365 const int i, const int j, const int k)
00366 {
00367 Real x1,x2,x3;
00368 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00369 #ifdef FARGO
00370 return pG->U[k][j][i].M1*pG->U[k][j][i].M3/pG->U[k][j][i].d;
00371 #else
00372 return pG->U[k][j][i].M1*(pG->U[k][j][i].M3/pG->U[k][j][i].d
00373 + qshear*Omega_0*x1);
00374 #endif
00375 }
00376
00377 Gasfun_t get_usr_expr(const char *expr)
00378 {
00379 if(strcmp(expr,"difdpar")==0) return expr_rhopardif;
00380 return NULL;
00381 }
00382
00383 VGFunout_t get_usr_out_fun(const char *name){
00384 return NULL;
00385 }
00386
00387 #ifdef PARTICLES
00388 PropFun_t get_usr_par_prop(const char *name)
00389 {
00390 if (strcmp(name,"limit")==0) return property_limit;
00391 if (strcmp(name,"trace")==0) return property_trace;
00392 if (strcmp(name,"type")==0) return property_type;
00393 return NULL;
00394 }
00395
00396
00397
00398
00399 void gasvshift(const Real x1, const Real x2, const Real x3,
00400 Real *u1, Real *u2, Real *u3)
00401 {
00402 return;
00403 }
00404
00405 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00406 const Real v1, const Real v2, const Real v3)
00407 {
00408 ft->x1 -= 2.0*etavk*Omega_0;
00409 return;
00410 }
00411 #endif
00412
00413 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00414 {
00415 return;
00416 }
00417
00418
00419
00420
00421
00422
00423
00424 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00425 {
00426 free(epsilon);
00427 free(wxNSH); free(wyNSH);
00428 return;
00429 }
00430
00431
00432
00433
00434
00435 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00436 {
00437 Real phi=0.0;
00438 #ifndef FARGO
00439 phi -= qshear*SQR(Omega_0*x1);
00440 #endif
00441 return phi;
00442 }
00443
00444
00445
00446 static int property_limit(const Grain *gr, const GrainAux *grsub)
00447 {
00448 if ((gr->pos == 1) && (gr->my_id<nlis))
00449 return 1;
00450 else
00451 return 0;
00452 }
00453
00454
00455
00456 static int property_trace(const Grain *gr, const GrainAux *grsub)
00457 {
00458 if ((gr->pos == 1) && (gr->my_id<ntrack))
00459 return 1;
00460 else
00461 return 0;
00462 }
00463
00464
00465
00466 static int property_type(const Grain *gr, const GrainAux *grsub)
00467 {
00468 if ((gr->pos == 1) && (gr->property == mytype))
00469 return 1;
00470 else
00471 return 0;
00472 }
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 void MultiNSH(int n, Real *tstop, Real *mratio, Real etavk,
00485 Real *uxNSH, Real *uyNSH, Real *wxNSH, Real *wyNSH)
00486 {
00487 int i,j;
00488 Real *Lambda1,**Lam1GamP1, **A, **B, **Tmp;
00489
00490 Lambda1 = (Real*)calloc_1d_array(n, sizeof(Real));
00491 Lam1GamP1=(Real**)calloc_2d_array(n, n, sizeof(Real));
00492 A = (Real**)calloc_2d_array(n, n, sizeof(Real));
00493 B = (Real**)calloc_2d_array(n, n, sizeof(Real));
00494 Tmp = (Real**)calloc_2d_array(n, n, sizeof(Real));
00495
00496
00497 for (i=0; i<n; i++){
00498 for (j=0; j<n; j++)
00499 Lam1GamP1[i][j] = mratio[j];
00500 Lam1GamP1[i][i] += 1.0;
00501 Lambda1[i] = 1.0/(tstop[i]+1.0e-16);
00502 for (j=0; j<n; j++)
00503 Lam1GamP1[i][j] *= Lambda1[i];
00504 }
00505
00506
00507 MatrixMult(Lam1GamP1, Lam1GamP1, n,n,n, Tmp);
00508 for (i=0; i<n; i++) Tmp[i][i] += 1.0;
00509 InverseMatrix(Tmp, n, B);
00510 for (i=0; i<n; i++)
00511 for (j=0; j<n; j++)
00512 B[i][j] *= Lambda1[j];
00513 MatrixMult(Lam1GamP1, B, n,n,n, A);
00514
00515
00516 *uxNSH = 0.0; *uyNSH = 0.0;
00517 for (i=0; i<n; i++){
00518 wxNSH[i] = 0.0;
00519 wyNSH[i] = 0.0;
00520 for (j=0; j<n; j++){
00521 wxNSH[i] -= B[i][j];
00522 wyNSH[i] -= A[i][j];
00523 }
00524 wxNSH[i] *= 2.0*etavk;
00525 wyNSH[i] *= etavk;
00526 *uxNSH -= mratio[i]*wxNSH[i];
00527 *uyNSH -= mratio[i]*wyNSH[i];
00528 wyNSH[i] += etavk;
00529 }
00530
00531 free(Lambda1);
00532 free_2d_array(A); free_2d_array(B);
00533 free_2d_array(Lam1GamP1); free_2d_array(Tmp);
00534
00535 return;
00536 }
00537
00538
00539
00540
00541 #define IM1 2147483563
00542 #define IM2 2147483399
00543 #define AM (1.0/IM1)
00544 #define IMM1 (IM1-1)
00545 #define IA1 40014
00546 #define IA2 40692
00547 #define IQ1 53668
00548 #define IQ2 52774
00549 #define IR1 12211
00550 #define IR2 3791
00551 #define NTAB 32
00552 #define NDIV (1+IMM1/NTAB)
00553 #define RNMX (1.0-DBL_EPSILON)
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568 double ran2(long int *idum)
00569 {
00570 int j;
00571 long int k;
00572 static long int idum2=123456789;
00573 static long int iy=0;
00574 static long int iv[NTAB];
00575 double temp;
00576
00577 if (*idum <= 0) {
00578 if (-(*idum) < 1) *idum=1;
00579 else *idum = -(*idum);
00580 idum2=(*idum);
00581 for (j=NTAB+7;j>=0;j--) {
00582 k=(*idum)/IQ1;
00583 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00584 if (*idum < 0) *idum += IM1;
00585 if (j < NTAB) iv[j] = *idum;
00586 }
00587 iy=iv[0];
00588 }
00589 k=(*idum)/IQ1;
00590 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00591 if (*idum < 0) *idum += IM1;
00592 k=idum2/IQ2;
00593 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00594 if (idum2 < 0) idum2 += IM2;
00595 j=(int)(iy/NDIV);
00596 iy=iv[j]-idum2;
00597 iv[j] = *idum;
00598 if (iy < 1) iy += IMM1;
00599 if ((temp=AM*iy) > RNMX) return RNMX;
00600 else return temp;
00601 }
00602
00603 #undef IM1
00604 #undef IM2
00605 #undef AM
00606 #undef IMM1
00607 #undef IA1
00608 #undef IA2
00609 #undef IQ1
00610 #undef IQ2
00611 #undef IR1
00612 #undef IR2
00613 #undef NTAB
00614 #undef NDIV
00615 #undef RNMX