Athena: particles/integrators_particle.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*===========================================================================*/
00003 /*! \file integrators_particle.c
00004  *  \brief Provide three kinds of particle integrators.
00005  *
00006  * PURPOSE: provide three kinds of particle integrators, namely, 2nd order
00007  *   explicit, 2nd order semi-implicit and 2nd order fully implicit.
00008  * 
00009  * CONTAINS PUBLIC FUNCTIONS:
00010  * - Integrate_Particles();
00011  * - int_par_exp   ()
00012  * - int_par_semimp()
00013  * - int_par_fulimp()
00014  * - feedback_predictor()
00015  * - feedback_corrector()
00016  *
00017  * PRIVATE FUNCTION PROTOTYPES:
00018  * - Delete_Ghost()   - delete ghost particles
00019  * - JudgeCrossing()  - judge if the particle cross the grid boundary
00020  * - Get_Drag()       - calculate the drag force
00021  * - Get_Force()      - calculate forces other than the drag
00022  * - Get_Term()       - calculate the termination particle velocity
00023  * - Get_ForceDiff()  - calculate the force difference between particle and gas
00024  * 
00025  * History:
00026  * - Written by Xuening Bai, Mar.2009                                         */
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         /* endif at the end of the file */
00039 
00040 /*==============================================================================
00041  * PRIVATE FUNCTION PROTOTYPES:
00042  *   Delete_Ghost()   - delete ghost particles
00043  *   JudgeCrossing()  - judge if the particle cross the grid boundary
00044  *   Get_Drag()       - calculate the drag force
00045  *   Get_Force()      - calculate forces other than the drag
00046  *   Get_Term()       - calculate the termination particle velocity
00047  *   Get_ForceDiff()  - calculate the force difference between particle and gas
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 /*=========================== PUBLIC FUNCTIONS ===============================*/
00061 /*----------------------------------------------------------------------------*/
00062 
00063 /*---------------------------------- Main Integrator -------------------------*/
00064 /*! \fn void Integrate_Particles(Grid *pG, Domain *pD)
00065  *  \brief Main particle integrator.
00066  *
00067  * Input: Grid which is already evolved in half a time step. Paricles unevolved.
00068  * Output: particle updated for one full time step; feedback array for corrector
00069  *         updated.
00070  * Note: This routine allows has the flexibility to freely choose the particle
00071  *       integrator.
00072  * Should use fully implicit integrator for tightly coupoled particles.
00073  * Otherwise the semi-implicit integrator performs better.
00074  */
00075 void Integrate_Particles(Grid *pG, Domain *pD)
00076 {
00077   Grain *curG, *curP, mygr;     /* pointer of the current working position */
00078   long p;                       /* particle index */
00079   Real dv1, dv2, dv3, ts, t1;   /* amount of velocity update, stopping time */
00080   Vector cell1;                 /* one over dx1, dx2, dx3 */
00081   Vector vterm;                 /* termination velocity */
00082 
00083 /* Initialization */
00084 #ifdef FEEDBACK
00085   feedback_clear(pG);   /* clean the feedback array */
00086 #endif /* FEEDBACK */
00087 
00088   curP = &(mygr);       /* temperory particle */
00089 
00090   /* cell1 is a shortcut expressions as well as dimension indicator */
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   /* delete all ghost particles */
00096   Delete_Ghost(pG);
00097 
00098   p = 0;
00099   while (p<pG->nparticle)
00100   {/* loop over all particles */
00101     curG = &(pG->particle[p]);
00102 
00103 /* Step 1: Calculate velocity update */
00104     switch(pG->grproperty[curG->property].integrator)
00105     {
00106       case 1: /* 2nd order explicit integrator */
00107         int_par_exp(pG,curG,cell1, &dv1,&dv2,&dv3, &ts);
00108         break;
00109 
00110       case 2: /* 2nd order semi-implicit integrator */
00111         int_par_semimp(pG,curG,cell1, &dv1,&dv2,&dv3, &ts);
00112         break;
00113 
00114       case 3: /* 2nd order fully implicit integrator */
00115         int_par_fulimp(pG,curG,cell1, &dv1,&dv2,&dv3, &ts);
00116         break;
00117 
00118       case 4: /* specific integrator in the strong coupling limit */
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 /* Step 2: particle update to curP */
00127 
00128     /* velocity update */
00129     curP->v1 = curG->v1 + dv1;
00130     curP->v2 = curG->v2 + dv2;
00131     curP->v3 = curG->v3 + dv3;
00132 
00133     /* position update */
00134     if (pG->Nx1 > 1)
00135       curP->x1 = curG->x1 + 0.5*pG->dt*(curG->v1 + curP->v1);
00136     else /* do not move if this dimension collapses */
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 /* do not move if this dimension collapses */
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 /* do not move if this dimension collapses */
00147       curP->x3 = curG->x3;
00148 
00149 #ifdef FARGO
00150     /* shift = -qshear * Omega_0 * x * dt */
00151     pG->parsub[p].shift = -0.5*qshear*Omega_0*(curG->x1+curP->x1)*pG->dt;
00152 #endif
00153 
00154     /* special treatment for integrator #4 */
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 /* Step 3: calculate feedback force to the gas */
00164 #ifdef FEEDBACK
00165     feedback_corrector(pG, curG, curP, cell1, dv1, dv2, dv3, ts);
00166 #endif /* FEEDBACK */
00167 
00168 /* Step 4: Final update of the particle */
00169     /* update particle status (crossing boundary or not) */
00170     JudgeCrossing(pG, curP->x1, curP->x2, curP->x3, curG);
00171 
00172     /* update the particle */
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   } /* end of the for loop */
00182 
00183 /* output the status */
00184   ath_pout(0, "In processor %d, there are %ld particles.\n",
00185                              pG->my_id, pG->nparticle);
00186 
00187   return;
00188 }
00189 
00190 /* ------------ 2nd order fully implicit particle integrator -----------------*/
00191 /*! \fn void int_par_fulimp(Grid *pG, Grain *curG, Vector cell1, 
00192  *                            Real *dv1, Real *dv2, Real *dv3, Real *ts)
00193  *  \brief 2nd order fully implicit particle integrator
00194  *
00195  * Input: 
00196  *   grid pointer (pG), grain pointer (curG), cell size indicator (cell1)
00197  * Output:
00198  *   dv1,dv2,dv3: velocity update
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;   /* first order new position at half a time step */
00204   Vector fd, fr;        /* drag force and other forces */
00205   Vector fc, fp, ft;    /* force at current & predicted position, total force */
00206   Real ts11, ts12;      /* 1/stopping time */
00207   Real b0,A,B,C,D,Det1; /* matrix elements and determinant */
00208 #ifdef SHEARING_BOX
00209   Real oh, oh2;         /* Omega_0*dt and its square */
00210 #endif
00211 
00212 /* step 1: predict of the particle position after one time step */
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); /* advection part */
00223 #endif
00224 #endif
00225 
00226 /* step 2: calculate the force at current position */
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 /* step 3: calculate the force at the predicted positoin */
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 /* step 4: calculate the velocity update */
00248   /* shortcut expressions */
00249   b0 = 1.0+pG->dt*ts11;
00250 
00251   /* Total force */
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) {/* 3D shearing sheet (x1,x2,x3)=(X,Y,Z) */
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 {         /* 2D shearing sheet (x1,x2,x3)=(X,Z,Y) */
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 /* SHEARING_BOX */
00274 
00275   /* calculate the inverse matrix elements */
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 /* FARGO */
00284   A = D - 2.0*oh2;
00285   C = -B;
00286 #endif /* FARGO */
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 /* SHEARING_BOX */
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 /* SHEARING_BOX */
00303 
00304   *ts = 0.5/ts11+0.5/ts12;
00305 
00306   return;
00307 }
00308 
00309 
00310 /*--------------- 2nd order semi-implicit particle integrator ----------------*/
00311 /*! \fn void int_par_semimp(Grid *pG, Grain *curG, Vector cell1, 
00312  *                            Real *dv1, Real *dv2, Real *dv3, Real *ts)
00313  *  \brief 2nd order semi-implicit particle integrator 
00314  *
00315  * Input: 
00316  *   grid pointer (pG), grain pointer (curG), cell size indicator (cell1)
00317  * Output:
00318  *   dv1,dv2,dv3: velocity update
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;    /* drag force and other forces, total force */
00324   Real ts1, b, b2;      /* other shortcut expressions */
00325   Real x1n, x2n, x3n;   /* first order new position at half a time step */
00326 #ifdef SHEARING_BOX
00327   Real b1, oh;          /* Omega_0*h */
00328 #endif
00329 
00330 /* step 1: predict of the particle position after half a time step */
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);/* advection part */
00341 #endif
00342 #endif
00343 
00344 /* Step 2: interpolation to get fluid density, velocity and the sound speed at\  * predicted position
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 /* step 3: calculate velocity update */
00356 
00357   /* shortcut expressions */
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 /* FARGO */
00366   b2 = b*b1;
00367 #else
00368   b2 = 1.0/b;
00369 #endif /* SHEARING BOX */
00370 
00371     /* velocity evolution */
00372 #ifdef SHEARING_BOX
00373   if (pG->Nx3>1)
00374   {/* 3D shearing sheet (x1,x2,x3)=(X,Y,Z) */
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 /* FARGO */
00381     *dv3 = pG->dt*2.0*ft.x3/b;
00382   }
00383   else
00384   {/* 2D shearing sheet (x1,x2,x3)=(X,Z,Y) */
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 /* SHEARING_BOX */
00398 
00399   *ts = 1.0/ts1;
00400 
00401   return;
00402 }
00403 
00404 
00405 /*------------------- 2nd order explicit particle integrator -----------------*/
00406 /*! \fn void int_par_exp(Grid *pG, Grain *curG, Vector cell1,
00407  *                         Real *dv1, Real *dv2, Real *dv3, Real *ts)
00408  *  \brief 2nd order explicit particle integrator 
00409  *
00410  * Input: 
00411  *   grid pointer (pG), grain pointer (curG), cell size indicator (cell1)
00412  * Output:
00413  *   dv1,dv2,dv3: velocity update
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;    /* drag force and other forces, total force */
00419   Real ts1;             /* 1/stopping time */
00420   Real x1n, x2n, x3n;   /* first order new position at half a time step */
00421   Real v1n, v2n, v3n;   /* first order new velocity at half a time step */
00422 
00423 /* step 1: predict of the particle position after half a time step */
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);/* advection part */
00437 #endif
00438 #endif
00439 
00440 /* step 2: calculate the force at current position */
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 /* step 3: calculate the force at the predicted positoin */
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 /* step 4: calculate velocity update */
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 /*------------------- 2nd order specific particle integrator -----------------*/
00475 /*! \fn void int_par_spec(Grid *pG, Grain *curG, Vector cell1,
00476  *                          Real *dv1, Real *dv2, Real *dv3, Real *ts)
00477  *  \brief 2nd order specific particle integrator;
00478  *  This integrator works ONLY in the strong coupling regime (t_stop<h)
00479  *
00480  * Input:
00481  *   grid pointer (pG), grain pointer (curG), cell size indicator (cell1)
00482  * Output:
00483  *   dv1,dv2,dv3: velocity update
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;         /* termination velocity */
00489   Real x1n, x2n, x3n;   /* predicted position at half a time step */
00490 
00491 /* step 1: predict of the particle position after half a time step */
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);/* advection part */
00505 #endif
00506 #endif
00507 
00508 /* step 2: get gas termination velocity */
00509   vterm = Get_Term(pG, curG->property, x1n, x2n, x3n, cell1, ts);
00510 
00511 /* step 3: calculate the velocity difference */
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 /*! \fn void feedback_predictor(Grid* pG)
00522  *  \brief Calculate the feedback of the drag force from the particle to the gas
00523  *
00524  * Serves for the predictor step. It deals with all the particles.
00525  * Input: pG: grid with particles
00526  * Output: pG: the array of drag forces exerted by the particle is updated
00527 */
00528 void feedback_predictor(Grid* pG)
00529 {
00530   int is,js,ks,i,j,k;
00531   long p;                   /* particle index */
00532   Real weight[3][3][3];     /* weight function */
00533   Real rho, cs, tstop;      /* density, sound speed, stopping time */
00534   Real u1, u2, u3;
00535   Real vd1, vd2, vd3, vd;   /* velocity difference between particle and gas */
00536   Real f1, f2, f3;          /* feedback force */
00537   Real m, ts1h;             /* grain mass, 0.5*dt/tstop */
00538   Vector cell1;             /* one over dx1, dx2, dx3 */
00539   Vector fb;                /* drag force, fluid velocity */
00540 #ifndef BAROTROPIC
00541   Real Elosspar;            /* energy dissipation rate due to drag */
00542 #endif
00543   Real stiffness;           /* stiffness parameter of feedback */
00544   Grain *cur;               /* pointer of the current working position */
00545 
00546   /* initialization */
00547   get_gasinfo(pG);              /* calculate gas information */
00548 
00549   for (k=klp; k<=kup; k++)
00550     for (j=jlp; j<=jup; j++)
00551       for (i=ilp; i<=iup; i++) {
00552         /* clean the feedback array */
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   /* convenient expressions */
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   /* loop over all particles to calculate the drag force */
00573   for (p=0; p<pG->nparticle; p++)
00574   {/* loop over all particle */
00575     cur = &(pG->particle[p]);
00576 
00577     /* interpolation to get fluid density and velocity */
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     { /* particle is in the grid */
00582 
00583       /* apply gas velocity shift due to pressure gradient */
00584       gasvshift(cur->x1, cur->x2, cur->x3, &u1, &u2, &u3);
00585       /* velocity difference */
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       /* calculate particle stopping time */
00592       tstop = get_ts(pG, cur->property, rho, cs, vd);
00593       ts1h = 0.5*pG->dt/tstop;
00594 
00595       /* Drag force density */
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       /* calculate feedback stiffness */
00602        stiffness = 2.0*m*ts1h;
00603 
00604 #ifndef BAROTROPIC
00605       Elosspar = fb.x1*vd1 + fb.x2*vd2 + fb.x3*vd3;
00606       /* distribute the drag force (density) to the grid */
00607       distrFB_pred(pG, weight, is, js, ks, fb, stiffness, Elosspar);
00608 #else
00609       /* distribute the drag force (density) to the grid */
00610       distrFB_pred(pG, weight, is, js, ks, fb, stiffness);
00611 #endif
00612     }
00613   }/* end of the for loop */
00614 
00615 /* normalize stiffness and correct for feedback */
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 //        stiffness = 1.0/MAX(1.0, pG->Coup[k][j][i].FBstiff);
00623 
00624 //        pG->Coup[k][j][i].fb1 *= stiffness;
00625 //        pG->Coup[k][j][i].fb2 *= stiffness;
00626 //        pG->Coup[k][j][i].fb3 *= stiffness;
00627       }
00628 
00629   return;
00630 }
00631 
00632 /*----------------------------------------------------------------------------*/
00633 /*! \fn void feedback_corrector(Grid *pG, Grain *gri, Grain *grf, Vector cell1,
00634  *                                Real dv1, Real dv2, Real dv3, Real ts)
00635  *  \brief  Calculate the feedback of the drag force from the particle 
00636  *          to the gas.
00637  *
00638  * Serves for the corrector step. It deals with one particle at a time.
00639  * Input: pG: grid with particles; gri,grf: initial and final particles;
00640  *        dv: velocity difference between gri and grf.
00641  * Output: pG: the array of drag forces exerted by the particle is updated
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;                        /* particle energy dissipation */
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   /* Force other than drag force */
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   /* energy dissipation */
00669   Elosspar = mgr*(SQR(fb.x1)+SQR(fb.x2)+SQR(fb.x3))*ts;
00670 
00671   /* Drag force density */
00672   fb.x1 = mgr*fb.x1;
00673   fb.x2 = mgr*fb.x2;
00674   fb.x3 = mgr*fb.x3;
00675 
00676   /* distribute the drag force (density) to the grid */
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 /* FEEDBACK */
00685 
00686 
00687 /*=========================== PRIVATE FUNCTIONS ==============================*/
00688 /*----------------------------------------------------------------------------*/
00689 /*! \fn void Delete_Ghost(Grid *pG)
00690  *  \brief Delete ghost particles */
00691 void Delete_Ghost(Grid *pG)
00692 {
00693   long p;
00694   Grain *gr;
00695 
00696   p = 0;
00697   while (p<pG->nparticle)
00698   {/* loop over all particles */
00699     gr = &(pG->particle[p]);
00700 
00701     if (gr->pos == 0)
00702     {/* gr is a ghost particle */
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 /*! \fn void JudgeCrossing(Grid *pG, Real x1, Real x2, Real x3, Grain *gr)
00716  *  \brief Judge if the particle is a crossing particle */
00717 void JudgeCrossing(Grid *pG, Real x1, Real x2, Real x3, Grain *gr)
00718 {
00719 #ifndef FARGO
00720     /* if it crosses the grid boundary, mark it as a crossing out particle */
00721     if ((x1>=x1upar) || (x1<x1lpar) || (x2>=x2upar) || (x2<x2lpar) ||
00722                                        (x3>=x3upar) || (x3<x3lpar) )
00723       gr->pos = 10;
00724 #else
00725     /* FARGO will naturally return the "crossing out" particles in the x2
00726         direction to the grid
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))) /* 2D problem */
00731       gr->pos = 10;
00732 #endif
00733     return;
00734 }
00735 
00736 /*--------------------------------------------------------------------------- */
00737 /*! \fn Vector Get_Drag(Grid *pG, int type, Real x1, Real x2, Real x3,
00738  *              Real v1, Real v2, Real v3, Vector cell1, Real *tstop1)
00739  *  \brief Calculate the drag force to the particles 
00740  *
00741  * Input:
00742  *   pG: grid;  type: particle type;    cell1: 1/dx1,1/dx2,1/dx3;
00743  *   x1,x2,x3,v1,v2,v3: particle position and velocity;
00744  * Output:
00745  *   tstop1: 1/stopping time;
00746  * Return:
00747  *   drag force;
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];         /* weight function */
00759   Vector fd;
00760 
00761   /* interpolation to get fluid density, velocity and the sound speed */
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   { /* particle in the grid */
00770 
00771     /* apply possible gas velocity shift (e.g., for fake gas velocity field) */
00772     gasvshift(x1, x2, x3, &u1, &u2, &u3);
00773 
00774     /* velocity difference */
00775     vd1 = v1-u1;
00776     vd2 = v2-u2;
00777     vd3 = v3-u3;
00778     vd = sqrt(SQR(vd1) + SQR(vd2) + SQR(vd3)); /* dimension independent */
00779 
00780     /* particle stopping time */
00781     tstop = get_ts(pG, type, rho, cs, vd);
00782 #ifdef FEEDBACK
00783 //    tstop *= MAX(1.0,stiffness);
00784 #endif
00785     ts1 = 1.0/tstop;
00786   }
00787   else
00788   { /* particle out of the grid, free motion, with warning sign */
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); /* warning! */
00792   }
00793 
00794   *tstop1 = ts1;
00795 
00796   /* Drag force */
00797   fd.x1 = -ts1*vd1;
00798   fd.x2 = -ts1*vd2;
00799   fd.x3 = -ts1*vd3;
00800 
00801   return fd;
00802 }
00803 
00804 /*--------------------------------------------------------------------------- */
00805 /*! \fn Vector Get_Force(Grid *pG, Real x1, Real x2, Real x3,
00806  *                         Real v1, Real v2, Real v3)
00807  *  \brief Calculate the forces to the particle other than the gas drag
00808  *
00809  * Input:
00810  *   pG: grid;
00811  *   x1,x2,x3,v1,v2,v3: particle position and velocity;
00812  * Return:
00813  *   forces;
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 /* User defined forces
00823  * Should be independent of velocity, or the integrators must be modified
00824  * Can also change the velocity to account for velocity difference.
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   {/* 3D shearing sheet (x1,x2,x3)=(X,Y,Z) */
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 /* FARGO */
00840   }
00841   else
00842   { /* 2D shearing sheet (x1,x2,x3)=(X,Z,Y) */
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 /* FARGO */
00850   }
00851 #endif /* SHEARING_BOX */
00852 
00853   return ft;
00854 }
00855 
00856 /*! \fn Vector Get_Term(Grid *pG, int type, Real x1, Real x2, Real x3, 
00857  *                      Vector cell1, Real *tstop)
00858  *  \brief Calculate the termination velocity of strongly coupled particles
00859  *
00860  * Used for the special integrator
00861  * Force difference include pressure gradient and momentum feedback
00862  * Input:
00863  *   pG: grid;  type: particle type;  x1,x2,x3: particle position;
00864  * Return:
00865  *   termination velocity, and the stopping time.
00866  */
00867 Vector Get_Term(Grid *pG, int type, Real x1, Real x2, Real x3, Vector cell1,
00868                                                        Real *tstop)
00869 {
00870   Vector vterm;             /* termination velocity */
00871   Vector ft;                /* force difference */
00872   Real rho, u1, u2, u3, cs; /* gas velocity */
00873   int is, js, ks;
00874 #ifdef FEEDBACK
00875   Real stiffness;
00876 #endif
00877   Real weight[3][3][3];     /* weight function */
00878 
00879   /* interpolation to get fluid density, velocity and the sound speed */
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   { /* position in the grid */
00888 
00889     /* apply possible gas velocity shift (e.g., for fake gas velocity field) */
00890     gasvshift(x1, x2, x3, &u1, &u2, &u3);
00891 
00892     /* particle stopping time */
00893     *tstop = get_ts(pG, type, rho, cs, 0.0);
00894 
00895     /* force difference */
00896     ft = Get_ForceDiff(pG, x1, x2, x3, u1, u2, u3);
00897 
00898     /* termination velocity */
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); /* warning! */
00908   }
00909 
00910   return vterm;
00911 }
00912 
00913 /*! \fn Vector Get_ForceDiff(Grid *pG, Real x1, Real x2, Real x3,
00914  *                             Real v1, Real v2, Real v3)
00915  *  \brief Calculate the force (density) difference between particles and gas
00916  *
00917  * Used ONLY for the special integrator.
00918  * THIS ROUTINE MUST BE EDITTED BY THE USER!
00919  * Force differences due to gas pressure gradient and momentum feedback are
00920  * automatically included. The user must provide other sources.
00921  *
00922  * Input:
00923  *   pG: grid;
00924  *   x1,x2,x3,v1,v2,v3: particle position and velocity;
00925  * Return:
00926  *   forces;
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   fd.x1 += x1;
00939   fd.x2 += x2;
00940   fd.x3 += x3;
00941 */
00942 
00943   return fd;
00944 }
00945 
00946 #endif /*PARTICLES*/