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
00028
00029
00030
00031
00032
00033
00034 #include <stdio.h>
00035 #include <stdlib.h>
00036 #include <math.h>
00037 #include "../defs.h"
00038 #include "../athena.h"
00039 #include "../prototypes.h"
00040 #include "prototypes.h"
00041 #include "particle.h"
00042 #include "../globals.h"
00043
00044 #ifdef PARTICLES
00045
00046
00047
00048
00049
00050
00051 int compare_gr(Grid *pG, Vector cell1, Grain gr1, Grain gr2);
00052 void quicksort_particle(Grid *pG, Vector cell1, long start, long end);
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 void getwei_linear(Grid *pG, Real x1, Real x2, Real x3, Vector cell1,
00085 Real weight[3][3][3], int *is, int *js, int *ks)
00086 {
00087 int i, j, k, i1, j1, k1;
00088 Real a, b, c;
00089 Real wei1[2], wei2[2], wei3[2];
00090
00091
00092
00093 if (cell1.x1 > 0.0) {
00094 i = celli(pG, x1, cell1.x1, &i1, &a);
00095 i1 = i+i1-1; *is = i1;
00096 wei1[1] = a - i1 - 0.5;
00097 wei1[0] = 1.0 - wei1[1];
00098 }
00099 else {
00100 *is = pG->is;
00101 wei1[1] = 0.0;
00102 wei1[0] = 1.0;
00103 }
00104
00105
00106 if (cell1.x2 > 0.0) {
00107 j = cellj(pG, x2, cell1.x2, &j1, &b);
00108 j1 = j+j1-1; *js = j1;
00109 wei2[1] = b - j1 - 0.5;
00110 wei2[0] = 1.0 - wei2[1];
00111 }
00112 else {
00113 *js = pG->js;
00114 wei2[1] = 0.0;
00115 wei2[0] = 1.0;
00116 }
00117
00118
00119 if (cell1.x3 > 0.0) {
00120 k = cellk(pG, x3, cell1.x3, &k1, &c);
00121 k1 = k+k1-1; *ks = k1;
00122 wei3[1] = c - k1 - 0.5;
00123 wei3[0] = 1.0 - wei3[1];
00124 }
00125 else {
00126 *ks = pG->ks;
00127 wei3[1] = 0.0;
00128 wei3[0] = 1.0;
00129 }
00130
00131
00132 for (k=0; k<2; k++)
00133 for (j=0; j<2; j++)
00134 for (i=0; i<2; i++)
00135 weight[k][j][i] = wei1[i] * wei2[j] * wei3[k];
00136
00137 return;
00138 }
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149 void getwei_TSC(Grid *pG, Real x1, Real x2, Real x3, Vector cell1,
00150 Real weight[3][3][3], int *is, int *js, int *ks)
00151 {
00152 int i, j, k, i1, j1, k1;
00153 Real a, b, c, d;
00154 Real wei1[3], wei2[3], wei3[3];
00155
00156
00157
00158 if (cell1.x1 > 0.0) {
00159 celli(pG, x1, cell1.x1, &i, &a);
00160 *is = i - 1;
00161 d = a - i;
00162 wei1[0] = 0.5*SQR(1.0-d);
00163 wei1[1] = 0.75-SQR(d-0.5);
00164 wei1[2] = 0.5*SQR(d);
00165 }
00166 else {
00167 *is = pG->is;
00168 wei1[1] = 0.0; wei1[2] = 0.0;
00169 wei1[0] = 1.0;
00170 }
00171
00172
00173 if (cell1.x2 > 0.0) {
00174 cellj(pG, x2, cell1.x2, &j, &b);
00175 *js = j - 1;
00176 d = b - j;
00177 wei2[0] = 0.5*SQR(1.0-d);
00178 wei2[1] = 0.75-SQR(d-0.5);
00179 wei2[2] = 0.5*SQR(d);
00180 }
00181 else {
00182 *js = pG->js;
00183 wei2[1] = 0.0; wei2[2] = 0.0;
00184 wei2[0] = 1.0;
00185 }
00186
00187
00188 if (cell1.x3 > 0.0) {
00189 cellk(pG, x3, cell1.x3, &k, &c);
00190 *ks = k - 1;
00191 d = c - k;
00192 wei3[0] = 0.5*SQR(1.0-d);
00193 wei3[1] = 0.75-SQR(d-0.5);
00194 wei3[2] = 0.5*SQR(d);
00195 }
00196 else {
00197 *ks = pG->ks;
00198 wei3[1] = 0.0; wei3[2] = 0.0;
00199 wei3[0] = 1.0;
00200 }
00201
00202
00203 for (k=0; k<3; k++)
00204 for (j=0; j<3; j++)
00205 for (i=0; i<3; i++)
00206 weight[k][j][i] = wei1[i] * wei2[j] * wei3[k];
00207
00208 return;
00209 }
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 void getwei_QP(Grid *pG, Real x1, Real x2, Real x3, Vector cell1,
00222 Real weight[3][3][3], int *is, int *js, int *ks)
00223 {
00224 int i, j, k, i1, j1, k1;
00225 Real a, b, c, d;
00226 Real wei1[3], wei2[3], wei3[3];
00227
00228
00229
00230 if (cell1.x1 > 0.0) {
00231 celli(pG, x1, cell1.x1, &i, &a);
00232 *is = i - 1;
00233 d = a - i;
00234 wei1[0] = 0.5*(0.5-d)*(1.5-d);
00235 wei1[1] = 1.0-SQR(d-0.5);
00236 wei1[2] = 0.5*(d-0.5)*(d+0.5);
00237 }
00238 else {
00239 *is = pG->is;
00240 wei1[1] = 0.0; wei1[2] = 0.0;
00241 wei1[0] = 1.0;
00242 }
00243
00244
00245 if (cell1.x2 > 0.0) {
00246 cellj(pG, x2, cell1.x2, &j, &b);
00247 *js = j - 1;
00248 d = b - j;
00249 wei2[0] = 0.5*(0.5-d)*(1.5-d);
00250 wei2[1] = 1.0-SQR(d-0.5);
00251 wei2[2] = 0.5*(d-0.5)*(d+0.5);
00252 }
00253 else {
00254 *js = pG->js;
00255 wei2[1] = 0.0; wei2[2] = 0.0;
00256 wei2[0] = 1.0;
00257 }
00258
00259
00260 if (cell1.x3 > 0.0) {
00261 cellk(pG, x3, cell1.x3, &k, &c);
00262 *ks = k - 1;
00263 d = c - k;
00264 wei3[0] = 0.5*(0.5-d)*(1.5-d);
00265 wei3[1] = 1.0-SQR(d-0.5);
00266 wei3[2] = 0.5*(d-0.5)*(d+0.5);
00267 }
00268 else {
00269 *ks = pG->ks;
00270 wei3[1] = 0.0; wei3[2] = 0.0;
00271 wei3[0] = 1.0;
00272 }
00273
00274
00275 for (k=0; k<3; k++)
00276 for (j=0; j<3; j++)
00277 for (i=0; i<3; i++)
00278 weight[k][j][i] = wei1[i] * wei2[j] * wei3[k];
00279
00280 return;
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296 int getvalues(Grid *pG, Real weight[3][3][3], int is, int js, int ks,
00297 #ifndef FEEDBACK
00298 Real *rho, Real *u1, Real *u2, Real *u3, Real *cs
00299 #else
00300 Real *rho, Real *u1, Real *u2, Real *u3, Real *cs, Real *stiff
00301 #endif
00302 ){
00303 int n0,i,j,k,i0,j0,k0,i1,j1,k1,i2,j2,k2;
00304 Real D, v1, v2, v3;
00305 GPCouple *pq;
00306 #ifndef ISOTHERMAL
00307 Real C = 0.0;
00308 #endif
00309 #ifdef FEEDBACK
00310 Real stiffness=0.0;
00311 #endif
00312 Real totwei, totwei1;
00313
00314
00315 D = 0.0; v1 = 0.0; v2 = 0.0; v3 = 0.0;
00316 totwei = 0.0; totwei1 = 1.0;
00317
00318
00319
00320 n0 = ncell-1;
00321 k1 = MAX(ks, klp); k2 = MIN(ks+n0, kup);
00322 j1 = MAX(js, jlp); j2 = MIN(js+n0, jup);
00323 i1 = MAX(is, ilp); i2 = MIN(is+n0, iup);
00324
00325 for (k=k1; k<=k2; k++) {
00326 k0=k-k1;
00327 for (j=j1; j<=j2; j++) {
00328 j0=j-j1;
00329 for (i=i1; i<=i2; i++) {
00330 i0=i-i1;
00331
00332 pq = &(pG->Coup[k][j][i]);
00333 D += weight[k0][j0][i0] * pq->grid_d;
00334 v1 += weight[k0][j0][i0] * pq->grid_v1;
00335 v2 += weight[k0][j0][i0] * pq->grid_v2;
00336 v3 += weight[k0][j0][i0] * pq->grid_v3;
00337 #ifndef ISOTHERMAL
00338 C += weight[k0][j0][i0] * pq->grid_cs;
00339 #endif
00340 #ifdef FEEDBACK
00341 stiffness += weight[k0][j0][i0] * pq->FBstiff;
00342 #endif
00343 totwei += weight[k0][j0][i0];
00344 }
00345 }
00346 }
00347 if (totwei < TINY_NUMBER)
00348 return -1;
00349
00350 totwei1 = 1.0/totwei;
00351 *rho = D*totwei1;
00352 *u1 = v1*totwei1; *u2 = v2*totwei1; *u3 = v3*totwei1;
00353 #ifdef ISOTHERMAL
00354 *cs = Iso_csound;
00355 #else
00356 *cs = C*totwei1;
00357 #endif
00358 #ifdef FEEDBACK
00359 *stiff = stiffness*totwei1;
00360 #endif
00361
00362 return 0;
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 Real get_ts_general(Grid *pG, int type, Real rho, Real cs, Real vd)
00384 {
00385 Real tstop;
00386 Real a, rhos;
00387 Real alam;
00388 Real rhoa;
00389 Real Re, CD;
00390
00391
00392 a = pG->grproperty[type].rad;
00393 rhos = pG->grproperty[type].rho;
00394 rhoa = grrhoa[type];
00395
00396
00397 alam = alamcoeff * a * rho;
00398
00399
00400 if (alam < 2.25) {
00401 tstop = rhoa/(rho*cs);
00402 }
00403 else {
00404 Re = 4.0*alam*vd/cs;
00405
00406 if (Re < 1.) CD = 24.0/Re;
00407 else if (Re < 800.0) CD = 24.0*exp(-0.6*log(Re));
00408 else CD = 0.44;
00409
00410 tstop = rhoa/(rho*vd*CD);
00411 }
00412
00413
00414 if (tstop < 1.0e-8*pG->dt)
00415 tstop = 1.0e-8*pG->dt;
00416
00417 return tstop;
00418 }
00419
00420
00421
00422
00423
00424
00425
00426 Real get_ts_epstein(Grid *pG, int type, Real rho, Real cs, Real vd)
00427 {
00428 Real tstop = grrhoa[type]/(rho*cs);
00429
00430
00431 if (tstop < 1.0e-8*pG->dt)
00432 tstop = 1.0e-8*pG->dt;
00433
00434 return tstop;
00435 }
00436
00437
00438
00439 Real get_ts_fixed(Grid *pG, int type, Real rho, Real cs, Real vd)
00440 {
00441 Real tstop = tstop0[type];
00442
00443
00444 if (tstop < 1.0e-8*pG->dt)
00445 tstop = 1.0e-8*pG->dt;
00446
00447 return tstop;
00448 }
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469 void get_gasinfo(Grid *pG)
00470 {
00471 int i,j,k;
00472 Real rho1;
00473 GPCouple *pq;
00474 #ifndef BAROTROPIC
00475 Real P;
00476 #endif
00477
00478
00479 for (k=klp; k<=kup; k++)
00480 for (j=jlp; j<=jup; j++)
00481 for (i=ilp; i<=iup; i++)
00482 {
00483 pq = &(pG->Coup[k][j][i]);
00484 rho1 = 1.0/(pG->U[k][j][i].d);
00485
00486 pq->grid_d = pG->U[k][j][i].d;
00487 pq->grid_v1 = pG->U[k][j][i].M1 * rho1;
00488 pq->grid_v2 = pG->U[k][j][i].M2 * rho1;
00489 pq->grid_v3 = pG->U[k][j][i].M3 * rho1;
00490
00491 #ifndef ISOTHERMAL
00492 #ifndef BAROTROPIC
00493
00494 P = pG->U[k][j][i].E - 0.5*pG->U[k][j][i].d*(SQR(pq->grid_v1) \
00495 + SQR(pq->grid_v2) + SQR(pq->grid_v3));
00496 #ifdef MHD
00497 P = P - 0.5*(SQR(pG->U[k][j][i].B1c)+SQR(pG->U[k][j][i].B2c)
00498 +SQR(pG->U[k][j][i].B3c));
00499 #endif
00500 P = MAX(Gamma_1*P, TINY_NUMBER);
00501 pq->grid_cs = sqrt(Gamma*P*rho1);
00502 #else
00503 ath_error("[get_gasinfo] can not calculate the sound speed!\n");
00504 #endif
00505 #endif
00506 }
00507
00508 return;
00509 }
00510
00511 #ifdef FEEDBACK
00512
00513
00514
00515 void feedback_clear(Grid *pG)
00516 {
00517 int i,j,k;
00518 GPCouple *pq;
00519
00520 for (k=klp; k<=kup; k++)
00521 for (j=jlp; j<=jup; j++)
00522 for (i=ilp; i<=iup; i++) {
00523 pq = &(pG->Coup[k][j][i]);
00524
00525 pq->fb1 = 0.0;
00526 pq->fb2 = 0.0;
00527 pq->fb3 = 0.0;
00528
00529 pq->Eloss = 0.0;
00530 }
00531
00532 return;
00533 }
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546 void distrFB_pred(Grid *pG, Real weight[3][3][3], int is, int js, int ks,
00547 #ifndef BAROTROPIC
00548 Vector fb, Real stiffness, Real Elosspar
00549 #else
00550 Vector fb, Real stiffness
00551 #endif
00552 ){
00553 int n0,i,j,k,i0,j0,k0,i1,j1,k1,i2,j2,k2;
00554 GPCouple *pq;
00555
00556
00557 n0 = ncell-1;
00558 k1 = MAX(ks, klp); k2 = MIN(ks+n0, kup);
00559 j1 = MAX(js, jlp); j2 = MIN(js+n0, jup);
00560 i1 = MAX(is, ilp); i2 = MIN(is+n0, iup);
00561 for (k=k1; k<=k2; k++) {
00562 k0 = k-k1;
00563 for (j=j1; j<=j2; j++) {
00564 j0 = j-j1;
00565 for (i=i1; i<=i2; i++) {
00566 i0 = i-i1;
00567 pq = &(pG->Coup[k][j][i]);
00568
00569 pq->fb1 += weight[k0][j0][i0] * fb.x1;
00570 pq->fb2 += weight[k0][j0][i0] * fb.x2;
00571 pq->fb3 += weight[k0][j0][i0] * fb.x3;
00572 #ifndef BAROTROPIC
00573 pq->Eloss += weight[k0][j0][i0] * Elosspar;
00574 #endif
00575 pq->FBstiff += weight[k0][j0][i0] * stiffness;
00576 }
00577 }
00578 }
00579
00580 return;
00581 }
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595 void distrFB_corr(Grid *pG, Real weight[3][3][3], int is, int js, int ks,
00596 Vector fb, Real Elosspar)
00597 {
00598 int n0,i,j,k,i0,j0,k0,i1,j1,k1,i2,j2,k2;
00599 GPCouple *pq;
00600
00601
00602 n0 = ncell-1;
00603 k1 = MAX(ks, klp); k2 = MIN(ks+n0, kup);
00604 j1 = MAX(js, jlp); j2 = MIN(js+n0, jup);
00605 i1 = MAX(is, ilp); i2 = MIN(is+n0, iup);
00606 for (k=k1; k<=k2; k++) {
00607 k0 = k-k1;
00608 for (j=j1; j<=j2; j++) {
00609 j0 = j-j1;
00610 for (i=i1; i<=i2; i++) {
00611 i0 = i-i1;
00612 pq = &(pG->Coup[k][j][i]);
00613
00614 pq->fb1 += weight[k0][j0][i0] * fb.x1;
00615 pq->fb2 += weight[k0][j0][i0] * fb.x2;
00616 pq->fb3 += weight[k0][j0][i0] * fb.x3;
00617
00618 pq->Eloss += weight[k0][j0][i0] * Elosspar;
00619 }
00620 }
00621 }
00622
00623 return;
00624 }
00625
00626 #endif
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645 void shuffle(Grid *pG)
00646 {
00647 Grain *cur, *allgr;
00648 Vector cell1;
00649
00650 if (pG->Nx1 > 1) cell1.x1 = 1.0/pG->dx1; else cell1.x1 = 0.0;
00651 if (pG->Nx2 > 1) cell1.x2 = 1.0/pG->dx2; else cell1.x2 = 0.0;
00652 if (pG->Nx3 > 1) cell1.x3 = 1.0/pG->dx3; else cell1.x3 = 0.0;
00653
00654
00655 ath_pout(0, "Resorting particles...\n");
00656
00657
00658 quicksort_particle(pG, cell1, 0, pG->nparticle-1);
00659
00660 return;
00661 }
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673 int compare_gr(Grid *pG, Vector cell1, Grain gr1, Grain gr2)
00674 {
00675 int i1,j1,k1, i2,j2,k2;
00676
00677 k1 = (int)((gr1.x3 - pG->x3_0) * cell1.x3);
00678 k2 = (int)((gr2.x3 - pG->x3_0) * cell1.x3);
00679 if (k1 < k2) return 1;
00680 if (k1 > k2) return 2;
00681
00682 j1 = (int)((gr1.x2 - pG->x2_0) * cell1.x2);
00683 j2 = (int)((gr2.x2 - pG->x2_0) * cell1.x2);
00684 if (j1 < j2) return 1;
00685 if (j1 > j2) return 2;
00686
00687 i1 = (int)((gr1.x1 - pG->x1_0) * cell1.x1);
00688 i2 = (int)((gr2.x1 - pG->x1_0) * cell1.x1);
00689 if (i1 < i2) return 1;
00690 if (i1 > i2) return 2;
00691
00692
00693 return 1;
00694 }
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706 void quicksort_particle(Grid *pG, Vector cell1, long start, long end)
00707 {
00708 long i, pivot;
00709 Grain gr;
00710 if (end <= start) return;
00711
00712
00713 pivot = (long)((start+end+1)/2);
00714
00715
00716 gr = pG->particle[pivot];
00717 pG->particle[pivot] = pG->particle[start];
00718 pG->particle[start] = gr;
00719
00720
00721 pivot = start;
00722 i = start + 1;
00723
00724
00725 while (i <= end) {
00726 if (compare_gr(pG, cell1, pG->particle[pivot], pG->particle[i]) == 2)
00727 {
00728 gr = pG->particle[pivot];
00729 pG->particle[pivot] = pG->particle[i];
00730 pG->particle[i] = pG->particle[pivot+1];
00731 pG->particle[pivot+1] = gr;
00732 pivot += 1;
00733 }
00734 i += 1;
00735 }
00736
00737
00738 quicksort_particle(pG, cell1, start, pivot-1);
00739 quicksort_particle(pG, cell1, pivot+1, end);
00740
00741 return;
00742 }
00743
00744 #endif