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