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