Athena: gravity/selfg.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file selfg.c
00004  *  \brief Contains functions to control solution of Poisson's equation for
00005  *   self-gravity.
00006  *
00007  * CONTAINS PUBLIC FUNCTIONS:
00008  * - selfg_fc()   - 2nd order flux-corrections for self-gravity terms
00009  * - selfg_init() - sets pointer to appropriate self-gravity function
00010  *============================================================================*/
00011 
00012 #include <math.h>
00013 #include <float.h>
00014 #include "../defs.h"
00015 #include "../athena.h"
00016 #include "../globals.h"
00017 #include "prototypes.h"
00018 #include "../prototypes.h"
00019 
00020 #ifdef SELF_GRAVITY
00021 /*----------------------------------------------------------------------------*/
00022 /*! \fn void selfg_fc(DomainS *pD)
00023  *  \brief Adds flux-correction to make the integration algorithms for the
00024  *   source terms in the momentum and energy equations second-order.  
00025  *
00026  *   This
00027  *   requires subtracting 1/2 the source terms computed with the old potential,
00028  *   and adding 1/2 the source terms computed with the new potential.
00029  *
00030  *   The source terms for momentum are computed using the divergence of the
00031  *   gravitational stress tensor to conserve momentum exactly.
00032  *   - dM/dt = -Div(G);   G = (gg - 0.5g^2)/4\piG;   g=-Grad(Phi);
00033  *
00034  *   The source terms for the energy are added using the mass fluxes at cell
00035  *   faces, to improve conservation.
00036  *   - S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
00037  */
00038 
00039 void selfg_fc(DomainS *pD)
00040 {
00041   GridS *pG = (pD->Grid);
00042   int i, is = pG->is, ie = pG->ie;
00043   int j, js = pG->js, je = pG->je;
00044   int k, ks = pG->ks, ke = pG->ke;
00045   int dim=0;
00046   Real dtodx1 = pG->dt/pG->dx1;
00047   Real dtodx2 = pG->dt/pG->dx2;
00048   Real dtodx3 = pG->dt/pG->dx3;
00049   Real phic,phil,phir,phil_old,phir_old,dphic,dphil,dphir;
00050   Real gxl,gxr,gyl,gyr,gzl,gzr;
00051   Real flx_m1r,flx_m1l,flx_m2r,flx_m2l,flx_m3r,flx_m3l;
00052   
00053 /* Calculate the dimensions  */
00054   if(pG->Nx[0] > 1) dim++;
00055   if(pG->Nx[1] > 1) dim++;
00056   if(pG->Nx[2] > 1) dim++;
00057 
00058 
00059 /* The divergence of the gravitational stress tensor depends on the dimensions
00060  * of the problem.
00061  */
00062 
00063   switch(dim){
00064 /*------------------------- 1D problem ---------------------------------------*/
00065   case 1:
00066 
00067 /* Step 1 for 1D.  Add fluxes and source terms due to (d/dx1) terms  */
00068 
00069     for (i=is; i<=ie; i++){
00070       phic = pG->Phi[ks][js][i];
00071       phil = 0.5*(pG->Phi[ks][js][i-1] + pG->Phi[ks][js][i  ]);
00072       phir = 0.5*(pG->Phi[ks][js][i  ] + pG->Phi[ks][js][i+1]);
00073       phil_old = 0.5*(pG->Phi_old[ks][js][i-1] + pG->Phi_old[ks][js][i  ]);
00074       phir_old = 0.5*(pG->Phi_old[ks][js][i  ] + pG->Phi_old[ks][js][i+1]);
00075 
00076       dphic = phic - pG->Phi_old[ks][js][i];
00077       dphil = phil - phil_old;
00078       dphir = phir - phir_old;
00079 
00080 /* momentum fluxes in x1.  gx centered at L and R x1-faces */
00081       gxl = (pG->Phi[ks][js][i-1] - pG->Phi[ks][js][i  ])/(pG->dx1);
00082       gxr = (pG->Phi[ks][js][i  ] - pG->Phi[ks][js][i+1])/(pG->dx1);
00083 
00084       flx_m1l = 0.5*(gxl*gxl)/four_pi_G + grav_mean_rho*phil;
00085       flx_m1r = 0.5*(gxr*gxr)/four_pi_G + grav_mean_rho*phir;
00086 
00087 /*  subtract off momentum fluxes from old potential */
00088       gxl = (pG->Phi_old[ks][js][i-1] - pG->Phi_old[ks][js][i  ])/(pG->dx1);
00089       gxr = (pG->Phi_old[ks][js][i  ] - pG->Phi_old[ks][js][i+1])/(pG->dx1);
00090 
00091       flx_m1l -= 0.5*(gxl*gxl)/four_pi_G + grav_mean_rho*phil_old;
00092       flx_m1r -= 0.5*(gxr*gxr)/four_pi_G + grav_mean_rho*phir_old;
00093 
00094 /* Update momenta and energy with d/dx1 terms  */
00095       pG->U[ks][js][i].M1 -= 0.5*dtodx1*(flx_m1r-flx_m1l);
00096 #ifndef ISOTHERMAL
00097       pG->U[ks][js][i].E -=
00098          0.5*dtodx1*(pG->x1MassFlux[ks][js][i  ]*(dphic - dphil) +
00099                      pG->x1MassFlux[ks][js][i+1]*(dphir - dphic));
00100 #endif
00101     }
00102     break;
00103 
00104 /*------------------------- 2D problem ---------------------------------------*/
00105   case 2:
00106 
00107 /* Step 1 for 2D.  Add fluxes and source terms due to (d/dx1) terms  */
00108 
00109     for (j=js; j<=je; j++){
00110       for (i=is; i<=ie; i++){
00111         phic = pG->Phi[ks][j][i];
00112         phil = 0.5*(pG->Phi[ks][j][i-1] + pG->Phi[ks][j][i  ]);
00113         phir = 0.5*(pG->Phi[ks][j][i  ] + pG->Phi[ks][j][i+1]);
00114         phil_old = 0.5*(pG->Phi_old[ks][j][i-1] + pG->Phi_old[ks][j][i  ]);
00115         phir_old = 0.5*(pG->Phi_old[ks][j][i  ] + pG->Phi_old[ks][j][i+1]);
00116 
00117         dphic = phic - pG->Phi_old[ks][j][i];
00118         dphil = phil - phil_old;
00119         dphir = phir - phir_old;
00120 
00121 /*  momentum fluxes in x1.  gx and gy centered at L and R x1-faces */
00122         gxl = (pG->Phi[ks][j][i-1] - pG->Phi[ks][j][i  ])/(pG->dx1);
00123         gxr = (pG->Phi[ks][j][i  ] - pG->Phi[ks][j][i+1])/(pG->dx1);
00124 
00125         gyl = (pG->Phi[ks][j-1][i-1] - pG->Phi[ks][j+1][i-1]) +
00126               (pG->Phi[ks][j-1][i  ] - pG->Phi[ks][j+1][i  ]);
00127         gyl *= (0.25/pG->dx2);
00128 
00129         gyr = (pG->Phi[ks][j-1][i  ] - pG->Phi[ks][j+1][i  ]) +
00130               (pG->Phi[ks][j-1][i+1] - pG->Phi[ks][j+1][i+1]);
00131         gyr *= (0.25/pG->dx2);
00132 
00133         flx_m1l = 0.5*(gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil;
00134         flx_m1r = 0.5*(gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir;
00135 
00136         flx_m2l = gxl*gyl/four_pi_G;
00137         flx_m2r = gxr*gyr/four_pi_G;
00138 
00139 /*  subtract off momentum fluxes from old potential */
00140         gxl = (pG->Phi_old[ks][j][i-1] - pG->Phi_old[ks][j][i  ])/(pG->dx1);
00141         gxr = (pG->Phi_old[ks][j][i  ] - pG->Phi_old[ks][j][i+1])/(pG->dx1);
00142 
00143         gyl = (pG->Phi_old[ks][j-1][i-1] - pG->Phi_old[ks][j+1][i-1]) +
00144               (pG->Phi_old[ks][j-1][i  ] - pG->Phi_old[ks][j+1][i  ]);
00145         gyl *= (0.25/pG->dx2);
00146 
00147         gyr = (pG->Phi_old[ks][j-1][i  ] - pG->Phi_old[ks][j+1][i  ]) +
00148               (pG->Phi_old[ks][j-1][i+1] - pG->Phi_old[ks][j+1][i+1]);
00149         gyr *= (0.25/pG->dx2);
00150 
00151         flx_m1l -= 0.5*(gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil_old;
00152         flx_m1r -= 0.5*(gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir_old;
00153 
00154         flx_m2l -= gxl*gyl/four_pi_G;
00155         flx_m2r -= gxr*gyr/four_pi_G;
00156 
00157 /* Update momenta and energy with d/dx1 terms  */
00158         pG->U[ks][j][i].M1 -= 0.5*dtodx1*(flx_m1r - flx_m1l);
00159         pG->U[ks][j][i].M2 -= 0.5*dtodx1*(flx_m2r - flx_m2l);
00160 #ifndef ISOTHERMAL
00161         pG->U[ks][j][i].E -=
00162            0.5*dtodx1*(pG->x1MassFlux[ks][j][i  ]*(dphic - dphil) +
00163                        pG->x1MassFlux[ks][j][i+1]*(dphir - dphic));
00164 #endif
00165       }
00166     }
00167 
00168 /* Step 2 for 2D.  Add fluxes and source terms due to (d/dx2) terms  */
00169 
00170     for (j=js; j<=je; j++){
00171       for (i=is; i<=ie; i++){
00172         phic = pG->Phi[ks][j][i];
00173         phil = 0.5*(pG->Phi[ks][j-1][i] + pG->Phi[ks][j  ][i]);
00174         phir = 0.5*(pG->Phi[ks][j  ][i] + pG->Phi[ks][j+1][i]);
00175         phil_old = 0.5*(pG->Phi_old[ks][j-1][i] + pG->Phi_old[ks][j  ][i]);
00176         phir_old = 0.5*(pG->Phi_old[ks][j  ][i] + pG->Phi_old[ks][j+1][i]);
00177 
00178         dphic = phic - pG->Phi[ks][j][i];
00179         dphil = phil - phil_old;
00180         dphir = phir - phir_old;
00181 
00182 /*  momentum fluxes in x2.  gx and gy centered at L and R x2-faces */
00183         gxl = (pG->Phi[ks][j-1][i-1] - pG->Phi[ks][j-1][i+1]) +
00184               (pG->Phi[ks][j  ][i-1] - pG->Phi[ks][j  ][i+1]);
00185         gxl *= (0.25/pG->dx1);
00186 
00187         gxr = (pG->Phi[ks][j  ][i-1] - pG->Phi[ks][j  ][i+1]) +
00188               (pG->Phi[ks][j+1][i-1] - pG->Phi[ks][j+1][i+1]);
00189         gxr *= (0.25/pG->dx1);
00190 
00191         gyl = (pG->Phi[ks][j-1][i] - pG->Phi[ks][j  ][i])/(pG->dx2);
00192         gyr = (pG->Phi[ks][j  ][i] - pG->Phi[ks][j+1][i])/(pG->dx2);
00193 
00194         flx_m1l = gyl*gxl/four_pi_G;
00195         flx_m1r = gyr*gxr/four_pi_G;
00196 
00197         flx_m2l = 0.5*(gyl*gyl-gxl*gxl)/four_pi_G + grav_mean_rho*phil;
00198         flx_m2r = 0.5*(gyr*gyr-gxr*gxr)/four_pi_G + grav_mean_rho*phir;
00199 
00200 /*  subtract off momentum fluxes from old potential */
00201         gxl = (pG->Phi_old[ks][j-1][i-1] - pG->Phi_old[ks][j-1][i+1]) +
00202               (pG->Phi_old[ks][j  ][i-1] - pG->Phi_old[ks][j  ][i+1]);
00203         gxl *= (0.25/pG->dx1);
00204 
00205         gxr = (pG->Phi_old[ks][j  ][i-1] - pG->Phi_old[ks][j  ][i+1]) +
00206               (pG->Phi_old[ks][j+1][i-1] - pG->Phi_old[ks][j+1][i+1]);
00207         gxr *= (0.25/pG->dx1);
00208 
00209         gyl = (pG->Phi_old[ks][j-1][i] - pG->Phi_old[ks][j  ][i])/(pG->dx2);
00210         gyr = (pG->Phi_old[ks][j  ][i] - pG->Phi_old[ks][j+1][i])/(pG->dx2);
00211 
00212         flx_m1l -= gyl*gxl/four_pi_G;
00213         flx_m1r -= gyr*gxr/four_pi_G;
00214 
00215         flx_m2l -= 0.5*(gyl*gyl-gxl*gxl)/four_pi_G + grav_mean_rho*phil_old;
00216         flx_m2r -= 0.5*(gyr*gyr-gxr*gxr)/four_pi_G + grav_mean_rho*phir_old;
00217 
00218 /* Update momenta and energy with d/dx2 terms  */
00219         pG->U[ks][j][i].M1 -= 0.5*dtodx2*(flx_m1r - flx_m1l);
00220         pG->U[ks][j][i].M2 -= 0.5*dtodx2*(flx_m2r - flx_m2l);
00221 #ifndef ISOTHERMAL
00222         pG->U[ks][j][i].E -=
00223            0.5*dtodx2*(pG->x2MassFlux[ks][j  ][i]*(dphic - dphil) +
00224                        pG->x2MassFlux[ks][j+1][i]*(dphir - dphic));
00225 #endif
00226       }
00227     }
00228 
00229     break;
00230 
00231 /*------------------------- 3D problem ---------------------------------------*/
00232   case 3:
00233 
00234 /* Step 1 for 3D.  Add fluxes and source terms due to (d/dx1) terms  */
00235 
00236     for (k=ks; k<=ke; k++){
00237     for (j=js; j<=je; j++){
00238       for (i=is; i<=ie; i++){
00239         phic = pG->Phi[k][j][i];
00240         phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j][i  ]);
00241         phir = 0.5*(pG->Phi[k][j][i  ] + pG->Phi[k][j][i+1]);
00242         phil_old = 0.5*(pG->Phi_old[k][j][i-1] + pG->Phi_old[k][j][i  ]);
00243         phir_old = 0.5*(pG->Phi_old[k][j][i  ] + pG->Phi_old[k][j][i+1]);
00244 
00245         dphic = phic - pG->Phi_old[k][j][i];
00246         dphil = phil - phil_old;
00247         dphir = phir - phir_old;
00248 
00249 /*  momentum fluxes in x1. gx, gy and gz centered at L and R x1-faces */
00250         gxl = (pG->Phi[k][j][i-1] - pG->Phi[k][j][i  ])/(pG->dx1);
00251         gxr = (pG->Phi[k][j][i  ] - pG->Phi[k][j][i+1])/(pG->dx1);
00252 
00253         gyl = (pG->Phi[k][j-1][i-1] - pG->Phi[k][j+1][i-1]) +
00254               (pG->Phi[k][j-1][i  ] - pG->Phi[k][j+1][i  ]);
00255         gyl *= (0.25/pG->dx2);
00256 
00257         gyr = (pG->Phi[k][j-1][i  ] - pG->Phi[k][j+1][i  ]) +
00258               (pG->Phi[k][j-1][i+1] - pG->Phi[k][j+1][i+1]);
00259         gyr *= (0.25/pG->dx2);
00260 
00261         gzl = (pG->Phi[k-1][j][i-1] - pG->Phi[k+1][j][i-1]) +
00262               (pG->Phi[k-1][j][i  ] - pG->Phi[k+1][j][i  ]);
00263         gzl *= (0.25/pG->dx3);
00264 
00265         gzr = (pG->Phi[k-1][j][i  ] - pG->Phi[k+1][j][i  ]) +
00266               (pG->Phi[k-1][j][i+1] - pG->Phi[k+1][j][i+1]);
00267         gzr *= (0.25/pG->dx3);
00268 
00269         flx_m1l = 0.5*(gxl*gxl-gyl*gyl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
00270         flx_m1r = 0.5*(gxr*gxr-gyr*gyr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
00271 
00272         flx_m2l = gxl*gyl/four_pi_G;
00273         flx_m2r = gxr*gyr/four_pi_G;
00274 
00275         flx_m3l = gxl*gzl/four_pi_G;
00276         flx_m3r = gxr*gzr/four_pi_G;
00277 
00278 /*  subtract off momentum fluxes from old potential */
00279         gxl = (pG->Phi_old[k][j][i-1] - pG->Phi_old[k][j][i  ])/(pG->dx1);
00280         gxr = (pG->Phi_old[k][j][i  ] - pG->Phi_old[k][j][i+1])/(pG->dx1);
00281 
00282         gyl = (pG->Phi_old[k][j-1][i-1] - pG->Phi_old[k][j+1][i-1]) +
00283               (pG->Phi_old[k][j-1][i  ] - pG->Phi_old[k][j+1][i  ]);
00284         gyl *= (0.25/pG->dx2);
00285 
00286         gyr = (pG->Phi_old[k][j-1][i  ] - pG->Phi_old[k][j+1][i  ]) +
00287               (pG->Phi_old[k][j-1][i+1] - pG->Phi_old[k][j+1][i+1]);
00288         gyr *= (0.25/pG->dx2);
00289 
00290         gzl = (pG->Phi_old[k-1][j][i-1] - pG->Phi_old[k+1][j][i-1]) +
00291               (pG->Phi_old[k-1][j][i  ] - pG->Phi_old[k+1][j][i  ]);
00292         gzl *= (0.25/pG->dx3);
00293 
00294         gzr = (pG->Phi_old[k-1][j][i  ] - pG->Phi_old[k+1][j][i  ]) +
00295               (pG->Phi_old[k-1][j][i+1] - pG->Phi_old[k+1][j][i+1]);
00296         gzr *= (0.25/pG->dx3);
00297 
00298         flx_m1l -= 0.5*(gxl*gxl-gyl*gyl-gzl*gzl)/four_pi_G 
00299                    + grav_mean_rho*phil_old;
00300         flx_m1r -= 0.5*(gxr*gxr-gyr*gyr-gzr*gzr)/four_pi_G
00301                    + grav_mean_rho*phir_old;
00302 
00303         flx_m2l -= gxl*gyl/four_pi_G;
00304         flx_m2r -= gxr*gyr/four_pi_G;
00305 
00306         flx_m3l -= gxl*gzl/four_pi_G;
00307         flx_m3r -= gxr*gzr/four_pi_G;
00308 /* Update momenta and energy with d/dx1 terms  */
00309         pG->U[k][j][i].M1 -= 0.5*dtodx1*(flx_m1r - flx_m1l);
00310         pG->U[k][j][i].M2 -= 0.5*dtodx1*(flx_m2r - flx_m2l);
00311         pG->U[k][j][i].M3 -= 0.5*dtodx1*(flx_m3r - flx_m3l);
00312 #ifdef ADIABATIC
00313         pG->U[k][j][i].E -= 0.5*dtodx1*
00314           (pG->x1MassFlux[k][j][i  ]*(dphic - dphil) +
00315            pG->x1MassFlux[k][j][i+1]*(dphir - dphic));
00316 #endif /* ADIABATIC */
00317       }
00318     }}
00319 
00320 /* Step 2 for 3D.  Add fluxes and source terms due to (d/dx2) terms  */
00321 
00322     for (k=ks; k<=ke; k++){
00323     for (j=js; j<=je; j++){
00324       for (i=is; i<=ie; i++){
00325         phic = pG->Phi[k][j][i];
00326         phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j  ][i]);
00327         phir = 0.5*(pG->Phi[k][j  ][i] + pG->Phi[k][j+1][i]);
00328         phil_old = 0.5*(pG->Phi_old[k][j-1][i] + pG->Phi_old[k][j  ][i]);
00329         phir_old = 0.5*(pG->Phi_old[k][j  ][i] + pG->Phi_old[k][j+1][i]);
00330 
00331         dphic = phic - pG->Phi_old[k][j][i];
00332         dphil = phil - phil_old;
00333         dphir = phir - phir_old;
00334 
00335 /* gx, gy and gz centered at L and R x2-faces */
00336         gxl = (pG->Phi[k][j-1][i-1] - pG->Phi[k][j-1][i+1]) +
00337               (pG->Phi[k][j  ][i-1] - pG->Phi[k][j  ][i+1]);
00338         gxl *= (0.25/pG->dx1);
00339 
00340         gxr = (pG->Phi[k][j  ][i-1] - pG->Phi[k][j  ][i+1]) +
00341               (pG->Phi[k][j+1][i-1] - pG->Phi[k][j+1][i+1]);
00342         gxr *= (0.25/pG->dx1);
00343 
00344         gyl = (pG->Phi[k][j-1][i] - pG->Phi[k][j  ][i])/(pG->dx2);
00345         gyr = (pG->Phi[k][j  ][i] - pG->Phi[k][j+1][i])/(pG->dx2);
00346 
00347         gzl = (pG->Phi[k-1][j-1][i] - pG->Phi[k+1][j-1][i]) +
00348               (pG->Phi[k-1][j  ][i] - pG->Phi[k+1][j  ][i]);
00349         gzl *= (0.25/pG->dx3);
00350 
00351         gzr = (pG->Phi[k-1][j  ][i] - pG->Phi[k+1][j  ][i]) +
00352               (pG->Phi[k-1][j+1][i] - pG->Phi[k+1][j+1][i]);
00353         gzr *= (0.25/pG->dx3);
00354 
00355         flx_m1l = gyl*gxl/four_pi_G;
00356         flx_m1r = gyr*gxr/four_pi_G;
00357 
00358         flx_m2l = 0.5*(gyl*gyl-gxl*gxl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
00359         flx_m2r = 0.5*(gyr*gyr-gxr*gxr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
00360 
00361         flx_m3l = gyl*gzl/four_pi_G;
00362         flx_m3r = gyr*gzr/four_pi_G;
00363 
00364 /*  subtract off momentum fluxes from old potential */
00365         gxl = (pG->Phi_old[k][j-1][i-1] - pG->Phi_old[k][j-1][i+1]) +
00366               (pG->Phi_old[k][j  ][i-1] - pG->Phi_old[k][j  ][i+1]);
00367         gxl *= (0.25/pG->dx1);
00368 
00369         gxr = (pG->Phi_old[k][j  ][i-1] - pG->Phi_old[k][j  ][i+1]) +
00370               (pG->Phi_old[k][j+1][i-1] - pG->Phi_old[k][j+1][i+1]);
00371         gxr *= (0.25/pG->dx1);
00372 
00373         gyl = (pG->Phi_old[k][j-1][i] - pG->Phi_old[k][j  ][i])/(pG->dx2);
00374         gyr = (pG->Phi_old[k][j  ][i] - pG->Phi_old[k][j+1][i])/(pG->dx2);
00375 
00376         gzl = (pG->Phi_old[k-1][j-1][i] - pG->Phi_old[k+1][j-1][i]) +
00377               (pG->Phi_old[k-1][j  ][i] - pG->Phi_old[k+1][j  ][i]);
00378         gzl *= (0.25/pG->dx3);
00379 
00380         gzr = (pG->Phi_old[k-1][j  ][i] - pG->Phi_old[k+1][j  ][i]) +
00381               (pG->Phi_old[k-1][j+1][i] - pG->Phi_old[k+1][j+1][i]);
00382         gzr *= (0.25/pG->dx3);
00383 
00384         flx_m1l -= gyl*gxl/four_pi_G;
00385         flx_m1r -= gyr*gxr/four_pi_G;
00386 
00387         flx_m2l -= 0.5*(gyl*gyl-gxl*gxl-gzl*gzl)/four_pi_G
00388                  + grav_mean_rho*phil_old;
00389         flx_m2r -= 0.5*(gyr*gyr-gxr*gxr-gzr*gzr)/four_pi_G 
00390                  + grav_mean_rho*phir_old;
00391 
00392         flx_m3l -= gyl*gzl/four_pi_G;
00393         flx_m3r -= gyr*gzr/four_pi_G;
00394 
00395 /* Update momenta and energy with d/dx2 terms  */
00396         pG->U[k][j][i].M1 -= 0.5*dtodx2*(flx_m1r - flx_m1l);
00397         pG->U[k][j][i].M2 -= 0.5*dtodx2*(flx_m2r - flx_m2l);
00398         pG->U[k][j][i].M3 -= 0.5*dtodx2*(flx_m3r - flx_m3l);
00399 #ifdef ADIABATIC
00400         pG->U[k][j][i].E -= 0.5*dtodx2*
00401           (pG->x2MassFlux[k][j  ][i]*(dphic - dphil) +
00402            pG->x2MassFlux[k][j+1][i]*(dphir - dphic));
00403 #endif /* ADIABATIC */
00404       }
00405     }}
00406 
00407 /* Step 3 for 3D.  Add fluxes and source terms due to (d/dx3) terms  */
00408 
00409     for (k=ks; k<=ke; k++){
00410     for (j=js; j<=je; j++){
00411       for (i=is; i<=ie; i++){
00412         phic = pG->Phi[k][j][i];
00413         phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k  ][j][i]);
00414         phir = 0.5*(pG->Phi[k  ][j][i] + pG->Phi[k+1][j][i]);
00415         phil_old = 0.5*(pG->Phi_old[k-1][j][i] + pG->Phi_old[k  ][j][i]);
00416         phir_old = 0.5*(pG->Phi_old[k  ][j][i] + pG->Phi_old[k+1][j][i]);
00417 
00418         dphic = phic - pG->Phi_old[k][j][i];
00419         dphil = phil - phil_old;
00420         dphir = phir - phir_old;
00421 
00422 /*  momentum fluxes in x3. gx, gy and gz centered at L and R x3-faces */
00423         gxl = (pG->Phi[k-1][j][i-1] - pG->Phi[k-1][j][i+1]) +
00424               (pG->Phi[k  ][j][i-1] - pG->Phi[k  ][j][i+1]);
00425         gxl *= (0.25/pG->dx1);
00426 
00427         gxr = (pG->Phi[k  ][j][i-1] - pG->Phi[k  ][j][i+1]) +
00428               (pG->Phi[k+1][j][i-1] - pG->Phi[k+1][j][i+1]);
00429         gxr *= (0.25/pG->dx1);
00430 
00431         gyl = (pG->Phi[k-1][j-1][i] - pG->Phi[k-1][j+1][i]) +
00432               (pG->Phi[k  ][j-1][i] - pG->Phi[k  ][j+1][i]);
00433         gyl *= (0.25/pG->dx2);
00434 
00435         gyr = (pG->Phi[k  ][j-1][i] - pG->Phi[k  ][j+1][i]) +
00436               (pG->Phi[k+1][j-1][i] - pG->Phi[k+1][j+1][i]);
00437         gyr *= (0.25/pG->dx2);
00438 
00439         gzl = (pG->Phi[k-1][j][i] - pG->Phi[k  ][j][i])/(pG->dx3);
00440         gzr = (pG->Phi[k  ][j][i] - pG->Phi[k+1][j][i])/(pG->dx3);
00441 
00442         flx_m1l = gzl*gxl/four_pi_G;
00443         flx_m1r = gzr*gxr/four_pi_G;
00444 
00445         flx_m2l = gzl*gyl/four_pi_G;
00446         flx_m2r = gzr*gyr/four_pi_G;
00447 
00448         flx_m3l = 0.5*(gzl*gzl-gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil;
00449         flx_m3r = 0.5*(gzr*gzr-gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir;
00450 
00451 /*  subtract off momentum fluxes from old potential */
00452         gxl = (pG->Phi_old[k-1][j][i-1] - pG->Phi_old[k-1][j][i+1]) +
00453               (pG->Phi_old[k  ][j][i-1] - pG->Phi_old[k  ][j][i+1]);
00454         gxl *= (0.25/pG->dx1);
00455 
00456         gxr = (pG->Phi_old[k  ][j][i-1] - pG->Phi_old[k  ][j][i+1]) +
00457               (pG->Phi_old[k+1][j][i-1] - pG->Phi_old[k+1][j][i+1]);
00458         gxr *= (0.25/pG->dx1);
00459 
00460         gyl = (pG->Phi_old[k-1][j-1][i] - pG->Phi_old[k-1][j+1][i]) +
00461               (pG->Phi_old[k  ][j-1][i] - pG->Phi_old[k  ][j+1][i]);
00462         gyl *= (0.25/pG->dx2);
00463 
00464         gyr = (pG->Phi_old[k  ][j-1][i] - pG->Phi_old[k  ][j+1][i]) +
00465               (pG->Phi_old[k+1][j-1][i] - pG->Phi_old[k+1][j+1][i]);
00466         gyr *= (0.25/pG->dx2);
00467 
00468         gzl = (pG->Phi_old[k-1][j][i] - pG->Phi_old[k  ][j][i])/(pG->dx3);
00469         gzr = (pG->Phi_old[k  ][j][i] - pG->Phi_old[k+1][j][i])/(pG->dx3);
00470 
00471         flx_m1l -= gzl*gxl/four_pi_G;
00472         flx_m1r -= gzr*gxr/four_pi_G;
00473 
00474         flx_m2l -= gzl*gyl/four_pi_G;
00475         flx_m2r -= gzr*gyr/four_pi_G;
00476 
00477         flx_m3l -= 0.5*(gzl*gzl-gxl*gxl-gyl*gyl)/four_pi_G 
00478                  + grav_mean_rho*phil_old;
00479         flx_m3r -= 0.5*(gzr*gzr-gxr*gxr-gyr*gyr)/four_pi_G
00480                  + grav_mean_rho*phir_old;
00481 
00482 /* Update momenta and energy with d/dx3 terms  */
00483         pG->U[k][j][i].M1 -= 0.5*dtodx3*(flx_m1r - flx_m1l);
00484         pG->U[k][j][i].M2 -= 0.5*dtodx3*(flx_m2r - flx_m2l);
00485         pG->U[k][j][i].M3 -= 0.5*dtodx3*(flx_m3r - flx_m3l);
00486 #ifdef ADIABATIC
00487         pG->U[k][j][i].E -= 0.5*dtodx3*
00488           (pG->x3MassFlux[k  ][j][i]*(dphic - dphil) +
00489            pG->x3MassFlux[k+1][j][i]*(dphir - dphic));
00490 #endif /* ADIABATIC */
00491       }
00492     }}
00493 
00494     break;
00495 
00496   } /* end of switch statement */
00497 
00498   return;
00499 }
00500 
00501 /*----------------------------------------------------------------------------*/
00502 /*! \fn VDFun_t selfg_init(MeshS *pM)
00503  *  \brief Initialize pointer to appropriate self-gravity f'n, 
00504  *   allocates memory for dPhi array used for flux correction.
00505  */
00506 
00507 VDFun_t selfg_init(MeshS *pM)
00508 {
00509   int dim = 0;
00510 
00511 /* Calculate the dimensions  */
00512   if(pM->Nx[0] > 1) dim++;
00513   if(pM->Nx[1] > 1) dim++;
00514   if(pM->Nx[2] > 1) dim++;
00515 
00516 /* test that user set values for constants */
00517   if (grav_mean_rho < 0.0)
00518     ath_error("[selfg_init] grav_mean_rho must be set >0 in prob generator\n");
00519   if (four_pi_G < 0.0)
00520     ath_error("[selfg_init] four_pi_G must be set >0 in prob generator\n");
00521 
00522 /* Return function pointer based on dimensions and algorithm */
00523 
00524   switch(dim){
00525 #ifdef SELF_GRAVITY_USING_MULTIGRID
00526   case 1:
00527     return selfg_multig_1d;
00528   case 2:
00529     return selfg_multig_2d;
00530   case 3:
00531     selfg_multig_3d_init(pM);
00532     return selfg_multig_3d;
00533 #endif
00534 
00535 /* for gravity using FFTs, also initialize plans and data for FFTW */
00536 #ifdef SELF_GRAVITY_USING_FFT
00537   case 1:
00538     return selfg_fft_1d;
00539   case 2:
00540     selfg_fft_2d_init(pM);
00541     return selfg_fft_2d;
00542   case 3:
00543     selfg_fft_3d_init(pM);
00544     return selfg_fft_3d;
00545 #endif
00546 /* for gravity using FFTs with open BC, also initialize plans and data for FFTW */
00547 #ifdef SELF_GRAVITY_USING_FFT_OBC
00548   case 1:
00549     ath_error("[selfg_init] FFT with open BC not defined for 1D \n");
00550   case 2:
00551     ath_error("[selfg_init] FFT with open BC not defined for 2D \n");
00552   case 3:
00553     selfg_fft_obc_3d_init(pM);
00554     return selfg_fft_obc_3d;
00555 #endif
00556 
00557 
00558 
00559   }
00560 
00561   return NULL;
00562 }
00563 #endif /* SELF_GRAVITY */