Go to the documentation of this file.00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
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
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
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
00054 if(pG->Nx[0] > 1) dim++;
00055 if(pG->Nx[1] > 1) dim++;
00056 if(pG->Nx[2] > 1) dim++;
00057
00058
00059
00060
00061
00062
00063 switch(dim){
00064
00065 case 1:
00066
00067
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
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
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
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
00105 case 2:
00106
00107
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
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
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
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
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
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
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
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
00232 case 3:
00233
00234
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
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
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
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
00317 }
00318 }}
00319
00320
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
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
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
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
00404 }
00405 }}
00406
00407
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
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
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
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
00491 }
00492 }}
00493
00494 break;
00495
00496 }
00497
00498 return;
00499 }
00500
00501
00502
00503
00504
00505
00506
00507 VDFun_t selfg_init(MeshS *pM)
00508 {
00509 int dim = 0;
00510
00511
00512 if(pM->Nx[0] > 1) dim++;
00513 if(pM->Nx[1] > 1) dim++;
00514 if(pM->Nx[2] > 1) dim++;
00515
00516
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
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
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
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