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 #include <stdio.h>
00029 #include <stdlib.h>
00030 #include <math.h>
00031 #include "../defs.h"
00032 #include "../athena.h"
00033 #include "../prototypes.h"
00034 #include "prototypes.h"
00035 #include "particle.h"
00036 #include "../globals.h"
00037
00038 #ifdef PARTICLES
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 void Delete_Ghost(Grid *pG);
00050 void JudgeCrossing(Grid *pG, Real x1, Real x2, Real x3, Grain *gr);
00051 Vector Get_Drag(Grid *pG, int type, Real x1, Real x2, Real x3,
00052 Real v1, Real v2, Real v3, Vector cell1, Real *tstop1);
00053 Vector Get_Force(Grid *pG, Real x1, Real x2, Real x3,
00054 Real v1, Real v2, Real v3);
00055 Vector Get_Term(Grid *pG, int type, Real x1, Real x2, Real x3, Vector cell1,
00056 Real *tstop);
00057 Vector Get_ForceDiff(Grid *pG, Real x1, Real x2, Real x3,
00058 Real v1, Real v2, Real v3);
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 void Integrate_Particles(Grid *pG, Domain *pD)
00076 {
00077 Grain *curG, *curP, mygr;
00078 long p;
00079 Real dv1, dv2, dv3, ts, t1;
00080 Vector cell1;
00081 Vector vterm;
00082
00083
00084 #ifdef FEEDBACK
00085 feedback_clear(pG);
00086 #endif
00087
00088 curP = &(mygr);
00089
00090
00091 if (pG->Nx1 > 1) cell1.x1 = 1.0/pG->dx1; else cell1.x1 = 0.0;
00092 if (pG->Nx2 > 1) cell1.x2 = 1.0/pG->dx2; else cell1.x2 = 0.0;
00093 if (pG->Nx3 > 1) cell1.x3 = 1.0/pG->dx3; else cell1.x3 = 0.0;
00094
00095
00096 Delete_Ghost(pG);
00097
00098 p = 0;
00099 while (p<pG->nparticle)
00100 {
00101 curG = &(pG->particle[p]);
00102
00103
00104 switch(pG->grproperty[curG->property].integrator)
00105 {
00106 case 1:
00107 int_par_exp(pG,curG,cell1, &dv1,&dv2,&dv3, &ts);
00108 break;
00109
00110 case 2:
00111 int_par_semimp(pG,curG,cell1, &dv1,&dv2,&dv3, &ts);
00112 break;
00113
00114 case 3:
00115 int_par_fulimp(pG,curG,cell1, &dv1,&dv2,&dv3, &ts);
00116 break;
00117
00118 case 4:
00119 int_par_spec(pG,curG,cell1, &dv1,&dv2,&dv3, &ts);
00120 break;
00121
00122 default:
00123 ath_error("[integrate_particle]: unknown integrator type!");
00124 }
00125
00126
00127
00128
00129 curP->v1 = curG->v1 + dv1;
00130 curP->v2 = curG->v2 + dv2;
00131 curP->v3 = curG->v3 + dv3;
00132
00133
00134 if (pG->Nx1 > 1)
00135 curP->x1 = curG->x1 + 0.5*pG->dt*(curG->v1 + curP->v1);
00136 else
00137 curP->x1 = curG->x1;
00138
00139 if (pG->Nx2 > 1)
00140 curP->x2 = curG->x2 + 0.5*pG->dt*(curG->v2 + curP->v2);
00141 else
00142 curP->x2 = curG->x2;
00143
00144 if (pG->Nx3 > 1)
00145 curP->x3 = curG->x3 + 0.5*pG->dt*(curG->v3 + curP->v3);
00146 else
00147 curP->x3 = curG->x3;
00148
00149 #ifdef FARGO
00150
00151 pG->parsub[p].shift = -0.5*qshear*Omega_0*(curG->x1+curP->x1)*pG->dt;
00152 #endif
00153
00154
00155 if (pG->grproperty[curG->property].integrator == 4)
00156 {
00157 vterm = Get_Term(pG,curG->property,curP->x1,curP->x2,curP->x3,cell1,&t1);
00158 curP->v1 = vterm.x1; curP->v2 = vterm.x2; curP->v2 = vterm.x2;
00159
00160 dv1=curP->v1-curG->v1; dv2=curP->v2-curG->v2; dv3=curP->v3-curG->v3;
00161 }
00162
00163
00164 #ifdef FEEDBACK
00165 feedback_corrector(pG, curG, curP, cell1, dv1, dv2, dv3, ts);
00166 #endif
00167
00168
00169
00170 JudgeCrossing(pG, curP->x1, curP->x2, curP->x3, curG);
00171
00172
00173 curG->x1 = curP->x1;
00174 curG->x2 = curP->x2;
00175 curG->x3 = curP->x3;
00176 curG->v1 = curP->v1;
00177 curG->v2 = curP->v2;
00178 curG->v3 = curP->v3;
00179 p++;
00180
00181 }
00182
00183
00184 ath_pout(0, "In processor %d, there are %ld particles.\n",
00185 pG->my_id, pG->nparticle);
00186
00187 return;
00188 }
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200 void int_par_fulimp(Grid *pG, Grain *curG, Vector cell1,
00201 Real *dv1, Real *dv2, Real *dv3, Real *ts)
00202 {
00203 Real x1n, x2n, x3n;
00204 Vector fd, fr;
00205 Vector fc, fp, ft;
00206 Real ts11, ts12;
00207 Real b0,A,B,C,D,Det1;
00208 #ifdef SHEARING_BOX
00209 Real oh, oh2;
00210 #endif
00211
00212
00213 if (pG->Nx1 > 1) x1n = curG->x1+curG->v1*pG->dt;
00214 else x1n = curG->x1;
00215 if (pG->Nx2 > 1) x2n = curG->x2+curG->v2*pG->dt;
00216 else x2n = curG->x2;
00217 if (pG->Nx3 > 1) x3n = curG->x3+curG->v3*pG->dt;
00218 else x3n = curG->x3;
00219
00220 #ifdef SHEARING_BOX
00221 #ifndef FARGO
00222 if (pG->Nx3 > 1) x2n -= 0.5*qshear*curG->v1*SQR(pG->dt);
00223 #endif
00224 #endif
00225
00226
00227 fd = Get_Drag(pG, curG->property, curG->x1, curG->x2, curG->x3,
00228 curG->v1, curG->v2, curG->v3, cell1, &ts11);
00229
00230 fr = Get_Force(pG, curG->x1, curG->x2, curG->x3,
00231 curG->v1, curG->v2, curG->v3);
00232
00233 fc.x1 = fd.x1+fr.x1;
00234 fc.x2 = fd.x2+fr.x2;
00235 fc.x3 = fd.x3+fr.x3;
00236
00237
00238 fd = Get_Drag(pG, curG->property, x1n, x2n, x3n,
00239 curG->v1, curG->v2, curG->v3, cell1, &ts12);
00240
00241 fr = Get_Force(pG, x1n, x2n, x3n, curG->v1, curG->v2, curG->v3);
00242
00243 fp.x1 = fd.x1+fr.x1;
00244 fp.x2 = fd.x2+fr.x2;
00245 fp.x3 = fd.x3+fr.x3;
00246
00247
00248
00249 b0 = 1.0+pG->dt*ts11;
00250
00251
00252 ft.x1 = 0.5*(fc.x1+b0*fp.x1);
00253 ft.x2 = 0.5*(fc.x2+b0*fp.x2);
00254 ft.x3 = 0.5*(fc.x3+b0*fp.x3);
00255
00256 #ifdef SHEARING_BOX
00257 oh = Omega_0*pG->dt;
00258 if (pG->Nx3 > 1) {
00259 ft.x1 += -oh*fp.x2;
00260 #ifdef FARGO
00261 ft.x2 += 0.5*(2.0-qshear)*oh*fp.x1;
00262 #else
00263 ft.x2 += oh*fp.x1;
00264 #endif
00265 } else {
00266 ft.x1 += -oh*fp.x3;
00267 #ifdef FARGO
00268 ft.x3 += 0.5*(2.0-qshear)*oh*fp.x1;
00269 #else
00270 ft.x3 += oh*fp.x1;
00271 #endif
00272 }
00273 #endif
00274
00275
00276 D = 1.0+0.5*pG->dt*(ts11 + ts12 + pG->dt*ts11*ts12);
00277 #ifdef SHEARING_BOX
00278 oh2 = SQR(oh);
00279 B = oh * (-2.0-(ts11+ts12)*pG->dt);
00280 #ifdef FARGO
00281 A = D - (2.0-qshear)*oh2;
00282 C = 0.5*(qshear-2.0)*B;
00283 #else
00284 A = D - 2.0*oh2;
00285 C = -B;
00286 #endif
00287 Det1 = 1.0/(SQR(A)-B*C);
00288 if (pG->Nx3>1) {
00289 *dv1 = pG->dt*Det1*(ft.x1*A-ft.x2*B);
00290 *dv2 = pG->dt*Det1*(-ft.x1*C+ft.x2*A);
00291 *dv3 = pG->dt*ft.x3/D;
00292 } else {
00293 *dv1 = pG->dt*Det1*(ft.x1*A-ft.x3*B);
00294 *dv3 = pG->dt*Det1*(-ft.x1*C+ft.x3*A);
00295 *dv2 = pG->dt*ft.x2/D;
00296 }
00297 #else
00298 D = 1.0/D;
00299 *dv1 = pG->dt*ft.x1*D;
00300 *dv2 = pG->dt*ft.x2*D;
00301 *dv3 = pG->dt*ft.x3*D;
00302 #endif
00303
00304 *ts = 0.5/ts11+0.5/ts12;
00305
00306 return;
00307 }
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320 void int_par_semimp(Grid *pG, Grain *curG, Vector cell1,
00321 Real *dv1, Real *dv2, Real *dv3, Real *ts)
00322 {
00323 Vector fd, fr, ft;
00324 Real ts1, b, b2;
00325 Real x1n, x2n, x3n;
00326 #ifdef SHEARING_BOX
00327 Real b1, oh;
00328 #endif
00329
00330
00331 if (pG->Nx1 > 1) x1n = curG->x1+0.5*curG->v1*pG->dt;
00332 else x1n = curG->x1;
00333 if (pG->Nx2 > 1) x2n = curG->x2+0.5*curG->v2*pG->dt;
00334 else x2n = curG->x2;
00335 if (pG->Nx3 > 1) x3n = curG->x3+0.5*curG->v3*pG->dt;
00336 else x3n = curG->x3;
00337
00338 #ifdef SHEARING_BOX
00339 #ifndef FARGO
00340 if (pG->Nx3 > 1) x2n -= 0.125*qshear*curG->v1*SQR(pG->dt);
00341 #endif
00342 #endif
00343
00344
00345
00346 fd = Get_Drag(pG, curG->property, x1n, x2n, x3n,
00347 curG->v1, curG->v2, curG->v3, cell1, &ts1);
00348
00349 fr = Get_Force(pG, x1n, x2n, x3n, curG->v1, curG->v2, curG->v3);
00350
00351 ft.x1 = fd.x1+fr.x1;
00352 ft.x2 = fd.x2+fr.x2;
00353 ft.x3 = fd.x3+fr.x3;
00354
00355
00356
00357
00358 b = pG->dt*ts1+2.0;
00359 #ifdef SHEARING_BOX
00360 oh = Omega_0*pG->dt;
00361 #ifdef FARGO
00362 b1 = 1.0/(SQR(b)+2.0*(2.0-qshear)*SQR(oh));
00363 #else
00364 b1 = 1.0/(SQR(b)+4.0*SQR(oh));
00365 #endif
00366 b2 = b*b1;
00367 #else
00368 b2 = 1.0/b;
00369 #endif
00370
00371
00372 #ifdef SHEARING_BOX
00373 if (pG->Nx3>1)
00374 {
00375 *dv1 = pG->dt*2.0*b2*ft.x1 + pG->dt*4.0*oh*b1*ft.x2;
00376 #ifdef FARGO
00377 *dv2 = pG->dt*2.0*b2*ft.x2 - 2.0*(2.0-qshear)*pG->dt*oh*b1*ft.x1;
00378 #else
00379 *dv2 = pG->dt*2.0*b2*ft.x2 - 4.0*pG->dt*oh*b1*ft.x1;
00380 #endif
00381 *dv3 = pG->dt*2.0*ft.x3/b;
00382 }
00383 else
00384 {
00385 *dv1 = pG->dt*2.0*b2*ft.x1 + pG->dt*4.0*oh*b1*ft.x3;
00386 *dv2 = pG->dt*2.0*ft.x2/b;
00387 #ifdef FARGO
00388 *dv3 = pG->dt*2.0*b2*ft.x3 - 2.0*(2.0-qshear)*pG->dt*oh*b1*ft.x1;
00389 #else
00390 *dv3 = pG->dt*2.0*b2*ft.x3 - 4.0*pG->dt*oh*b1*ft.x1;
00391 #endif
00392 }
00393 #else
00394 *dv1 = pG->dt*2.0*b2*ft.x1;
00395 *dv2 = pG->dt*2.0*b2*ft.x2;
00396 *dv3 = pG->dt*2.0*b2*ft.x3;
00397 #endif
00398
00399 *ts = 1.0/ts1;
00400
00401 return;
00402 }
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415 void int_par_exp(Grid *pG, Grain *curG, Vector cell1,
00416 Real *dv1, Real *dv2, Real *dv3, Real *ts)
00417 {
00418 Vector fd, fr, ft;
00419 Real ts1;
00420 Real x1n, x2n, x3n;
00421 Real v1n, v2n, v3n;
00422
00423
00424 if (pG->Nx1 > 1)
00425 x1n = curG->x1+0.5*curG->v1*pG->dt;
00426 else x1n = curG->x1;
00427 if (pG->Nx2 > 1)
00428 x2n = curG->x2+0.5*curG->v2*pG->dt;
00429 else x2n = curG->x2;
00430 if (pG->Nx3 > 1)
00431 x3n = curG->x3+0.5*curG->v3*pG->dt;
00432 else x3n = curG->x3;
00433
00434 #ifdef SHEARING_BOX
00435 #ifndef FARGO
00436 if (pG->Nx3 > 1) x2n -= 0.125*qshear*curG->v1*SQR(pG->dt);
00437 #endif
00438 #endif
00439
00440
00441 fd = Get_Drag(pG, curG->property, curG->x1, curG->x2, curG->x3,
00442 curG->v1, curG->v2, curG->v3, cell1, &ts1);
00443
00444 fr = Get_Force(pG, curG->x1, curG->x2, curG->x3,
00445 curG->v1, curG->v2, curG->v3);
00446
00447 ft.x1 = fd.x1+fr.x1;
00448 ft.x2 = fd.x2+fr.x2;
00449 ft.x3 = fd.x3+fr.x3;
00450
00451 v1n = curG->v1 + 0.5*ft.x1*pG->dt;
00452 v2n = curG->v2 + 0.5*ft.x2*pG->dt;
00453 v3n = curG->v3 + 0.5*ft.x3*pG->dt;
00454
00455
00456 fd = Get_Drag(pG, curG->property, x1n, x2n, x3n, v1n, v2n, v3n, cell1, &ts1);
00457
00458 fr = Get_Force(pG, x1n, x2n, x3n, v1n, v2n, v3n);
00459
00460 ft.x1 = fd.x1+fr.x1;
00461 ft.x2 = fd.x2+fr.x2;
00462 ft.x3 = fd.x3+fr.x3;
00463
00464
00465 *dv1 = ft.x1*pG->dt;
00466 *dv2 = ft.x2*pG->dt;
00467 *dv3 = ft.x3*pG->dt;
00468
00469 *ts = 1.0/ts1;
00470
00471 return;
00472 }
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485 void int_par_spec(Grid *pG, Grain *curG, Vector cell1,
00486 Real *dv1, Real *dv2, Real *dv3, Real *ts)
00487 {
00488 Vector vterm;
00489 Real x1n, x2n, x3n;
00490
00491
00492 if (pG->Nx1 > 1)
00493 x1n = curG->x1+0.5*curG->v1*pG->dt;
00494 else x1n = curG->x1;
00495 if (pG->Nx2 > 1)
00496 x2n = curG->x2+0.5*curG->v2*pG->dt;
00497 else x2n = curG->x2;
00498 if (pG->Nx3 > 1)
00499 x3n = curG->x3+0.5*curG->v3*pG->dt;
00500 else x3n = curG->x3;
00501
00502 #ifdef SHEARING_BOX
00503 #ifndef FARGO
00504 if (pG->Nx3 > 1) x2n -= 0.125*qshear*curG->v1*SQR(pG->dt);
00505 #endif
00506 #endif
00507
00508
00509 vterm = Get_Term(pG, curG->property, x1n, x2n, x3n, cell1, ts);
00510
00511
00512 *dv1 = 2.0*(vterm.x1 - curG->v1);
00513 *dv2 = 2.0*(vterm.x2 - curG->v2);
00514 *dv3 = 2.0*(vterm.x3 - curG->v3);
00515
00516 return;
00517 }
00518
00519 #ifdef FEEDBACK
00520
00521
00522
00523
00524
00525
00526
00527
00528 void feedback_predictor(Grid* pG)
00529 {
00530 int is,js,ks,i,j,k;
00531 long p;
00532 Real weight[3][3][3];
00533 Real rho, cs, tstop;
00534 Real u1, u2, u3;
00535 Real vd1, vd2, vd3, vd;
00536 Real f1, f2, f3;
00537 Real m, ts1h;
00538 Vector cell1;
00539 Vector fb;
00540 #ifndef BAROTROPIC
00541 Real Elosspar;
00542 #endif
00543 Real stiffness;
00544 Grain *cur;
00545
00546
00547 get_gasinfo(pG);
00548
00549 for (k=klp; k<=kup; k++)
00550 for (j=jlp; j<=jup; j++)
00551 for (i=ilp; i<=iup; i++) {
00552
00553 pG->Coup[k][j][i].fb1 = 0.0;
00554 pG->Coup[k][j][i].fb2 = 0.0;
00555 pG->Coup[k][j][i].fb3 = 0.0;
00556 #ifndef BAROTROPIC
00557 pG->Coup[k][j][i].Eloss = 0.0;
00558 #endif
00559 pG->Coup[k][j][i].FBstiff = 0.0;
00560 }
00561
00562
00563 if (pG->Nx1 > 1) cell1.x1 = 1.0/pG->dx1;
00564 else cell1.x1 = 0.0;
00565
00566 if (pG->Nx2 > 1) cell1.x2 = 1.0/pG->dx2;
00567 else cell1.x2 = 0.0;
00568
00569 if (pG->Nx3 > 1) cell1.x3 = 1.0/pG->dx3;
00570 else cell1.x3 = 0.0;
00571
00572
00573 for (p=0; p<pG->nparticle; p++)
00574 {
00575 cur = &(pG->particle[p]);
00576
00577
00578 getweight(pG, cur->x1, cur->x2, cur->x3, cell1, weight, &is, &js, &ks);
00579 if (getvalues(pG, weight, is, js, ks,
00580 &rho, &u1, &u2, &u3, &cs, &stiffness) == 0)
00581 {
00582
00583
00584 gasvshift(cur->x1, cur->x2, cur->x3, &u1, &u2, &u3);
00585
00586 vd1 = u1-cur->v1;
00587 vd2 = u2-cur->v2;
00588 vd3 = u3-cur->v3;
00589 vd = sqrt(vd1*vd1 + vd2*vd2 + vd3*vd3);
00590
00591
00592 tstop = get_ts(pG, cur->property, rho, cs, vd);
00593 ts1h = 0.5*pG->dt/tstop;
00594
00595
00596 m = pG->grproperty[cur->property].m;
00597 fb.x1 = m * vd1 * ts1h;
00598 fb.x2 = m * vd2 * ts1h;
00599 fb.x3 = m * vd3 * ts1h;
00600
00601
00602 stiffness = 2.0*m*ts1h;
00603
00604 #ifndef BAROTROPIC
00605 Elosspar = fb.x1*vd1 + fb.x2*vd2 + fb.x3*vd3;
00606
00607 distrFB_pred(pG, weight, is, js, ks, fb, stiffness, Elosspar);
00608 #else
00609
00610 distrFB_pred(pG, weight, is, js, ks, fb, stiffness);
00611 #endif
00612 }
00613 }
00614
00615
00616 for (k=klp; k<=kup; k++)
00617 for (j=jlp; j<=jup; j++)
00618 for (i=ilp; i<=iup; i++)
00619 {
00620 pG->Coup[k][j][i].FBstiff /= pG->U[k][j][i].d;
00621
00622
00623
00624
00625
00626
00627 }
00628
00629 return;
00630 }
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 void feedback_corrector(Grid *pG, Grain *gri, Grain *grf, Vector cell1,
00644 Real dv1, Real dv2, Real dv3, Real ts)
00645 {
00646 int is, js, ks;
00647 Real x1, x2, x3, v1, v2, v3;
00648 Real mgr;
00649 Real weight[3][3][3];
00650 Real Elosspar;
00651 Vector fb;
00652
00653 mgr = pG->grproperty[gri->property].m;
00654 x1 = 0.5*(gri->x1+grf->x1);
00655 x2 = 0.5*(gri->x2+grf->x2);
00656 x3 = 0.5*(gri->x3+grf->x3);
00657 v1 = 0.5*(gri->v1+grf->v1);
00658 v2 = 0.5*(gri->v2+grf->v2);
00659 v3 = 0.5*(gri->v3+grf->v3);
00660
00661
00662 fb = Get_Force(pG, x1, x2, x3, v1, v2, v3);
00663
00664 fb.x1 = dv1 - pG->dt*fb.x1;
00665 fb.x2 = dv2 - pG->dt*fb.x2;
00666 fb.x3 = dv3 - pG->dt*fb.x3;
00667
00668
00669 Elosspar = mgr*(SQR(fb.x1)+SQR(fb.x2)+SQR(fb.x3))*ts;
00670
00671
00672 fb.x1 = mgr*fb.x1;
00673 fb.x2 = mgr*fb.x2;
00674 fb.x3 = mgr*fb.x3;
00675
00676
00677 getweight(pG, x1, x2, x3, cell1, weight, &is, &js, &ks);
00678 distrFB_corr(pG, weight, is, js, ks, fb, Elosspar);
00679
00680 return;
00681
00682 }
00683
00684 #endif
00685
00686
00687
00688
00689
00690
00691 void Delete_Ghost(Grid *pG)
00692 {
00693 long p;
00694 Grain *gr;
00695
00696 p = 0;
00697 while (p<pG->nparticle)
00698 {
00699 gr = &(pG->particle[p]);
00700
00701 if (gr->pos == 0)
00702 {
00703 pG->nparticle -= 1;
00704 pG->grproperty[gr->property].num -= 1;
00705 pG->particle[p] = pG->particle[pG->nparticle];
00706 }
00707 else
00708 p++;
00709 }
00710
00711 return;
00712 }
00713
00714
00715
00716
00717 void JudgeCrossing(Grid *pG, Real x1, Real x2, Real x3, Grain *gr)
00718 {
00719 #ifndef FARGO
00720
00721 if ((x1>=x1upar) || (x1<x1lpar) || (x2>=x2upar) || (x2<x2lpar) ||
00722 (x3>=x3upar) || (x3<x3lpar) )
00723 gr->pos = 10;
00724 #else
00725
00726
00727
00728 if ((x1>=x1upar) || (x1<x1lpar) || (x3>=x3upar) || (x3<x3lpar))
00729 gr->pos = 10;
00730 else if ((pG->Nx3==1) && ((x2>=x2upar) || (x2<x2lpar)))
00731 gr->pos = 10;
00732 #endif
00733 return;
00734 }
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749 Vector Get_Drag(Grid *pG, int type, Real x1, Real x2, Real x3,
00750 Real v1, Real v2, Real v3, Vector cell1, Real *tstop1)
00751 {
00752 int is, js, ks;
00753 Real rho, u1, u2, u3, cs;
00754 Real vd1, vd2, vd3, vd, tstop, ts1;
00755 #ifdef FEEDBACK
00756 Real stiffness;
00757 #endif
00758 Real weight[3][3][3];
00759 Vector fd;
00760
00761
00762 getweight(pG, x1, x2, x3, cell1, weight, &is, &js, &ks);
00763
00764 #ifndef FEEDBACK
00765 if (getvalues(pG, weight, is, js, ks, &rho, &u1, &u2, &u3, &cs) == 0)
00766 #else
00767 if (getvalues(pG, weight, is, js, ks, &rho,&u1,&u2,&u3,&cs, &stiffness) == 0)
00768 #endif
00769 {
00770
00771
00772 gasvshift(x1, x2, x3, &u1, &u2, &u3);
00773
00774
00775 vd1 = v1-u1;
00776 vd2 = v2-u2;
00777 vd3 = v3-u3;
00778 vd = sqrt(SQR(vd1) + SQR(vd2) + SQR(vd3));
00779
00780
00781 tstop = get_ts(pG, type, rho, cs, vd);
00782 #ifdef FEEDBACK
00783
00784 #endif
00785 ts1 = 1.0/tstop;
00786 }
00787 else
00788 {
00789 vd1 = 0.0; vd2 = 0.0; vd3 = 0.0; ts1 = 0.0;
00790 ath_perr(0, "Particle move out of grid %d with position (%f,%f,%f)!\n",
00791 pG->my_id,x1,x2,x3);
00792 }
00793
00794 *tstop1 = ts1;
00795
00796
00797 fd.x1 = -ts1*vd1;
00798 fd.x2 = -ts1*vd2;
00799 fd.x3 = -ts1*vd3;
00800
00801 return fd;
00802 }
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815 Vector Get_Force(Grid *pG, Real x1, Real x2, Real x3,
00816 Real v1, Real v2, Real v3)
00817 {
00818 Vector ft;
00819
00820 ft.x1 = ft.x2 = ft.x3 = 0.0;
00821
00822
00823
00824
00825
00826 Userforce_particle(&ft, x1, x2, x3, v1, v2, v3);
00827
00828 #ifdef SHEARING_BOX
00829 Real omg2 = SQR(Omega_0);
00830
00831 if (pG->Nx3 > 1)
00832 {
00833 #ifdef FARGO
00834 ft.x1 += 2.0*v2*Omega_0;
00835 ft.x2 += (qshear-2.0)*v1*Omega_0;
00836 #else
00837 ft.x1 += 2.0*(qshear*omg2*x1 + v2*Omega_0);
00838 ft.x2 += -2.0*v1*Omega_0;
00839 #endif
00840 }
00841 else
00842 {
00843 #ifdef FARGO
00844 ft.x1 += 2.0*v3*Omega_0;
00845 ft.x3 += (qshear-2.0)*v1*Omega_0;
00846 #else
00847 ft.x1 += 2.0*(qshear*omg2*x1 + v3*Omega_0);
00848 ft.x3 += -2.0*v1*Omega_0;
00849 #endif
00850 }
00851 #endif
00852
00853 return ft;
00854 }
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867 Vector Get_Term(Grid *pG, int type, Real x1, Real x2, Real x3, Vector cell1,
00868 Real *tstop)
00869 {
00870 Vector vterm;
00871 Vector ft;
00872 Real rho, u1, u2, u3, cs;
00873 int is, js, ks;
00874 #ifdef FEEDBACK
00875 Real stiffness;
00876 #endif
00877 Real weight[3][3][3];
00878
00879
00880 getweight(pG, x1, x2, x3, cell1, weight, &is, &js, &ks);
00881
00882 #ifndef FEEDBACK
00883 if (getvalues(pG, weight, is, js, ks, &rho, &u1, &u2, &u3, &cs) == 0)
00884 #else
00885 if (getvalues(pG, weight, is, js, ks, &rho,&u1,&u2,&u3,&cs, &stiffness) == 0)
00886 #endif
00887 {
00888
00889
00890 gasvshift(x1, x2, x3, &u1, &u2, &u3);
00891
00892
00893 *tstop = get_ts(pG, type, rho, cs, 0.0);
00894
00895
00896 ft = Get_ForceDiff(pG, x1, x2, x3, u1, u2, u3);
00897
00898
00899 vterm.x1 = u1 + *tstop*ft.x1;
00900 vterm.x2 = u2 + *tstop*ft.x2;
00901 vterm.x3 = u3 + *tstop*ft.x3;
00902 }
00903 else
00904 {
00905 vterm.x1 = 0.0; vterm.x2 = 0.0; vterm.x3 = 0.0;
00906 ath_perr(0, "[get_term]: Position (%f,%f,%f) is out of grid!\n",
00907 pG->my_id,x1,x2,x3);
00908 }
00909
00910 return vterm;
00911 }
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928 Vector Get_ForceDiff(Grid *pG, Real x1, Real x2, Real x3,
00929 Real v1, Real v2, Real v3)
00930 {
00931 Vector fd;
00932
00933 fd.x1 = 0.0; fd.x2 = 0.0; fd.x3 = 0.0;
00934
00935 Userforce_particle(&fd, x1, x2, x3, v1, v2, v3);
00936
00937
00938
00939
00940
00941
00942
00943 return fd;
00944 }
00945
00946 #endif