00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include <math.h>
00038 #include <stdio.h>
00039 #include <stdlib.h>
00040 #include <string.h>
00041 #include "../defs.h"
00042 #include "../athena.h"
00043 #include "../globals.h"
00044 #include "prototypes.h"
00045 #include "../prototypes.h"
00046 #ifdef PARTICLES
00047 #include "../particles/particle.h"
00048 #endif
00049
00050 #if defined(CTU_INTEGRATOR)
00051 #ifdef SPECIAL_RELATIVITY
00052 #error : The CTU integrator cannot be used for special relativity.
00053 #endif
00054
00055
00056 static Cons1DS ***Ul_x1Face=NULL, ***Ur_x1Face=NULL;
00057 static Cons1DS ***Ul_x2Face=NULL, ***Ur_x2Face=NULL;
00058 static Cons1DS ***Ul_x3Face=NULL, ***Ur_x3Face=NULL;
00059 Cons1DS ***x1Flux=NULL, ***x2Flux=NULL, ***x3Flux=NULL;
00060
00061
00062 #ifdef MHD
00063 static Real ***B1_x1Face=NULL, ***B2_x2Face=NULL, ***B3_x3Face=NULL;
00064 Real ***emf1=NULL, ***emf2=NULL, ***emf3=NULL;
00065 static Real ***emf1_cc=NULL, ***emf2_cc=NULL, ***emf3_cc=NULL;
00066 #endif
00067
00068
00069 static Real *Bxc=NULL, *Bxi=NULL;
00070 static Prim1DS *W=NULL, *Wl=NULL, *Wr=NULL;
00071 static Cons1DS *U1d=NULL;
00072
00073
00074 static Real ***dhalf = NULL, ***phalf=NULL;
00075
00076
00077 extern Real etah;
00078 #ifdef H_CORRECTION
00079 static Real ***eta1=NULL, ***eta2=NULL, ***eta3=NULL;
00080 #endif
00081
00082
00083 #ifdef SHEARING_BOX
00084 static Real **remapEyiib=NULL, **remapEyoib=NULL;
00085 #endif
00086
00087
00088 #ifdef CYLINDRICAL
00089 static Real ***geom_src=NULL;
00090 #endif
00091
00092
00093
00094
00095
00096
00097
00098
00099 #ifdef MHD
00100 static void integrate_emf1_corner(const GridS *pG);
00101 static void integrate_emf2_corner(const GridS *pG);
00102 static void integrate_emf3_corner(const GridS *pG);
00103 #endif
00104
00105
00106
00107
00108
00109
00110 void integrate_3d_ctu(DomainS *pD)
00111 {
00112 GridS *pG=(pD->Grid);
00113 Real dtodx1=pG->dt/pG->dx1, dtodx2=pG->dt/pG->dx2, dtodx3=pG->dt/pG->dx3;
00114 Real hdt = 0.5*pG->dt, dx2=pG->dx2;
00115 Real q1 = 0.5*dtodx1, q2 = 0.5*dtodx2, q3 = 0.5*dtodx3;
00116 int i,il,iu, is = pG->is, ie = pG->ie;
00117 int j,jl,ju, js = pG->js, je = pG->je;
00118 int k,kl,ku, ks = pG->ks, ke = pG->ke;
00119 Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic,M1h,M2h,M3h,Bx=0.0;
00120 #ifndef BAROTROPIC
00121 Real coolfl,coolfr,coolf,Eh=0.0;
00122 #endif
00123 #ifdef MHD
00124 Real MHD_src_By,MHD_src_Bz,mdb1,mdb2,mdb3;
00125 Real db1,db2,db3,l1,l2,l3,B1,B2,B3,V1,V2,V3;
00126 Real B1ch,B2ch,B3ch;
00127 #endif
00128 #if defined(MHD) || defined(SELF_GRAVITY)
00129 Real dx1i=1.0/pG->dx1, dx2i=1.0/pG->dx2, dx3i=1.0/pG->dx3;
00130 #endif
00131 #ifdef H_CORRECTION
00132 Real cfr,cfl,lambdar,lambdal;
00133 #endif
00134 #if (NSCALARS > 0)
00135 int n;
00136 #endif
00137 #ifdef SELF_GRAVITY
00138 Real gxl,gxr,gyl,gyr,gzl,gzr,flx_m1l,flx_m1r,flx_m2l,flx_m2r,flx_m3l,flx_m3r;
00139 #endif
00140 #ifdef FEEDBACK
00141 Real dt1 = 1.0/pG->dt;
00142 #endif
00143 #ifdef SHEARING_BOX
00144 int my_iproc,my_jproc,my_kproc;
00145 Real M1n, dM2n;
00146 Real M1e, dM2e;
00147 Real flx1_dM2, frx1_dM2, flx2_dM2, frx2_dM2, flx3_dM2, frx3_dM2;
00148 Real fact, qom, om_dt = Omega_0*pG->dt;
00149 #endif
00150 #ifdef STATIC_MESH_REFINEMENT
00151 int ncg,npg,dim;
00152 int ii,ics,ice,jj,jcs,jce,kk,kcs,kce,ips,ipe,jps,jpe,kps,kpe;
00153 #endif
00154
00155
00156 #ifdef CYLINDRICAL
00157 #ifndef ISOTHERMAL
00158 Real Pavgh;
00159 #endif
00160 Real g,gl,gr,rinv;
00161 Real geom_src_d, geom_src_Vx, geom_src_Vy, geom_src_P, geom_src_By, geom_src_Bz;
00162 const Real *r=pG->r, *ri=pG->ri;
00163 #ifdef FARGO
00164 Real Om, qshear, Mrn, Mpn, Mre, Mpe, Mrav, Mpav;
00165 #endif
00166 #endif
00167 Real lsf=1.0, rsf=1.0;
00168
00169
00170 #ifdef PARTICLES
00171 Real d1;
00172 il = is - 3;
00173 iu = ie + 3;
00174 jl = js - 3;
00175 ju = je + 3;
00176 kl = ks - 3;
00177 ku = ke + 3;
00178 #else
00179 il = is - 2;
00180 iu = ie + 2;
00181 jl = js - 2;
00182 ju = je + 2;
00183 kl = ks - 2;
00184 ku = ke + 2;
00185 #endif
00186
00187
00188 etah = 0.0;
00189
00190
00191 #ifdef FEEDBACK
00192 feedback_predictor(pG);
00193 #endif
00194
00195
00196
00197
00198
00199
00200
00201
00202 for (k=kl; k<=ku; k++) {
00203 for (j=jl; j<=ju; j++) {
00204 for (i=is-nghost; i<=ie+nghost; i++) {
00205 U1d[i].d = pG->U[k][j][i].d;
00206 U1d[i].Mx = pG->U[k][j][i].M1;
00207 U1d[i].My = pG->U[k][j][i].M2;
00208 U1d[i].Mz = pG->U[k][j][i].M3;
00209 #ifndef BAROTROPIC
00210 U1d[i].E = pG->U[k][j][i].E;
00211 #endif
00212 #ifdef MHD
00213 U1d[i].By = pG->U[k][j][i].B2c;
00214 U1d[i].Bz = pG->U[k][j][i].B3c;
00215 Bxc[i] = pG->U[k][j][i].B1c;
00216 Bxi[i] = pG->B1i[k][j][i];
00217 B1_x1Face[k][j][i] = pG->B1i[k][j][i];
00218 #endif
00219 #if (NSCALARS > 0)
00220 for (n=0; n<NSCALARS; n++) U1d[i].s[n] = pG->U[k][j][i].s[n];
00221 #endif
00222 }
00223
00224
00225
00226
00227
00228 for (i=is-nghost; i<=ie+nghost; i++) {
00229 W[i] = Cons1D_to_Prim1D(&U1d[i],&Bxc[i]);
00230
00231
00232
00233
00234 #ifdef CYLINDRICAL
00235 geom_src[k][j][i] = W[i].d*SQR(W[i].Vy);
00236 #ifdef MHD
00237 geom_src[k][j][i] += 0.5*(SQR(Bxc[i]) - SQR(W[i].By) + SQR(W[i].Bz));
00238 #endif
00239 #ifdef ISOTHERMAL
00240 geom_src[k][j][i] += Iso_csound2*W[i].d;
00241 #else
00242 geom_src[k][j][i] += W[i].P;
00243 #endif
00244 geom_src[k][j][i] /= x1vc(pG,i);
00245 #endif
00246 }
00247
00248 lr_states(pG,W,Bxc,pG->dt,pG->dx1,il+1,iu-1,Wl,Wr,1);
00249
00250 #ifdef MHD
00251 for (i=il+1; i<=iu; i++) {
00252
00253 #ifdef CYLINDRICAL
00254 rsf = ri[i]/r[i-1]; lsf = ri[i-1]/r[i-1];
00255 dx2i = 1.0/(r[i-1]*pG->dx2);
00256 #endif
00257 db1 = (rsf*pG->B1i[k ][j ][i ] - lsf*pG->B1i[k][j][i-1])*dx1i;
00258 db2 = ( pG->B2i[k ][j+1][i-1] - pG->B2i[k][j][i-1])*dx2i;
00259 db3 = ( pG->B3i[k+1][j ][i-1] - pG->B3i[k][j][i-1])*dx3i;
00260
00261 if(db1 >= 0.0){
00262 l3 = db1 < -db3 ? db1 : -db3;
00263 l3 = l3 > 0.0 ? l3 : 0.0;
00264
00265 l2 = db1 < -db2 ? db1 : -db2;
00266 l2 = l2 > 0.0 ? l2 : 0.0;
00267 }
00268 else{
00269 l3 = db1 > -db3 ? db1 : -db3;
00270 l3 = l3 < 0.0 ? l3 : 0.0;
00271
00272 l2 = db1 > -db2 ? db1 : -db2;
00273 l2 = l2 < 0.0 ? l2 : 0.0;
00274 }
00275
00276 MHD_src_By = (pG->U[k][j][i-1].M2/pG->U[k][j][i-1].d)*l2;
00277 MHD_src_Bz = (pG->U[k][j][i-1].M3/pG->U[k][j][i-1].d)*l3;
00278
00279 Wl[i].By += hdt*MHD_src_By;
00280 Wl[i].Bz += hdt*MHD_src_Bz;
00281
00282
00283 #ifdef CYLINDRICAL
00284 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
00285 dx2i = 1.0/(r[i]*pG->dx2);
00286 #endif
00287 db1 = (rsf*pG->B1i[k ][j ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
00288 db2 = ( pG->B2i[k ][j+1][i ] - pG->B2i[k][j][i])*dx2i;
00289 db3 = ( pG->B3i[k+1][j ][i ] - pG->B3i[k][j][i])*dx3i;
00290
00291 if(db1 >= 0.0){
00292 l3 = db1 < -db3 ? db1 : -db3;
00293 l3 = l3 > 0.0 ? l3 : 0.0;
00294
00295 l2 = db1 < -db2 ? db1 : -db2;
00296 l2 = l2 > 0.0 ? l2 : 0.0;
00297 }
00298 else{
00299 l3 = db1 > -db3 ? db1 : -db3;
00300 l3 = l3 < 0.0 ? l3 : 0.0;
00301
00302 l2 = db1 > -db2 ? db1 : -db2;
00303 l2 = l2 < 0.0 ? l2 : 0.0;
00304 }
00305
00306 MHD_src_By = (pG->U[k][j][i].M2/pG->U[k][j][i].d)*l2;
00307 MHD_src_Bz = (pG->U[k][j][i].M3/pG->U[k][j][i].d)*l3;
00308
00309 Wr[i].By += hdt*MHD_src_By;
00310 Wr[i].Bz += hdt*MHD_src_Bz;
00311 }
00312 #endif
00313
00314
00315
00316
00317
00318 if (StaticGravPot != NULL){
00319 for (i=il+1; i<=iu; i++) {
00320 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00321 #ifdef CYLINDRICAL
00322 gl = (*x1GravAcc)(x1vc(pG,i-1),x2,x3);
00323 gr = (*x1GravAcc)(x1vc(pG,i),x2,x3);
00324 #ifdef FARGO
00325
00326 gl = gl - x1vc(pG,i-1)*SQR((*OrbitalProfile)(x1vc(pG,i-1)));
00327 gr = gr - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
00328 #endif
00329
00330
00331 Wl[i].Vx -= hdt*gl;
00332 Wr[i].Vx -= hdt*gr;
00333 #else
00334 phicr = (*StaticGravPot)( x1 ,x2,x3);
00335 phicl = (*StaticGravPot)((x1- pG->dx1),x2,x3);
00336 phifc = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00337
00338 Wl[i].Vx -= dtodx1*(phifc - phicl);
00339 Wr[i].Vx -= dtodx1*(phicr - phifc);
00340 #endif
00341 }
00342 }
00343
00344
00345
00346
00347
00348 #ifdef SELF_GRAVITY
00349 for (i=il+1; i<=iu; i++) {
00350 Wl[i].Vx -= q1*(pG->Phi[k][j][i] - pG->Phi[k][j][i-1]);
00351 Wr[i].Vx -= q1*(pG->Phi[k][j][i] - pG->Phi[k][j][i-1]);
00352 }
00353 #endif
00354
00355
00356
00357
00358
00359 #ifndef BAROTROPIC
00360 if (CoolingFunc != NULL){
00361 for (i=il+1; i<=iu; i++) {
00362 coolfl = (*CoolingFunc)(Wl[i].d,Wl[i].P,(0.5*pG->dt));
00363 coolfr = (*CoolingFunc)(Wr[i].d,Wr[i].P,(0.5*pG->dt));
00364
00365 Wl[i].P -= 0.5*pG->dt*Gamma_1*coolfl;
00366 Wr[i].P -= 0.5*pG->dt*Gamma_1*coolfr;
00367 }
00368 }
00369 #endif
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379 #ifdef SHEARING_BOX
00380 if (ShearingBoxPot != NULL){
00381 for (i=il+1; i<=iu; i++) {
00382 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00383 phicr = (*ShearingBoxPot)( x1 ,x2,x3);
00384 phicl = (*ShearingBoxPot)((x1- pG->dx1),x2,x3);
00385 phifc = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
00386
00387 Wl[i].Vx -= dtodx1*(phifc - phicl);
00388 Wr[i].Vx -= dtodx1*(phicr - phifc);
00389 }
00390 }
00391
00392 for (i=il+1; i<=iu; i++) {
00393 Wl[i].Vx += pG->dt*Omega_0*W[i-1].Vy;
00394 Wr[i].Vx += pG->dt*Omega_0*W[i].Vy;
00395 #ifdef FARGO
00396 Wl[i].Vy += hdt*(qshear-2.)*Omega_0*W[i-1].Vx;
00397 Wr[i].Vy += hdt*(qshear-2.)*Omega_0*W[i].Vx;
00398 #else
00399 Wl[i].Vy -= pG->dt*Omega_0*W[i-1].Vx;
00400 Wr[i].Vy -= pG->dt*Omega_0*W[i].Vx;
00401 #endif
00402 }
00403 #endif
00404
00405 #if defined(CYLINDRICAL) && defined(FARGO)
00406 for (i=il+1; i<=iu; i++) {
00407 Om = (*OrbitalProfile)(x1vc(pG,i-1));
00408 qshear = (*ShearProfile)(x1vc(pG,i-1));
00409 Wl[i].Vx += (pG->dt)*Om*W[i-1].Vy;
00410 Wl[i].Vy += hdt*(qshear - 2.0)*Om*W[i-1].Vx;
00411
00412 Om = (*OrbitalProfile)(x1vc(pG,i));
00413 qshear = (*ShearProfile)(x1vc(pG,i));
00414 Wr[i].Vx += (pG->dt)*Om*W[i].Vy;
00415 Wr[i].Vy += hdt*(qshear - 2.0)*Om*W[i].Vx;
00416 }
00417 #endif
00418
00419
00420
00421
00422
00423 #ifdef FEEDBACK
00424 for (i=il+1; i<=iu; i++) {
00425 d1 = 1.0/W[i-1].d;
00426 Wl[i].Vx -= pG->Coup[k][j][i-1].fb1*d1;
00427 Wl[i].Vy -= pG->Coup[k][j][i-1].fb2*d1;
00428 Wl[i].Vz -= pG->Coup[k][j][i-1].fb3*d1;
00429
00430 d1 = 1.0/W[i].d;
00431 Wr[i].Vx -= pG->Coup[k][j][i].fb1*d1;
00432 Wr[i].Vy -= pG->Coup[k][j][i].fb2*d1;
00433 Wr[i].Vz -= pG->Coup[k][j][i].fb3*d1;
00434
00435 #ifndef BAROTROPIC
00436 Wl[i].P += pG->Coup[k][j][i-1].Eloss*Gamma_1;
00437 Wr[i].P += pG->Coup[k][j][i].Eloss*Gamma_1;
00438 #endif
00439
00440 }
00441 #endif
00442
00443
00444
00445
00446
00447 #ifdef CYLINDRICAL
00448 for (i=is-1; i<=ie+2; i++) {
00449
00450
00451 rinv = 1.0/x1vc(pG,i-1);
00452 geom_src_d = -W[i-1].d*W[i-1].Vx*rinv;
00453 geom_src_Vx = SQR(W[i-1].Vy);
00454 geom_src_Vy = -W[i-1].Vx*W[i-1].Vy;
00455 #ifdef MHD
00456 geom_src_Vx -= SQR(W[i-1].By)/W[i-1].d;
00457 geom_src_Vy += Bxc[i-1]*W[i-1].By/W[i-1].d;
00458 geom_src_By = -W[i-1].Vy*Bxc[i-1]*rinv;
00459 geom_src_Bz = -W[i-1].Vx*W[i-1].Bz*rinv;
00460 #endif
00461 geom_src_Vx *= rinv;
00462 geom_src_Vy *= rinv;
00463 #ifndef ISOTHERMAL
00464 geom_src_P = -Gamma*W[i-1].P*W[i-1].Vx*rinv;
00465 #endif
00466
00467
00468 Wl[i].d += hdt*geom_src_d;
00469 Wl[i].Vx += hdt*geom_src_Vx;
00470 Wl[i].Vy += hdt*geom_src_Vy;
00471 #ifdef MHD
00472 Wl[i].By += hdt*geom_src_By;
00473 Wl[i].Bz += hdt*geom_src_Bz;
00474 #endif
00475 #ifndef ISOTHERMAL
00476 Wl[i].P += hdt*geom_src_P;
00477 #endif
00478
00479
00480
00481 rinv = 1.0/x1vc(pG,i);
00482 geom_src_d = -W[i].d*W[i].Vx*rinv;
00483 geom_src_Vx = SQR(W[i].Vy);
00484 geom_src_Vy = -W[i].Vx*W[i].Vy;
00485 #ifdef MHD
00486 geom_src_Vx -= SQR(W[i].By)/W[i].d;
00487 geom_src_Vy += Bxc[i]*W[i].By/W[i].d;
00488 geom_src_By = -W[i].Vy*Bxc[i]*rinv;
00489 geom_src_Bz = -W[i].Vx*W[i].Bz*rinv;
00490 #endif
00491 geom_src_Vx *= rinv;
00492 geom_src_Vy *= rinv;
00493 #ifndef ISOTHERMAL
00494 geom_src_P = -Gamma*W[i].P*W[i].Vx*rinv;
00495 #endif
00496
00497
00498 Wr[i].d += hdt*geom_src_d;
00499 Wr[i].Vx += hdt*geom_src_Vx;
00500 Wr[i].Vy += hdt*geom_src_Vy;
00501 #ifdef MHD
00502 Wr[i].By += hdt*geom_src_By;
00503 Wr[i].Bz += hdt*geom_src_Bz;
00504 #endif
00505 #ifndef ISOTHERMAL
00506 Wr[i].P += hdt*geom_src_P;
00507 #endif
00508 }
00509 #endif
00510
00511
00512
00513
00514
00515 for (i=il+1; i<=iu; i++) {
00516 Ul_x1Face[k][j][i] = Prim1D_to_Cons1D(&Wl[i],&Bxi[i]);
00517 Ur_x1Face[k][j][i] = Prim1D_to_Cons1D(&Wr[i],&Bxi[i]);
00518
00519 #ifdef MHD
00520 Bx = B1_x1Face[k][j][i];
00521 #endif
00522 fluxes(Ul_x1Face[k][j][i],Ur_x1Face[k][j][i],Wl[i],Wr[i],Bx,
00523 &x1Flux[k][j][i]);
00524 }
00525 }
00526 }
00527
00528
00529
00530
00531
00532
00533
00534
00535 for (k=kl; k<=ku; k++) {
00536 for (i=il; i<=iu; i++) {
00537 #ifdef CYLINDRICAL
00538 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
00539 dx2 = r[i]*pG->dx2;
00540 dtodx2 = pG->dt/dx2;
00541 #endif
00542 for (j=js-nghost; j<=je+nghost; j++) {
00543 U1d[j].d = pG->U[k][j][i].d;
00544 U1d[j].Mx = pG->U[k][j][i].M2;
00545 U1d[j].My = pG->U[k][j][i].M3;
00546 U1d[j].Mz = pG->U[k][j][i].M1;
00547 #ifndef BAROTROPIC
00548 U1d[j].E = pG->U[k][j][i].E;
00549 #endif
00550 #ifdef MHD
00551 U1d[j].By = pG->U[k][j][i].B3c;
00552 U1d[j].Bz = pG->U[k][j][i].B1c;
00553 Bxc[j] = pG->U[k][j][i].B2c;
00554 Bxi[j] = pG->B2i[k][j][i];
00555 B2_x2Face[k][j][i] = pG->B2i[k][j][i];
00556 #endif
00557 #if (NSCALARS > 0)
00558 for (n=0; n<NSCALARS; n++) U1d[j].s[n] = pG->U[k][j][i].s[n];
00559 #endif
00560 }
00561
00562
00563
00564
00565
00566 for (j=js-nghost; j<=je+nghost; j++) {
00567 W[j] = Cons1D_to_Prim1D(&U1d[j],&Bxc[j]);
00568 }
00569
00570 lr_states(pG,W,Bxc,pG->dt,dx2,jl+1,ju-1,Wl,Wr,2);
00571
00572 #ifdef MHD
00573 #ifdef CYLINDRICAL
00574 dx2i = 1.0/dx2;
00575 #endif
00576 for (j=jl+1; j<=ju; j++) {
00577
00578 db1 = (rsf*pG->B1i[k ][j-1][i+1] - lsf*pG->B1i[k][j-1][i])*dx1i;
00579 db2 = ( pG->B2i[k ][j ][i ] - pG->B2i[k][j-1][i])*dx2i;
00580 db3 = ( pG->B3i[k+1][j-1][i ] - pG->B3i[k][j-1][i])*dx3i;
00581
00582 if(db2 >= 0.0){
00583 l1 = db2 < -db1 ? db2 : -db1;
00584 l1 = l1 > 0.0 ? l1 : 0.0;
00585
00586 l3 = db2 < -db3 ? db2 : -db3;
00587 l3 = l3 > 0.0 ? l3 : 0.0;
00588 }
00589 else{
00590 l1 = db2 > -db1 ? db2 : -db1;
00591 l1 = l1 < 0.0 ? l1 : 0.0;
00592
00593 l3 = db2 > -db3 ? db2 : -db3;
00594 l3 = l3 < 0.0 ? l3 : 0.0;
00595 }
00596
00597 MHD_src_By = (pG->U[k][j-1][i].M3/pG->U[k][j-1][i].d)*l3;
00598 MHD_src_Bz = (pG->U[k][j-1][i].M1/pG->U[k][j-1][i].d)*l1;
00599
00600 Wl[j].By += hdt*MHD_src_By;
00601 Wl[j].Bz += hdt*MHD_src_Bz;
00602
00603
00604 db1 = (rsf*pG->B1i[k ][j ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
00605 db2 = ( pG->B2i[k ][j+1][i ] - pG->B2i[k][j][i])*dx2i;
00606 db3 = ( pG->B3i[k+1][j ][i ] - pG->B3i[k][j][i])*dx3i;
00607
00608 if(db2 >= 0.0){
00609 l1 = db2 < -db1 ? db2 : -db1;
00610 l1 = l1 > 0.0 ? l1 : 0.0;
00611
00612 l3 = db2 < -db3 ? db2 : -db3;
00613 l3 = l3 > 0.0 ? l3 : 0.0;
00614 }
00615 else{
00616 l1 = db2 > -db1 ? db2 : -db1;
00617 l1 = l1 < 0.0 ? l1 : 0.0;
00618
00619 l3 = db2 > -db3 ? db2 : -db3;
00620 l3 = l3 < 0.0 ? l3 : 0.0;
00621 }
00622
00623 MHD_src_By = (pG->U[k][j][i].M3/pG->U[k][j][i].d)*l3;
00624 MHD_src_Bz = (pG->U[k][j][i].M1/pG->U[k][j][i].d)*l1;
00625
00626 Wr[j].By += hdt*MHD_src_By;
00627 Wr[j].Bz += hdt*MHD_src_Bz;
00628 }
00629 #endif
00630
00631
00632
00633
00634
00635 if (StaticGravPot != NULL){
00636 for (j=jl+1; j<=ju; j++) {
00637 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00638 phicr = (*StaticGravPot)(x1, x2 ,x3);
00639 phicl = (*StaticGravPot)(x1,(x2- pG->dx2),x3);
00640 phifc = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00641
00642 Wl[j].Vx -= dtodx2*(phifc - phicl);
00643 Wr[j].Vx -= dtodx2*(phicr - phifc);
00644 }
00645 }
00646
00647
00648
00649
00650
00651 #ifdef SELF_GRAVITY
00652 for (j=jl+1; j<=ju; j++) {
00653 Wl[j].Vx -= q2*(pG->Phi[k][j][i] - pG->Phi[k][j-1][i]);
00654 Wr[j].Vx -= q2*(pG->Phi[k][j][i] - pG->Phi[k][j-1][i]);
00655 }
00656 #endif
00657
00658
00659
00660
00661
00662 #ifndef BAROTROPIC
00663 if (CoolingFunc != NULL){
00664 for (j=jl+1; j<=ju; j++) {
00665 coolfl = (*CoolingFunc)(Wl[j].d,Wl[j].P,(0.5*pG->dt));
00666 coolfr = (*CoolingFunc)(Wr[j].d,Wr[j].P,(0.5*pG->dt));
00667
00668 Wl[j].P -= 0.5*pG->dt*Gamma_1*coolfl;
00669 Wr[j].P -= 0.5*pG->dt*Gamma_1*coolfr;
00670 }
00671 }
00672 #endif
00673
00674
00675
00676
00677
00678 #ifdef FEEDBACK
00679 for (j=jl+1; j<=ju; j++) {
00680 d1 = 1.0/W[j-1].d;
00681 Wl[j].Vx -= pG->Coup[k][j-1][i].fb2*d1;
00682 Wl[j].Vy -= pG->Coup[k][j-1][i].fb3*d1;
00683 Wl[j].Vz -= pG->Coup[k][j-1][i].fb1*d1;
00684
00685 d1 = 1.0/W[j].d;
00686 Wr[j].Vx -= pG->Coup[k][j][i].fb2*d1;
00687 Wr[j].Vy -= pG->Coup[k][j][i].fb3*d1;
00688 Wr[j].Vz -= pG->Coup[k][j][i].fb1*d1;
00689
00690 #ifndef BAROTROPIC
00691 Wl[i].P += pG->Coup[k][j-1][i].Eloss*Gamma_1;
00692 Wr[i].P += pG->Coup[k][j][i].Eloss*Gamma_1;
00693 #endif
00694
00695 }
00696 #endif
00697
00698
00699
00700
00701
00702
00703 for (j=jl+1; j<=ju; j++) {
00704 Ul_x2Face[k][j][i] = Prim1D_to_Cons1D(&Wl[j],&Bxi[j]);
00705 Ur_x2Face[k][j][i] = Prim1D_to_Cons1D(&Wr[j],&Bxi[j]);
00706
00707 #ifdef MHD
00708 Bx = B2_x2Face[k][j][i];
00709 #endif
00710 fluxes(Ul_x2Face[k][j][i],Ur_x2Face[k][j][i],Wl[j],Wr[j],Bx,
00711 &x2Flux[k][j][i]);
00712 }
00713 }
00714 }
00715
00716
00717
00718
00719
00720
00721
00722
00723 for (j=jl; j<=ju; j++) {
00724 for (i=il; i<=iu; i++) {
00725 for (k=ks-nghost; k<=ke+nghost; k++) {
00726 U1d[k].d = pG->U[k][j][i].d;
00727 U1d[k].Mx = pG->U[k][j][i].M3;
00728 U1d[k].My = pG->U[k][j][i].M1;
00729 U1d[k].Mz = pG->U[k][j][i].M2;
00730 #ifndef BAROTROPIC
00731 U1d[k].E = pG->U[k][j][i].E;
00732 #endif
00733 #ifdef MHD
00734 U1d[k].By = pG->U[k][j][i].B1c;
00735 U1d[k].Bz = pG->U[k][j][i].B2c;
00736 Bxc[k] = pG->U[k][j][i].B3c;
00737 Bxi[k] = pG->B3i[k][j][i];
00738 B3_x3Face[k][j][i] = pG->B3i[k][j][i];
00739 #endif
00740 #if (NSCALARS > 0)
00741 for (n=0; n<NSCALARS; n++) U1d[k].s[n] = pG->U[k][j][i].s[n];
00742 #endif
00743 }
00744
00745
00746
00747
00748
00749 for (k=ks-nghost; k<=ke+nghost; k++) {
00750 W[k] = Cons1D_to_Prim1D(&U1d[k],&Bxc[k]);
00751 }
00752
00753 lr_states(pG,W,Bxc,pG->dt,pG->dx3,kl+1,ku-1,Wl,Wr,3);
00754
00755 #ifdef MHD
00756 #ifdef CYLINDRICAL
00757 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
00758 dx2i = 1.0/(r[i]*pG->dx2);
00759 #endif
00760 for (k=kl+1; k<=ku; k++) {
00761
00762 db1 = (rsf*pG->B1i[k-1][j ][i+1] - lsf*pG->B1i[k-1][j][i])*dx1i;
00763 db2 = ( pG->B2i[k-1][j+1][i ] - pG->B2i[k-1][j][i])*dx2i;
00764 db3 = ( pG->B3i[k ][j ][i ] - pG->B3i[k-1][j][i])*dx3i;
00765
00766 if(db3 >= 0.0){
00767 l1 = db3 < -db1 ? db3 : -db1;
00768 l1 = l1 > 0.0 ? l1 : 0.0;
00769
00770 l2 = db3 < -db2 ? db3 : -db2;
00771 l2 = l2 > 0.0 ? l2 : 0.0;
00772 }
00773 else{
00774 l1 = db3 > -db1 ? db3 : -db1;
00775 l1 = l1 < 0.0 ? l1 : 0.0;
00776
00777 l2 = db3 > -db2 ? db3 : -db2;
00778 l2 = l2 < 0.0 ? l2 : 0.0;
00779 }
00780
00781 MHD_src_By = (pG->U[k-1][j][i].M1/pG->U[k-1][j][i].d)*l1;
00782 MHD_src_Bz = (pG->U[k-1][j][i].M2/pG->U[k-1][j][i].d)*l2;
00783
00784 Wl[k].By += hdt*MHD_src_By;
00785 Wl[k].Bz += hdt*MHD_src_Bz;
00786
00787
00788 db1 = (rsf*pG->B1i[k][j][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
00789 db2 = ( pG->B2i[k][j+1][i] - pG->B2i[k][j][i])*dx2i;
00790 db3 = ( pG->B3i[k+1][j][i] - pG->B3i[k][j][i])*dx3i;
00791
00792 if(db3 >= 0.0){
00793 l1 = db3 < -db1 ? db3 : -db1;
00794 l1 = l1 > 0.0 ? l1 : 0.0;
00795
00796 l2 = db3 < -db2 ? db3 : -db2;
00797 l2 = l2 > 0.0 ? l2 : 0.0;
00798 }
00799 else{
00800 l1 = db3 > -db1 ? db3 : -db1;
00801 l1 = l1 < 0.0 ? l1 : 0.0;
00802
00803 l2 = db3 > -db2 ? db3 : -db2;
00804 l2 = l2 < 0.0 ? l2 : 0.0;
00805 }
00806
00807 MHD_src_By = (pG->U[k][j][i].M1/pG->U[k][j][i].d)*l1;
00808 MHD_src_Bz = (pG->U[k][j][i].M2/pG->U[k][j][i].d)*l2;
00809
00810 Wr[k].By += hdt*MHD_src_By;
00811 Wr[k].Bz += hdt*MHD_src_Bz;
00812 }
00813 #endif
00814
00815
00816
00817
00818
00819 if (StaticGravPot != NULL){
00820 for (k=kl+1; k<=ku; k++) {
00821 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00822 phicr = (*StaticGravPot)(x1,x2, x3 );
00823 phicl = (*StaticGravPot)(x1,x2,(x3- pG->dx3));
00824 phifc = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
00825
00826 Wl[k].Vx -= dtodx3*(phifc - phicl);
00827 Wr[k].Vx -= dtodx3*(phicr - phifc);
00828 }
00829 }
00830
00831
00832
00833
00834
00835 #ifdef SELF_GRAVITY
00836 for (k=kl+1; k<=ku; k++) {
00837 Wl[k].Vx -= q3*(pG->Phi[k][j][i] - pG->Phi[k-1][j][i]);
00838 Wr[k].Vx -= q3*(pG->Phi[k][j][i] - pG->Phi[k-1][j][i]);
00839 }
00840 #endif
00841
00842
00843
00844
00845
00846 #ifndef BAROTROPIC
00847 if (CoolingFunc != NULL){
00848 for (k=kl+1; k<=ku; k++) {
00849 coolfl = (*CoolingFunc)(Wl[k].d,Wl[k].P,(0.5*pG->dt));
00850 coolfr = (*CoolingFunc)(Wr[k].d,Wr[k].P,(0.5*pG->dt));
00851
00852 Wl[k].P -= 0.5*pG->dt*Gamma_1*coolfl;
00853 Wr[k].P -= 0.5*pG->dt*Gamma_1*coolfr;
00854 }
00855 }
00856 #endif
00857
00858
00859
00860
00861
00862 #ifdef FEEDBACK
00863 for (k=kl+1; k<=ku; k++) {
00864 d1 = 1.0/W[k-1].d;
00865 Wl[k].Vx -= pG->Coup[k-1][j][i].fb3*d1;
00866 Wl[k].Vy -= pG->Coup[k-1][j][i].fb1*d1;
00867 Wl[k].Vz -= pG->Coup[k-1][j][i].fb2*d1;
00868
00869 d1 = 1.0/W[k].d;
00870 Wr[k].Vx -= pG->Coup[k][j][i].fb3*d1;
00871 Wr[k].Vy -= pG->Coup[k][j][i].fb1*d1;
00872 Wr[k].Vz -= pG->Coup[k][j][i].fb2*d1;
00873
00874 #ifndef BAROTROPIC
00875 Wl[i].P += pG->Coup[k-1][j][i].Eloss*Gamma_1;
00876 Wr[i].P += pG->Coup[k][j][i].Eloss*Gamma_1;
00877 #endif
00878 }
00879 #endif
00880
00881
00882
00883
00884
00885
00886 for (k=kl+1; k<=ku; k++) {
00887 Ul_x3Face[k][j][i] = Prim1D_to_Cons1D(&Wl[k],&Bxi[k]);
00888 Ur_x3Face[k][j][i] = Prim1D_to_Cons1D(&Wr[k],&Bxi[k]);
00889
00890 #ifdef MHD
00891 Bx = B3_x3Face[k][j][i];
00892 #endif
00893 fluxes(Ul_x3Face[k][j][i],Ur_x3Face[k][j][i],Wl[k],Wr[k],Bx,
00894 &x3Flux[k][j][i]);
00895 }
00896 }
00897 }
00898
00899
00900
00901
00902
00903
00904
00905
00906 #ifdef MHD
00907
00908 for (k=kl; k<=ku; k++) {
00909 for (j=jl; j<=ju; j++) {
00910 for (i=il; i<=iu; i++) {
00911 emf1_cc[k][j][i] = (pG->U[k][j][i].B2c*pG->U[k][j][i].M3 -
00912 pG->U[k][j][i].B3c*pG->U[k][j][i].M2)
00913 /pG->U[k][j][i].d;
00914 emf2_cc[k][j][i] = (pG->U[k][j][i].B3c*pG->U[k][j][i].M1 -
00915 pG->U[k][j][i].B1c*pG->U[k][j][i].M3)
00916 /pG->U[k][j][i].d;
00917 emf3_cc[k][j][i] = (pG->U[k][j][i].B1c*pG->U[k][j][i].M2 -
00918 pG->U[k][j][i].B2c*pG->U[k][j][i].M1)
00919 /pG->U[k][j][i].d;
00920 }
00921 }
00922 }
00923 integrate_emf1_corner(pG);
00924 integrate_emf2_corner(pG);
00925 integrate_emf3_corner(pG);
00926
00927
00928
00929
00930
00931 for (k=kl+1; k<=ku-1; k++) {
00932 for (j=jl+1; j<=ju-1; j++) {
00933 for (i=il+1; i<=iu-1; i++) {
00934 #ifdef CYLINDRICAL
00935 q2 = hdt/(ri[i]*pG->dx2);
00936 #endif
00937 B1_x1Face[k][j][i] += q3*(emf2[k+1][j ][i ] - emf2[k][j][i]) -
00938 q2*(emf3[k ][j+1][i ] - emf3[k][j][i]);
00939 B2_x2Face[k][j][i] += q1*(emf3[k ][j ][i+1] - emf3[k][j][i]) -
00940 q3*(emf1[k+1][j ][i ] - emf1[k][j][i]);
00941 #ifdef CYLINDRICAL
00942 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
00943 q2 = hdt/(r[i]*pG->dx2);
00944 #endif
00945 B3_x3Face[k][j][i] += q2*( emf1[k ][j+1][i ] - emf1[k][j][i]) -
00946 q1*(rsf*emf2[k ][j ][i+1] - lsf*emf2[k][j][i]);
00947 }
00948 #ifdef CYLINDRICAL
00949 q2 = hdt/(ri[iu]*pG->dx2);
00950 #endif
00951 B1_x1Face[k][j][iu] += q3*(emf2[k+1][j ][iu]-emf2[k][j][iu]) -
00952 q2*(emf3[k ][j+1][iu]-emf3[k][j][iu]);
00953 }
00954 for (i=il+1; i<=iu-1; i++) {
00955 B2_x2Face[k][ju][i] += q1*(emf3[k ][ju][i+1]-emf3[k][ju][i]) -
00956 q3*(emf1[k+1][ju][i ]-emf1[k][ju][i]);
00957 }
00958 }
00959 for (j=jl+1; j<=ju-1; j++) {
00960 for (i=il+1; i<=iu-1; i++) {
00961 #ifdef CYLINDRICAL
00962 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
00963 q2 = hdt/(r[i]*pG->dx2);
00964 #endif
00965 B3_x3Face[ku][j][i] += q2*( emf1[ku][j+1][i ] - emf1[ku][j][i]) -
00966 q1*(rsf*emf2[ku][j ][i+1] - lsf*emf2[ku][j][i]);
00967 }
00968 }
00969 #endif
00970
00971
00972
00973
00974
00975
00976
00977
00978 for (k=kl+1; k<=ku-1; k++) {
00979 for (j=jl+1; j<=ju-1; j++) {
00980 for (i=il+1; i<=iu; i++) {
00981 #ifdef CYLINDRICAL
00982 q2 = hdt/(r[i-1]*pG->dx2);
00983 #endif
00984 Ul_x1Face[k][j][i].d -=q2*(x2Flux[k][j+1][i-1].d -x2Flux[k][j][i-1].d );
00985 Ul_x1Face[k][j][i].Mx-=q2*(x2Flux[k][j+1][i-1].Mz-x2Flux[k][j][i-1].Mz);
00986 Ul_x1Face[k][j][i].My-=q2*(x2Flux[k][j+1][i-1].Mx-x2Flux[k][j][i-1].Mx);
00987 Ul_x1Face[k][j][i].Mz-=q2*(x2Flux[k][j+1][i-1].My-x2Flux[k][j][i-1].My);
00988 #ifndef BAROTROPIC
00989 Ul_x1Face[k][j][i].E -=q2*(x2Flux[k][j+1][i-1].E -x2Flux[k][j][i-1].E );
00990 #endif
00991 #ifdef MHD
00992
00993 Ul_x1Face[k][j][i].Bz+=q2*0.5*
00994 ((emf1[k ][j+1][i-1] - emf1[k ][j][i-1]) +
00995 (emf1[k+1][j+1][i-1] - emf1[k+1][j][i-1]));
00996 #endif
00997
00998 #ifdef CYLINDRICAL
00999 q2 = hdt/(r[i]*pG->dx2);
01000 #endif
01001 Ur_x1Face[k][j][i].d -=q2*(x2Flux[k][j+1][i ].d -x2Flux[k][j][i ].d );
01002 Ur_x1Face[k][j][i].Mx-=q2*(x2Flux[k][j+1][i ].Mz-x2Flux[k][j][i ].Mz);
01003 Ur_x1Face[k][j][i].My-=q2*(x2Flux[k][j+1][i ].Mx-x2Flux[k][j][i ].Mx);
01004 Ur_x1Face[k][j][i].Mz-=q2*(x2Flux[k][j+1][i ].My-x2Flux[k][j][i ].My);
01005 #ifndef BAROTROPIC
01006 Ur_x1Face[k][j][i].E -=q2*(x2Flux[k][j+1][i ].E -x2Flux[k][j][i ].E );
01007 #endif
01008 #ifdef MHD
01009
01010 Ur_x1Face[k][j][i].Bz+=q2*0.5*
01011 ((emf1[k ][j+1][i] - emf1[k ][j][i]) +
01012 (emf1[k+1][j+1][i] - emf1[k+1][j][i]));
01013 #endif
01014 #if (NSCALARS > 0)
01015 for (n=0; n<NSCALARS; n++) {
01016 Ul_x1Face[k][j][i].s[n] -=
01017 q2*(x2Flux[k][j+1][i-1].s[n] - x2Flux[k][j][i-1].s[n]);
01018 Ur_x1Face[k][j][i].s[n] -=
01019 q2*(x2Flux[k][j+1][i ].s[n] - x2Flux[k][j][i ].s[n]);
01020 }
01021 #endif
01022
01023
01024
01025
01026
01027
01028 Ul_x1Face[k][j][i].d -=q3*(x3Flux[k+1][j][i-1].d -x3Flux[k][j][i-1].d );
01029 Ul_x1Face[k][j][i].Mx-=q3*(x3Flux[k+1][j][i-1].My-x3Flux[k][j][i-1].My);
01030 Ul_x1Face[k][j][i].My-=q3*(x3Flux[k+1][j][i-1].Mz-x3Flux[k][j][i-1].Mz);
01031 Ul_x1Face[k][j][i].Mz-=q3*(x3Flux[k+1][j][i-1].Mx-x3Flux[k][j][i-1].Mx);
01032 #ifndef BAROTROPIC
01033 Ul_x1Face[k][j][i].E -=q3*(x3Flux[k+1][j][i-1].E -x3Flux[k][j][i-1].E );
01034 #endif
01035 #ifdef MHD
01036
01037 Ul_x1Face[k][j][i].By-=q3*0.5*
01038 ((emf1[k+1][j ][i-1] - emf1[k][j ][i-1]) +
01039 (emf1[k+1][j+1][i-1] - emf1[k][j+1][i-1]));
01040 #endif
01041
01042 Ur_x1Face[k][j][i].d -=q3*(x3Flux[k+1][j][i ].d -x3Flux[k][j][i ].d );
01043 Ur_x1Face[k][j][i].Mx-=q3*(x3Flux[k+1][j][i ].My-x3Flux[k][j][i ].My);
01044 Ur_x1Face[k][j][i].My-=q3*(x3Flux[k+1][j][i ].Mz-x3Flux[k][j][i ].Mz);
01045 Ur_x1Face[k][j][i].Mz-=q3*(x3Flux[k+1][j][i ].Mx-x3Flux[k][j][i ].Mx);
01046 #ifndef BAROTROPIC
01047 Ur_x1Face[k][j][i].E -=q3*(x3Flux[k+1][j][i ].E -x3Flux[k][j][i ].E );
01048 #endif
01049 #ifdef MHD
01050
01051 Ur_x1Face[k][j][i].By-=q3*0.5*
01052 ((emf1[k+1][j ][i] - emf1[k][j ][i]) +
01053 (emf1[k+1][j+1][i] - emf1[k][j+1][i]));
01054 #endif
01055 #if (NSCALARS > 0)
01056 for (n=0; n<NSCALARS; n++) {
01057 Ul_x1Face[k][j][i].s[n] -=
01058 q3*(x3Flux[k+1][j][i-1].s[n] - x3Flux[k][j][i-1].s[n]);
01059 Ur_x1Face[k][j][i].s[n] -=
01060 q3*(x3Flux[k+1][j][i ].s[n] - x3Flux[k][j][i ].s[n]);
01061 }
01062 #endif
01063 }
01064 }
01065 }
01066
01067
01068
01069
01070
01071
01072 #ifdef MHD
01073 for (k=kl+1; k<=ku-1; k++) {
01074 for (j=jl+1; j<=ju-1; j++) {
01075 for (i=il+1; i<=iu; i++) {
01076 #ifdef CYLINDRICAL
01077 rsf = ri[i]/r[i-1]; lsf = ri[i-1]/r[i-1];
01078 dx2i = 1.0/(r[i-1]*pG->dx2);
01079 #endif
01080 db1 = (rsf*pG->B1i[k ][j ][i ] - lsf*pG->B1i[k][j][i-1])*dx1i;
01081 db2 = ( pG->B2i[k ][j+1][i-1] - pG->B2i[k][j][i-1])*dx2i;
01082 db3 = ( pG->B3i[k+1][j ][i-1] - pG->B3i[k][j][i-1])*dx3i;
01083 B1 = pG->U[k][j][i-1].B1c;
01084 B2 = pG->U[k][j][i-1].B2c;
01085 B3 = pG->U[k][j][i-1].B3c;
01086 V2 = pG->U[k][j][i-1].M2/pG->U[k][j][i-1].d;
01087 V3 = pG->U[k][j][i-1].M3/pG->U[k][j][i-1].d;
01088
01089
01090 if(db1 > 0.0 && db2 < 0.0){
01091 mdb2 = db2 > -db1 ? db2 : -db1;
01092 }
01093 else if(db1 < 0.0 && db2 > 0.0){
01094 mdb2 = db2 < -db1 ? db2 : -db1;
01095 }
01096 else mdb2 = 0.0;
01097
01098
01099 if(db1 > 0.0 && db3 < 0.0){
01100 mdb3 = db3 > -db1 ? db3 : -db1;
01101 }
01102 else if(db1 < 0.0 && db3 > 0.0){
01103 mdb3 = db3 < -db1 ? db3 : -db1;
01104 }
01105 else mdb3 = 0.0;
01106
01107 Ul_x1Face[k][j][i].Mx += hdt*B1*db1;
01108 Ul_x1Face[k][j][i].My += hdt*B2*db1;
01109 Ul_x1Face[k][j][i].Mz += hdt*B3*db1;
01110 Ul_x1Face[k][j][i].By += hdt*V2*(-mdb3);
01111 Ul_x1Face[k][j][i].Bz += hdt*V3*(-mdb2);
01112 #ifndef BAROTROPIC
01113 Ul_x1Face[k][j][i].E += hdt*(B2*V2*(-mdb3) + B3*V3*(-mdb2) );
01114 #endif
01115
01116 #ifdef CYLINDRICAL
01117 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01118 dx2i = 1.0/(r[i]*pG->dx2);
01119 #endif
01120 db1 = (rsf*pG->B1i[k ][j ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
01121 db2 = ( pG->B2i[k ][j+1][i ] - pG->B2i[k][j][i])*dx2i;
01122 db3 = ( pG->B3i[k+1][j ][i ] - pG->B3i[k][j][i])*dx3i;
01123 B1 = pG->U[k][j][i].B1c;
01124 B2 = pG->U[k][j][i].B2c;
01125 B3 = pG->U[k][j][i].B3c;
01126 V2 = pG->U[k][j][i].M2/pG->U[k][j][i].d;
01127 V3 = pG->U[k][j][i].M3/pG->U[k][j][i].d;
01128
01129
01130 if(db1 > 0.0 && db2 < 0.0){
01131 mdb2 = db2 > -db1 ? db2 : -db1;
01132 }
01133 else if(db1 < 0.0 && db2 > 0.0){
01134 mdb2 = db2 < -db1 ? db2 : -db1;
01135 }
01136 else mdb2 = 0.0;
01137
01138
01139 if(db1 > 0.0 && db3 < 0.0){
01140 mdb3 = db3 > -db1 ? db3 : -db1;
01141 }
01142 else if(db1 < 0.0 && db3 > 0.0){
01143 mdb3 = db3 < -db1 ? db3 : -db1;
01144 }
01145 else mdb3 = 0.0;
01146
01147 Ur_x1Face[k][j][i].Mx += hdt*B1*db1;
01148 Ur_x1Face[k][j][i].My += hdt*B2*db1;
01149 Ur_x1Face[k][j][i].Mz += hdt*B3*db1;
01150 Ur_x1Face[k][j][i].By += hdt*V2*(-mdb3);
01151 Ur_x1Face[k][j][i].Bz += hdt*V3*(-mdb2);
01152 #ifndef BAROTROPIC
01153 Ur_x1Face[k][j][i].E += hdt*(B2*V2*(-mdb3) + B3*V3*(-mdb2) );
01154 #endif
01155 }
01156 }
01157 }
01158 #endif
01159
01160
01161
01162
01163
01164
01165
01166
01167 if (StaticGravPot != NULL){
01168 for (k=kl+1; k<=ku-1; k++) {
01169 for (j=jl+1; j<=ju-1; j++) {
01170 for (i=il+1; i<=iu; i++) {
01171 cc_pos(pG,i,j,k,&x1,&x2,&x3);
01172 phic = (*StaticGravPot)(x1, x2 ,x3);
01173 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01174 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01175
01176
01177 #ifdef CYLINDRICAL
01178 q2 = hdt/(r[i]*pG->dx2);
01179 #endif
01180 Ur_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i].d;
01181 #ifndef BAROTROPIC
01182 Ur_x1Face[k][j][i].E -= q2*(x2Flux[k][j ][i ].d*(phic - phil)
01183 + x2Flux[k][j+1][i ].d*(phir - phic));
01184 #endif
01185
01186 phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
01187 phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
01188
01189 Ur_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i].d;
01190 #ifndef BAROTROPIC
01191 Ur_x1Face[k][j][i].E -= q3*(x3Flux[k ][j][i ].d*(phic - phil)
01192 + x3Flux[k+1][j][i ].d*(phir - phic));
01193 #endif
01194
01195
01196 phic = (*StaticGravPot)((x1-pG->dx1), x2 ,x3);
01197 phir = (*StaticGravPot)((x1-pG->dx1),(x2+0.5*pG->dx2),x3);
01198 phil = (*StaticGravPot)((x1-pG->dx1),(x2-0.5*pG->dx2),x3);
01199
01200 #ifdef CYLINDRICAL
01201 q2 = hdt/(r[i-1]*pG->dx2);
01202 #endif
01203 Ul_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i-1].d;
01204 #ifndef BAROTROPIC
01205 Ul_x1Face[k][j][i].E -= q2*(x2Flux[k][j ][i-1].d*(phic - phil)
01206 + x2Flux[k][j+1][i-1].d*(phir - phic));
01207 #endif
01208
01209 phir = (*StaticGravPot)((x1-pG->dx1),x2,(x3+0.5*pG->dx3));
01210 phil = (*StaticGravPot)((x1-pG->dx1),x2,(x3-0.5*pG->dx3));
01211
01212 Ul_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i-1].d;
01213 #ifndef BAROTROPIC
01214 Ul_x1Face[k][j][i].E -= q3*(x3Flux[k ][j][i-1].d*(phic - phil)
01215 + x3Flux[k+1][j][i-1].d*(phir - phic));
01216 #endif
01217 }
01218 }
01219 }}
01220
01221
01222
01223
01224
01225
01226 #ifdef SELF_GRAVITY
01227 for (k=kl+1; k<=ku-1; k++) {
01228 for (j=jl+1; j<=ju-1; j++) {
01229 for (i=il+1; i<=iu; i++) {
01230 phic = pG->Phi[k][j][i];
01231 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
01232 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
01233
01234
01235 Ur_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i].d;
01236 #ifndef BAROTROPIC
01237 Ur_x1Face[k][j][i].E -= q2*(x2Flux[k][j ][i ].d*(phic - phil)
01238 + x2Flux[k][j+1][i ].d*(phir - phic));
01239 #endif
01240
01241 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
01242 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
01243
01244 Ur_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i].d;
01245 #ifndef BAROTROPIC
01246 Ur_x1Face[k][j][i].E -= q3*(x3Flux[k ][j][i ].d*(phic - phil)
01247 + x3Flux[k+1][j][i ].d*(phir - phic));
01248 #endif
01249
01250
01251 phic = pG->Phi[k][j][i-1];
01252 phir = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j+1][i-1]);
01253 phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j-1][i-1]);
01254
01255 Ul_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i-1].d;
01256 #ifndef BAROTROPIC
01257 Ul_x1Face[k][j][i].E -= q2*(x2Flux[k][j ][i-1].d*(phic - phil)
01258 + x2Flux[k][j+1][i-1].d*(phir - phic));
01259 #endif
01260
01261 phir = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k+1][j][i-1]);
01262 phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k-1][j][i-1]);
01263
01264 Ul_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i-1].d;
01265 #ifndef BAROTROPIC
01266 Ul_x1Face[k][j][i].E -= q3*(x3Flux[k ][j][i-1].d*(phic - phil)
01267 + x3Flux[k+1][j][i-1].d*(phir - phic));
01268 #endif
01269 }
01270 }
01271 }
01272 #endif
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282 for (k=kl+1; k<=ku-1; k++) {
01283 for (j=jl+1; j<=ju; j++) {
01284 for (i=il+1; i<=iu-1; i++) {
01285 #ifdef CYLINDRICAL
01286 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01287 #endif
01288 Ul_x2Face[k][j][i].d -=q1*(rsf*x1Flux[k][j-1][i+1].d -lsf*x1Flux[k][j-1][i].d );
01289 Ul_x2Face[k][j][i].Mx-=q1*(SQR(rsf)*x1Flux[k][j-1][i+1].My-SQR(lsf)*x1Flux[k][j-1][i].My);
01290 Ul_x2Face[k][j][i].My-=q1*(rsf*x1Flux[k][j-1][i+1].Mz-lsf*x1Flux[k][j-1][i].Mz);
01291 Ul_x2Face[k][j][i].Mz-=q1*(rsf*x1Flux[k][j-1][i+1].Mx-lsf*x1Flux[k][j-1][i].Mx);
01292 #ifndef BAROTROPIC
01293 Ul_x2Face[k][j][i].E -=q1*(rsf*x1Flux[k][j-1][i+1].E -lsf*x1Flux[k][j-1][i].E );
01294 #endif
01295 #ifdef MHD
01296
01297 Ul_x2Face[k][j][i].By-=q1*0.5*
01298 ((rsf*emf2[k ][j-1][i+1] - lsf*emf2[k ][j-1][i]) +
01299 (rsf*emf2[k+1][j-1][i+1] - lsf*emf2[k+1][j-1][i]));
01300 #endif
01301
01302 Ur_x2Face[k][j][i].d -=q1*(rsf*x1Flux[k][j ][i+1].d -lsf*x1Flux[k][j ][i].d );
01303 Ur_x2Face[k][j][i].Mx-=q1*(SQR(rsf)*x1Flux[k][j ][i+1].My-SQR(lsf)*x1Flux[k][j ][i].My);
01304 Ur_x2Face[k][j][i].My-=q1*(rsf*x1Flux[k][j ][i+1].Mz-lsf*x1Flux[k][j ][i].Mz);
01305 Ur_x2Face[k][j][i].Mz-=q1*(rsf*x1Flux[k][j ][i+1].Mx-lsf*x1Flux[k][j ][i].Mx);
01306 #ifndef BAROTROPIC
01307 Ur_x2Face[k][j][i].E -=q1*(rsf*x1Flux[k][j ][i+1].E -lsf*x1Flux[k][j ][i].E );
01308 #endif
01309 #ifdef MHD
01310
01311 Ur_x2Face[k][j][i].By-=q1*0.5*
01312 ((rsf*emf2[k ][j][i+1] - lsf*emf2[k ][j][i]) +
01313 (rsf*emf2[k+1][j][i+1] - lsf*emf2[k+1][j][i]));
01314 #endif
01315 #if (NSCALARS > 0)
01316 for (n=0; n<NSCALARS; n++) {
01317 Ul_x2Face[k][j][i].s[n] -=
01318 q1*(rsf*x1Flux[k][j-1][i+1].s[n] - lsf*x1Flux[k][j-1][i].s[n]);
01319 Ur_x2Face[k][j][i].s[n] -=
01320 q1*(rsf*x1Flux[k][j ][i+1].s[n] - lsf*x1Flux[k][j ][i].s[n]);
01321 }
01322 #endif
01323
01324
01325
01326
01327
01328
01329 Ul_x2Face[k][j][i].d -=q3*(x3Flux[k+1][j-1][i].d -x3Flux[k][j-1][i].d );
01330 Ul_x2Face[k][j][i].Mx-=q3*(x3Flux[k+1][j-1][i].Mz-x3Flux[k][j-1][i].Mz);
01331 Ul_x2Face[k][j][i].My-=q3*(x3Flux[k+1][j-1][i].Mx-x3Flux[k][j-1][i].Mx);
01332 Ul_x2Face[k][j][i].Mz-=q3*(x3Flux[k+1][j-1][i].My-x3Flux[k][j-1][i].My);
01333 #ifndef BAROTROPIC
01334 Ul_x2Face[k][j][i].E -=q3*(x3Flux[k+1][j-1][i].E -x3Flux[k][j-1][i].E );
01335 #endif
01336 #ifdef MHD
01337
01338 Ul_x2Face[k][j][i].Bz+=q3*0.5*
01339 (lsf*(emf2[k+1][j-1][i ] - emf2[k][j-1][i ]) +
01340 rsf*(emf2[k+1][j-1][i+1] - emf2[k][j-1][i+1]));
01341 #endif
01342
01343 Ur_x2Face[k][j][i].d -=q3*(x3Flux[k+1][j ][i].d -x3Flux[k][j ][i].d );
01344 Ur_x2Face[k][j][i].Mx-=q3*(x3Flux[k+1][j ][i].Mz-x3Flux[k][j ][i].Mz);
01345 Ur_x2Face[k][j][i].My-=q3*(x3Flux[k+1][j ][i].Mx-x3Flux[k][j ][i].Mx);
01346 Ur_x2Face[k][j][i].Mz-=q3*(x3Flux[k+1][j ][i].My-x3Flux[k][j ][i].My);
01347 #ifndef BAROTROPIC
01348 Ur_x2Face[k][j][i].E -=q3*(x3Flux[k+1][j ][i].E -x3Flux[k][j ][i].E );
01349 #endif
01350 #ifdef MHD
01351
01352 Ur_x2Face[k][j][i].Bz+=q3*0.5*
01353 (lsf*(emf2[k+1][j][i ] - emf2[k][j][i ]) +
01354 rsf*(emf2[k+1][j][i+1] - emf2[k][j][i+1]));
01355 #endif
01356 #if (NSCALARS > 0)
01357 for (n=0; n<NSCALARS; n++) {
01358 Ul_x2Face[k][j][i].s[n] -=
01359 q3*(x3Flux[k+1][j-1][i].s[n] - x3Flux[k][j-1][i].s[n]);
01360 Ur_x2Face[k][j][i].s[n] -=
01361 q3*(x3Flux[k+1][j ][i].s[n] - x3Flux[k][j ][i].s[n]);
01362 }
01363 #endif
01364 }
01365 }
01366 }
01367
01368
01369
01370
01371
01372
01373 #ifdef MHD
01374 for (k=kl+1; k<=ku-1; k++) {
01375 for (j=jl+1; j<=ju; j++) {
01376 for (i=il+1; i<=iu-1; i++) {
01377 #ifdef CYLINDRICAL
01378 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01379 #endif
01380 db1 = (rsf*pG->B1i[k ][j-1][i+1] - lsf*pG->B1i[k][j-1][i])*dx1i;
01381 db2 = ( pG->B2i[k ][j ][i ] - pG->B2i[k][j-1][i])*dx2i;
01382 db3 = ( pG->B3i[k+1][j-1][i ] - pG->B3i[k][j-1][i])*dx3i;
01383 B1 = pG->U[k][j-1][i].B1c;
01384 B2 = pG->U[k][j-1][i].B2c;
01385 B3 = pG->U[k][j-1][i].B3c;
01386 V1 = pG->U[k][j-1][i].M1/pG->U[k][j-1][i].d;
01387 V3 = pG->U[k][j-1][i].M3/pG->U[k][j-1][i].d;
01388
01389
01390 if(db2 > 0.0 && db1 < 0.0){
01391 mdb1 = db1 > -db2 ? db1 : -db2;
01392 }
01393 else if(db2 < 0.0 && db1 > 0.0){
01394 mdb1 = db1 < -db2 ? db1 : -db2;
01395 }
01396 else mdb1 = 0.0;
01397
01398
01399 if(db2 > 0.0 && db3 < 0.0){
01400 mdb3 = db3 > -db2 ? db3 : -db2;
01401 }
01402 else if(db2 < 0.0 && db3 > 0.0){
01403 mdb3 = db3 < -db2 ? db3 : -db2;
01404 }
01405 else mdb3 = 0.0;
01406
01407 Ul_x2Face[k][j][i].Mz += hdt*B1*db2;
01408 Ul_x2Face[k][j][i].Mx += hdt*B2*db2;
01409 Ul_x2Face[k][j][i].My += hdt*B3*db2;
01410 Ul_x2Face[k][j][i].By += hdt*V3*(-mdb1);
01411 Ul_x2Face[k][j][i].Bz += hdt*V1*(-mdb3);
01412 #ifndef BAROTROPIC
01413 Ul_x2Face[k][j][i].E += hdt*(B3*V3*(-mdb1) + B1*V1*(-mdb3) );
01414 #endif
01415
01416 db1 = (rsf*pG->B1i[k ][j ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
01417 db2 = ( pG->B2i[k ][j+1][i ] - pG->B2i[k][j][i])*dx2i;
01418 db3 = ( pG->B3i[k+1][j ][i ] - pG->B3i[k][j][i])*dx3i;
01419 B1 = pG->U[k][j][i].B1c;
01420 B2 = pG->U[k][j][i].B2c;
01421 B3 = pG->U[k][j][i].B3c;
01422 V1 = pG->U[k][j][i].M1/pG->U[k][j][i].d;
01423 V3 = pG->U[k][j][i].M3/pG->U[k][j][i].d;
01424
01425
01426 if(db2 > 0.0 && db1 < 0.0){
01427 mdb1 = db1 > -db2 ? db1 : -db2;
01428 }
01429 else if(db2 < 0.0 && db1 > 0.0){
01430 mdb1 = db1 < -db2 ? db1 : -db2;
01431 }
01432 else mdb1 = 0.0;
01433
01434
01435 if(db2 > 0.0 && db3 < 0.0){
01436 mdb3 = db3 > -db2 ? db3 : -db2;
01437 }
01438 else if(db2 < 0.0 && db3 > 0.0){
01439 mdb3 = db3 < -db2 ? db3 : -db2;
01440 }
01441 else mdb3 = 0.0;
01442
01443 Ur_x2Face[k][j][i].Mz += hdt*B1*db2;
01444 Ur_x2Face[k][j][i].Mx += hdt*B2*db2;
01445 Ur_x2Face[k][j][i].My += hdt*B3*db2;
01446 Ur_x2Face[k][j][i].By += hdt*V3*(-mdb1);
01447 Ur_x2Face[k][j][i].Bz += hdt*V1*(-mdb3);
01448 #ifndef BAROTROPIC
01449 Ur_x2Face[k][j][i].E += hdt*(B3*V3*(-mdb1) + B1*V1*(-mdb3) );
01450 #endif
01451 }
01452 }
01453 }
01454 #endif
01455
01456
01457
01458
01459
01460
01461
01462
01463 if (StaticGravPot != NULL){
01464 for (k=kl+1; k<=ku-1; k++) {
01465 for (j=jl+1; j<=ju; j++) {
01466 for (i=il+1; i<=iu-1; i++) {
01467 cc_pos(pG,i,j,k,&x1,&x2,&x3);
01468 phic = (*StaticGravPot)((x1 ),x2,x3);
01469 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
01470 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
01471
01472
01473 #ifdef CYLINDRICAL
01474 g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
01475 #ifdef FARGO
01476 g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
01477 #endif
01478 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01479 Ur_x2Face[k][j][i].Mz -= hdt*pG->U[k][j][i].d*g;
01480 #else
01481 Ur_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j][i].d;
01482 #endif
01483 #ifndef BAROTROPIC
01484 Ur_x2Face[k][j][i].E -= q1*(lsf*x1Flux[k][j ][i ].d*(phic - phil)
01485 + rsf*x1Flux[k][j ][i+1].d*(phir - phic));
01486 #endif
01487
01488 phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
01489 phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
01490
01491 Ur_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j][i].d;
01492 #ifndef BAROTROPIC
01493 Ur_x2Face[k][j][i].E -= q3*(x3Flux[k ][j ][i].d*(phic - phil)
01494 + x3Flux[k+1][j ][i].d*(phir - phic));
01495 #endif
01496
01497
01498 phic = (*StaticGravPot)((x1 ),(x2-pG->dx2),x3);
01499 phir = (*StaticGravPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
01500 phil = (*StaticGravPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
01501
01502 #ifdef CYLINDRICAL
01503 g = (*x1GravAcc)(x1vc(pG,i),(x2-pG->dx2),x3);
01504 #ifdef FARGO
01505 g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
01506 #endif
01507
01508 Ul_x2Face[k][j][i].Mz -= hdt*pG->U[k][j-1][i].d*g;
01509 #else
01510 Ul_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j-1][i].d;
01511 #endif
01512 #ifndef BAROTROPIC
01513 Ul_x2Face[k][j][i].E -= q1*(lsf*x1Flux[k][j-1][i ].d*(phic - phil)
01514 + rsf*x1Flux[k][j-1][i+1].d*(phir - phic));
01515 #endif
01516 phir = (*StaticGravPot)(x1,(x2-pG->dx2),(x3+0.5*pG->dx3));
01517 phil = (*StaticGravPot)(x1,(x2-pG->dx2),(x3-0.5*pG->dx3));
01518
01519 Ul_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j-1][i].d;
01520 #ifndef BAROTROPIC
01521 Ul_x2Face[k][j][i].E -= q3*(x3Flux[k ][j-1][i].d*(phic - phil)
01522 + x3Flux[k+1][j-1][i].d*(phir - phic));
01523 #endif
01524 }
01525 }
01526 }}
01527
01528
01529
01530
01531
01532
01533 #ifdef SELF_GRAVITY
01534 for (k=kl+1; k<=ku-1; k++) {
01535 for (j=jl+1; j<=ju; j++) {
01536 for (i=il+1; i<=iu-1; i++) {
01537 phic = pG->Phi[k][j][i];
01538 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
01539 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
01540
01541
01542 Ur_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j][i].d;
01543 #ifndef BAROTROPIC
01544 Ur_x2Face[k][j][i].E -= q1*(x1Flux[k][j][i ].d*(phic - phil)
01545 + x1Flux[k][j][i+1].d*(phir - phic));
01546 #endif
01547
01548 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
01549 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
01550
01551 Ur_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j][i].d;
01552 #ifndef BAROTROPIC
01553 Ur_x2Face[k][j][i].E -= q3*(x3Flux[k ][j][i].d*(phic - phil)
01554 + x3Flux[k+1][j][i].d*(phir - phic));
01555 #endif
01556
01557
01558 phic = pG->Phi[k][j-1][i];
01559 phir = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j-1][i+1]);
01560 phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j-1][i-1]);
01561
01562 Ul_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j-1][i].d;
01563 #ifndef BAROTROPIC
01564 Ul_x2Face[k][j][i].E -= q1*(x1Flux[k][j-1][i ].d*(phic - phil)
01565 + x1Flux[k][j-1][i+1].d*(phir - phic));
01566 #endif
01567 phir = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k+1][j-1][i]);
01568 phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k-1][j-1][i]);
01569
01570 Ul_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j-1][i].d;
01571 #ifndef BAROTROPIC
01572 Ul_x2Face[k][j][i].E -= q3*(x3Flux[k ][j-1][i].d*(phic - phil)
01573 + x3Flux[k+1][j-1][i].d*(phir - phic));
01574 #endif
01575 }
01576 }
01577 }
01578 #endif
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588 #ifdef SHEARING_BOX
01589 if (ShearingBoxPot != NULL){
01590 for (k=kl+1; k<=ku-1; k++) {
01591 for (j=jl+1; j<=ju; j++) {
01592 for (i=il+1; i<=iu-1; i++) {
01593 cc_pos(pG,i,j,k,&x1,&x2,&x3);
01594 phic = (*ShearingBoxPot)((x1 ),x2,x3);
01595 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
01596 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
01597
01598
01599 Ur_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j][i].d;
01600 #ifndef BAROTROPIC
01601 Ur_x2Face[k][j][i].E -= q1*(x1Flux[k][j ][i ].d*(phic - phil)
01602 + x1Flux[k][j ][i+1].d*(phir - phic));
01603 #endif
01604
01605 phir = (*ShearingBoxPot)(x1,x2,(x3+0.5*pG->dx3));
01606 phil = (*ShearingBoxPot)(x1,x2,(x3-0.5*pG->dx3));
01607
01608 Ur_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j][i].d;
01609 #ifndef BAROTROPIC
01610 Ur_x2Face[k][j][i].E -= q3*(x3Flux[k ][j ][i].d*(phic - phil)
01611 + x3Flux[k+1][j ][i].d*(phir - phic));
01612 #endif
01613
01614
01615 phic = (*ShearingBoxPot)((x1 ),(x2-pG->dx2),x3);
01616 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
01617 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
01618
01619 Ul_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j-1][i].d;
01620 #ifndef BAROTROPIC
01621 Ul_x2Face[k][j][i].E -= q1*(x1Flux[k][j-1][i ].d*(phic - phil)
01622 + x1Flux[k][j-1][i+1].d*(phir - phic));
01623 #endif
01624 phir = (*ShearingBoxPot)(x1,(x2-pG->dx2),(x3+0.5*pG->dx3));
01625 phil = (*ShearingBoxPot)(x1,(x2-pG->dx2),(x3-0.5*pG->dx3));
01626
01627 Ul_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j-1][i].d;
01628 #ifndef BAROTROPIC
01629 Ul_x2Face[k][j][i].E -= q3*(x3Flux[k ][j-1][i].d*(phic - phil)
01630 + x3Flux[k+1][j-1][i].d*(phir - phic));
01631 #endif
01632 }
01633 }
01634 }}
01635
01636 for (k=kl+1; k<=ku-1; k++) {
01637 for (j=jl+1; j<=ju; j++) {
01638 for (i=il+1; i<=iu-1; i++) {
01639 Ur_x2Face[k][j][i].Mz += pG->dt*Omega_0*pG->U[k][j][i].M2;
01640 Ul_x2Face[k][j][i].Mz += pG->dt*Omega_0*pG->U[k][j-1][i].M2;
01641 #ifdef FARGO
01642 Ur_x2Face[k][j][i].Mx += hdt*(qshear-2.)*Omega_0*pG->U[k][j][i].M1;
01643 Ul_x2Face[k][j][i].Mx += hdt*(qshear-2.)*Omega_0*pG->U[k][j-1][i].M1;
01644 #else
01645 Ur_x2Face[k][j][i].Mx -= pG->dt*Omega_0*pG->U[k][j][i].M1;
01646 Ul_x2Face[k][j][i].Mx -= pG->dt*Omega_0*pG->U[k][j-1][i].M1;
01647 #endif
01648 }
01649 }
01650 }
01651 #endif
01652
01653 #if defined(CYLINDRICAL) && defined(FARGO)
01654 for (k=kl+1; k<=ku-1; k++) {
01655 for (j=jl+1; j<=ju; j++) {
01656 for (i=il+1; i<=iu-1; i++) {
01657 Om = (*OrbitalProfile)(x1vc(pG,i));
01658 qshear = (*ShearProfile)(x1vc(pG,i));
01659 Ur_x2Face[k][j][i].Mz += pG->dt*Om*pG->U[k][j][i].M2;
01660 Ur_x2Face[k][j][i].Mx += hdt*(qshear-2.0)*Om*pG->U[k][j][i].M1;
01661
01662 Ul_x2Face[k][j][i].Mz += pG->dt*Om*pG->U[k][j-1][i].M2;
01663 Ul_x2Face[k][j][i].Mx += hdt*(qshear-2.0)*Om*pG->U[k][j-1][i].M1;
01664 }
01665 }
01666 }
01667 #endif
01668
01669
01670
01671
01672
01673 #ifdef CYLINDRICAL
01674 for (k=kl+1; k<=ku-1; k++) {
01675 for (j=jl+1; j<=ju; j++) {
01676 for (i=il+1; i<=iu-1; i++) {
01677 Ur_x2Face[k][j][i].Mz += hdt*geom_src[k][j ][i];
01678 Ul_x2Face[k][j][i].Mz += hdt*geom_src[k][j-1][i];
01679 }
01680 }
01681 }
01682 #endif
01683
01684
01685
01686
01687
01688
01689
01690
01691 for (k=kl+1; k<=ku; k++) {
01692 for (j=jl+1; j<=ju-1; j++) {
01693 for (i=il+1; i<=iu-1; i++) {
01694 #ifdef CYLINDRICAL
01695 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01696 q2 = hdt/(r[i]*pG->dx2);
01697 #endif
01698 Ul_x3Face[k][j][i].d -=q1*(rsf*x1Flux[k-1][j][i+1].d -lsf*x1Flux[k-1][j][i].d );
01699 Ul_x3Face[k][j][i].Mx-=q1*(rsf*x1Flux[k-1][j][i+1].Mz-lsf*x1Flux[k-1][j][i].Mz);
01700 Ul_x3Face[k][j][i].My-=q1*(rsf*x1Flux[k-1][j][i+1].Mx-lsf*x1Flux[k-1][j][i].Mx);
01701 Ul_x3Face[k][j][i].Mz-=q1*(SQR(rsf)*x1Flux[k-1][j][i+1].My-SQR(lsf)*x1Flux[k-1][j][i].My);
01702 #ifndef BAROTROPIC
01703 Ul_x3Face[k][j][i].E -=q1*(rsf*x1Flux[k-1][j][i+1].E -lsf*x1Flux[k-1][j][i].E );
01704 #endif
01705 #ifdef MHD
01706
01707 Ul_x3Face[k][j][i].Bz+=q1*0.5*
01708 ((emf3[k-1][j ][i+1] - emf3[k-1][j ][i]) +
01709 (emf3[k-1][j+1][i+1] - emf3[k-1][j+1][i]));
01710 #endif
01711
01712 Ur_x3Face[k][j][i].d -=q1*(rsf*x1Flux[k ][j][i+1].d -lsf*x1Flux[k ][j][i].d );
01713 Ur_x3Face[k][j][i].Mx-=q1*(rsf*x1Flux[k ][j][i+1].Mz-lsf*x1Flux[k ][j][i].Mz);
01714 Ur_x3Face[k][j][i].My-=q1*(rsf*x1Flux[k ][j][i+1].Mx-lsf*x1Flux[k ][j][i].Mx);
01715 Ur_x3Face[k][j][i].Mz-=q1*(SQR(rsf)*x1Flux[k ][j][i+1].My-SQR(lsf)*x1Flux[k ][j][i].My);
01716 #ifndef BAROTROPIC
01717 Ur_x3Face[k][j][i].E -=q1*(rsf*x1Flux[k ][j][i+1].E -lsf*x1Flux[k ][j][i].E );
01718 #endif
01719 #ifdef MHD
01720
01721 Ur_x3Face[k][j][i].Bz+=q1*0.5*
01722 ((emf3[k][j ][i+1] - emf3[k][j ][i]) +
01723 (emf3[k][j+1][i+1] - emf3[k][j+1][i]));
01724 #endif
01725 #if (NSCALARS > 0)
01726 for (n=0; n<NSCALARS; n++) {
01727 Ul_x3Face[k][j][i].s[n] -=
01728 q1*(rsf*x1Flux[k-1][j][i+1].s[n] - lsf*x1Flux[k-1][j][i].s[n]);
01729 Ur_x3Face[k][j][i].s[n] -=
01730 q1*(rsf*x1Flux[k ][j][i+1].s[n] - lsf*x1Flux[k ][j][i].s[n]);
01731 }
01732 #endif
01733
01734
01735
01736
01737
01738
01739 Ul_x3Face[k][j][i].d -=q2*(x2Flux[k-1][j+1][i].d -x2Flux[k-1][j][i].d );
01740 Ul_x3Face[k][j][i].Mx-=q2*(x2Flux[k-1][j+1][i].My-x2Flux[k-1][j][i].My);
01741 Ul_x3Face[k][j][i].My-=q2*(x2Flux[k-1][j+1][i].Mz-x2Flux[k-1][j][i].Mz);
01742 Ul_x3Face[k][j][i].Mz-=q2*(x2Flux[k-1][j+1][i].Mx-x2Flux[k-1][j][i].Mx);
01743 #ifndef BAROTROPIC
01744 Ul_x3Face[k][j][i].E -=q2*(x2Flux[k-1][j+1][i].E -x2Flux[k-1][j][i].E );
01745 #endif
01746 #ifdef MHD
01747
01748 Ul_x3Face[k][j][i].By-=q2*0.5*
01749 ((emf3[k-1][j+1][i ] - emf3[k-1][j][i ]) +
01750 (emf3[k-1][j+1][i+1] - emf3[k-1][j][i+1]));
01751 #endif
01752
01753 Ur_x3Face[k][j][i].d -=q2*(x2Flux[k ][j+1][i].d -x2Flux[k ][j][i].d );
01754 Ur_x3Face[k][j][i].Mx-=q2*(x2Flux[k ][j+1][i].My-x2Flux[k ][j][i].My);
01755 Ur_x3Face[k][j][i].My-=q2*(x2Flux[k ][j+1][i].Mz-x2Flux[k ][j][i].Mz);
01756 Ur_x3Face[k][j][i].Mz-=q2*(x2Flux[k ][j+1][i].Mx-x2Flux[k ][j][i].Mx);
01757 #ifndef BAROTROPIC
01758 Ur_x3Face[k][j][i].E -=q2*(x2Flux[k ][j+1][i].E -x2Flux[k ][j][i].E );
01759 #endif
01760 #ifdef MHD
01761
01762 Ur_x3Face[k][j][i].By-=q2*0.5*
01763 ((emf3[k][j+1][i ] - emf3[k][j][i ]) +
01764 (emf3[k][j+1][i+1] - emf3[k][j][i+1]));
01765 #endif
01766 #if (NSCALARS > 0)
01767 for (n=0; n<NSCALARS; n++) {
01768 Ul_x3Face[k][j][i].s[n] -=
01769 q2*(x2Flux[k-1][j+1][i].s[n] - x2Flux[k-1][j][i].s[n]);
01770 Ur_x3Face[k][j][i].s[n] -=
01771 q2*(x2Flux[k ][j+1][i].s[n] - x2Flux[k ][j][i].s[n]);
01772 }
01773 #endif
01774 }
01775 }
01776 }
01777
01778
01779
01780
01781
01782
01783 #ifdef MHD
01784 for (k=kl+1; k<=ku; k++) {
01785 for (j=jl+1; j<=ju-1; j++) {
01786 for (i=il+1; i<=iu-1; i++) {
01787 #ifdef CYLINDRICAL
01788 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01789 #endif
01790 db1 = (rsf*pG->B1i[k-1][j ][i+1] - lsf*pG->B1i[k-1][j][i])*dx1i;
01791 db2 = ( pG->B2i[k-1][j+1][i ] - pG->B2i[k-1][j][i])*dx2i;
01792 db3 = ( pG->B3i[k ][j ][i ] - pG->B3i[k-1][j][i])*dx3i;
01793 B1 = pG->U[k-1][j][i].B1c;
01794 B2 = pG->U[k-1][j][i].B2c;
01795 B3 = pG->U[k-1][j][i].B3c;
01796 V1 = pG->U[k-1][j][i].M1/pG->U[k-1][j][i].d;
01797 V2 = pG->U[k-1][j][i].M2/pG->U[k-1][j][i].d;
01798
01799
01800 if(db3 > 0.0 && db1 < 0.0){
01801 mdb1 = db1 > -db3 ? db1 : -db3;
01802 }
01803 else if(db3 < 0.0 && db1 > 0.0){
01804 mdb1 = db1 < -db3 ? db1 : -db3;
01805 }
01806 else mdb1 = 0.0;
01807
01808
01809 if(db3 > 0.0 && db2 < 0.0){
01810 mdb2 = db2 > -db3 ? db2 : -db3;
01811 }
01812 else if(db3 < 0.0 && db2 > 0.0){
01813 mdb2 = db2 < -db3 ? db2 : -db3;
01814 }
01815 else mdb2 = 0.0;
01816
01817 Ul_x3Face[k][j][i].My += hdt*B1*db3;
01818 Ul_x3Face[k][j][i].Mz += hdt*B2*db3;
01819 Ul_x3Face[k][j][i].Mx += hdt*B3*db3;
01820 Ul_x3Face[k][j][i].By += hdt*V1*(-mdb2);
01821 Ul_x3Face[k][j][i].Bz += hdt*V2*(-mdb1);
01822 #ifndef BAROTROPIC
01823 Ul_x3Face[k][j][i].E += hdt*(B1*V1*(-mdb2) + B2*V2*(-mdb1) );
01824 #endif
01825
01826 db1 = (rsf*pG->B1i[k ][j ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
01827 db2 = ( pG->B2i[k ][j+1][i ] - pG->B2i[k][j][i])*dx2i;
01828 db3 = ( pG->B3i[k+1][j ][i ] - pG->B3i[k][j][i])*dx3i;
01829 B1 = pG->U[k][j][i].B1c;
01830 B2 = pG->U[k][j][i].B2c;
01831 B3 = pG->U[k][j][i].B3c;
01832 V1 = pG->U[k][j][i].M1/pG->U[k][j][i].d;
01833 V2 = pG->U[k][j][i].M2/pG->U[k][j][i].d;
01834
01835
01836 if(db3 > 0.0 && db1 < 0.0){
01837 mdb1 = db1 > -db3 ? db1 : -db3;
01838 }
01839 else if(db3 < 0.0 && db1 > 0.0){
01840 mdb1 = db1 < -db3 ? db1 : -db3;
01841 }
01842 else mdb1 = 0.0;
01843
01844
01845 if(db3 > 0.0 && db2 < 0.0){
01846 mdb2 = db2 > -db3 ? db2 : -db3;
01847 }
01848 else if(db3 < 0.0 && db2 > 0.0){
01849 mdb2 = db2 < -db3 ? db2 : -db3;
01850 }
01851 else mdb2 = 0.0;
01852
01853 Ur_x3Face[k][j][i].My += hdt*B1*db3;
01854 Ur_x3Face[k][j][i].Mz += hdt*B2*db3;
01855 Ur_x3Face[k][j][i].Mx += hdt*B3*db3;
01856 Ur_x3Face[k][j][i].By += hdt*V1*(-mdb2);
01857 Ur_x3Face[k][j][i].Bz += hdt*V2*(-mdb1);
01858 #ifndef BAROTROPIC
01859 Ur_x3Face[k][j][i].E += hdt*(B1*V1*(-mdb2) + B2*V2*(-mdb1) );
01860 #endif
01861 }
01862 }
01863 }
01864 #endif
01865
01866
01867
01868
01869
01870
01871
01872
01873 if (StaticGravPot != NULL){
01874 for (k=kl+1; k<=ku; k++) {
01875 for (j=jl+1; j<=ju-1; j++) {
01876 for (i=il+1; i<=iu-1; i++) {
01877 cc_pos(pG,i,j,k,&x1,&x2,&x3);
01878 phic = (*StaticGravPot)((x1 ),x2,x3);
01879 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
01880 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
01881
01882
01883 #ifdef CYLINDRICAL
01884 g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
01885 #ifdef FARGO
01886 g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
01887 #endif
01888 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01889 q2 = hdt/(r[i]*pG->dx2);
01890 Ur_x3Face[k][j][i].My -= hdt*pG->U[k][j][i].d*g;
01891 #else
01892 Ur_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k][j][i].d;
01893 #endif
01894 #ifndef BAROTROPIC
01895 Ur_x3Face[k][j][i].E -= q1*(lsf*x1Flux[k ][j][i ].d*(phic - phil)
01896 + rsf*x1Flux[k ][j][i+1].d*(phir - phic));
01897 #endif
01898
01899 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01900 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01901
01902 Ur_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k][j][i].d;
01903 #ifndef BAROTROPIC
01904 Ur_x3Face[k][j][i].E -= q2*(x2Flux[k ][j ][i].d*(phic - phil)
01905 + x2Flux[k ][j+1][i].d*(phir - phic));
01906 #endif
01907
01908
01909 phic = (*StaticGravPot)((x1 ),x2,(x3-pG->dx3));
01910 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,(x3-pG->dx3));
01911 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,(x3-pG->dx3));
01912
01913 #ifdef CYLINDRICAL
01914 g = (*x1GravAcc)(x1vc(pG,i),x2,(x3-pG->dx3));
01915 #ifdef FARGO
01916 g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
01917 #endif
01918 Ul_x3Face[k][j][i].My -= hdt*pG->U[k-1][j][i].d*g;
01919 #else
01920 Ul_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k-1][j][i].d;
01921 #endif
01922 #ifndef BAROTROPIC
01923 Ul_x3Face[k][j][i].E -= q1*(lsf*x1Flux[k-1][j][i ].d*(phic - phil)
01924 + rsf*x1Flux[k-1][j][i+1].d*(phir - phic));
01925 #endif
01926
01927 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),(x3-pG->dx3));
01928 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),(x3-pG->dx3));
01929
01930 Ul_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k-1][j][i].d;
01931 #ifndef BAROTROPIC
01932 Ul_x3Face[k][j][i].E -= q2*(x2Flux[k-1][j ][i].d*(phic - phil)
01933 + x2Flux[k-1][j+1][i].d*(phir - phic));
01934 #endif
01935 }
01936 }
01937 }}
01938
01939
01940
01941
01942
01943
01944 #ifdef SELF_GRAVITY
01945 for (k=kl+1; k<=ku; k++) {
01946 for (j=jl+1; j<=ju-1; j++) {
01947 for (i=il+1; i<=iu-1; i++) {
01948 phic = pG->Phi[k][j][i];
01949 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
01950 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
01951
01952
01953 Ur_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k][j][i].d;
01954 #ifndef BAROTROPIC
01955 Ur_x3Face[k][j][i].E -= q1*(x1Flux[k][j][i ].d*(phic - phil)
01956 + x1Flux[k][j][i+1].d*(phir - phic));
01957 #endif
01958
01959 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
01960 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
01961
01962 Ur_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k][j][i].d;
01963 #ifndef BAROTROPIC
01964 Ur_x3Face[k][j][i].E -= q2*(x2Flux[k][j ][i].d*(phic - phil)
01965 + x2Flux[k][j+1][i].d*(phir - phic));
01966 #endif
01967
01968
01969 phic = pG->Phi[k-1][j][i];
01970 phir = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j][i+1]);
01971 phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j][i-1]);
01972
01973 Ul_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k-1][j][i].d;
01974 #ifndef BAROTROPIC
01975 Ul_x3Face[k][j][i].E -= q1*(x1Flux[k-1][j][i ].d*(phic - phil)
01976 + x1Flux[k-1][j][i+1].d*(phir - phic));
01977 #endif
01978
01979 phir = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j+1][i]);
01980 phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j-1][i]);
01981
01982 Ul_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k-1][j][i].d;
01983 #ifndef BAROTROPIC
01984 Ul_x3Face[k][j][i].E -= q2*(x2Flux[k-1][j ][i].d*(phic - phil)
01985 + x2Flux[k-1][j+1][i].d*(phir - phic));
01986 #endif
01987 }
01988 }
01989 }
01990 #endif
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000 #ifdef SHEARING_BOX
02001 if (ShearingBoxPot != NULL){
02002 for (k=kl+1; k<=ku; k++) {
02003 for (j=jl+1; j<=ju-1; j++) {
02004 for (i=il+1; i<=iu-1; i++) {
02005 cc_pos(pG,i,j,k,&x1,&x2,&x3);
02006 phic = (*ShearingBoxPot)((x1 ),x2,x3);
02007 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
02008 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
02009
02010
02011 Ur_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k][j][i].d;
02012 #ifndef BAROTROPIC
02013 Ur_x3Face[k][j][i].E -= q1*(x1Flux[k ][j][i ].d*(phic - phil)
02014 + x1Flux[k ][j][i+1].d*(phir - phic));
02015 #endif
02016
02017 phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
02018 phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
02019
02020 Ur_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k][j][i].d;
02021 #ifndef BAROTROPIC
02022 Ur_x3Face[k][j][i].E -= q2*(x2Flux[k ][j ][i].d*(phic - phil)
02023 + x2Flux[k ][j+1][i].d*(phir - phic));
02024 #endif
02025
02026
02027 phic = (*ShearingBoxPot)((x1 ),x2,(x3-pG->dx3));
02028 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,(x3-pG->dx3));
02029 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,(x3-pG->dx3));
02030
02031 Ul_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k-1][j][i].d;
02032 #ifndef BAROTROPIC
02033 Ul_x3Face[k][j][i].E -= q1*(x1Flux[k-1][j][i ].d*(phic - phil)
02034 + x1Flux[k-1][j][i+1].d*(phir - phic));
02035 #endif
02036
02037 phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),(x3-pG->dx3));
02038 phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),(x3-pG->dx3));
02039
02040 Ul_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k-1][j][i].d;
02041 #ifndef BAROTROPIC
02042 Ul_x3Face[k][j][i].E -= q2*(x2Flux[k-1][j ][i].d*(phic - phil)
02043 + x2Flux[k-1][j+1][i].d*(phir - phic));
02044 #endif
02045 }
02046 }
02047 }}
02048
02049 for (k=kl+1; k<=ku; k++) {
02050 for (j=jl+1; j<=ju-1; j++) {
02051 for (i=il+1; i<=iu-1; i++) {
02052 Ur_x3Face[k][j][i].My += pG->dt*Omega_0*pG->U[k][j][i].M2;
02053 Ul_x3Face[k][j][i].My += pG->dt*Omega_0*pG->U[k-1][j][i].M2;
02054 #ifdef FARGO
02055 Ur_x3Face[k][j][i].Mz += hdt*(qshear-2.)*Omega_0*pG->U[k][j][i].M1;
02056 Ul_x3Face[k][j][i].Mz += hdt*(qshear-2.)*Omega_0*pG->U[k-1][j][i].M1;
02057 #else
02058 Ur_x3Face[k][j][i].Mz -= pG->dt*Omega_0*pG->U[k][j][i].M1;
02059 Ul_x3Face[k][j][i].Mz -= pG->dt*Omega_0*pG->U[k-1][j][i].M1;
02060 #endif
02061 }
02062 }
02063 }
02064 #endif
02065
02066 #if defined(CYLINDRICAL) && defined(FARGO)
02067 for (k=kl+1; k<=ku; k++) {
02068 for (j=jl+1; j<=ju-1; j++) {
02069 for (i=il+1; i<=iu-1; i++) {
02070 Om = (*OrbitalProfile)(x1vc(pG,i));
02071 qshear = (*ShearProfile)(x1vc(pG,i));
02072
02073 Ur_x3Face[k][j][i].My += pG->dt*Om*pG->U[k][j][i].M2;
02074 Ur_x3Face[k][j][i].Mz += hdt*(qshear-2.0)*Om*pG->U[k][j][i].M1;
02075
02076 Ul_x3Face[k][j][i].My += pG->dt*Om*pG->U[k-1][j][i].M2;
02077 Ul_x3Face[k][j][i].Mz += hdt*(qshear-2.0)*Om*pG->U[k-1][j][i].M1;
02078 }
02079 }
02080 }
02081 #endif
02082
02083
02084
02085
02086 #ifdef CYLINDRICAL
02087 for (k=kl+1; k<=ku; k++) {
02088 for (j=jl+1; j<=ju-1; j++) {
02089 for (i=il+1; i<=iu-1; i++) {
02090 Ul_x3Face[k][j][i].My += hdt*geom_src[k-1][j][i];
02091 Ur_x3Face[k][j][i].My += hdt*geom_src[k ][j][i];
02092 }
02093 }
02094 }
02095 #endif
02096
02097
02098
02099
02100
02101
02102 #ifndef MHD
02103 #ifndef PARTICLES
02104 if ((StaticGravPot != NULL) || (CoolingFunc != NULL))
02105 #endif
02106 #endif
02107 {
02108 for (k=kl+1; k<=ku-1; k++) {
02109 for (j=jl+1; j<=ju-1; j++) {
02110 for (i=il+1; i<=iu-1; i++) {
02111 #ifdef CYLINDRICAL
02112 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
02113 q2 = hdt/(r[i]*pG->dx2);
02114 #endif
02115 dhalf[k][j][i] = pG->U[k][j][i].d
02116 - q1*(rsf*x1Flux[k ][j ][i+1].d - lsf*x1Flux[k][j][i].d)
02117 - q2*( x2Flux[k ][j+1][i ].d - x2Flux[k][j][i].d)
02118 - q3*( x3Flux[k+1][j ][i ].d - x3Flux[k][j][i].d);
02119 #ifdef PARTICLES
02120 pG->Coup[k][j][i].grid_d = dhalf[k][j][i];
02121 #endif
02122 }
02123 }
02124 }
02125 }
02126
02127
02128
02129
02130
02131 #ifndef MHD
02132 #ifndef PARTICLES
02133 if (CoolingFunc != NULL)
02134 #endif
02135 #endif
02136 {
02137 for (k=kl+1; k<=ku-1; k++) {
02138 for (j=jl+1; j<=ju-1; j++) {
02139 for (i=il+1; i<=iu-1; i++) {
02140 #ifdef CYLINDRICAL
02141 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
02142 q2 = hdt/(r[i]/pG->dx2);
02143 #endif
02144 M1h = pG->U[k][j][i].M1
02145 - q1*(rsf*x1Flux[k ][j ][i+1].Mx - lsf*x1Flux[k][j][i].Mx)
02146 - q2*( x2Flux[k ][j+1][i ].Mz - x2Flux[k][j][i].Mz)
02147 - q3*( x3Flux[k+1][j ][i ].My - x3Flux[k][j][i].My);
02148
02149 M2h = pG->U[k][j][i].M2
02150 - q1*(SQR(rsf)*x1Flux[k ][j ][i+1].My - SQR(lsf)*x1Flux[k][j][i].My)
02151 - q2*( x2Flux[k ][j+1][i ].Mx - x2Flux[k][j][i].Mx)
02152 - q3*( x3Flux[k+1][j ][i ].Mz - x3Flux[k][j][i].Mz);
02153
02154 M3h = pG->U[k][j][i].M3
02155 - q1*(lsf*x1Flux[k ][j ][i+1].Mz - rsf*x1Flux[k][j][i].Mz)
02156 - q2*( x2Flux[k ][j+1][i ].My - x2Flux[k][j][i].My)
02157 - q3*( x3Flux[k+1][j ][i ].Mx - x3Flux[k][j][i].Mx);
02158
02159 #ifndef BAROTROPIC
02160 Eh = pG->U[k][j][i].E
02161 - q1*(rsf*x1Flux[k ][j ][i+1].E - lsf*x1Flux[k][j][i].E)
02162 - q2*( x2Flux[k ][j+1][i ].E - x2Flux[k][j][i].E)
02163 - q3*( x3Flux[k+1][j ][i ].E - x3Flux[k][j][i].E);
02164 #endif
02165
02166
02167 if (StaticGravPot != NULL){
02168 cc_pos(pG,i,j,k,&x1,&x2,&x3);
02169 #ifdef CYLINDRICAL
02170 g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
02171 #ifdef FARGO
02172 g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
02173 #endif
02174 M1h -= hdt*pG->U[k][j][i].d*g;
02175 #else
02176 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
02177 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
02178 M1h -= q1*(phir-phil)*pG->U[k][j][i].d;
02179 #endif
02180
02181 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
02182 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
02183 M2h -= q2*(phir-phil)*pG->U[k][j][i].d;
02184
02185 phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
02186 phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
02187 M3h -= q3*(phir-phil)*pG->U[k][j][i].d;
02188 }
02189
02190
02191 #ifdef SELF_GRAVITY
02192 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
02193 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
02194 M1h -= q1*(phir-phil)*pG->U[k][j][i].d;
02195
02196 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
02197 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
02198 M2h -= q2*(phir-phil)*pG->U[k][j][i].d;
02199
02200 phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
02201 phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
02202 M3h -= q3*(phir-phil)*pG->U[k][j][i].d;
02203 #endif
02204
02205
02206 #ifdef SHEARING_BOX
02207 if (ShearingBoxPot != NULL){
02208 cc_pos(pG,i,j,k,&x1,&x2,&x3);
02209 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
02210 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
02211 M1h -= q1*(phir-phil)*pG->U[k][j][i].d;
02212
02213 phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
02214 phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
02215 M2h -= q2*(phir-phil)*pG->U[k][j][i].d;
02216
02217 phir = (*ShearingBoxPot)(x1,x2,(x3+0.5*pG->dx3));
02218 phil = (*ShearingBoxPot)(x1,x2,(x3-0.5*pG->dx3));
02219 M3h -= q3*(phir-phil)*pG->U[k][j][i].d;
02220 }
02221
02222 M1h += pG->dt*Omega_0*pG->U[k][j][i].M2;
02223 #ifdef FARGO
02224 M2h += hdt*(qshear-2.)*Omega_0*pG->U[k][j][i].M1;
02225 #else
02226 M2h -= pG->dt*Omega_0*pG->U[k][j][i].M1;
02227 #endif
02228 #endif
02229 #if defined(CYLINDRICAL) && defined(FARGO)
02230 Om = (*OrbitalProfile)(x1vc(pG,i));
02231 qshear = (*ShearProfile)(x1vc(pG,i));
02232 M1h += hdt*2.0*Om*pG->U[k][j][i].M2;
02233 M2h += hdt*Om*(qshear-2.0)*pG->U[k][j][i].M1;
02234 #endif
02235
02236
02237 #ifdef FEEDBACK
02238 M1h -= pG->Coup[k][j][i].fb1;
02239 M2h -= pG->Coup[k][j][i].fb2;
02240 M3h -= pG->Coup[k][j][i].fb3;
02241 #endif
02242
02243
02244 #ifdef CYLINDRICAL
02245 M1h += hdt*geom_src[k][j][i];
02246 #endif
02247
02248 #ifndef BAROTROPIC
02249 phalf[k][j][i] = Eh - 0.5*(M1h*M1h + M2h*M2h + M3h*M3h)/dhalf[k][j][i];
02250 #endif
02251
02252 #ifdef MHD
02253 B1ch = 0.5*(lsf*B1_x1Face[k][j][i] + rsf*B1_x1Face[k ][j ][i+1]);
02254 B2ch = 0.5*( B2_x2Face[k][j][i] + B2_x2Face[k ][j+1][i ]);
02255 B3ch = 0.5*( B3_x3Face[k][j][i] + B3_x3Face[k+1][j ][i ]);
02256 emf1_cc[k][j][i] = (B2ch*M3h - B3ch*M2h)/dhalf[k][j][i];
02257 emf2_cc[k][j][i] = (B3ch*M1h - B1ch*M3h)/dhalf[k][j][i];
02258 emf3_cc[k][j][i] = (B1ch*M2h - B2ch*M1h)/dhalf[k][j][i];
02259 #ifndef BAROTROPIC
02260 phalf[k][j][i] -= 0.5*(B1ch*B1ch + B2ch*B2ch + B3ch*B3ch);
02261 #endif
02262 #endif
02263
02264 #ifndef BAROTROPIC
02265 phalf[k][j][i] *= Gamma_1;
02266 #endif
02267
02268 #ifdef PARTICLES
02269 d1 = 1.0/dhalf[k][j][i];
02270 pG->Coup[k][j][i].grid_v1 = M1h*d1;
02271 pG->Coup[k][j][i].grid_v2 = M2h*d1;
02272 pG->Coup[k][j][i].grid_v3 = M3h*d1;
02273 #ifndef BAROTROPIC
02274 pG->Coup[k][j][i].grid_cs = sqrt(Gamma*phalf[k][j][i]*d1);
02275 #endif
02276 #endif
02277
02278 }
02279 }
02280 }
02281 }
02282
02283
02284
02285 #ifdef PARTICLES
02286 Integrate_Particles(pG,pD);
02287 #ifdef FEEDBACK
02288 exchange_feedback(pG, pD);
02289 #endif
02290 #endif
02291
02292
02293
02294
02295
02296
02297
02298
02299 #ifdef H_CORRECTION
02300 for (k=ks-1; k<=ke+1; k++) {
02301 for (j=js-1; j<=je+1; j++) {
02302 for (i=is-1; i<=ie+2; i++) {
02303 #ifdef MHD
02304 Bx = B1_x1Face[k][j][i];
02305 #endif
02306 cfr = cfast(&(Ur_x1Face[k][j][i]),&Bx);
02307 cfl = cfast(&(Ul_x1Face[k][j][i]),&Bx);
02308 lambdar = Ur_x1Face[k][j][i].Mx/Ur_x1Face[k][j][i].d + cfr;
02309 lambdal = Ul_x1Face[k][j][i].Mx/Ul_x1Face[k][j][i].d - cfl;
02310 eta1[k][j][i] = 0.5*fabs(lambdar - lambdal);
02311 }
02312 }
02313 }
02314
02315 for (k=ks-1; k<=ke+1; k++) {
02316 for (j=js-1; j<=je+2; j++) {
02317 for (i=is-1; i<=ie+1; i++) {
02318 #ifdef MHD
02319 Bx = B2_x2Face[k][j][i];
02320 #endif
02321 cfr = cfast(&(Ur_x2Face[k][j][i]),&Bx);
02322 cfl = cfast(&(Ul_x2Face[k][j][i]),&Bx);
02323 lambdar = Ur_x2Face[k][j][i].Mx/Ur_x2Face[k][j][i].d + cfr;
02324 lambdal = Ul_x2Face[k][j][i].Mx/Ul_x2Face[k][j][i].d - cfl;
02325 eta2[k][j][i] = 0.5*fabs(lambdar - lambdal);
02326 }
02327 }
02328 }
02329
02330 for (k=ks-1; k<=ke+2; k++) {
02331 for (j=js-1; j<=je+1; j++) {
02332 for (i=is-1; i<=ie+1; i++) {
02333 #ifdef MHD
02334 Bx = B3_x3Face[k][j][i];
02335 #endif
02336 cfr = cfast(&(Ur_x3Face[k][j][i]),&Bx);
02337 cfl = cfast(&(Ul_x3Face[k][j][i]),&Bx);
02338 lambdar = Ur_x3Face[k][j][i].Mx/Ur_x3Face[k][j][i].d + cfr;
02339 lambdal = Ul_x3Face[k][j][i].Mx/Ul_x3Face[k][j][i].d - cfl;
02340 eta3[k][j][i] = 0.5*fabs(lambdar - lambdal);
02341 }
02342 }
02343 }
02344 #endif
02345
02346
02347
02348
02349
02350 for (k=ks-1; k<=ke+1; k++) {
02351 for (j=js-1; j<=je+1; j++) {
02352 for (i=is; i<=ie+1; i++) {
02353 #ifdef H_CORRECTION
02354 etah = MAX(eta2[k][j][i-1],eta2[k][j][i]);
02355 etah = MAX(etah,eta2[k][j+1][i-1]);
02356 etah = MAX(etah,eta2[k][j+1][i ]);
02357
02358 etah = MAX(etah,eta3[k ][j][i-1]);
02359 etah = MAX(etah,eta3[k ][j][i ]);
02360 etah = MAX(etah,eta3[k+1][j][i-1]);
02361 etah = MAX(etah,eta3[k+1][j][i ]);
02362
02363 etah = MAX(etah,eta1[k ][j][i ]);
02364 #endif
02365 #ifdef MHD
02366 Bx = B1_x1Face[k][j][i];
02367 #endif
02368 Wl[i] = Cons1D_to_Prim1D(&Ul_x1Face[k][j][i],&Bx);
02369 Wr[i] = Cons1D_to_Prim1D(&Ur_x1Face[k][j][i],&Bx);
02370
02371 fluxes(Ul_x1Face[k][j][i],Ur_x1Face[k][j][i],Wl[i],Wr[i],Bx,
02372 &x1Flux[k][j][i]);
02373 }
02374 }
02375 }
02376
02377
02378
02379
02380
02381 for (k=ks-1; k<=ke+1; k++) {
02382 for (j=js; j<=je+1; j++) {
02383 for (i=is-1; i<=ie+1; i++) {
02384 #ifdef H_CORRECTION
02385 etah = MAX(eta1[k][j-1][i],eta1[k][j][i]);
02386 etah = MAX(etah,eta1[k][j-1][i+1]);
02387 etah = MAX(etah,eta1[k][j ][i+1]);
02388
02389 etah = MAX(etah,eta3[k ][j-1][i]);
02390 etah = MAX(etah,eta3[k ][j ][i]);
02391 etah = MAX(etah,eta3[k+1][j-1][i]);
02392 etah = MAX(etah,eta3[k+1][j ][i]);
02393
02394 etah = MAX(etah,eta2[k ][j ][i]);
02395 #endif
02396 #ifdef MHD
02397 Bx = B2_x2Face[k][j][i];
02398 #endif
02399 Wl[i] = Cons1D_to_Prim1D(&Ul_x2Face[k][j][i],&Bx);
02400 Wr[i] = Cons1D_to_Prim1D(&Ur_x2Face[k][j][i],&Bx);
02401
02402 fluxes(Ul_x2Face[k][j][i],Ur_x2Face[k][j][i],Wl[i],Wr[i],Bx,
02403 &x2Flux[k][j][i]);
02404 }
02405 }
02406 }
02407
02408
02409
02410
02411
02412 for (k=ks; k<=ke+1; k++) {
02413 for (j=js-1; j<=je+1; j++) {
02414 for (i=is-1; i<=ie+1; i++) {
02415 #ifdef H_CORRECTION
02416 etah = MAX(eta1[k-1][j][i],eta1[k][j][i]);
02417 etah = MAX(etah,eta1[k-1][j][i+1]);
02418 etah = MAX(etah,eta1[k][j ][i+1]);
02419
02420 etah = MAX(etah,eta2[k-1][j ][i]);
02421 etah = MAX(etah,eta2[k ][j ][i]);
02422 etah = MAX(etah,eta2[k-1][j+1][i]);
02423 etah = MAX(etah,eta2[k ][j+1][i]);
02424
02425 etah = MAX(etah,eta3[k ][j ][i]);
02426 #endif
02427 #ifdef MHD
02428 Bx = B3_x3Face[k][j][i];
02429 #endif
02430 Wl[i] = Cons1D_to_Prim1D(&Ul_x3Face[k][j][i],&Bx);
02431 Wr[i] = Cons1D_to_Prim1D(&Ur_x3Face[k][j][i],&Bx);
02432
02433 fluxes(Ul_x3Face[k][j][i],Ur_x3Face[k][j][i],Wl[i],Wr[i],Bx,
02434 &x3Flux[k][j][i]);
02435 }
02436 }
02437 }
02438
02439
02440
02441
02442
02443
02444
02445 #ifdef MHD
02446 integrate_emf1_corner(pG);
02447 integrate_emf2_corner(pG);
02448 integrate_emf3_corner(pG);
02449
02450
02451 #ifdef SHEARING_BOX
02452 get_myGridIndex(pD, myID_Comm_world, &my_iproc, &my_jproc, &my_kproc);
02453
02454
02455
02456 if (my_iproc == 0) {
02457 RemapEy_ix1(pD, emf2, remapEyiib);
02458 }
02459 if (my_iproc == (pD->NGrid[0]-1)) {
02460 RemapEy_ox1(pD, emf2, remapEyoib);
02461 }
02462
02463
02464
02465 if (my_iproc == 0) {
02466 for(k=ks; k<=ke+1; k++) {
02467 for(j=js; j<=je; j++){
02468 emf2[k][j][is] = 0.5*(emf2[k][j][is] + remapEyiib[k][j]);
02469 }
02470 }
02471 }
02472
02473 if (my_iproc == (pD->NGrid[0]-1)) {
02474 for(k=ks; k<=ke+1; k++) {
02475 for(j=js; j<=je; j++){
02476 emf2[k][j][ie+1] = 0.5*(emf2[k][j][ie+1] + remapEyoib[k][j]);
02477 }
02478 }
02479 }
02480 #endif
02481
02482
02483
02484
02485
02486 for (k=ks; k<=ke; k++) {
02487 for (j=js; j<=je; j++) {
02488 for (i=is; i<=ie; i++) {
02489 #ifdef CYLINDRICAL
02490 dtodx2 = pG->dt/(ri[i]*pG->dx2);
02491 #endif
02492 pG->B1i[k][j][i] += dtodx3*(emf2[k+1][j ][i ] - emf2[k][j][i]) -
02493 dtodx2*(emf3[k ][j+1][i ] - emf3[k][j][i]);
02494 pG->B2i[k][j][i] += dtodx1*(emf3[k ][j ][i+1] - emf3[k][j][i]) -
02495 dtodx3*(emf1[k+1][j ][i ] - emf1[k][j][i]);
02496 #ifdef CYLINDRICAL
02497 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
02498 dtodx2 = pG->dt/(r[i]*pG->dx2);
02499 #endif
02500 pG->B3i[k][j][i] += dtodx2*( emf1[k ][j+1][i ] - emf1[k][j][i]) -
02501 dtodx1*(rsf*emf2[k ][j ][i+1] - lsf*emf2[k][j][i]);
02502 }
02503 #ifdef CYLINDRICAL
02504 dtodx2 = pG->dt/(ri[ie+1]*pG->dx2);
02505 #endif
02506 pG->B1i[k][j][ie+1] +=
02507 dtodx3*(emf2[k+1][j ][ie+1] - emf2[k][j][ie+1]) -
02508 dtodx2*(emf3[k ][j+1][ie+1] - emf3[k][j][ie+1]);
02509 }
02510 for (i=is; i<=ie; i++) {
02511 pG->B2i[k][je+1][i] +=
02512 dtodx1*(emf3[k ][je+1][i+1] - emf3[k][je+1][i]) -
02513 dtodx3*(emf1[k+1][je+1][i ] - emf1[k][je+1][i]);
02514 }
02515 }
02516 for (j=js; j<=je; j++) {
02517 for (i=is; i<=ie; i++) {
02518 #ifdef CYLINDRICAL
02519 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
02520 dtodx2 = pG->dt/(r[i]*pG->dx2);
02521 #endif
02522 pG->B3i[ke+1][j][i] +=
02523 dtodx2*( emf1[ke+1][j+1][i ] - emf1[ke+1][j][i]) -
02524 dtodx1*(rsf*emf2[ke+1][j ][i+1] - lsf*emf2[ke+1][j][i]);
02525 }
02526 }
02527 #endif
02528
02529
02530
02531
02532
02533
02534 #if defined(CYLINDRICAL) && !defined(FARGO)
02535 for (k=ks; k<=ke; k++) {
02536 for (j=js; j<=je; j++) {
02537 for (i=is; i<=ie; i++) {
02538 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
02539 q2 = hdt/(r[i]*pG->dx2);
02540
02541
02542 dhalf[k][j][i] = pG->U[k][j][i].d
02543 - q1*(rsf*x1Flux[k ][j ][i+1].d - lsf*x1Flux[k][j][i].d)
02544 - q2*( x2Flux[k ][j+1][i ].d - x2Flux[k][j][i].d)
02545 - q3*( x3Flux[k+1][j ][i ].d - x3Flux[k][j][i].d);
02546
02547
02548 M2h = pG->U[k][j][i].M2
02549 - q1*(SQR(rsf)*x1Flux[k ][j ][i+1].My - SQR(lsf)*x1Flux[k][j][i].My)
02550 - q2*( x2Flux[k ][j+1][i ].Mx - x2Flux[k][j][i].Mx)
02551 - q3*( x3Flux[k+1][j ][i ].Mz - x3Flux[k][j][i].Mz);
02552
02553
02554 if (StaticGravPot != NULL){
02555 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
02556 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
02557 M2h -= q2*(phir-phil)*pG->U[k][j][i].d;
02558 }
02559
02560
02561 geom_src[k][j][i] = SQR(M2h)/dhalf[k][j][i];
02562 #ifdef MHD
02563 B2ch = 0.5*(B2_x2Face[k][j][i] + B2_x2Face[k][j+1][i]);
02564 geom_src[k][j][i] -= SQR(B2ch);
02565 #endif
02566 #ifdef ISOTHERMAL
02567 geom_src[k][j][i] += Iso_csound2*dhalf[k][j][i];
02568 #ifdef MHD
02569 B1ch = 0.5*(lsf*B1_x1Face[k][j][i] + rsf*B1_x1Face[k ][j ][i+1]);
02570 B3ch = 0.5*( B3_x3Face[k][j][i] + B3_x3Face[k+1][j ][i ]);
02571 geom_src[k][j][i] += 0.5*(SQR(B1ch)+SQR(B2ch)+SQR(B3ch));
02572 #endif
02573 #else
02574 Pavgh = 0.5*(lsf*x1Flux[k][j][i ].Pflux + rsf*x1Flux[k][j][i+1].Pflux);
02575 geom_src[k][j][i] += Pavgh;
02576 #endif
02577 geom_src[k][j][i] /= x1vc(pG,i);
02578
02579
02580 pG->U[k][j][i].M1 += pG->dt*geom_src[k][j][i];
02581 }
02582 }
02583 }
02584 #endif
02585
02586
02587 #if defined(CYLINDRICAL) && defined(FARGO)
02588 for (k=ks; k<=ke; k++) {
02589 for (j=js; j<=je; j++) {
02590 for (i=is; i<=ie; i++) {
02591 cc_pos(pG,i,j,k,&x1,&x2,&x3);
02592 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
02593 q2 = hdt/(r[i]*pG->dx2);
02594
02595
02596 dhalf[k][j][i] = pG->U[k][j][i].d
02597 - q1*(rsf*x1Flux[k ][j ][i+1].d - lsf*x1Flux[k][j][i].d)
02598 - q2*( x2Flux[k ][j+1][i ].d - x2Flux[k][j][i].d)
02599 - q3*( x3Flux[k+1][j ][i ].d - x3Flux[k][j][i].d);
02600
02601 Mrn = pG->U[k][j][i].M1;
02602 Mpn = pG->U[k][j][i].M2;
02603
02604 Om = (*OrbitalProfile)(x1vc(pG,i));
02605 qshear = (*ShearProfile)(x1vc(pG,i));
02606 g = (*x1GravAcc)(x1vc(pG,i),x2,x3) - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
02607
02608 Mre = Mrn
02609 - (dtodx1/r[i])*(ri[i+1]*x1Flux[k][j][i+1].Mx - ri[i]*x1Flux[k][j][i].Mx)
02610 - (dtodx2/r[i])*( x2Flux[k][j+1][i].Mz - x2Flux[k][j][i].Mz)
02611 - (dtodx3)*( x3Flux[k+1][j][i].My - x3Flux[k][j][i].My);
02612 Mre += pG->dt*( 2.0*Om*Mpn + geom_src[k][j][i] - pG->U[k][j][i].d*g);
02613
02614 Mpe = Mpn + pG->dt*Om*(qshear-2.0)*Mrn
02615 - (dtodx1/SQR(r[i]))*(SQR(ri[i+1])*x1Flux[k ][j ][i+1].My - SQR(ri[i])*x1Flux[k][j][i].My)
02616 - (dtodx2/r[i])*( x2Flux[k ][j+1][i ].Mx - x2Flux[k][j][i].Mx)
02617 - (dtodx3)*( x3Flux[k+1][j ][i ].Mz - x3Flux[k][j][i].Mz);
02618
02619 Mrav = 0.5*(Mrn+Mre);
02620 Mpav = 0.5*(Mpn+Mpe);
02621
02622
02623 geom_src[k][j][i] = SQR(Mpav)/dhalf[k][j][i];
02624 #ifdef MHD
02625 B1ch = 0.5*(ri[i]*B1_x1Face[k][j][i] + ri[i+1]*B1_x1Face[k ][j ][i+1])/r[i];
02626 B2ch = 0.5*(B2_x2Face[k][j][i] + B2_x2Face[k][j+1][i]);
02627 B3ch = 0.5*(B3_x3Face[k][j][i] + B3_x3Face[k+1][j ][i ]);
02628 geom_src[k][j][i] += 0.5*( SQR(B1ch) - SQR(B2ch) + SQR(B3ch));
02629 #endif
02630 #ifdef ISOTHERMAL
02631 geom_src[k][j][i] += Iso_csound2*dhalf[k][j][i];
02632 #endif
02633 geom_src[k][j][i] /= x1vc(pG,i);
02634
02635
02636 pG->U[k][j][i].M1 += pG->dt*( 2.0*Om*Mpav + geom_src[k][j][i]);
02637 pG->U[k][j][i].M2 += pG->dt*( Om*(qshear-2.0)*Mrav);
02638
02639 }
02640 }
02641 }
02642
02643 #endif
02644
02645
02646
02647
02648
02649
02650
02651
02652
02653 #ifdef SHEARING_BOX
02654 fact = om_dt/(2. + (2.-qshear)*om_dt*om_dt);
02655 qom = qshear*Omega_0;
02656 for(k=ks; k<=ke; k++) {
02657 for(j=js; j<=je; j++) {
02658 for(i=is; i<=ie; i++) {
02659 cc_pos(pG,i,j,k,&x1,&x2,&x3);
02660
02661
02662 M1n = pG->U[k][j][i].M1;
02663 #ifdef FARGO
02664 dM2n = pG->U[k][j][i].M2;
02665 #else
02666 dM2n = pG->U[k][j][i].M2 + qom*x1*pG->U[k][j][i].d;
02667 #endif
02668
02669
02670 frx1_dM2 = x1Flux[k][j][i+1].My;
02671 flx1_dM2 = x1Flux[k][j][i ].My;
02672 frx2_dM2 = x2Flux[k][j+1][i].Mx;
02673 flx2_dM2 = x2Flux[k][j ][i].Mx;
02674 frx3_dM2 = x3Flux[k+1][j][i].Mz;
02675 flx3_dM2 = x3Flux[k ][j][i].Mz;
02676 #ifndef FARGO
02677 frx1_dM2 += qom*(x1+0.5*pG->dx1)*x1Flux[k][j][i+1].d;
02678 flx1_dM2 += qom*(x1-0.5*pG->dx1)*x1Flux[k][j][i ].d;
02679 frx2_dM2 += qom*(x1 )*x2Flux[k][j+1][i].d;
02680 flx2_dM2 += qom*(x1 )*x2Flux[k][j ][i].d;
02681 frx3_dM2 += qom*(x1 )*x3Flux[k+1][j][i].d;
02682 flx3_dM2 += qom*(x1 )*x3Flux[k ][j][i].d;
02683 #endif
02684
02685
02686 M1e = M1n - q1*(x1Flux[k][j][i+1].Mx - x1Flux[k][j][i].Mx)
02687 - q2*(x2Flux[k][j+1][i].Mz - x2Flux[k][j][i].Mz)
02688 - q3*(x3Flux[k+1][j][i].My - x3Flux[k][j][i].My);
02689
02690 dM2e = dM2n - q1*(frx1_dM2 - flx1_dM2)
02691 - q2*(frx2_dM2 - flx2_dM2)
02692 - q3*(frx3_dM2 - flx3_dM2);
02693
02694 #ifdef FEEDBACK
02695 M1e -= 0.5*pG->Coup[k][j][i].fb1;
02696 dM2e -= 0.5*pG->Coup[k][j][i].fb2;
02697 #endif
02698
02699
02700
02701
02702
02703 pG->U[k][j][i].M1 += (4.0*dM2e + 2.0*(qshear-2.)*om_dt*M1e)*fact;
02704 pG->U[k][j][i].M2 += 2.0*(qshear-2.)*(M1e + om_dt*dM2e)*fact;
02705
02706 #ifndef FARGO
02707 pG->U[k][j][i].M2 -= 0.5*qshear*om_dt*
02708 (x1Flux[k][j][i].d + x1Flux[k][j][i+1].d);
02709 #endif
02710
02711
02712
02713
02714 phic = (*ShearingBoxPot)((x1 ),x2,x3);
02715 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
02716 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
02717 #ifndef BAROTROPIC
02718 pG->U[k][j][i].E -= dtodx1*(x1Flux[k][j][i ].d*(phic - phil) +
02719 x1Flux[k][j][i+1].d*(phir - phic));
02720 #endif
02721
02722 phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
02723 phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
02724 #ifndef BAROTROPIC
02725 pG->U[k][j][i].E -= dtodx2*(x2Flux[k][j ][i].d*(phic - phil) +
02726 x2Flux[k][j+1][i].d*(phir - phic));
02727 #endif
02728
02729 phir = (*ShearingBoxPot)(x1,x2,(x3+0.5*pG->dx3));
02730 phil = (*ShearingBoxPot)(x1,x2,(x3-0.5*pG->dx3));
02731 #ifndef BAROTROPIC
02732 pG->U[k][j][i].E -= dtodx3*(x3Flux[k ][j][i].d*(phic - phil) +
02733 x3Flux[k+1][j][i].d*(phir - phic));
02734 #endif
02735 }
02736 }
02737 }
02738
02739 #endif
02740
02741 if (StaticGravPot != NULL){
02742 for (k=ks; k<=ke; k++) {
02743 for (j=js; j<=je; j++) {
02744 for (i=is; i<=ie; i++) {
02745 cc_pos(pG,i,j,k,&x1,&x2,&x3);
02746 phic = (*StaticGravPot)((x1 ),x2,x3);
02747 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
02748 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
02749
02750 #ifdef CYLINDRICAL
02751 g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
02752 #ifdef FARGO
02753 g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
02754 #endif
02755 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
02756 dtodx2 = pG->dt/(r[i]*pG->dx2);
02757 pG->U[k][j][i].M1 -= pG->dt*dhalf[k][j][i]*g;
02758 #else
02759 pG->U[k][j][i].M1 -= dtodx1*(phir-phil)*dhalf[k][j][i];
02760 #endif
02761 #ifndef BAROTROPIC
02762 pG->U[k][j][i].E -= dtodx1*(lsf*x1Flux[k][j][i ].d*(phic - phil) +
02763 rsf*x1Flux[k][j][i+1].d*(phir - phic));
02764 #endif
02765 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
02766 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
02767 pG->U[k][j][i].M2 -= dtodx2*(phir-phil)*dhalf[k][j][i];
02768 #ifndef BAROTROPIC
02769 pG->U[k][j][i].E -= dtodx2*(x2Flux[k][j ][i].d*(phic - phil) +
02770 x2Flux[k][j+1][i].d*(phir - phic));
02771 #endif
02772 phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
02773 phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
02774 pG->U[k][j][i].M3 -= dtodx3*(phir-phil)*dhalf[k][j][i];
02775 #ifndef BAROTROPIC
02776 pG->U[k][j][i].E -= dtodx3*(x3Flux[k ][j][i].d*(phic - phil) +
02777 x3Flux[k+1][j][i].d*(phir - phic));
02778 #endif
02779 }
02780 }
02781 }
02782 }
02783
02784
02785
02786
02787
02788
02789
02790 #ifdef SELF_GRAVITY
02791
02792
02793 for (k=ks; k<=ke; k++){
02794 for (j=js; j<=je; j++){
02795 for (i=is; i<=ie; i++){
02796 phic = pG->Phi[k][j][i];
02797 phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j][i ]);
02798 phir = 0.5*(pG->Phi[k][j][i ] + pG->Phi[k][j][i+1]);
02799
02800
02801 gxl = (pG->Phi[k][j][i-1] - pG->Phi[k][j][i ])*(dx1i);
02802 gxr = (pG->Phi[k][j][i ] - pG->Phi[k][j][i+1])*(dx1i);
02803
02804 gyl = 0.25*((pG->Phi[k][j-1][i-1] - pG->Phi[k][j+1][i-1]) +
02805 (pG->Phi[k][j-1][i ] - pG->Phi[k][j+1][i ]) )*(dx2i);
02806 gyr = 0.25*((pG->Phi[k][j-1][i ] - pG->Phi[k][j+1][i ]) +
02807 (pG->Phi[k][j-1][i+1] - pG->Phi[k][j+1][i+1]) )*(dx2i);
02808
02809 gzl = 0.25*((pG->Phi[k-1][j][i-1] - pG->Phi[k+1][j][i-1]) +
02810 (pG->Phi[k-1][j][i ] - pG->Phi[k+1][j][i ]) )*(dx3i);
02811 gzr = 0.25*((pG->Phi[k-1][j][i ] - pG->Phi[k+1][j][i ]) +
02812 (pG->Phi[k-1][j][i+1] - pG->Phi[k+1][j][i+1]) )*(dx3i);
02813
02814
02815 flx_m1l = 0.5*(gxl*gxl-gyl*gyl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
02816 flx_m1r = 0.5*(gxr*gxr-gyr*gyr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
02817
02818 flx_m2l = gxl*gyl/four_pi_G;
02819 flx_m2r = gxr*gyr/four_pi_G;
02820
02821 flx_m3l = gxl*gzl/four_pi_G;
02822 flx_m3r = gxr*gzr/four_pi_G;
02823
02824
02825 pG->U[k][j][i].M1 -= dtodx1*(flx_m1r - flx_m1l);
02826 pG->U[k][j][i].M2 -= dtodx1*(flx_m2r - flx_m2l);
02827 pG->U[k][j][i].M3 -= dtodx1*(flx_m3r - flx_m3l);
02828 #ifndef BAROTROPIC
02829 pG->U[k][j][i].E -= dtodx1*(x1Flux[k][j][i ].d*(phic - phil) +
02830 x1Flux[k][j][i+1].d*(phir - phic));
02831 #endif
02832 }
02833 }
02834 }
02835
02836
02837
02838 for (k=ks; k<=ke; k++){
02839 for (j=js; j<=je; j++){
02840 for (i=is; i<=ie; i++){
02841 phic = pG->Phi[k][j][i];
02842 phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j ][i]);
02843 phir = 0.5*(pG->Phi[k][j ][i] + pG->Phi[k][j+1][i]);
02844
02845
02846 gxl = 0.25*((pG->Phi[k][j-1][i-1] - pG->Phi[k][j-1][i+1]) +
02847 (pG->Phi[k][j ][i-1] - pG->Phi[k][j ][i+1]) )*(dx1i);
02848 gxr = 0.25*((pG->Phi[k][j ][i-1] - pG->Phi[k][j ][i+1]) +
02849 (pG->Phi[k][j+1][i-1] - pG->Phi[k][j+1][i+1]) )*(dx1i);
02850
02851 gyl = (pG->Phi[k][j-1][i] - pG->Phi[k][j ][i])*(dx2i);
02852 gyr = (pG->Phi[k][j ][i] - pG->Phi[k][j+1][i])*(dx2i);
02853
02854 gzl = 0.25*((pG->Phi[k-1][j-1][i] - pG->Phi[k+1][j-1][i]) +
02855 (pG->Phi[k-1][j ][i] - pG->Phi[k+1][j ][i]) )*(dx3i);
02856 gzr = 0.25*((pG->Phi[k-1][j ][i] - pG->Phi[k+1][j ][i]) +
02857 (pG->Phi[k-1][j+1][i] - pG->Phi[k+1][j+1][i]) )*(dx3i);
02858
02859
02860 flx_m1l = gyl*gxl/four_pi_G;
02861 flx_m1r = gyr*gxr/four_pi_G;
02862
02863 flx_m2l = 0.5*(gyl*gyl-gxl*gxl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
02864 flx_m2r = 0.5*(gyr*gyr-gxr*gxr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
02865
02866 flx_m3l = gyl*gzl/four_pi_G;
02867 flx_m3r = gyr*gzr/four_pi_G;
02868
02869
02870 pG->U[k][j][i].M1 -= dtodx2*(flx_m1r - flx_m1l);
02871 pG->U[k][j][i].M2 -= dtodx2*(flx_m2r - flx_m2l);
02872 pG->U[k][j][i].M3 -= dtodx2*(flx_m3r - flx_m3l);
02873 #ifndef BAROTROPIC
02874 pG->U[k][j][i].E -= dtodx2*(x2Flux[k][j ][i].d*(phic - phil) +
02875 x2Flux[k][j+1][i].d*(phir - phic));
02876 #endif
02877 }
02878 }
02879 }
02880
02881
02882
02883 for (k=ks; k<=ke; k++){
02884 for (j=js; j<=je; j++){
02885 for (i=is; i<=ie; i++){
02886 phic = pG->Phi[k][j][i];
02887 phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k ][j][i]);
02888 phir = 0.5*(pG->Phi[k ][j][i] + pG->Phi[k+1][j][i]);
02889
02890
02891 gxl = 0.25*((pG->Phi[k-1][j][i-1] - pG->Phi[k-1][j][i+1]) +
02892 (pG->Phi[k ][j][i-1] - pG->Phi[k ][j][i+1]) )*(dx1i);
02893 gxr = 0.25*((pG->Phi[k ][j][i-1] - pG->Phi[k ][j][i+1]) +
02894 (pG->Phi[k+1][j][i-1] - pG->Phi[k+1][j][i+1]) )*(dx1i);
02895
02896 gyl = 0.25*((pG->Phi[k-1][j-1][i] - pG->Phi[k-1][j+1][i]) +
02897 (pG->Phi[k ][j-1][i] - pG->Phi[k ][j+1][i]) )*(dx2i);
02898 gyr = 0.25*((pG->Phi[k ][j-1][i] - pG->Phi[k ][j+1][i]) +
02899 (pG->Phi[k+1][j-1][i] - pG->Phi[k+1][j+1][i]) )*(dx2i);
02900
02901 gzl = (pG->Phi[k-1][j][i] - pG->Phi[k ][j][i])*(dx3i);
02902 gzr = (pG->Phi[k ][j][i] - pG->Phi[k+1][j][i])*(dx3i);
02903
02904
02905 flx_m1l = gzl*gxl/four_pi_G;
02906 flx_m1r = gzr*gxr/four_pi_G;
02907
02908 flx_m2l = gzl*gyl/four_pi_G;
02909 flx_m2r = gzr*gyr/four_pi_G;
02910
02911 flx_m3l = 0.5*(gzl*gzl-gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil;
02912 flx_m3r = 0.5*(gzr*gzr-gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir;
02913
02914
02915 pG->U[k][j][i].M1 -= dtodx3*(flx_m1r - flx_m1l);
02916 pG->U[k][j][i].M2 -= dtodx3*(flx_m2r - flx_m2l);
02917 pG->U[k][j][i].M3 -= dtodx3*(flx_m3r - flx_m3l);
02918 #ifndef BAROTROPIC
02919 pG->U[k][j][i].E -= dtodx3*(x3Flux[k ][j][i].d*(phic - phil) +
02920 x3Flux[k+1][j][i].d*(phir - phic));
02921 #endif
02922 }
02923 }
02924 }
02925
02926
02927
02928 for (k=ks; k<=ke+1; k++) {
02929 for (j=js; j<=je+1; j++) {
02930 for (i=is; i<=ie+1; i++) {
02931 pG->x1MassFlux[k][j][i] = x1Flux[k][j][i].d;
02932 pG->x2MassFlux[k][j][i] = x2Flux[k][j][i].d;
02933 pG->x3MassFlux[k][j][i] = x3Flux[k][j][i].d;
02934 }
02935 }
02936 }
02937 #endif
02938
02939
02940
02941
02942
02943 #ifndef BAROTROPIC
02944 if (CoolingFunc != NULL){
02945 for (k=ks; k<=ke; k++){
02946 for (j=js; j<=je; j++){
02947 for (i=is; i<=ie; i++){
02948 coolf = (*CoolingFunc)(dhalf[k][j][i],phalf[k][j][i],pG->dt);
02949 pG->U[k][j][i].E -= pG->dt*coolf;
02950 }
02951 }
02952 }
02953 }
02954 #endif
02955
02956
02957
02958
02959
02960 #ifdef FEEDBACK
02961 for (k=ks; k<=ke; k++)
02962 for (j=js; j<=je; j++)
02963 for (i=is; i<=ie; i++) {
02964 pG->U[k][j][i].M1 -= pG->Coup[k][j][i].fb1;
02965 pG->U[k][j][i].M2 -= pG->Coup[k][j][i].fb2;
02966 pG->U[k][j][i].M3 -= pG->Coup[k][j][i].fb3;
02967 #ifndef BAROTROPIC
02968 pG->U[k][j][i].E += pG->Coup[k][j][i].Eloss;
02969 pG->Coup[k][j][i].Eloss *= dt1;
02970 pG->Eloss[k][j][i] *= dt1;
02971 #endif
02972 }
02973 #endif
02974
02975
02976
02977
02978
02979
02980
02981 for (k=ks; k<=ke; k++) {
02982 for (j=js; j<=je; j++) {
02983 for (i=is; i<=ie; i++) {
02984 #ifdef CYLINDRICAL
02985 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
02986 #endif
02987 pG->U[k][j][i].d -= dtodx1*(rsf*x1Flux[k][j][i+1].d -lsf*x1Flux[k][j][i].d );
02988 pG->U[k][j][i].M1 -= dtodx1*(rsf*x1Flux[k][j][i+1].Mx-lsf*x1Flux[k][j][i].Mx);
02989 pG->U[k][j][i].M2 -= dtodx1*(SQR(rsf)*x1Flux[k][j][i+1].My-SQR(lsf)*x1Flux[k][j][i].My);
02990 pG->U[k][j][i].M3 -= dtodx1*(rsf*x1Flux[k][j][i+1].Mz-lsf*x1Flux[k][j][i].Mz);
02991 #ifndef BAROTROPIC
02992 pG->U[k][j][i].E -= dtodx1*(rsf*x1Flux[k][j][i+1].E -lsf*x1Flux[k][j][i].E );
02993 #endif
02994 #if (NSCALARS > 0)
02995 for (n=0; n<NSCALARS; n++)
02996 pG->U[k][j][i].s[n] -= dtodx1*(rsf*x1Flux[k][j][i+1].s[n]
02997 - lsf*x1Flux[k][j][i ].s[n]);
02998 #endif
02999 }
03000 }
03001 }
03002
03003
03004
03005
03006
03007 for (k=ks; k<=ke; k++) {
03008 for (j=js; j<=je; j++) {
03009 for (i=is; i<=ie; i++) {
03010 #ifdef CYLINDRICAL
03011 dtodx2 = pG->dt/(r[i]*pG->dx2);
03012 #endif
03013 pG->U[k][j][i].d -= dtodx2*(x2Flux[k][j+1][i].d -x2Flux[k][j][i].d );
03014 pG->U[k][j][i].M1 -= dtodx2*(x2Flux[k][j+1][i].Mz-x2Flux[k][j][i].Mz);
03015 pG->U[k][j][i].M2 -= dtodx2*(x2Flux[k][j+1][i].Mx-x2Flux[k][j][i].Mx);
03016 pG->U[k][j][i].M3 -= dtodx2*(x2Flux[k][j+1][i].My-x2Flux[k][j][i].My);
03017 #ifndef BAROTROPIC
03018 pG->U[k][j][i].E -=dtodx2*(x2Flux[k][j+1][i].E -x2Flux[k][j][i].E );
03019 #endif
03020 #if (NSCALARS > 0)
03021 for (n=0; n<NSCALARS; n++)
03022 pG->U[k][j][i].s[n] -= dtodx2*(x2Flux[k][j+1][i].s[n]
03023 - x2Flux[k][j ][i].s[n]);
03024 #endif
03025 }
03026 }
03027 }
03028
03029
03030
03031
03032
03033 for (k=ks; k<=ke; k++) {
03034 for (j=js; j<=je; j++) {
03035 for (i=is; i<=ie; i++) {
03036 pG->U[k][j][i].d -= dtodx3*(x3Flux[k+1][j][i].d -x3Flux[k][j][i].d );
03037 pG->U[k][j][i].M1 -= dtodx3*(x3Flux[k+1][j][i].My-x3Flux[k][j][i].My);
03038 pG->U[k][j][i].M2 -= dtodx3*(x3Flux[k+1][j][i].Mz-x3Flux[k][j][i].Mz);
03039 pG->U[k][j][i].M3 -= dtodx3*(x3Flux[k+1][j][i].Mx-x3Flux[k][j][i].Mx);
03040 #ifndef BAROTROPIC
03041 pG->U[k][j][i].E -= dtodx3*(x3Flux[k+1][j][i].E -x3Flux[k][j][i].E );
03042 #endif
03043 #if (NSCALARS > 0)
03044 for (n=0; n<NSCALARS; n++)
03045 pG->U[k][j][i].s[n] -= dtodx3*(x3Flux[k+1][j][i].s[n]
03046 - x3Flux[k ][j][i].s[n]);
03047 #endif
03048 }
03049 }
03050 }
03051
03052
03053
03054
03055
03056 #ifdef MHD
03057 for (k=ks; k<=ke; k++) {
03058 for (j=js; j<=je; j++) {
03059 for (i=is; i<=ie; i++) {
03060 #ifdef CYLINDRICAL
03061 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
03062 #endif
03063 pG->U[k][j][i].B1c = 0.5*(lsf*pG->B1i[k][j][i] + rsf*pG->B1i[k][j][i+1]);
03064 pG->U[k][j][i].B2c = 0.5*( pG->B2i[k][j][i] + pG->B2i[k][j+1][i]);
03065 pG->U[k][j][i].B3c = 0.5*( pG->B3i[k][j][i] + pG->B3i[k+1][j][i]);
03066 }
03067 }
03068 }
03069 #endif
03070
03071 #ifdef STATIC_MESH_REFINEMENT
03072
03073
03074
03075
03076 for (ncg=0; ncg<pG->NCGrid; ncg++) {
03077
03078
03079
03080 for (dim=0; dim<2; dim++){
03081 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
03082
03083 if (dim==0) i = pG->CGrid[ncg].ijks[0];
03084 if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
03085 jcs = pG->CGrid[ncg].ijks[1];
03086 jce = pG->CGrid[ncg].ijke[1];
03087 kcs = pG->CGrid[ncg].ijks[2];
03088 kce = pG->CGrid[ncg].ijke[2];
03089
03090 for (k=kcs, kk=0; k<=kce; k++, kk++){
03091 for (j=jcs, jj=0; j<=jce; j++, jj++){
03092 pG->CGrid[ncg].myFlx[dim][kk][jj].d = x1Flux[k][j][i].d;
03093 pG->CGrid[ncg].myFlx[dim][kk][jj].M1 = x1Flux[k][j][i].Mx;
03094 pG->CGrid[ncg].myFlx[dim][kk][jj].M2 = x1Flux[k][j][i].My;
03095 pG->CGrid[ncg].myFlx[dim][kk][jj].M3 = x1Flux[k][j][i].Mz;
03096 #ifndef BAROTROPIC
03097 pG->CGrid[ncg].myFlx[dim][kk][jj].E = x1Flux[k][j][i].E;
03098 #endif
03099 #ifdef MHD
03100 pG->CGrid[ncg].myFlx[dim][kk][jj].B1c = 0.0;
03101 pG->CGrid[ncg].myFlx[dim][kk][jj].B2c = x1Flux[k][j][i].By;
03102 pG->CGrid[ncg].myFlx[dim][kk][jj].B3c = x1Flux[k][j][i].Bz;
03103 #endif
03104 #if (NSCALARS > 0)
03105 for (n=0; n<NSCALARS; n++)
03106 pG->CGrid[ncg].myFlx[dim][kk][jj].s[n] = x1Flux[k][j][i].s[n];
03107 #endif
03108 }
03109 }
03110 #ifdef MHD
03111 for (k=kcs, kk=0; k<=kce+1; k++, kk++){
03112 for (j=jcs, jj=0; j<=jce; j++, jj++){
03113 pG->CGrid[ncg].myEMF2[dim][kk][jj] = emf2[k][j][i];
03114 }
03115 }
03116 for (k=kcs, kk=0; k<=kce; k++, kk++){
03117 for (j=jcs, jj=0; j<=jce+1; j++, jj++){
03118 pG->CGrid[ncg].myEMF3[dim][kk][jj] = emf3[k][j][i];
03119 }
03120 }
03121 #endif
03122 }
03123 }
03124
03125
03126
03127 for (dim=2; dim<4; dim++){
03128 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
03129
03130 ics = pG->CGrid[ncg].ijks[0];
03131 ice = pG->CGrid[ncg].ijke[0];
03132 if (dim==2) j = pG->CGrid[ncg].ijks[1];
03133 if (dim==3) j = pG->CGrid[ncg].ijke[1] + 1;
03134 kcs = pG->CGrid[ncg].ijks[2];
03135 kce = pG->CGrid[ncg].ijke[2];
03136
03137 for (k=kcs, kk=0; k<=kce; k++, kk++){
03138 for (i=ics, ii=0; i<=ice; i++, ii++){
03139 pG->CGrid[ncg].myFlx[dim][kk][ii].d = x2Flux[k][j][i].d;
03140 pG->CGrid[ncg].myFlx[dim][kk][ii].M1 = x2Flux[k][j][i].Mz;
03141 pG->CGrid[ncg].myFlx[dim][kk][ii].M2 = x2Flux[k][j][i].Mx;
03142 pG->CGrid[ncg].myFlx[dim][kk][ii].M3 = x2Flux[k][j][i].My;
03143 #ifndef BAROTROPIC
03144 pG->CGrid[ncg].myFlx[dim][kk][ii].E = x2Flux[k][j][i].E;
03145 #endif
03146 #ifdef MHD
03147 pG->CGrid[ncg].myFlx[dim][kk][ii].B1c = x2Flux[k][j][i].Bz;
03148 pG->CGrid[ncg].myFlx[dim][kk][ii].B2c = 0.0;
03149 pG->CGrid[ncg].myFlx[dim][kk][ii].B3c = x2Flux[k][j][i].By;
03150 #endif
03151 #if (NSCALARS > 0)
03152 for (n=0; n<NSCALARS; n++)
03153 pG->CGrid[ncg].myFlx[dim][kk][ii].s[n] = x2Flux[k][j][i].s[n];
03154 #endif
03155 }
03156 }
03157 #ifdef MHD
03158 for (k=kcs, kk=0; k<=kce+1; k++, kk++){
03159 for (i=ics, ii=0; i<=ice; i++, ii++){
03160 pG->CGrid[ncg].myEMF1[dim][kk][ii] = emf1[k][j][i];
03161 }
03162 }
03163 for (k=kcs, kk=0; k<=kce; k++, kk++){
03164 for (i=ics, ii=0; i<=ice+1; i++, ii++){
03165 pG->CGrid[ncg].myEMF3[dim][kk][ii] = emf3[k][j][i];
03166 }
03167 }
03168 #endif
03169 }
03170 }
03171
03172
03173
03174 for (dim=4; dim<6; dim++){
03175 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
03176
03177 ics = pG->CGrid[ncg].ijks[0];
03178 ice = pG->CGrid[ncg].ijke[0];
03179 jcs = pG->CGrid[ncg].ijks[1];
03180 jce = pG->CGrid[ncg].ijke[1];
03181 if (dim==4) k = pG->CGrid[ncg].ijks[2];
03182 if (dim==5) k = pG->CGrid[ncg].ijke[2] + 1;
03183
03184 for (j=jcs, jj=0; j<=jce; j++, jj++){
03185 for (i=ics, ii=0; i<=ice; i++, ii++){
03186 pG->CGrid[ncg].myFlx[dim][jj][ii].d = x3Flux[k][j][i].d;
03187 pG->CGrid[ncg].myFlx[dim][jj][ii].M1 = x3Flux[k][j][i].My;
03188 pG->CGrid[ncg].myFlx[dim][jj][ii].M2 = x3Flux[k][j][i].Mz;
03189 pG->CGrid[ncg].myFlx[dim][jj][ii].M3 = x3Flux[k][j][i].Mx;
03190 #ifndef BAROTROPIC
03191 pG->CGrid[ncg].myFlx[dim][jj][ii].E = x3Flux[k][j][i].E;
03192 #endif
03193 #ifdef MHD
03194 pG->CGrid[ncg].myFlx[dim][jj][ii].B1c = x3Flux[k][j][i].By;
03195 pG->CGrid[ncg].myFlx[dim][jj][ii].B2c = x3Flux[k][j][i].Bz;
03196 pG->CGrid[ncg].myFlx[dim][jj][ii].B3c = 0.0;
03197 #endif
03198 #if (NSCALARS > 0)
03199 for (n=0; n<NSCALARS; n++)
03200 pG->CGrid[ncg].myFlx[dim][jj][ii].s[n] = x3Flux[k][j][i].s[n];
03201 #endif
03202 }
03203 }
03204 #ifdef MHD
03205 for (j=jcs, jj=0; j<=jce+1; j++, jj++){
03206 for (i=ics, ii=0; i<=ice; i++, ii++){
03207 pG->CGrid[ncg].myEMF1[dim][jj][ii] = emf1[k][j][i];
03208 }
03209 }
03210 for (j=jcs, jj=0; j<=jce; j++, jj++){
03211 for (i=ics, ii=0; i<=ice+1; i++, ii++){
03212 pG->CGrid[ncg].myEMF2[dim][jj][ii] = emf2[k][j][i];
03213 }
03214 }
03215 #endif
03216 }
03217 }
03218 }
03219
03220
03221
03222 for (npg=0; npg<pG->NPGrid; npg++) {
03223
03224
03225
03226 for (dim=0; dim<2; dim++){
03227 if (pG->PGrid[npg].myFlx[dim] != NULL) {
03228
03229 if (dim==0) i = pG->PGrid[npg].ijks[0];
03230 if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
03231 jps = pG->PGrid[npg].ijks[1];
03232 jpe = pG->PGrid[npg].ijke[1];
03233 kps = pG->PGrid[npg].ijks[2];
03234 kpe = pG->PGrid[npg].ijke[2];
03235
03236 for (k=kps, kk=0; k<=kpe; k++, kk++){
03237 for (j=jps, jj=0; j<=jpe; j++, jj++){
03238 pG->PGrid[npg].myFlx[dim][kk][jj].d = x1Flux[k][j][i].d;
03239 pG->PGrid[npg].myFlx[dim][kk][jj].M1 = x1Flux[k][j][i].Mx;
03240 pG->PGrid[npg].myFlx[dim][kk][jj].M2 = x1Flux[k][j][i].My;
03241 pG->PGrid[npg].myFlx[dim][kk][jj].M3 = x1Flux[k][j][i].Mz;
03242 #ifndef BAROTROPIC
03243 pG->PGrid[npg].myFlx[dim][kk][jj].E = x1Flux[k][j][i].E;
03244 #endif
03245 #ifdef MHD
03246 pG->PGrid[npg].myFlx[dim][kk][jj].B1c = 0.0;
03247 pG->PGrid[npg].myFlx[dim][kk][jj].B2c = x1Flux[k][j][i].By;
03248 pG->PGrid[npg].myFlx[dim][kk][jj].B3c = x1Flux[k][j][i].Bz;
03249 #endif
03250 #if (NSCALARS > 0)
03251 for (n=0; n<NSCALARS; n++)
03252 pG->PGrid[npg].myFlx[dim][kk][jj].s[n] = x1Flux[k][j][i].s[n];
03253 #endif
03254 }
03255 }
03256 #ifdef MHD
03257 for (k=kps, kk=0; k<=kpe+1; k++, kk++){
03258 for (j=jps, jj=0; j<=jpe; j++, jj++){
03259 pG->PGrid[npg].myEMF2[dim][kk][jj] = emf2[k][j][i];
03260 }
03261 }
03262 for (k=kps, kk=0; k<=kpe; k++, kk++){
03263 for (j=jps, jj=0; j<=jpe+1; j++, jj++){
03264 pG->PGrid[npg].myEMF3[dim][kk][jj] = emf3[k][j][i];
03265 }
03266 }
03267 #endif
03268 }
03269 }
03270
03271
03272
03273 for (dim=2; dim<4; dim++){
03274 if (pG->PGrid[npg].myFlx[dim] != NULL) {
03275
03276 ips = pG->PGrid[npg].ijks[0];
03277 ipe = pG->PGrid[npg].ijke[0];
03278 if (dim==2) j = pG->PGrid[npg].ijks[1];
03279 if (dim==3) j = pG->PGrid[npg].ijke[1] + 1;
03280 kps = pG->PGrid[npg].ijks[2];
03281 kpe = pG->PGrid[npg].ijke[2];
03282
03283 for (k=kps, kk=0; k<=kpe; k++, kk++){
03284 for (i=ips, ii=0; i<=ipe; i++, ii++){
03285 pG->PGrid[npg].myFlx[dim][kk][ii].d = x2Flux[k][j][i].d;
03286 pG->PGrid[npg].myFlx[dim][kk][ii].M1 = x2Flux[k][j][i].Mz;
03287 pG->PGrid[npg].myFlx[dim][kk][ii].M2 = x2Flux[k][j][i].Mx;
03288 pG->PGrid[npg].myFlx[dim][kk][ii].M3 = x2Flux[k][j][i].My;
03289 #ifndef BAROTROPIC
03290 pG->PGrid[npg].myFlx[dim][kk][ii].E = x2Flux[k][j][i].E;
03291 #endif
03292 #ifdef MHD
03293 pG->PGrid[npg].myFlx[dim][kk][ii].B1c = x2Flux[k][j][i].Bz;
03294 pG->PGrid[npg].myFlx[dim][kk][ii].B2c = 0.0;
03295 pG->PGrid[npg].myFlx[dim][kk][ii].B3c = x2Flux[k][j][i].By;
03296 #endif
03297 #if (NSCALARS > 0)
03298 for (n=0; n<NSCALARS; n++)
03299 pG->PGrid[npg].myFlx[dim][kk][ii].s[n] = x2Flux[k][j][i].s[n];
03300 #endif
03301 }
03302 }
03303 #ifdef MHD
03304 for (k=kps, kk=0; k<=kpe+1; k++, kk++){
03305 for (i=ips, ii=0; i<=ipe; i++, ii++){
03306 pG->PGrid[npg].myEMF1[dim][kk][ii] = emf1[k][j][i];
03307 }
03308 }
03309 for (k=kps, kk=0; k<=kpe; k++, kk++){
03310 for (i=ips, ii=0; i<=ipe+1; i++, ii++){
03311 pG->PGrid[npg].myEMF3[dim][kk][ii] = emf3[k][j][i];
03312 }
03313 }
03314 #endif
03315 }
03316 }
03317
03318
03319
03320 for (dim=4; dim<6; dim++){
03321 if (pG->PGrid[npg].myFlx[dim] != NULL) {
03322
03323 ips = pG->PGrid[npg].ijks[0];
03324 ipe = pG->PGrid[npg].ijke[0];
03325 jps = pG->PGrid[npg].ijks[1];
03326 jpe = pG->PGrid[npg].ijke[1];
03327 if (dim==4) k = pG->PGrid[npg].ijks[2];
03328 if (dim==5) k = pG->PGrid[npg].ijke[2] + 1;
03329
03330 for (j=jps, jj=0; j<=jpe; j++, jj++){
03331 for (i=ips, ii=0; i<=ipe; i++, ii++){
03332 pG->PGrid[npg].myFlx[dim][jj][ii].d = x3Flux[k][j][i].d;
03333 pG->PGrid[npg].myFlx[dim][jj][ii].M1 = x3Flux[k][j][i].My;
03334 pG->PGrid[npg].myFlx[dim][jj][ii].M2 = x3Flux[k][j][i].Mz;
03335 pG->PGrid[npg].myFlx[dim][jj][ii].M3 = x3Flux[k][j][i].Mx;
03336 #ifndef BAROTROPIC
03337 pG->PGrid[npg].myFlx[dim][jj][ii].E = x3Flux[k][j][i].E;
03338 #endif
03339 #ifdef MHD
03340 pG->PGrid[npg].myFlx[dim][jj][ii].B1c = x3Flux[k][j][i].By;
03341 pG->PGrid[npg].myFlx[dim][jj][ii].B2c = x3Flux[k][j][i].Bz;
03342 pG->PGrid[npg].myFlx[dim][jj][ii].B3c = 0.0;
03343 #endif
03344 #if (NSCALARS > 0)
03345 for (n=0; n<NSCALARS; n++)
03346 pG->PGrid[npg].myFlx[dim][jj][ii].s[n] = x3Flux[k][j][i].s[n];
03347 #endif
03348 }
03349 }
03350 #ifdef MHD
03351 for (j=jps, jj=0; j<=jpe+1; j++, jj++){
03352 for (i=ips, ii=0; i<=ipe; i++, ii++){
03353 pG->PGrid[npg].myEMF1[dim][jj][ii] = emf1[k][j][i];
03354 }
03355 }
03356 for (j=jps, jj=0; j<=jpe; j++, jj++){
03357 for (i=ips, ii=0; i<=ipe+1; i++, ii++){
03358 pG->PGrid[npg].myEMF2[dim][jj][ii] = emf2[k][j][i];
03359 }
03360 }
03361 #endif
03362 }
03363 }
03364 }
03365
03366 #endif
03367 return;
03368 }
03369
03370
03371
03372
03373
03374 void integrate_init_3d(MeshS *pM)
03375 {
03376 int nmax,size1=0,size2=0,size3=0,nl,nd;
03377
03378
03379 for (nl=0; nl<(pM->NLevels); nl++){
03380 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
03381 if (pM->Domain[nl][nd].Grid != NULL) {
03382 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
03383 size1 = pM->Domain[nl][nd].Grid->Nx[0];
03384 }
03385 if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
03386 size2 = pM->Domain[nl][nd].Grid->Nx[1];
03387 }
03388 if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
03389 size3 = pM->Domain[nl][nd].Grid->Nx[2];
03390 }
03391 }
03392 }
03393 }
03394
03395 size1 = size1 + 2*nghost;
03396 size2 = size2 + 2*nghost;
03397 size3 = size3 + 2*nghost;
03398 nmax = MAX((MAX(size1,size2)),size3);
03399
03400 #ifdef MHD
03401 if ((emf1 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03402 goto on_error;
03403 if ((emf2 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03404 goto on_error;
03405 if ((emf3 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03406 goto on_error;
03407
03408 if ((emf1_cc=(Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03409 goto on_error;
03410 if ((emf2_cc=(Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03411 goto on_error;
03412 if ((emf3_cc=(Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03413 goto on_error;
03414 #endif
03415 #ifdef H_CORRECTION
03416 if ((eta1 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03417 goto on_error;
03418 if ((eta2 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03419 goto on_error;
03420 if ((eta3 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03421 goto on_error;
03422 #endif
03423
03424 if ((Bxc = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
03425 if ((Bxi = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
03426
03427 #ifdef MHD
03428 if ((B1_x1Face = (Real***)calloc_3d_array(size3,size2,size1, sizeof(Real)))
03429 == NULL) goto on_error;
03430 if ((B2_x2Face = (Real***)calloc_3d_array(size3,size2,size1, sizeof(Real)))
03431 == NULL) goto on_error;
03432 if ((B3_x3Face = (Real***)calloc_3d_array(size3,size2,size1, sizeof(Real)))
03433 == NULL) goto on_error;
03434 #endif
03435
03436 if ((U1d=(Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
03437 if ((W =(Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
03438 if ((Wl =(Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
03439 if ((Wr =(Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
03440
03441 if ((Ul_x1Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03442 == NULL) goto on_error;
03443 if ((Ur_x1Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03444 == NULL) goto on_error;
03445 if ((Ul_x2Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03446 == NULL) goto on_error;
03447 if ((Ur_x2Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03448 == NULL) goto on_error;
03449 if ((Ul_x3Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03450 == NULL) goto on_error;
03451 if ((Ur_x3Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03452 == NULL) goto on_error;
03453 if ((x1Flux =(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03454 == NULL) goto on_error;
03455 if ((x2Flux =(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03456 == NULL) goto on_error;
03457 if ((x3Flux =(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03458 == NULL) goto on_error;
03459
03460 #ifdef CYLINDRICAL
03461 #ifndef MHD
03462 #ifndef PARTICLES
03463 if((StaticGravPot != NULL) || (CoolingFunc != NULL))
03464 #endif
03465 #endif
03466 #endif
03467 {
03468 if ((dhalf = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03469 goto on_error;
03470 if ((phalf = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03471 goto on_error;
03472 }
03473
03474 #ifdef SHEARING_BOX
03475 if ((remapEyiib = (Real**)calloc_2d_array(size3,size2, sizeof(Real))) == NULL)
03476 goto on_error;
03477 if ((remapEyoib = (Real**)calloc_2d_array(size3,size2, sizeof(Real))) == NULL)
03478 goto on_error;
03479 #endif
03480
03481
03482 #ifdef CYLINDRICAL
03483 if ((geom_src = (Real***)calloc_3d_array(size3, size2, size1, sizeof(Real))) == NULL)
03484 goto on_error;
03485 #endif
03486
03487 return;
03488
03489 on_error:
03490 integrate_destruct();
03491 ath_error("[integrate_init]: malloc returned a NULL pointer\n");
03492 }
03493
03494
03495
03496
03497
03498 void integrate_destruct_3d(void)
03499 {
03500
03501 #ifdef MHD
03502 if (emf1 != NULL) free_3d_array(emf1);
03503 if (emf2 != NULL) free_3d_array(emf2);
03504 if (emf3 != NULL) free_3d_array(emf3);
03505 if (emf1_cc != NULL) free_3d_array(emf1_cc);
03506 if (emf2_cc != NULL) free_3d_array(emf2_cc);
03507 if (emf3_cc != NULL) free_3d_array(emf3_cc);
03508 #endif
03509 #ifdef H_CORRECTION
03510 if (eta1 != NULL) free_3d_array(eta1);
03511 if (eta2 != NULL) free_3d_array(eta2);
03512 if (eta3 != NULL) free_3d_array(eta3);
03513 #endif
03514
03515 if (Bxc != NULL) free(Bxc);
03516 if (Bxi != NULL) free(Bxi);
03517 #ifdef MHD
03518 if (B1_x1Face != NULL) free_3d_array(B1_x1Face);
03519 if (B2_x2Face != NULL) free_3d_array(B2_x2Face);
03520 if (B3_x3Face != NULL) free_3d_array(B3_x3Face);
03521 #endif
03522
03523 if (U1d != NULL) free(U1d);
03524 if (W != NULL) free(W);
03525 if (Wl != NULL) free(Wl);
03526 if (Wr != NULL) free(Wr);
03527
03528 if (Ul_x1Face != NULL) free_3d_array(Ul_x1Face);
03529 if (Ur_x1Face != NULL) free_3d_array(Ur_x1Face);
03530 if (Ul_x2Face != NULL) free_3d_array(Ul_x2Face);
03531 if (Ur_x2Face != NULL) free_3d_array(Ur_x2Face);
03532 if (Ul_x3Face != NULL) free_3d_array(Ul_x3Face);
03533 if (Ur_x3Face != NULL) free_3d_array(Ur_x3Face);
03534 if (x1Flux != NULL) free_3d_array(x1Flux);
03535 if (x2Flux != NULL) free_3d_array(x2Flux);
03536 if (x3Flux != NULL) free_3d_array(x3Flux);
03537 if (dhalf != NULL) free_3d_array(dhalf);
03538 if (phalf != NULL) free_3d_array(phalf);
03539 #ifdef SHEARING_BOX
03540 if (remapEyiib != NULL) free_2d_array(remapEyiib);
03541 if (remapEyoib != NULL) free_2d_array(remapEyoib);
03542 #endif
03543
03544
03545 #ifdef CYLINDRICAL
03546 if (geom_src != NULL) free_3d_array(geom_src);
03547 #endif
03548
03549 return;
03550 }
03551
03552
03553
03554
03555
03556
03557
03558
03559
03560
03561
03562
03563
03564
03565
03566 #ifdef MHD
03567 static void integrate_emf1_corner(const GridS *pG)
03568 {
03569 int i, is = pG->is, ie = pG->ie;
03570 int j, js = pG->js, je = pG->je;
03571 int k, ks = pG->ks, ke = pG->ke;
03572 Real de1_l2, de1_r2, de1_l3, de1_r3;
03573
03574 for (k=ks-1; k<=ke+2; k++) {
03575 for (j=js-1; j<=je+2; j++) {
03576 for (i=is-2; i<=ie+2; i++) {
03577
03578
03579 if (x2Flux[k-1][j][i].d > 0.0)
03580 de1_l3 = x3Flux[k][j-1][i].Bz - emf1_cc[k-1][j-1][i];
03581 else if (x2Flux[k-1][j][i].d < 0.0)
03582 de1_l3 = x3Flux[k][j][i].Bz - emf1_cc[k-1][j][i];
03583 else {
03584 de1_l3 = 0.5*(x3Flux[k][j-1][i].Bz - emf1_cc[k-1][j-1][i] +
03585 x3Flux[k][j ][i].Bz - emf1_cc[k-1][j ][i] );
03586 }
03587
03588 if (x2Flux[k][j][i].d > 0.0)
03589 de1_r3 = x3Flux[k][j-1][i].Bz - emf1_cc[k][j-1][i];
03590 else if (x2Flux[k][j][i].d < 0.0)
03591 de1_r3 = x3Flux[k][j][i].Bz - emf1_cc[k][j][i];
03592 else {
03593 de1_r3 = 0.5*(x3Flux[k][j-1][i].Bz - emf1_cc[k][j-1][i] +
03594 x3Flux[k][j ][i].Bz - emf1_cc[k][j ][i] );
03595 }
03596
03597 if (x3Flux[k][j-1][i].d > 0.0)
03598 de1_l2 = -x2Flux[k-1][j][i].By - emf1_cc[k-1][j-1][i];
03599 else if (x3Flux[k][j-1][i].d < 0.0)
03600 de1_l2 = -x2Flux[k][j][i].By - emf1_cc[k][j-1][i];
03601 else {
03602 de1_l2 = 0.5*(-x2Flux[k-1][j][i].By - emf1_cc[k-1][j-1][i]
03603 -x2Flux[k ][j][i].By - emf1_cc[k ][j-1][i] );
03604 }
03605
03606 if (x3Flux[k][j][i].d > 0.0)
03607 de1_r2 = -x2Flux[k-1][j][i].By - emf1_cc[k-1][j][i];
03608 else if (x3Flux[k][j][i].d < 0.0)
03609 de1_r2 = -x2Flux[k][j][i].By - emf1_cc[k][j][i];
03610 else {
03611 de1_r2 = 0.5*(-x2Flux[k-1][j][i].By - emf1_cc[k-1][j][i]
03612 -x2Flux[k ][j][i].By - emf1_cc[k ][j][i] );
03613 }
03614
03615 emf1[k][j][i] = 0.25*( x3Flux[k][j][i].Bz + x3Flux[k][j-1][i].Bz
03616 - x2Flux[k][j][i].By - x2Flux[k-1][j][i].By
03617 + de1_l2 + de1_r2 + de1_l3 + de1_r3);
03618 }
03619 }
03620 }
03621
03622 return;
03623 }
03624
03625
03626
03627
03628
03629
03630
03631
03632
03633
03634
03635
03636 static void integrate_emf2_corner(const GridS *pG)
03637 {
03638 int i, is = pG->is, ie = pG->ie;
03639 int j, js = pG->js, je = pG->je;
03640 int k, ks = pG->ks, ke = pG->ke;
03641 Real de2_l1, de2_r1, de2_l3, de2_r3;
03642
03643 for (k=ks-1; k<=ke+2; k++) {
03644 for (j=js-2; j<=je+2; j++) {
03645 for (i=is-1; i<=ie+2; i++) {
03646
03647
03648 if (x1Flux[k-1][j][i].d > 0.0)
03649 de2_l3 = -x3Flux[k][j][i-1].By - emf2_cc[k-1][j][i-1];
03650 else if (x1Flux[k-1][j][i].d < 0.0)
03651 de2_l3 = -x3Flux[k][j][i].By - emf2_cc[k-1][j][i];
03652 else {
03653 de2_l3 = 0.5*(-x3Flux[k][j][i-1].By - emf2_cc[k-1][j][i-1]
03654 -x3Flux[k][j][i ].By - emf2_cc[k-1][j][i ] );
03655 }
03656
03657 if (x1Flux[k][j][i].d > 0.0)
03658 de2_r3 = -x3Flux[k][j][i-1].By - emf2_cc[k][j][i-1];
03659 else if (x1Flux[k][j][i].d < 0.0)
03660 de2_r3 = -x3Flux[k][j][i].By - emf2_cc[k][j][i];
03661 else {
03662 de2_r3 = 0.5*(-x3Flux[k][j][i-1].By - emf2_cc[k][j][i-1]
03663 -x3Flux[k][j][i ].By - emf2_cc[k][j][i ] );
03664 }
03665
03666 if (x3Flux[k][j][i-1].d > 0.0)
03667 de2_l1 = x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i-1];
03668 else if (x3Flux[k][j][i-1].d < 0.0)
03669 de2_l1 = x1Flux[k][j][i].Bz - emf2_cc[k][j][i-1];
03670 else {
03671 de2_l1 = 0.5*(x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i-1] +
03672 x1Flux[k ][j][i].Bz - emf2_cc[k ][j][i-1] );
03673 }
03674
03675 if (x3Flux[k][j][i].d > 0.0)
03676 de2_r1 = x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i];
03677 else if (x3Flux[k][j][i].d < 0.0)
03678 de2_r1 = x1Flux[k][j][i].Bz - emf2_cc[k][j][i];
03679 else {
03680 de2_r1 = 0.5*(x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i] +
03681 x1Flux[k ][j][i].Bz - emf2_cc[k ][j][i] );
03682 }
03683
03684 emf2[k][j][i] = 0.25*( x1Flux[k][j][i].Bz + x1Flux[k-1][j][i ].Bz
03685 - x3Flux[k][j][i].By - x3Flux[k ][j][i-1].By
03686 + de2_l1 + de2_r1 + de2_l3 + de2_r3);
03687 }
03688 }
03689 }
03690
03691 return;
03692 }
03693
03694
03695
03696
03697
03698
03699
03700
03701
03702
03703
03704
03705 static void integrate_emf3_corner(const GridS *pG)
03706 {
03707 int i, is = pG->is, ie = pG->ie;
03708 int j, js = pG->js, je = pG->je;
03709 int k, ks = pG->ks, ke = pG->ke;
03710 Real de3_l1, de3_r1, de3_l2, de3_r2;
03711 Real rsf=1.0,lsf=1.0;
03712
03713 for (k=ks-2; k<=ke+2; k++) {
03714 for (j=js-1; j<=je+2; j++) {
03715 for (i=is-1; i<=ie+2; i++) {
03716
03717
03718 #ifdef CYLINDRICAL
03719 rsf = pG->ri[i]/pG->r[i]; lsf = pG->ri[i]/pG->r[i-1];
03720 #endif
03721 if (x1Flux[k][j-1][i].d > 0.0)
03722 de3_l2 = (x2Flux[k][j][i-1].Bz - emf3_cc[k][j-1][i-1])*lsf;
03723 else if (x1Flux[k][j-1][i].d < 0.0)
03724 de3_l2 = (x2Flux[k][j][i].Bz - emf3_cc[k][j-1][i])*rsf;
03725 else {
03726 de3_l2 = 0.5*((x2Flux[k][j][i-1].Bz - emf3_cc[k][j-1][i-1])*lsf +
03727 (x2Flux[k][j][i ].Bz - emf3_cc[k][j-1][i ])*rsf );
03728 }
03729
03730 if (x1Flux[k][j][i].d > 0.0)
03731 de3_r2 = (x2Flux[k][j][i-1].Bz - emf3_cc[k][j][i-1])*lsf;
03732 else if (x1Flux[k][j][i].d < 0.0)
03733 de3_r2 = (x2Flux[k][j][i].Bz - emf3_cc[k][j][i])*rsf;
03734 else {
03735 de3_r2 = 0.5*((x2Flux[k][j][i-1].Bz - emf3_cc[k][j][i-1])*lsf +
03736 (x2Flux[k][j][i ].Bz - emf3_cc[k][j][i ])*rsf );
03737 }
03738
03739 if (x2Flux[k][j][i-1].d > 0.0)
03740 de3_l1 = -x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i-1];
03741 else if (x2Flux[k][j][i-1].d < 0.0)
03742 de3_l1 = -x1Flux[k][j][i].By - emf3_cc[k][j][i-1];
03743 else {
03744 de3_l1 = 0.5*(-x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i-1]
03745 -x1Flux[k][j ][i].By - emf3_cc[k][j ][i-1] );
03746 }
03747
03748 if (x2Flux[k][j][i].d > 0.0)
03749 de3_r1 = -x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i];
03750 else if (x2Flux[k][j][i].d < 0.0)
03751 de3_r1 = -x1Flux[k][j][i].By - emf3_cc[k][j][i];
03752 else {
03753 de3_r1 = 0.5*(-x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i]
03754 -x1Flux[k][j ][i].By - emf3_cc[k][j ][i] );
03755 }
03756
03757 emf3[k][j][i] = 0.25*( x2Flux[k][j ][i-1].Bz + x2Flux[k][j][i].Bz
03758 - x1Flux[k][j-1][i ].By - x1Flux[k][j][i].By
03759 + de3_l1 + de3_r1 + de3_l2 + de3_r2);
03760 }
03761 }
03762 }
03763
03764 return;
03765 }
03766 #endif
03767
03768 #endif