Athena: particles/utils_particle.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file utils_particle.c
00004  *  \brief Contains most of the utilities for the particle code.
00005  *
00006  * PURPOSE: Contains most of the utilities for the particle code: all the
00007  *   interpolation functions, stopping time calculation functions, shuffle
00008  *   algorithms. Also contained are the default (and trivial) gas velocity
00009  *   shift function. The get_gasinfo(Grid *pG) routine is used for test
00010  *   purposes only.
00011  * 
00012  * CONTAINS PUBLIC FUNCTIONS:
00013  * 
00014  * - getwei_linear()
00015  * - getwei_TSC   ()
00016  * - getwei_QP    ()
00017  * - getvalues()
00018  * - get_ts_epstein()
00019  * - get_ts_general()
00020  * - get_ts_fixed  ()
00021  * - get_gasinfo()
00022  * - feedback_clear()
00023  * - distrFB      ()
00024  * - void shuffle()
00025  * - void gasvshift_zero()
00026  * 
00027  * PRIVATE FUNCTION PROTOTYPES:
00028  * - compare_gr()         - compare the location of the two particles
00029  * - quicksort_particle() - sort the particles using the quicksort
00030  *
00031  * History:
00032  * - Written by Xuening Bai, Apr. 2009                                        */
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         /* endif at the end of the file */
00045 
00046 /*==============================================================================
00047  * PRIVATE FUNCTION PROTOTYPES:
00048  *   compare_gr()         - compare the location of the two particles
00049  *   quicksort_particle() - sort the particles using the quicksort
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 /*============================== ALL FUNCTIONS ===============================*/
00056 /*------------------------------------------------------------------------------
00057  * interpolation functions;
00058  * stopping time functions;
00059  * feedback related functions;
00060  * shuffle related functions;
00061  *----------------------------------------------------------------------------*/
00062 
00063 
00064 /*============================================================================*/
00065 /*-----------------------------INTERPOLATION----------------------------------
00066  *
00067  * getwei_linear()
00068  * getwei_TSC()
00069  * getwei_QP ()
00070  * getvalues();
00071  */
00072 /*============================================================================*/
00073 
00074 /*----------------------------------------------------------------------------*/
00075 /*! \fn void getwei_linear(Grid *pG, Real x1, Real x2, Real x3, Vector cell1,
00076  *                           Real weight[3][3][3], int *is, int *js, int *ks)
00077  *  \brief Get weight using linear interpolation
00078  *
00079  * Input: pG: grid; x1,x2,x3: global coordinate; cell1: 1 over dx1,dx2,dx3
00080  * Output: weight: weight function; is,js,ks: starting cell indices in the grid.
00081  * Note: this interpolation works in any 1-3 dimensions.
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;         /* grid coordinate for the position (x1,x2,x3) */
00089   Real wei1[2], wei2[2], wei3[2];/* weight function in x1,x2,x3 directions */
00090 
00091   /* find cell locations and calculate 1D weight */
00092   /* x1 direction */
00093   if (cell1.x1 > 0.0) {
00094     i = celli(pG, x1, cell1.x1, &i1, &a);       /* x1 index */
00095     i1 = i+i1-1;        *is = i1;               /* starting x1 index */
00096     wei1[1] = a - i1 - 0.5;                     /* one direction weight */
00097     wei1[0] = 1.0 - wei1[1];                    /* 0: left; 1: right */
00098   }
00099   else { /* x1 dimension collapses */
00100     *is = pG->is;
00101     wei1[1] = 0.0;
00102     wei1[0] = 1.0;
00103   }
00104 
00105   /* x2 direction */
00106   if (cell1.x2 > 0.0) {
00107     j = cellj(pG, x2, cell1.x2, &j1, &b);       /* x2 index */
00108     j1 = j+j1-1;        *js = j1;               /* starting x2 index */
00109     wei2[1] = b - j1 - 0.5;                     /* one direction weight */
00110     wei2[0] = 1.0 - wei2[1];                    /* 0: left; 1: right */
00111   }
00112   else { /* x2 dimension collapses */
00113     *js = pG->js;
00114     wei2[1] = 0.0;
00115     wei2[0] = 1.0;
00116   }
00117 
00118   /* x3 direction */
00119   if (cell1.x3 > 0.0) {
00120     k = cellk(pG, x3, cell1.x3, &k1, &c);       /* x3 index */
00121     k1 = k+k1-1;        *ks = k1;               /* starting x3 index */
00122     wei3[1] = c - k1 - 0.5;                     /* one direction weight */
00123     wei3[0] = 1.0 - wei3[1];                    /* 0: left; 1: right */
00124   }
00125   else { /* x3 dimension collapses */
00126     *ks = pG->ks;
00127     wei3[1] = 0.0;
00128     wei3[0] = 1.0;
00129   }
00130 
00131   /* calculate 3D weight */
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 /*! \fn void getwei_TSC(Grid *pG, Real x1, Real x2, Real x3, Vector cell1,
00142  *                        Real weight[3][3][3], int *is, int *js, int *ks)
00143  *
00144  *  \brief Get weight using Triangular Shaped Cloud (TSC) interpolation 
00145  * Input: pG: grid; x1,x2,x3: global coordinate; cell1: 1 over dx1,dx2,dx3
00146  * Output: weight: weight function; is,js,ks: starting cell indices in the grid.
00147  * Note: this interpolation works in any 1-3 dimensions.
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;      /* grid coordinate for the position (x1,x2,x3) */
00154   Real wei1[3], wei2[3], wei3[3];/* weight function in x1,x2,x3 directions */
00155 
00156   /* find cell locations and calculate 1D weight */
00157   /* x1 direction */
00158   if (cell1.x1 > 0.0) {
00159     celli(pG, x1, cell1.x1, &i, &a);            /* x1 index */
00160     *is = i - 1;                                /* starting x1 index, wei[0] */
00161     d = a - i;
00162     wei1[0] = 0.5*SQR(1.0-d);                   /* 0: left; 2: right */
00163     wei1[1] = 0.75-SQR(d-0.5);                  /* one direction weight */
00164     wei1[2] = 0.5*SQR(d);
00165   }
00166   else { /* x1 dimension collapses */
00167     *is = pG->is;
00168     wei1[1] = 0.0;      wei1[2] = 0.0;
00169     wei1[0] = 1.0;
00170   }
00171 
00172   /* x2 direction */
00173   if (cell1.x2 > 0.0) {
00174     cellj(pG, x2, cell1.x2, &j, &b);            /* x2 index */
00175     *js = j - 1;                                /* starting x2 index */
00176     d = b - j;
00177     wei2[0] = 0.5*SQR(1.0-d);                   /* 0: left; 2: right */
00178     wei2[1] = 0.75-SQR(d-0.5);                  /* one direction weight */
00179     wei2[2] = 0.5*SQR(d);
00180   }
00181   else { /* x2 dimension collapses */
00182     *js = pG->js;
00183     wei2[1] = 0.0;      wei2[2] = 0.0;
00184     wei2[0] = 1.0;
00185   }
00186 
00187   /* x3 direction */
00188   if (cell1.x3 > 0.0) {
00189     cellk(pG, x3, cell1.x3, &k, &c);            /* x3 index */
00190     *ks = k - 1;                                /* starting x3 index */
00191     d = c - k;
00192     wei3[0] = 0.5*SQR(1.0-d);                   /* 0: left; 2: right */
00193     wei3[1] = 0.75-SQR(d-0.5);                  /* one direction weight */
00194     wei3[2] = 0.5*SQR(d);
00195   }
00196   else { /* x3 dimension collapses */
00197     *ks = pG->ks;
00198     wei3[1] = 0.0;      wei3[2] = 0.0;
00199     wei3[0] = 1.0;
00200   }
00201 
00202   /* calculate 3D weight */
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 /*! \fn void getwei_QP(Grid *pG, Real x1, Real x2, Real x3, Vector cell1,
00214  *                       Real weight[3][3][3], int *is, int *js, int *ks)
00215  *  \brief Get weight using quadratic polynomial interpolation 
00216  *
00217  * Input: pG: grid; x1,x2,x3: global coordinate; cell1: 1 over dx1,dx2,dx3
00218  * Output: weight: weight function; is,js,ks: starting cell indices in the grid.
00219  * Note: this interpolation works in any 1-3 dimensions.
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;      /* grid coordinate for the position (x1,x2,x3) */
00226   Real wei1[3], wei2[3], wei3[3];/* weight function in x1,x2,x3 directions */
00227 
00228   /* find cell locations and calculate 1D weight */
00229   /* x1 direction */
00230   if (cell1.x1 > 0.0) {
00231     celli(pG, x1, cell1.x1, &i, &a);            /* x1 index */
00232     *is = i - 1;                                /* starting x1 index, wei[0] */
00233     d = a - i;
00234     wei1[0] = 0.5*(0.5-d)*(1.5-d);              /* 0: left; 2: right */
00235     wei1[1] = 1.0-SQR(d-0.5);                   /* one direction weight */
00236     wei1[2] = 0.5*(d-0.5)*(d+0.5);
00237   }
00238   else { /* x1 dimension collapses */
00239     *is = pG->is;
00240     wei1[1] = 0.0;      wei1[2] = 0.0;
00241     wei1[0] = 1.0;
00242   }
00243 
00244   /* x2 direction */
00245   if (cell1.x2 > 0.0) {
00246     cellj(pG, x2, cell1.x2, &j, &b);            /* x2 index */
00247     *js = j - 1;                                /* starting x2 index */
00248     d = b - j;
00249     wei2[0] = 0.5*(0.5-d)*(1.5-d);              /* 0: left; 2: right */
00250     wei2[1] = 1.0-SQR(d-0.5);                   /* one direction weight */
00251     wei2[2] = 0.5*(d-0.5)*(d+0.5);
00252   }
00253   else { /* x2 dimension collapses */
00254     *js = pG->js;
00255     wei2[1] = 0.0;      wei2[2] = 0.0;
00256     wei2[0] = 1.0;
00257   }
00258 
00259   /* x3 direction */
00260   if (cell1.x3 > 0.0) {
00261     cellk(pG, x3, cell1.x3, &k, &c);            /* x3 index */
00262     *ks = k - 1;                                /* starting x3 index */
00263     d = c - k;
00264     wei3[0] = 0.5*(0.5-d)*(1.5-d);              /* 0: left; 2: right */
00265     wei3[1] = 1.0-SQR(d-0.5);                   /* one direction weight */
00266     wei3[2] = 0.5*(d-0.5)*(d+0.5);
00267   }
00268   else { /* x3 dimension collapses */
00269     *ks = pG->ks;
00270     wei3[1] = 0.0;      wei3[2] = 0.0;
00271     wei3[0] = 1.0;
00272   }
00273 
00274   /* calculate 3D weight */
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 /*! \fn int getvalues()
00285  *  \brief Get interpolated value using the weight
00286  *
00287  * Input:
00288  * - pG: grid; weight: weight function;
00289  * - is,js,ks: starting cell indices in the grid.
00290  * Output:
00291  * - interpolated values of density, velocity and sound speed of the fluid
00292  *
00293  * Return: 0: normal exit; -1: particle lie out of the grid, cannot interpolate
00294  * Note: this interpolation works in any 1-3 dimensions.
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;           /* density and velocity of the fluid */
00305   GPCouple *pq;
00306 #ifndef ISOTHERMAL
00307   Real C = 0.0;                 /* fluid sound speed */
00308 #endif
00309 #ifdef FEEDBACK
00310   Real stiffness=0.0;
00311 #endif
00312   Real totwei, totwei1;         /* total weight (in case of edge cells) */
00313 
00314   /* linear interpolation */
00315   D = 0.0; v1 = 0.0; v2 = 0.0; v3 = 0.0;
00316   totwei = 0.0;         totwei1 = 1.0;
00317 
00318   /* Interpolate density, velocity and sound speed */
00319   /* Note: in lower dimensions only wei[0] is non-zero */
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) /* particle lies out of the grid, warning! */
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 /* ISOTHERMAL */
00358 #ifdef FEEDBACK
00359   *stiff = stiffness*totwei1;
00360 #endif
00361 
00362   return 0;
00363 }
00364 
00365 
00366 /*============================================================================*/
00367 /*-----------------------------STOPPING TIME----------------------------------
00368  *
00369  * get_ts_general()
00370  * get_ts_epstein()
00371  * get_ts_fixed()
00372  */
00373 /*============================================================================*/
00374 
00375 /*----------------------------------------------------------------------------*/
00376 /*! \fn Real get_ts_general(Grid *pG, int type, Real rho, Real cs, Real vd)
00377  *  \brief Calculate the stopping time for the most general case
00378  *
00379  * The relavent scale to calculate is: 
00380  * - 1. a/lambda_m == alam
00381  * - 2. rho_s*a in normalized unit == rhoa
00382  */
00383 Real get_ts_general(Grid *pG, int type, Real rho, Real cs, Real vd)
00384 {
00385   Real tstop;   /* stopping time */
00386   Real a, rhos; /* primitive properties: size, solid density in cgs unit */
00387   Real alam;    /* a/lambda: particle size/gas mean free path */
00388   Real rhoa;    /* rhoa: rho_s*a in normalized unit (density, velocity, time) */
00389   Real Re, CD;  /* Reynolds number and drag coefficient */
00390 
00391   /* particle properties */
00392   a = pG->grproperty[type].rad;
00393   rhos = pG->grproperty[type].rho;
00394   rhoa = grrhoa[type];
00395 
00396   /* calculate particle size/gas mean free path */
00397   alam = alamcoeff * a * rho;  /* alamcoeff is defined in global.h */
00398 
00399   /* calculate the stopping time */
00400   if (alam < 2.25) {            /* Epstein regime */
00401     tstop = rhoa/(rho*cs);
00402   }
00403   else {
00404     Re = 4.0*alam*vd/cs;        /* the Reynolds number */
00405 
00406     if (Re < 1.) CD = 24.0/Re;  /* Stokes regime */
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   } /* endif */
00412 
00413   /* avoid zero stopping time */
00414   if (tstop < 1.0e-8*pG->dt)
00415     tstop = 1.0e-8*pG->dt;
00416 
00417   return tstop;
00418 }
00419 
00420 /*----------------------------------------------------------------------------*/
00421 /*! \fn Real get_ts_epstein(Grid *pG, int type, Real rho, Real cs, Real vd)
00422  *  \brief Calculate the stopping time in the Epstein regime  
00423  *
00424  * Note grrhoa == rho_s*a in normalized unit
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   /* avoid zero stopping time */
00431   if (tstop < 1.0e-8*pG->dt)
00432     tstop = 1.0e-8*pG->dt;
00433 
00434   return tstop;
00435 }
00436 
00437 /*! \fn Real get_ts_fixed(Grid *pG, int type, Real rho, Real cs, Real vd)
00438  *  \brief Return the fixed stopping time */
00439 Real get_ts_fixed(Grid *pG, int type, Real rho, Real cs, Real vd)
00440 {
00441   Real tstop = tstop0[type];
00442 
00443   /* avoid zero stopping time */
00444   if (tstop < 1.0e-8*pG->dt)
00445     tstop = 1.0e-8*pG->dt;
00446 
00447   return tstop;
00448 }
00449 
00450 
00451 /*============================================================================*/
00452 /*--------------------------------FEEDBACK------------------------------------
00453  *
00454  * get_gasinfo()
00455  * feedback_clear()
00456  * distrFB_????()
00457  */
00458 /*============================================================================*/
00459 
00460 /*----------------------------------------------------------------------------*/
00461 /*! \fn void get_gasinfo(Grid *pG)
00462  *  \brief Calculate the gas information from conserved variables for 
00463  *  feedback_predictor
00464  *
00465  * Input: pG: grid (not evolved yet).
00466  * Output: calculate 3D array grid_v/grid_cs in the grid structure.
00467  *         Calculated are gas velocity and sound speed.
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   /* get gas information */
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         /* E = P/(gamma-1) + rho*v^2/2 + B^2/2 */
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 /* MHD */
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 /* BAROTROPIC */
00505 #endif /* ISOTHERMAL */
00506       }
00507 
00508   return;
00509 }
00510 
00511 #ifdef FEEDBACK
00512 /*----------------------------------------------------------------------------*/
00513 /*! \fn void feedback_clear(Grid *pG)
00514  *  \brief clean the feedback array */
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 /*! \fn void distrFB_pred()
00537  *  \brief Distribute the feedback force to grid cells for the predict step
00538  *
00539  * Input: 
00540  * - pG: grid;   weight: weight function; 
00541  * - is,js,ks: starting cell indices in the grid.
00542  * - f1, f2, f3: feedback force from one particle.
00543  * Output:
00544  * - pG: feedback array is updated.
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   /* distribute feedback force */
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 /*! \fn void distrFB_corr(Grid *pG, Real weight[3][3][3], int is, int js,int ks,
00585  *                     Vector fb, Real Elosspar)
00586  *  \brief Distribute the feedback force to grid cells for the correct step
00587  *
00588  * Input:
00589  * - pG: grid;   weight: weight function;
00590  * - is,js,ks: starting cell indices in the grid.
00591  * - f1, f2, f3: feedback force from one particle.
00592  * Output:
00593  * - pG: feedback array is updated.
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   /* distribute feedback force */
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 /* FEEDBACK */
00627 
00628 /*============================================================================*/
00629 /*---------------------------------SHUFFLE------------------------------------
00630  *
00631  * shuffle()
00632  * compare_gr()
00633  * quicksort_particle()
00634  */
00635 /*============================================================================*/
00636 
00637 /*----------------------------------------------------------------------------*/
00638 /*! \fn void shuffle(Grid *pG)
00639  *  \brief Shuffle the particles
00640  *
00641  * Input: pG: grid with particles;
00642  * Output: pG: particles in the linked list are rearranged by the order of their
00643  *         locations that are consistent with grid cell storage.
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   /* output status */
00655   ath_pout(0, "Resorting particles...\n");
00656 
00657   /* sort the particles according to their positions */
00658   quicksort_particle(pG, cell1, 0, pG->nparticle-1);
00659 
00660   return;
00661 }
00662 
00663 /*----------------------------------------------------------------------------*/
00664 /*! \fn int compare_gr(Grid *pG, Vector cell1, Grain gr1, Grain gr2)
00665  *  \brief Compare the order of two particles according to their positions in 
00666  *  the grid
00667  *
00668  * Input: pG: grid; 
00669  * -      cell1: 1/dx1,1/dx2,1/dx3, or 0 if that dimension collapses.
00670  * -      gr1,gr2: pointers of the two particles to be compared.
00671  * Output: pointer of the particle that should be put in front of the other.
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);   /* x3 index of gr1 */
00678   k2 = (int)((gr2.x3 - pG->x3_0) * cell1.x3);   /* x3 index of gr2 */
00679   if (k1 < k2) return 1;
00680   if (k1 > k2) return 2;
00681 
00682   j1 = (int)((gr1.x2 - pG->x2_0) * cell1.x2);   /* x2 index of gr1 */
00683   j2 = (int)((gr2.x2 - pG->x2_0) * cell1.x2);   /* x2 index of gr2 */
00684   if (j1 < j2) return 1;
00685   if (j1 > j2) return 2;
00686 
00687   i1 = (int)((gr1.x1 - pG->x1_0) * cell1.x1);   /* x1 index of gr1 */
00688   i2 = (int)((gr2.x1 - pG->x1_0) * cell1.x1);   /* x1 index of gr2 */
00689   if (i1 < i2) return 1;
00690   if (i1 > i2) return 2;
00691 
00692   /* if they have equal indices, arbitrarily choose gr1 */
00693   return 1;
00694 }
00695 
00696 /*----------------------------------------------------------------------------*/
00697 /*! \fn void quicksort_particle(Grid *pG, Vector cell1, long start, long end)
00698  *  \brief Quick sort algorithm to shuffle the particles
00699  *
00700  * Input: pG, cell1: for compare_gr subroutine only. See above.
00701  *        head, rear: head and rear of the linked list.
00702  *          They do not contain data, or equal the pivot in the recursion.
00703  *        length: length of the linked list (does not contain head or rear).
00704  * Output: *head: linked list with shuffling finished.
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;     /* automatically sorted already */
00711 
00712   /* location of the pivot at half chain length */
00713   pivot = (long)((start+end+1)/2);
00714 
00715   /* move the pivot to the start */
00716   gr = pG->particle[pivot];
00717   pG->particle[pivot] = pG->particle[start];
00718   pG->particle[start] = gr;
00719 
00720   /* initial configuration */
00721   pivot = start;
00722   i = start + 1;
00723 
00724   /* move the particles that are "smaller" than the pivot before it */
00725   while (i <= end) {
00726     if (compare_gr(pG, cell1, pG->particle[pivot], pG->particle[i]) == 2)
00727     {/* the ith particle is smaller, move it before the pivot */
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   /* recursively call this routine to complete sorting */
00738   quicksort_particle(pG, cell1, start, pivot-1);
00739   quicksort_particle(pG, cell1, pivot+1, end);
00740 
00741   return;
00742 }
00743 
00744 #endif /*PARTICLES*/