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
00027 #include <float.h>
00028 #include <math.h>
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 #include <time.h>
00033 #include "defs.h"
00034 #include "athena.h"
00035 #include "globals.h"
00036 #include "prototypes.h"
00037 #include "particles/particle.h"
00038
00039 #ifndef SHEARING_BOX
00040 #error : The streaming3d problem requires shearing-box to be enabled.
00041 #endif
00042
00043 #ifndef PARTICLES
00044 #error : The streaming3d problem requires particles to be enabled.
00045 #endif
00046
00047 #ifndef ISOTHERMAL
00048 #error : The streaming3d problem requires isothermal equation of state.
00049 #endif
00050
00051
00052
00053 Real rho0, mratio, etavk, uxNSH, uyNSH, wxNSH, wyNSH;
00054
00055 int Npar,Npar3,downsamp,Nx;
00056
00057 Real Reux,Imux,Reuy,Imuy,Reuz,Imuz,Rewx,Imwx,Rewy,Imwy,Rewz,Imwz,Rerho,Imrho,omg,s;
00058
00059 Real amp, kx, kz;
00060 int ipert, nwave;
00061
00062 char name[50];
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 double ran2(long int *idum);
00073 void OutputModeAmplitude(Grid *pGrid, Domain *pDomain, Output *pOut);
00074 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3);
00075 static Real pert_even(Real fR, Real fI, Real x, Real z, Real t);
00076 static Real pert_odd(Real fR, Real fI, Real x, Real z, Real t);
00077 static int property_mybin(const Grain *gr, const GrainAux *grsub);
00078 extern Real expr_V3(const Grid *pG, const int i, const int j, const int k);
00079 extern Real expr_V1par(const Grid *pG, const int i, const int j, const int k);
00080 extern Real expr_V2par(const Grid *pG, const int i, const int j, const int k);
00081 extern Real expr_V3par(const Grid *pG, const int i, const int j, const int k);
00082
00083
00084
00085
00086
00087
00088 void problem(Grid *pGrid, Domain *pDomain)
00089 {
00090 int i,j,k,ip,jp,kp,interp;
00091 long p;
00092 Real x1,x2,x3,t,x1l,x1u,x2l,x2u,x3l,x3u,x1p,x2p,x3p,paramp,factor2,reduct;
00093 Real x1min,x1max,x2min,x2max,x3min,x3max,Lx,Ly,Lz;
00094 Real rhog,cs,u1,u2,u3,w1,w2,w3,Kxn,Kzn,denorm1,interval;
00095 long int iseed = -1;
00096
00097 if (par_geti("grid","Nx3") == 1) {
00098 ath_error("[streaming3d]: streaming3D only works for 3D problem.\n");
00099 }
00100
00101
00102 x1min = par_getd("grid","x1min");
00103 x1max = par_getd("grid","x1max");
00104 Lx = x1max - x1min;
00105
00106 x2min = par_getd("grid","x2min");
00107 x2max = par_getd("grid","x2max");
00108 Ly = x2max - x2min;
00109
00110 x3min = par_getd("grid","x3min");
00111 x3max = par_getd("grid","x3max");
00112 Lz = x3max - x3min;
00113
00114 Nx = par_getd("grid","Nx1");
00115
00116
00117 rho0 = 1.0;
00118 Omega_0 = par_getd("problem","omega");
00119 qshear = par_getd_def("problem","qshear",1.5);
00120 amp = par_getd_def("problem","amp",0.0);
00121 ipert = par_geti_def("problem","ipert", 1);
00122 etavk = par_getd_def("problem","etavk",0.05);
00123
00124
00125 if (par_geti("particle","partypes") != 1)
00126 ath_error("[streaming3d]: This test only allows ONE particle species!\n");
00127
00128 Npar = (int)(pow(par_geti("particle","parnumcell"),1.0/3.0));
00129 Npar3 = Npar*SQR(Npar);
00130
00131 pGrid->nparticle = Npar3*pGrid->Nx1*pGrid->Nx2*pGrid->Nx3;
00132 pGrid->grproperty[0].num = pGrid->nparticle;
00133
00134 if (pGrid->nparticle+2 > pGrid->arrsize)
00135 particle_realloc(pGrid, pGrid->nparticle+2);
00136
00137
00138
00139
00140 downsamp = par_geti_def("problem","downsamp",Npar3);
00141
00142
00143 tstop0[0] = par_getd("problem","tstop");
00144 if (par_geti("particle","tsmode") != 3)
00145 ath_error("[streaming3d]: This test works only for fixed stopping time!\n");
00146
00147
00148 #ifdef FEEDBACK
00149 mratio = par_getd_def("problem","mratio",0.0);
00150 if (mratio < 0.0)
00151 ath_error("[streaming3d]: mratio must be positive!\n");
00152 #else
00153 mratio = 0.0;
00154 if ((ipert == 1) || (ipert == 2))
00155 ath_error("[streaming3d]: linear growth test requires FEEDBACK to be\
00156 turned on!\n");
00157 #endif
00158 amp = amp*mratio;
00159
00160 if (ipert == 1) {
00161 Kxn = 30.0;
00162 Kzn = 30.0;
00163 Reux = -0.1691398*amp;
00164 Imux = 0.0361553*amp;
00165 Reuy = 0.1336704*amp;
00166 Imuy = 0.0591695*amp;
00167 Reuz = 0.1691389*amp;
00168 Imuz = -0.0361555*amp;
00169 Rerho = 0.0000224*amp;
00170 Imrho = 0.0000212*amp;
00171 Rewx = -0.1398623*amp;
00172 Imwx = 0.0372951*amp;
00173 Rewy = 0.1305628*amp;
00174 Imwy = 0.0640574*amp;
00175 Rewz = 0.1639549*amp;
00176 Imwz = -0.0233277*amp;
00177 omg = -0.3480127*Omega_0;
00178 s = 0.4190204*Omega_0;
00179 mratio= 3.0;
00180 tstop0[0]= 0.1/Omega_0;
00181 etavk = 0.05;
00182 }
00183
00184 if (ipert == 2) {
00185 Kxn = 6.0;
00186 Kzn = 6.0;
00187 Reux = -0.0174121*amp;
00188 Imux = -0.2770347*amp;
00189 Reuy = 0.2767976*amp;
00190 Imuy = -0.0187568*amp;
00191 Reuz = 0.0174130*amp;
00192 Imuz = 0.2770423*amp;
00193 Rerho = -0.0000067*amp;
00194 Imrho = -0.0000691*amp;
00195 Rewx = 0.0462916*amp;
00196 Imwx = -0.2743072*amp;
00197 Rewy = 0.2739304*amp;
00198 Imwy = 0.0039293*amp;
00199 Rewz = 0.0083263*amp;
00200 Imwz = 0.2768866*amp;
00201 omg = 0.4998786*Omega_0;
00202 s = 0.0154764*Omega_0;
00203 mratio= 0.2;
00204 tstop0[0]= 0.1/Omega_0;
00205 etavk = 0.05;
00206 }
00207
00208 #ifdef FEEDBACK
00209 pGrid->grproperty[0].m = rho0*mratio/Npar3;
00210 #endif
00211
00212
00213 if ((ipert == 1) || (ipert == 2))
00214 {
00215 if (Lx != Lz)
00216 ath_error("[streaming3d]: Linear test requires Lx1=Lx2!\n");
00217 nwave = par_geti_def("problem","nwave",1);
00218 kx = 2.0*PI/Lx*nwave;
00219 kz = 2.0*PI/Lz*nwave;
00220
00221 cs = Kxn/kx/etavk*Omega_0;
00222 if (cs != Iso_csound)
00223 ath_pout(0, "[streaming3d]: Iso_csound=%f is modified to cs=%f.\n",
00224 Iso_csound, cs);
00225 Iso_csound = cs;
00226 Iso_csound2 = SQR(Iso_csound);
00227
00228 interp = par_geti("particle","interp");
00229 if (interp == 3) {
00230 paramp = amp*kx*pGrid->dx1/sin(kx*pGrid->dx1);
00231 factor2 = 0.5*tan(kx*pGrid->dx1)/(kx*pGrid->dx1);
00232 reduct = 1.0;
00233 }
00234 else if (interp == 2) {
00235 paramp = 4.0*amp*kx*pGrid->dx1/sin(kx*pGrid->dx1)/
00236 (3.0+cos(kz*pGrid->dx2));
00237 factor2 = 0.5*tan(kx*pGrid->dx1)/(kx*pGrid->dx1);
00238 reduct = 1.0/(1.0-0.25*SQR(kx*pGrid->dx1)); reduct=1.0;
00239 }
00240 else {
00241 paramp = amp;
00242 factor2 = 0.5;
00243 reduct = 1.0;
00244 }
00245 }
00246 etavk = etavk * Iso_csound;
00247
00248
00249 denorm1 = 1.0/(SQR(1.0+mratio)+SQR(tstop0[0]*Omega_0));
00250
00251 wxNSH = -2.0*tstop0[0]*Omega_0*denorm1*etavk;
00252 wyNSH = -(1.0+mratio)*denorm1*etavk;
00253
00254 uxNSH = -mratio*wxNSH;
00255 uyNSH = -mratio*wyNSH;
00256
00257 wyNSH += etavk;
00258
00259 ath_pout(0,"etavk=%f, Iso_csound=%f\n",etavk,Iso_csound);
00260
00261
00262 t = 0.0;
00263
00264 for (k=pGrid->ks; k<=pGrid->ke; k++) {
00265 for (j=pGrid->js; j<=pGrid->je; j++) {
00266 for (i=pGrid->is; i<=pGrid->ie; i++) {
00267 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00268
00269 if ((ipert == 1) || (ipert == 2)) {
00270 rhog = rho0 * (1.0+reduct*pert_even(Rerho,Imrho,x1,x3,t));
00271 u1 = etavk * reduct * pert_even(Reux,Imux,x1,x3,t);
00272 u3 = etavk * reduct * pert_odd (Reuz,Imuz,x1,x3,t);
00273 u2 = etavk * reduct * pert_even(Reuy,Imuy,x1,x3,t);
00274 } else {
00275 rhog = rho0; u1 = u2 = u3 = 0.0;
00276 }
00277
00278 pGrid->U[k][j][i].d = rhog;
00279
00280 pGrid->U[k][j][i].M1 = rhog * (uxNSH+u1);
00281 pGrid->U[k][j][i].M2 = rhog * (uyNSH+u2);
00282
00283 pGrid->U[k][j][i].M3 = rhog * u3;
00284 #ifndef FARGO
00285 pGrid->U[k][j][i].M2 -= qshear*rhog*Omega*x1;
00286 #endif
00287
00288 }}}
00289
00290
00291 p = 0;
00292 Lx = pGrid->Nx1*pGrid->dx1;
00293 x1min = pGrid->x1_0 + (pGrid->is + pGrid->idisp)*pGrid->dx1;
00294
00295 Ly = pGrid->Nx2*pGrid->dx2;
00296 x2min = pGrid->x2_0 + (pGrid->js + pGrid->jdisp)*pGrid->dx2;
00297
00298 Lz = pGrid->Nx3*pGrid->dx3;
00299 x3min = pGrid->x3_0 + (pGrid->ks + pGrid->kdisp)*pGrid->dx3;
00300
00301 for (k=pGrid->ks; k<=pGrid->ke; k++)
00302 {
00303 x3l = pGrid->x3_0 + (k+pGrid->kdisp)*pGrid->dx3;
00304 x3u = pGrid->x3_0 + ((k+pGrid->kdisp)+1.0)*pGrid->dx3;
00305
00306 for (j=pGrid->js; j<=pGrid->je; j++)
00307 {
00308 x2l = pGrid->x2_0 + (j+pGrid->jdisp)*pGrid->dx2;
00309 x2u = pGrid->x2_0 + ((j+pGrid->jdisp)+1.0)*pGrid->dx2;
00310
00311 for (i=pGrid->is; i<=pGrid->ie; i++)
00312 {
00313 x1l = pGrid->x1_0 + (i + pGrid->idisp)*pGrid->dx1;
00314 x1u = pGrid->x1_0 + ((i + pGrid->idisp) + 1.0)*pGrid->dx1;
00315
00316 for (ip=0;ip<Npar;ip++)
00317 {
00318
00319 if (ipert != 3)
00320 x1p = x1l+pGrid->dx1/Npar*(ip+0.5);
00321
00322
00323
00324 for (jp=0;jp<Npar;jp++)
00325 {
00326 if (ipert != 3)
00327 x2p = x2l+pGrid->dx2/Npar*(jp+0.5);
00328
00329
00330
00331 for (kp=0;kp<Npar;kp++)
00332 {
00333 if (ipert == 3){
00334
00335 x1p = x1min + Lx*ran2(&iseed);
00336 x2p = x2min + Ly*ran2(&iseed);
00337 x3p = x3min + Lz*ran2(&iseed);
00338 }
00339 else
00340 x3p = x3l+pGrid->dx3/Npar*(kp+0.5);
00341
00342 pGrid->particle[p].property = 0;
00343 pGrid->particle[p].x1 = x1p;
00344 pGrid->particle[p].x2 = x2p;
00345 pGrid->particle[p].x3 = x3p;
00346
00347 if ((ipert == 1) || (ipert == 2)) {
00348 pGrid->particle[p].x1 += paramp*cos(kz*x3p)*(-sin(kx*x1p)
00349 +factor2*paramp*sin(2.0*kx*x1p))/kx;
00350
00351
00352 w1 = etavk * pert_even(Rewx,Imwx,pGrid->particle[p].x1,x3p,t);
00353 w3 = etavk * pert_odd (Rewz,Imwz,pGrid->particle[p].x1,x3p,t);
00354 w2 = etavk * pert_even(Rewy,Imwy,pGrid->particle[p].x1,x3p,t);
00355 } else {
00356 w1 = w2 = w3 = 0.0;
00357 }
00358
00359 pGrid->particle[p].v1 = wxNSH+w1;
00360 pGrid->particle[p].v2 = wyNSH+w2;
00361
00362 pGrid->particle[p].v3 = w3;
00363 #ifndef FARGO
00364 pGrid->particle[p].v2 -= qshear*Omega*x1p;
00365 #endif
00366 pGrid->particle[p].pos = 1;
00367 pGrid->particle[p].my_id = p;
00368 #ifdef MPI_PARALLEL
00369 pGrid->particle[p].init_id = pGrid->my_id;
00370 #endif
00371 p += 1;
00372 }
00373 }
00374 }
00375 }
00376 }
00377 }
00378
00379
00380 StaticGravPot = ShearingBoxPot;
00381
00382 if (pGrid->my_id == 0) {
00383
00384 sprintf(name, "%s_%d_%d.dat", pGrid->outfilename,Nx,ipert);
00385 FILE *fid = fopen(name,"w");
00386 fclose(fid);
00387 #ifdef MPI_PARALLEL
00388 sprintf(name, "../%s_%d_%d.dat", pGrid->outfilename,Nx,ipert);
00389 #endif
00390 }
00391
00392
00393 interval = par_getd_def("problem","interval",0.01/Omega_0);
00394 data_output_enroll(pGrid->time,interval,0,OutputModeAmplitude,
00395 NULL,NULL,0,0.0,0.0,0,0,1,property_all);
00396
00397 return;
00398 }
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00411 {
00412 fwrite(&rho0, sizeof(Real),1,fp); fwrite(&mratio, sizeof(Real),1,fp);
00413 fwrite(&etavk, sizeof(Real),1,fp);
00414 fwrite(&uxNSH, sizeof(Real),1,fp); fwrite(&uyNSH, sizeof(Real),1,fp);
00415 fwrite(&wxNSH, sizeof(Real),1,fp); fwrite(&wyNSH, sizeof(Real),1,fp);
00416 if ((ipert==1) || (ipert==2)) {
00417 fwrite(&Reux, sizeof(Real),1,fp); fwrite(&Imux, sizeof(Real),1,fp);
00418 fwrite(&Reuy, sizeof(Real),1,fp); fwrite(&Imuy, sizeof(Real),1,fp);
00419 fwrite(&Reuz, sizeof(Real),1,fp); fwrite(&Imuz, sizeof(Real),1,fp);
00420 fwrite(&Rewx, sizeof(Real),1,fp); fwrite(&Imwx, sizeof(Real),1,fp);
00421 fwrite(&Rewy, sizeof(Real),1,fp); fwrite(&Imwy, sizeof(Real),1,fp);
00422 fwrite(&Rewz, sizeof(Real),1,fp); fwrite(&Imwz, sizeof(Real),1,fp);
00423 fwrite(&Rerho, sizeof(Real),1,fp);fwrite(&Imrho, sizeof(Real),1,fp);
00424 fwrite(&omg, sizeof(Real),1,fp); fwrite(&s, sizeof(Real),1,fp);
00425 fwrite(&Iso_csound, sizeof(Real),1,fp);
00426 fwrite(&kx, sizeof(Real),1,fp); fwrite(&kz, sizeof(Real),1,fp);
00427 }
00428 return;
00429 }
00430
00431 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00432 {
00433 Real interval;
00434
00435 StaticGravPot = ShearingBoxPot;
00436
00437 Omega_0 = par_getd("problem","omega");
00438 qshear = par_getd_def("problem","qshear",1.5);
00439 amp = par_getd_def("problem","amp",0.0);
00440 ipert = par_geti_def("problem","ipert", 1);
00441
00442 Nx = par_getd("grid","Nx1");
00443 Npar = (int)(sqrt(par_geti("particle","parnumcell")));
00444 Npar3 = SQR(Npar)*Npar;
00445 downsamp = par_geti_def("problem","downsamp",Npar3);
00446
00447 fread(&rho0, sizeof(Real),1,fp); fread(&mratio, sizeof(Real),1,fp);
00448 fread(&etavk, sizeof(Real),1,fp);
00449 fread(&uxNSH, sizeof(Real),1,fp); fread(&uyNSH, sizeof(Real),1,fp);
00450 fread(&wxNSH, sizeof(Real),1,fp); fread(&wyNSH, sizeof(Real),1,fp);
00451 if ((ipert==1) || (ipert==2)) {
00452 fread(&Reux, sizeof(Real),1,fp); fread(&Imux, sizeof(Real),1,fp);
00453 fread(&Reuy, sizeof(Real),1,fp); fread(&Imuy, sizeof(Real),1,fp);
00454 fread(&Reuz, sizeof(Real),1,fp); fread(&Imuz, sizeof(Real),1,fp);
00455 fread(&Rewx, sizeof(Real),1,fp); fread(&Imwx, sizeof(Real),1,fp);
00456 fread(&Rewy, sizeof(Real),1,fp); fread(&Imwy, sizeof(Real),1,fp);
00457 fread(&Rewz, sizeof(Real),1,fp); fread(&Imwz, sizeof(Real),1,fp);
00458 fread(&Rerho, sizeof(Real),1,fp);fread(&Imrho, sizeof(Real),1,fp);
00459 fread(&omg, sizeof(Real),1,fp); fread(&s, sizeof(Real),1,fp);
00460 fread(&Iso_csound, sizeof(Real),1,fp); Iso_csound2 = SQR(Iso_csound);
00461 fread(&kx, sizeof(Real),1,fp); fread(&kz, sizeof(Real),1,fp);
00462 }
00463
00464 if (pG->my_id == 0)
00465 #ifdef MPI_PARALLEL
00466 sprintf(name, "../%s_%d_%d.dat", pG->outfilename,Nx,ipert);
00467 #else
00468 sprintf(name, "%s_%d_%d.dat", pG->outfilename,Nx,ipert);
00469 #endif
00470
00471
00472 interval = par_getd_def("problem","interval",0.01/Omega_0);
00473 data_output_enroll(pG->time,interval,0,OutputModeAmplitude,
00474 NULL,NULL,0,0.0,0.0,0,0,1,property_all);
00475
00476 return;
00477 }
00478
00479
00480
00481
00482 static Real expr_rhodif(const Grid *pG, const int i, const int j, const int k)
00483 {
00484 Real x1,x2,x3;
00485 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00486 return pG->U[k][j][i].d - rho0;
00487 }
00488
00489
00490
00491
00492 static Real expr_dVx(const Grid *pG, const int i, const int j, const int k)
00493 {
00494 Real x1,x2,x3;
00495 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00496 return pG->U[k][j][i].M1/pG->U[k][j][i].d - uxNSH;
00497 }
00498
00499
00500
00501
00502 static Real expr_dVy(const Grid *pG, const int i, const int j, const int k)
00503 {
00504 Real x1,x2,x3;
00505 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00506 #ifdef FARGO
00507 return pG->U[k][j][i].M2/pG->U[k][j][i].d - uyNSH;
00508 #else
00509 return (pG->U[k][j][i].M2/pG->U[k][j][i].d -uyNSH + qshear*Omega_0*x1);
00510 #endif
00511 }
00512
00513
00514
00515
00516 static Real expr_rhopardif(const Grid *pG, const int i, const int j, const int k)
00517 {
00518 Real x1,x2,x3;
00519 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00520 return pG->Coup[k][j][i].grid_d - rho0*mratio;
00521 }
00522
00523
00524
00525
00526 static Real expr_dVxpar(const Grid *pG, const int i, const int j, const int k)
00527 {
00528 Real x1,x2,x3;
00529 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00530 return expr_V1par(pG,i,j,k) - wxNSH;
00531 }
00532
00533
00534
00535
00536 static Real expr_dVypar(const Grid *pG, const int i, const int j, const int k)
00537 {
00538 Real x1,x2,x3;
00539 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00540 #ifdef FARGO
00541 return expr_V2par(pG,i,j,k) - wyNSH;
00542 #else
00543 return expr_V2par(pG,i,j,k) - wyNSH + qshear*Omega_0*x1);
00544 #endif
00545 }
00546
00547 Gasfun_t get_usr_expr(const char *expr)
00548 {
00549 if(strcmp(expr,"difd")==0) return expr_rhodif;
00550 if(strcmp(expr,"dVx")==0) return expr_dVx;
00551 if(strcmp(expr,"dVy")==0) return expr_dVy;
00552 if(strcmp(expr,"difdpar")==0) return expr_rhopardif;
00553 if(strcmp(expr,"dVxpar")==0) return expr_dVxpar;
00554 if(strcmp(expr,"dVypar")==0) return expr_dVypar;
00555 return NULL;
00556 }
00557
00558 VGFunout_t get_usr_out_fun(const char *name){
00559 return NULL;
00560 }
00561
00562 #ifdef PARTICLES
00563 PropFun_t get_usr_par_prop(const char *name)
00564 {
00565 if (strcmp(name,"downsamp")==0) return property_mybin;
00566 return NULL;
00567 }
00568
00569
00570
00571
00572 void gasvshift(const Real x1, const Real x2, const Real x3,
00573 Real *u1, Real *u2, Real *u3)
00574 {
00575 return;
00576 }
00577
00578 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00579 const Real v1, const Real v2, const Real v3)
00580 {
00581 ft->x1 -= 2.0*etavk*Omega_0;
00582 return;
00583 }
00584 #endif
00585
00586 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00587 {
00588 return;
00589 }
00590
00591
00592
00593
00594
00595
00596
00597 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00598 {
00599 return;
00600 }
00601
00602
00603
00604
00605
00606
00607 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00608 {
00609 Real phi=0.0;
00610 #ifndef FARGO
00611 phi -= qshear*SQR(Omega_0*x1);
00612 #endif
00613 return phi;
00614 }
00615
00616
00617
00618 static Real pert_even(Real fR, Real fI, Real x, Real z, Real t)
00619 {
00620 return (fR*cos(kx*x-omg*t)-fI*sin(kx*x-omg*t))*cos(kz*z)*exp(s*t);
00621 }
00622
00623
00624
00625 static Real pert_odd(Real fR, Real fI, Real x, Real z, Real t)
00626 {
00627 return -(fR*sin(kx*x-omg*t)+fI*cos(kx*x-omg*t))*sin(kz*z)*exp(s*t);
00628 }
00629
00630
00631
00632 static int property_mybin(const Grain *gr, const GrainAux *grsub)
00633 {
00634 long a,b,c,d,e,ds,sp;
00635
00636 sp = MAX(downsamp/Npar3,1);
00637 ds = Npar3*sp;
00638
00639 a = gr->my_id/ds;
00640 b = gr->my_id - a*ds;
00641
00642 c = gr->my_id/(Npar3*Nx);
00643 d = c/sp;
00644 e = c-sp*d;
00645
00646 if ((e == 0) && (b == 0) && (gr->pos == 1))
00647 return 1;
00648 else
00649 return 0;
00650 }
00651
00652
00653
00654
00655 void OutputModeAmplitude(Grid *pGrid, Domain *pDomain, Output *pOut)
00656 {
00657 FILE *fid;
00658 int i,j,k;
00659 Real dm,dparm,uxm,uym,uzm,wxm,wym,wzm;
00660
00661 dm=0.0; dparm=0.0; uxm=0.0; uym=0.0; uzm=0.0; wxm=0.0; wym=0.0; wzm=0.0;
00662
00663 for (k=pGrid->ks; k<=pGrid->ke; k++) {
00664 for (j=pGrid->js; j<=pGrid->je; j++) {
00665 for (i=pGrid->is; i<=pGrid->ie; i++) {
00666 dm = MAX(dm,fabs(expr_rhodif(pGrid,i,j,k)));
00667 dparm=MAX(dparm,fabs(expr_rhopardif(pGrid,i,j,k)));
00668 uxm = MAX(uxm,fabs(expr_dVx(pGrid,i,j,k)));
00669 wxm = MAX(wxm,fabs(expr_dVxpar(pGrid,i,j,k)));
00670 uym = MAX(uym,fabs(expr_dVy(pGrid,i,j,k)));
00671 wym = MAX(wym,fabs(expr_dVypar(pGrid,i,j,k)));
00672 uzm = MAX(uzm,fabs(expr_V3(pGrid,i,j,k)));
00673 wzm = MAX(wzm,fabs(expr_V3par(pGrid,i,j,k)));
00674 }}}
00675
00676 #ifdef MPI_PARALLEL
00677 Real sendbuf[8], recvbuf[8];
00678 int err;
00679 sendbuf[0]=dm; sendbuf[1]=dparm;
00680 sendbuf[2]=uxm; sendbuf[3]=wxm;
00681 sendbuf[4]=uym; sendbuf[5]=wym;
00682 sendbuf[6]=uzm; sendbuf[7]=wzm;
00683
00684 err = MPI_Reduce(sendbuf,recvbuf,8,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
00685 if(err) ath_error("[streaming3d]: MPI_Reduce returned error code %d\n",err);
00686
00687 if (pGrid->my_id == 0) {
00688 dm=recvbuf[0]; dparm=recvbuf[1];
00689 uxm=recvbuf[2]; wxm=recvbuf[3];
00690 uym=recvbuf[4]; wym=recvbuf[5];
00691 uzm=recvbuf[6]; wzm=recvbuf[7];
00692 }
00693 #endif
00694
00695 if (pGrid->my_id == 0) {
00696 fid = fopen(name,"a+");
00697 fprintf(fid,"%e %e %e %e %e %e %e %e %e\n",pGrid->time,dm,dparm,uxm,wxm,uym,wym,uzm,wzm);
00698 fclose(fid);
00699 }
00700
00701 return;
00702 }
00703
00704
00705
00706
00707
00708 #define IM1 2147483563
00709 #define IM2 2147483399
00710 #define AM (1.0/IM1)
00711 #define IMM1 (IM1-1)
00712 #define IA1 40014
00713 #define IA2 40692
00714 #define IQ1 53668
00715 #define IQ2 52774
00716 #define IR1 12211
00717 #define IR2 3791
00718 #define NTAB 32
00719 #define NDIV (1+IMM1/NTAB)
00720 #define RNMX (1.0-DBL_EPSILON)
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735 double ran2(long int *idum)
00736 {
00737 int j;
00738 long int k;
00739 static long int idum2=123456789;
00740 static long int iy=0;
00741 static long int iv[NTAB];
00742 double temp;
00743
00744 if (*idum <= 0) {
00745 if (-(*idum) < 1) *idum=1;
00746 else *idum = -(*idum);
00747 idum2=(*idum);
00748 for (j=NTAB+7;j>=0;j--) {
00749 k=(*idum)/IQ1;
00750 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00751 if (*idum < 0) *idum += IM1;
00752 if (j < NTAB) iv[j] = *idum;
00753 }
00754 iy=iv[0];
00755 }
00756 k=(*idum)/IQ1;
00757 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00758 if (*idum < 0) *idum += IM1;
00759 k=idum2/IQ2;
00760 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00761 if (idum2 < 0) idum2 += IM2;
00762 j=(int)(iy/NDIV);
00763 iy=iv[j]-idum2;
00764 iv[j] = *idum;
00765 if (iy < 1) iy += IMM1;
00766 if ((temp=AM*iy) > RNMX) return RNMX;
00767 else return temp;
00768 }
00769
00770 #undef IM1
00771 #undef IM2
00772 #undef AM
00773 #undef IMM1
00774 #undef IA1
00775 #undef IA2
00776 #undef IQ1
00777 #undef IQ2
00778 #undef IR1
00779 #undef IR2
00780 #undef NTAB
00781 #undef NDIV
00782 #undef RNMX