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 #include <math.h>
00032 #include <stdio.h>
00033 #include <stdlib.h>
00034 #include <string.h>
00035 #include "../defs.h"
00036 #include "../athena.h"
00037 #include "../globals.h"
00038 #include "prototypes.h"
00039 #include "../prototypes.h"
00040 #ifdef PARTICLES
00041 #include "../particles/particle.h"
00042 #endif
00043
00044 #if defined(CTU_INTEGRATOR)
00045 #ifdef SPECIAL_RELATIVITY
00046 #error : The CTU integrator cannot be used for special relativity.
00047 #endif
00048
00049
00050 static Cons1DS **Ul_x1Face=NULL, **Ur_x1Face=NULL;
00051 static Cons1DS **Ul_x2Face=NULL, **Ur_x2Face=NULL;
00052 static Cons1DS **x1Flux=NULL, **x2Flux=NULL;
00053
00054
00055 #ifdef MHD
00056 static Real **B1_x1Face=NULL, **B2_x2Face=NULL;
00057 static Real **emf3=NULL, **emf3_cc=NULL;
00058 #endif
00059
00060
00061 static Real *Bxc=NULL, *Bxi=NULL;
00062 static Prim1DS *W=NULL, *Wl=NULL, *Wr=NULL;
00063 static Cons1DS *U1d=NULL, *Ul=NULL, *Ur=NULL;
00064
00065
00066 static Real **dhalf = NULL,**phalf = NULL;
00067
00068
00069 extern Real etah;
00070 #ifdef H_CORRECTION
00071 static Real **eta1=NULL, **eta2=NULL;
00072 #endif
00073
00074
00075 #ifdef CYLINDRICAL
00076 static Real **geom_src=NULL;
00077 #endif
00078
00079
00080
00081
00082
00083
00084 #ifdef MHD
00085 static void integrate_emf3_corner(GridS *pG);
00086 #endif
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 void integrate_2d_ctu(DomainS *pD)
00098 {
00099 GridS *pG=(pD->Grid);
00100 Real dtodx1 = pG->dt/pG->dx1, dtodx2 = pG->dt/pG->dx2;
00101 Real hdtodx1 = 0.5*dtodx1, hdtodx2 = 0.5*dtodx2;
00102 Real hdt = 0.5*pG->dt, dx2 = pG->dx2;
00103 int i,il,iu,is=pG->is, ie=pG->ie;
00104 int j,jl,ju,js=pG->js, je=pG->je;
00105 int ks=pG->ks;
00106 Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic,M1h,M2h,M3h,Bx=0.0;
00107 #ifndef BAROTROPIC
00108 Real coolfl,coolfr,coolf,Eh=0.0;
00109 #endif
00110 #ifdef MHD
00111 Real MHD_src,dbx,dby,B1,B2,B3,V3;
00112 Real B1ch, B2ch, B3ch;
00113 #endif
00114 #ifdef H_CORRECTION
00115 Real cfr,cfl,lambdar,lambdal;
00116 #endif
00117 #if (NSCALARS > 0)
00118 int n;
00119 #endif
00120 #ifdef SELF_GRAVITY
00121 Real gxl,gxr,gyl,gyr,flux_m1l,flux_m1r,flux_m2l,flux_m2r;
00122 #endif
00123 #ifdef FEEDBACK
00124 Real dt1 = 1.0/pG->dt;
00125 #endif
00126 #ifdef SHEARING_BOX
00127
00128 Real Vphi, Mphi;
00129 Real M1n, dM2n, dM3n;
00130 Real M1e, dM2e, dM3e;
00131 Real flx1_dM2, frx1_dM2, flx2_dM2, frx2_dM2;
00132 Real flx1_dM3, frx1_dM3, flx2_dM3, frx2_dM3;
00133 Real fact, qom, om_dt = Omega_0*pG->dt;
00134 #endif
00135 #ifdef STATIC_MESH_REFINEMENT
00136 int ncg,npg,dim;
00137 int ii,ics,ice,jj,jcs,jce,ips,ipe,jps,jpe;
00138 #endif
00139
00140
00141 #ifdef CYLINDRICAL
00142 #ifndef ISOTHERMAL
00143 Real Pavgh;
00144 #endif
00145 Real g,gl,gr,rinv;
00146 Real geom_src_d,geom_src_Vx,geom_src_Vy,geom_src_P,geom_src_By,geom_src_Bz;
00147 const Real *r=pG->r, *ri=pG->ri;
00148 #endif
00149 Real lsf=1.0, rsf=1.0;
00150
00151
00152 #ifdef PARTICLES
00153 Real d1;
00154 il = is - 3;
00155 iu = ie + 3;
00156 jl = js - 3;
00157 ju = je + 3;
00158 #else
00159 il = is - 2;
00160 iu = ie + 2;
00161 jl = js - 2;
00162 ju = je + 2;
00163 #endif
00164
00165
00166 etah = 0.0;
00167
00168
00169 #ifdef FEEDBACK
00170 feedback_predictor(pG);
00171 #endif
00172
00173
00174
00175
00176
00177
00178
00179
00180 for (j=jl; j<=ju; j++) {
00181 for (i=is-nghost; i<=ie+nghost; i++) {
00182 U1d[i].d = pG->U[ks][j][i].d;
00183 U1d[i].Mx = pG->U[ks][j][i].M1;
00184 U1d[i].My = pG->U[ks][j][i].M2;
00185 U1d[i].Mz = pG->U[ks][j][i].M3;
00186 #ifndef BAROTROPIC
00187 U1d[i].E = pG->U[ks][j][i].E;
00188 #endif
00189 #ifdef MHD
00190 U1d[i].By = pG->U[ks][j][i].B2c;
00191 U1d[i].Bz = pG->U[ks][j][i].B3c;
00192 Bxc[i] = pG->U[ks][j][i].B1c;
00193 Bxi[i] = pG->B1i[ks][j][i];
00194 B1_x1Face[j][i] = pG->B1i[ks][j][i];
00195 #endif
00196 #if (NSCALARS > 0)
00197 for (n=0; n<NSCALARS; n++) U1d[i].s[n] = pG->U[ks][j][i].s[n];
00198 #endif
00199 }
00200
00201
00202
00203
00204
00205 for (i=is-nghost; i<=ie+nghost; i++) {
00206 W[i] = Cons1D_to_Prim1D(&U1d[i],&Bxc[i]);
00207
00208
00209
00210
00211 #ifdef CYLINDRICAL
00212 geom_src[j][i] = W[i].d*SQR(W[i].Vy);
00213 #ifdef MHD
00214 geom_src[j][i] += 0.5*(SQR(Bxc[i]) - SQR(W[i].By) + SQR(W[i].Bz));
00215 #endif
00216 #ifdef ISOTHERMAL
00217 geom_src[j][i] += Iso_csound2*W[i].d;
00218 #else
00219 geom_src[j][i] += W[i].P;
00220 #endif
00221 geom_src[j][i] /= x1vc(pG,i);
00222 #endif
00223 }
00224
00225 lr_states(pG,W,Bxc,pG->dt,pG->dx1,il+1,iu-1,Wl,Wr,1);
00226
00227 #ifdef MHD
00228 for (i=il+1; i<=iu; i++) {
00229 #ifdef CYLINDRICAL
00230 rsf = ri[i]/r[i-1]; lsf = ri[i-1]/r[i-1];
00231 #endif
00232 MHD_src = (pG->U[ks][j][i-1].M2/pG->U[ks][j][i-1].d)*
00233 (rsf*pG->B1i[ks][j][i] - lsf*pG->B1i[ks][j][i-1])/pG->dx1;
00234 Wl[i].By += hdt*MHD_src;
00235
00236 #ifdef CYLINDRICAL
00237 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
00238 #endif
00239 MHD_src = (pG->U[ks][j][i].M2/pG->U[ks][j][i].d)*
00240 (rsf*pG->B1i[ks][j][i+1] - lsf*pG->B1i[ks][j][i])/pG->dx1;
00241 Wr[i].By += hdt*MHD_src;
00242 }
00243 #endif
00244
00245
00246
00247
00248
00249 if (StaticGravPot != NULL){
00250 for (i=il+1; i<=iu; i++) {
00251 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00252 #ifdef CYLINDRICAL
00253 gl = (*x1GravAcc)(x1vc(pG,i-1),x2,x3);
00254 gr = (*x1GravAcc)(x1vc(pG,i),x2,x3);
00255
00256
00257 Wl[i].Vx -= hdt*gl;
00258 Wr[i].Vx -= hdt*gr;
00259 #else
00260 phicr = (*StaticGravPot)( x1 ,x2,x3);
00261 phicl = (*StaticGravPot)((x1- pG->dx1),x2,x3);
00262 phifc = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00263
00264 Wl[i].Vx -= dtodx1*(phifc - phicl);
00265 Wr[i].Vx -= dtodx1*(phicr - phifc);
00266 #endif
00267 }
00268 }
00269
00270
00271
00272
00273
00274 #ifdef SELF_GRAVITY
00275 for (i=il+1; i<=iu; i++) {
00276 Wl[i].Vx -= hdtodx1*(pG->Phi[ks][j][i] - pG->Phi[ks][j][i-1]);
00277 Wr[i].Vx -= hdtodx1*(pG->Phi[ks][j][i] - pG->Phi[ks][j][i-1]);
00278 }
00279 #endif
00280
00281
00282
00283
00284
00285 #ifndef BAROTROPIC
00286 if (CoolingFunc != NULL){
00287 for (i=il+1; i<=iu; i++) {
00288 coolfl = (*CoolingFunc)(Wl[i].d,Wl[i].P,(0.5*pG->dt));
00289 coolfr = (*CoolingFunc)(Wr[i].d,Wr[i].P,(0.5*pG->dt));
00290
00291 Wl[i].P -= 0.5*pG->dt*Gamma_1*coolfl;
00292 Wr[i].P -= 0.5*pG->dt*Gamma_1*coolfr;
00293 }
00294 }
00295 #endif
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 #ifdef SHEARING_BOX
00307 if (ShearingBoxPot != NULL){
00308 for (i=il+1; i<=iu; i++) {
00309 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00310 phicr = (*ShearingBoxPot)( x1 ,x2,x3);
00311 phicl = (*ShearingBoxPot)((x1- pG->dx1),x2,x3);
00312 phifc = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
00313
00314 Wl[i].Vx -= dtodx1*(phifc - phicl);
00315 Wr[i].Vx -= dtodx1*(phicr - phifc);
00316 }
00317 }
00318
00319 if (ShBoxCoord == xz){
00320 for (i=il+1; i<=iu; i++) {
00321 Wl[i].Vx += pG->dt*Omega_0*W[i-1].Vz;
00322 Wr[i].Vx += pG->dt*Omega_0*W[i].Vz;
00323 #ifdef FARGO
00324 Wl[i].Vz += hdt*(qshear-2.)*Omega_0*W[i-1].Vx;
00325 Wr[i].Vz += hdt*(qshear-2.)*Omega_0*W[i].Vx;
00326 #else
00327 Wl[i].Vz -= pG->dt*Omega_0*W[i-1].Vx;
00328 Wr[i].Vz -= pG->dt*Omega_0*W[i].Vx;
00329 #endif
00330 }
00331 }
00332
00333 if (ShBoxCoord == xy) {
00334 for (i=il+1; i<=iu; i++) {
00335 Wl[i].Vx += pG->dt*Omega_0*W[i-1].Vy;
00336 Wr[i].Vx += pG->dt*Omega_0*W[i].Vy;
00337 #ifdef FARGO
00338 Wl[i].Vy += hdt*(qshear-2.)*Omega_0*W[i-1].Vx;
00339 Wr[i].Vy += hdt*(qshear-2.)*Omega_0*W[i].Vx;
00340 #else
00341 Wl[i].Vy -= pG->dt*Omega_0*W[i-1].Vx;
00342 Wr[i].Vy -= pG->dt*Omega_0*W[i].Vx;
00343 #endif
00344 }
00345 }
00346 #endif
00347
00348
00349
00350
00351
00352 #ifdef FEEDBACK
00353 for (i=il+1; i<=iu; i++) {
00354 d1 = 1.0/W[i-1].d;
00355 Wl[i].Vx -= pG->Coup[ks][j][i-1].fb1*d1;
00356 Wl[i].Vy -= pG->Coup[ks][j][i-1].fb2*d1;
00357 Wl[i].Vz -= pG->Coup[ks][j][i-1].fb3*d1;
00358
00359 d1 = 1.0/W[i].d;
00360 Wr[i].Vx -= pG->Coup[ks][j][i].fb1*d1;
00361 Wr[i].Vy -= pG->Coup[ks][j][i].fb2*d1;
00362 Wr[i].Vz -= pG->Coup[ks][j][i].fb3*d1;
00363
00364 #ifndef BAROTROPIC
00365 Wl[i].P += pG->Coup[ks][j][i-1].Eloss*Gamma_1;
00366 Wr[i].P += pG->Coup[ks][j][i].Eloss*Gamma_1;
00367 #endif
00368 }
00369 #endif
00370
00371
00372
00373
00374
00375 #ifdef CYLINDRICAL
00376 for (i=il+1; i<=iu; i++) {
00377
00378
00379 rinv = 1.0/x1vc(pG,i-1);
00380 geom_src_d = -W[i-1].d*W[i-1].Vx*rinv;
00381 geom_src_Vx = SQR(W[i-1].Vy);
00382 geom_src_Vy = -W[i-1].Vx*W[i-1].Vy;
00383 #ifdef MHD
00384 geom_src_Vx -= SQR(W[i-1].By)/W[i-1].d;
00385 geom_src_Vy += Bxc[i-1]*W[i-1].By/W[i-1].d;
00386 geom_src_By = -W[i-1].Vy*Bxc[i-1]*rinv;
00387 geom_src_Bz = -W[i-1].Vx*W[i-1].Bz*rinv;
00388 #endif
00389 geom_src_Vx *= rinv;
00390 geom_src_Vy *= rinv;
00391 #ifndef ISOTHERMAL
00392 geom_src_P = -Gamma*W[i-1].P*W[i-1].Vx*rinv;
00393 #endif
00394
00395
00396 Wl[i].d += hdt*geom_src_d;
00397 Wl[i].Vx += hdt*geom_src_Vx;
00398 Wl[i].Vy += hdt*geom_src_Vy;
00399 #ifdef MHD
00400 Wl[i].By += hdt*geom_src_By;
00401 Wl[i].Bz += hdt*geom_src_Bz;
00402 #endif
00403 #ifndef ISOTHERMAL
00404 Wl[i].P += hdt*geom_src_P;
00405 #endif
00406
00407
00408 rinv = 1.0/x1vc(pG,i);
00409 geom_src_d = -W[i].d*W[i].Vx*rinv;
00410 geom_src_Vx = SQR(W[i].Vy);
00411 geom_src_Vy = -W[i].Vx*W[i].Vy;
00412 #ifdef MHD
00413 geom_src_Vx -= SQR(W[i].By)/W[i].d;
00414 geom_src_Vy += Bxc[i]*W[i].By/W[i].d;
00415 geom_src_By = -W[i].Vy*Bxc[i]*rinv;
00416 geom_src_Bz = -W[i].Vx*W[i].Bz*rinv;
00417 #endif
00418 geom_src_Vx *= rinv;
00419 geom_src_Vy *= rinv;
00420 #ifndef ISOTHERMAL
00421 geom_src_P = -Gamma*W[i].P*W[i].Vx*rinv;
00422 #endif
00423
00424
00425 Wr[i].d += hdt*geom_src_d;
00426 Wr[i].Vx += hdt*geom_src_Vx;
00427 Wr[i].Vy += hdt*geom_src_Vy;
00428 #ifdef MHD
00429 Wr[i].By += hdt*geom_src_By;
00430 Wr[i].Bz += hdt*geom_src_Bz;
00431 #endif
00432 #ifndef ISOTHERMAL
00433 Wr[i].P += hdt*geom_src_P;
00434 #endif
00435 }
00436 #endif
00437
00438
00439
00440
00441
00442 for (i=il+1; i<=iu; i++) {
00443 Ul_x1Face[j][i] = Prim1D_to_Cons1D(&Wl[i],&Bxi[i]);
00444 Ur_x1Face[j][i] = Prim1D_to_Cons1D(&Wr[i],&Bxi[i]);
00445
00446 #ifdef MHD
00447 Bx = B1_x1Face[j][i];
00448 #endif
00449 fluxes(Ul_x1Face[j][i],Ur_x1Face[j][i],Wl[i],Wr[i],Bx,&x1Flux[j][i]);
00450 }
00451 }
00452
00453
00454
00455
00456
00457
00458
00459
00460 for (i=il; i<=iu; i++) {
00461 #ifdef CYLINDRICAL
00462 dx2 = r[i]*pG->dx2;
00463 dtodx2 = pG->dt/dx2;
00464 hdtodx2 = 0.5*dtodx2;
00465 #endif
00466 for (j=js-nghost; j<=je+nghost; j++) {
00467 U1d[j].d = pG->U[ks][j][i].d;
00468 U1d[j].Mx = pG->U[ks][j][i].M2;
00469 U1d[j].My = pG->U[ks][j][i].M3;
00470 U1d[j].Mz = pG->U[ks][j][i].M1;
00471 #ifndef BAROTROPIC
00472 U1d[j].E = pG->U[ks][j][i].E;
00473 #endif
00474 #ifdef MHD
00475 U1d[j].By = pG->U[ks][j][i].B3c;
00476 U1d[j].Bz = pG->U[ks][j][i].B1c;
00477 Bxc[j] = pG->U[ks][j][i].B2c;
00478 Bxi[j] = pG->B2i[ks][j][i];
00479 B2_x2Face[j][i] = pG->B2i[ks][j][i];
00480 #endif
00481 #if (NSCALARS > 0)
00482 for (n=0; n<NSCALARS; n++) U1d[j].s[n] = pG->U[ks][j][i].s[n];
00483 #endif
00484 }
00485
00486
00487
00488
00489
00490 for (j=js-nghost; j<=je+nghost; j++) {
00491 W[j] = Cons1D_to_Prim1D(&U1d[j],&Bxc[j]);
00492 }
00493
00494 lr_states(pG,W,Bxc,pG->dt,dx2,jl+1,ju-1,Wl,Wr,2);
00495
00496 #ifdef MHD
00497 for (j=jl+1; j<=ju; j++) {
00498 MHD_src = (pG->U[ks][j-1][i].M1/pG->U[ks][j-1][i].d)*
00499 (pG->B2i[ks][j][i] - pG->B2i[ks][j-1][i])/dx2;
00500 Wl[j].Bz += hdt*MHD_src;
00501
00502 MHD_src = (pG->U[ks][j][i].M1/pG->U[ks][j][i].d)*
00503 (pG->B2i[ks][j+1][i] - pG->B2i[ks][j][i])/dx2;
00504 Wr[j].Bz += hdt*MHD_src;
00505 }
00506 #endif
00507
00508
00509
00510
00511
00512 if (StaticGravPot != NULL){
00513 for (j=jl+1; j<=ju; j++) {
00514 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00515 phicr = (*StaticGravPot)(x1, x2 ,x3);
00516 phicl = (*StaticGravPot)(x1,(x2- pG->dx2),x3);
00517 phifc = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00518
00519 Wl[j].Vx -= dtodx2*(phifc - phicl);
00520 Wr[j].Vx -= dtodx2*(phicr - phifc);
00521 }
00522 }
00523
00524
00525
00526
00527
00528 #ifdef SELF_GRAVITY
00529 for (j=jl+1; j<=ju; j++) {
00530 Wl[j].Vx -= hdtodx2*(pG->Phi[ks][j][i] - pG->Phi[ks][j-1][i]);
00531 Wr[j].Vx -= hdtodx2*(pG->Phi[ks][j][i] - pG->Phi[ks][j-1][i]);
00532 }
00533 #endif
00534
00535
00536
00537
00538
00539 #ifndef BAROTROPIC
00540 if (CoolingFunc != NULL){
00541 for (j=jl+1; j<=ju; j++) {
00542 coolfl = (*CoolingFunc)(Wl[j].d,Wl[j].P,(0.5*pG->dt));
00543 coolfr = (*CoolingFunc)(Wr[j].d,Wr[j].P,(0.5*pG->dt));
00544
00545 Wl[j].P -= 0.5*pG->dt*Gamma_1*coolfl;
00546 Wr[j].P -= 0.5*pG->dt*Gamma_1*coolfr;
00547 }
00548 }
00549 #endif
00550
00551
00552
00553
00554
00555 #ifdef FEEDBACK
00556 for (j=jl+1; j<=ju; j++) {
00557 d1 = 1.0/W[j-1].d;
00558 Wl[j].Vx -= pG->Coup[ks][j-1][i].fb2*d1;
00559 Wl[j].Vy -= pG->Coup[ks][j-1][i].fb3*d1;
00560 Wl[j].Vz -= pG->Coup[ks][j-1][i].fb1*d1;
00561
00562 d1 = 1.0/W[j].d;
00563 Wr[j].Vx -= pG->Coup[ks][j][i].fb2*d1;
00564 Wr[j].Vy -= pG->Coup[ks][j][i].fb3*d1;
00565 Wr[j].Vz -= pG->Coup[ks][j][i].fb1*d1;
00566
00567 #ifndef BAROTROPIC
00568 Wl[i].P += pG->Coup[ks][j-1][i].Eloss*Gamma_1;
00569 Wr[i].P += pG->Coup[ks][j][i].Eloss*Gamma_1;
00570 #endif
00571 }
00572 #endif
00573
00574
00575
00576
00577
00578 for (j=jl+1; j<=ju; j++) {
00579 Ul_x2Face[j][i] = Prim1D_to_Cons1D(&Wl[j],&Bxi[j]);
00580 Ur_x2Face[j][i] = Prim1D_to_Cons1D(&Wr[j],&Bxi[j]);
00581 #ifdef MHD
00582 Bx = B2_x2Face[j][i];
00583 #endif
00584 fluxes(Ul_x2Face[j][i],Ur_x2Face[j][i],Wl[j],Wr[j],Bx,&x2Flux[j][i]);
00585 }
00586 }
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597 #ifdef MHD
00598 for (j=jl; j<=ju; j++) {
00599 for (i=il; i<=iu; i++) {
00600 emf3_cc[j][i] =
00601 (pG->U[ks][j][i].B1c*pG->U[ks][j][i].M2 -
00602 pG->U[ks][j][i].B2c*pG->U[ks][j][i].M1 )/pG->U[ks][j][i].d;
00603 }
00604 }
00605 integrate_emf3_corner(pG);
00606
00607
00608
00609
00610
00611 for (j=jl+1; j<=ju-1; j++) {
00612 for (i=il+1; i<=iu-1; i++) {
00613 #ifdef CYLINDRICAL
00614 hdtodx2 = hdt/(ri[i]*pG->dx2);
00615 #endif
00616 B1_x1Face[j][i] -= hdtodx2*(emf3[j+1][i ] - emf3[j][i]);
00617 B2_x2Face[j][i] += hdtodx1*(emf3[j ][i+1] - emf3[j][i]);
00618 }
00619 #ifdef CYLINDRICAL
00620 hdtodx2 = hdt/(ri[iu]*pG->dx2);
00621 #endif
00622 B1_x1Face[j][iu] -= hdtodx2*(emf3[j+1][iu] - emf3[j][iu]);
00623 }
00624 for (i=il+1; i<=iu-1; i++) {
00625 B2_x2Face[ju][i] += hdtodx1*(emf3[ju][i+1] - emf3[ju][i]);
00626 }
00627 #endif
00628
00629
00630
00631
00632
00633
00634
00635
00636 for (j=jl+1; j<=ju-1; j++) {
00637 for (i=il+1; i<=iu; i++) {
00638 #ifdef CYLINDRICAL
00639 hdtodx2 = hdt/(r[i-1]*pG->dx2);
00640 #endif
00641 Ul_x1Face[j][i].d -= hdtodx2*(x2Flux[j+1][i-1].d - x2Flux[j][i-1].d );
00642 Ul_x1Face[j][i].Mx -= hdtodx2*(x2Flux[j+1][i-1].Mz - x2Flux[j][i-1].Mz);
00643 Ul_x1Face[j][i].My -= hdtodx2*(x2Flux[j+1][i-1].Mx - x2Flux[j][i-1].Mx);
00644 Ul_x1Face[j][i].Mz -= hdtodx2*(x2Flux[j+1][i-1].My - x2Flux[j][i-1].My);
00645 #ifndef BAROTROPIC
00646 Ul_x1Face[j][i].E -= hdtodx2*(x2Flux[j+1][i-1].E - x2Flux[j][i-1].E );
00647 #endif
00648 #ifdef MHD
00649 Ul_x1Face[j][i].Bz -= hdtodx2*(x2Flux[j+1][i-1].By - x2Flux[j][i-1].By);
00650 #endif
00651 #if (NSCALARS > 0)
00652 for (n=0; n<NSCALARS; n++) {
00653 Ul_x1Face[j][i].s[n]-=hdtodx2*(x2Flux[j+1][i-1].s[n]-x2Flux[j][i-1].s[n]);
00654 }
00655 #endif
00656
00657 #ifdef CYLINDRICAL
00658 hdtodx2 = hdt/(r[i]*pG->dx2);
00659 #endif
00660 Ur_x1Face[j][i].d -= hdtodx2*(x2Flux[j+1][i ].d - x2Flux[j][i ].d );
00661 Ur_x1Face[j][i].Mx -= hdtodx2*(x2Flux[j+1][i ].Mz - x2Flux[j][i ].Mz);
00662 Ur_x1Face[j][i].My -= hdtodx2*(x2Flux[j+1][i ].Mx - x2Flux[j][i ].Mx);
00663 Ur_x1Face[j][i].Mz -= hdtodx2*(x2Flux[j+1][i ].My - x2Flux[j][i ].My);
00664 #ifndef BAROTROPIC
00665 Ur_x1Face[j][i].E -= hdtodx2*(x2Flux[j+1][i ].E - x2Flux[j][i ].E );
00666 #endif
00667 #ifdef MHD
00668 Ur_x1Face[j][i].Bz -= hdtodx2*(x2Flux[j+1][i ].By - x2Flux[j][i ].By);
00669 #endif
00670 #if (NSCALARS > 0)
00671 for (n=0; n<NSCALARS; n++) {
00672 Ur_x1Face[j][i].s[n]-=hdtodx2*(x2Flux[j+1][i ].s[n]-x2Flux[j][i ].s[n]);
00673 }
00674 #endif
00675 }
00676 }
00677
00678
00679
00680
00681
00682
00683
00684 #ifdef MHD
00685 for (j=jl+1; j<=ju-1; j++) {
00686 for (i=il+1; i<=iu; i++) {
00687 #ifdef CYLINDRICAL
00688 dbx = (ri[i]*pG->B1i[ks][j][i] - ri[i-1]*pG->B1i[ks][j][i-1])/r[i-1];
00689 #else
00690 dbx = pG->B1i[ks][j][i] - pG->B1i[ks][j][i-1];
00691 #endif
00692 B1 = pG->U[ks][j][i-1].B1c;
00693 B2 = pG->U[ks][j][i-1].B2c;
00694 B3 = pG->U[ks][j][i-1].B3c;
00695 V3 = pG->U[ks][j][i-1].M3/pG->U[ks][j][i-1].d;
00696
00697 Ul_x1Face[j][i].Mx += hdtodx1*B1*dbx;
00698 Ul_x1Face[j][i].My += hdtodx1*B2*dbx;
00699 Ul_x1Face[j][i].Mz += hdtodx1*B3*dbx;
00700 Ul_x1Face[j][i].Bz += hdtodx1*V3*dbx;
00701 #ifndef BAROTROPIC
00702 Ul_x1Face[j][i].E += hdtodx1*B3*V3*dbx;
00703 #endif
00704
00705 #ifdef CYLINDRICAL
00706 dbx = (ri[i+1]*pG->B1i[ks][j][i+1] - ri[i]*pG->B1i[ks][j][i])/r[i];
00707 #else
00708 dbx = pG->B1i[ks][j][i+1] - pG->B1i[ks][j][i];
00709 #endif
00710 B1 = pG->U[ks][j][i].B1c;
00711 B2 = pG->U[ks][j][i].B2c;
00712 B3 = pG->U[ks][j][i].B3c;
00713 V3 = pG->U[ks][j][i].M3/pG->U[ks][j][i].d;
00714
00715 Ur_x1Face[j][i].Mx += hdtodx1*B1*dbx;
00716 Ur_x1Face[j][i].My += hdtodx1*B2*dbx;
00717 Ur_x1Face[j][i].Mz += hdtodx1*B3*dbx;
00718 Ur_x1Face[j][i].Bz += hdtodx1*V3*dbx;
00719 #ifndef BAROTROPIC
00720 Ur_x1Face[j][i].E += hdtodx1*B3*V3*dbx;
00721 #endif
00722 }
00723 }
00724 #endif
00725
00726
00727
00728
00729
00730
00731
00732
00733 if (StaticGravPot != NULL){
00734 for (j=jl+1; j<=ju-1; j++) {
00735 for (i=il+1; i<=iu; i++) {
00736 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00737 phic = (*StaticGravPot)(x1, x2 ,x3);
00738 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
00739 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00740
00741 #ifdef CYLINDRICAL
00742 hdtodx2 = hdt/(r[i]*pG->dx2);
00743 #endif
00744 Ur_x1Face[j][i].My -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
00745 #ifndef BAROTROPIC
00746 Ur_x1Face[j][i].E -= hdtodx2*(x2Flux[j ][i ].d*(phic - phil) +
00747 x2Flux[j+1][i ].d*(phir - phic));
00748 #endif
00749
00750 phic = (*StaticGravPot)((x1-pG->dx1), x2 ,x3);
00751 phir = (*StaticGravPot)((x1-pG->dx1),(x2+0.5*pG->dx2),x3);
00752 phil = (*StaticGravPot)((x1-pG->dx1),(x2-0.5*pG->dx2),x3);
00753
00754 #ifdef CYLINDRICAL
00755 hdtodx2 = hdt/(r[i-1]*pG->dx2);
00756 #endif
00757 Ul_x1Face[j][i].My -= hdtodx2*(phir-phil)*pG->U[ks][j][i-1].d;
00758 #ifndef BAROTROPIC
00759 Ul_x1Face[j][i].E -= hdtodx2*(x2Flux[j ][i-1].d*(phic - phil) +
00760 x2Flux[j+1][i-1].d*(phir - phic));
00761 #endif
00762 }
00763 }
00764 }
00765
00766
00767
00768
00769
00770
00771 #ifdef SELF_GRAVITY
00772 for (j=jl+1; j<=ju-1; j++) {
00773 for (i=il+1; i<=iu; i++) {
00774 phic = pG->Phi[ks][j][i];
00775 phir = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j+1][i]);
00776 phil = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j-1][i]);
00777
00778 Ur_x1Face[j][i].My -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
00779 #ifndef BAROTROPIC
00780 Ur_x1Face[j][i].E -= hdtodx2*(x2Flux[j ][i ].d*(phic - phil) +
00781 x2Flux[j+1][i ].d*(phir - phic));
00782 #endif
00783
00784 phic = pG->Phi[ks][j][i-1];
00785 phir = 0.5*(pG->Phi[ks][j][i-1] + pG->Phi[ks][j+1][i-1]);
00786 phil = 0.5*(pG->Phi[ks][j][i-1] + pG->Phi[ks][j-1][i-1]);
00787
00788 Ul_x1Face[j][i].My -= hdtodx2*(phir-phil)*pG->U[ks][j][i-1].d;
00789 #ifndef BAROTROPIC
00790 Ul_x1Face[j][i].E -= hdtodx2*(x2Flux[j ][i-1].d*(phic - phil) +
00791 x2Flux[j+1][i-1].d*(phir - phic));
00792 #endif
00793 }
00794 }
00795 #endif
00796
00797
00798
00799
00800
00801
00802
00803
00804 for (j=jl+1; j<=ju; j++) {
00805 for (i=il+1; i<=iu-1; i++) {
00806 #ifdef CYLINDRICAL
00807 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
00808 #endif
00809 Ul_x2Face[j][i].d -= hdtodx1*(rsf*x1Flux[j-1][i+1].d - lsf*x1Flux[j-1][i].d );
00810 Ul_x2Face[j][i].Mx -= hdtodx1*(SQR(rsf)*x1Flux[j-1][i+1].My - SQR(lsf)*x1Flux[j-1][i].My);
00811 Ul_x2Face[j][i].My -= hdtodx1*(rsf*x1Flux[j-1][i+1].Mz - lsf*x1Flux[j-1][i].Mz);
00812 Ul_x2Face[j][i].Mz -= hdtodx1*(rsf*x1Flux[j-1][i+1].Mx - lsf*x1Flux[j-1][i].Mx);
00813 #ifndef BAROTROPIC
00814 Ul_x2Face[j][i].E -= hdtodx1*(rsf*x1Flux[j-1][i+1].E - lsf*x1Flux[j-1][i].E );
00815 #endif
00816 #ifdef MHD
00817 Ul_x2Face[j][i].By -= hdtodx1*(rsf*x1Flux[j-1][i+1].Bz - lsf*x1Flux[j-1][i].Bz);
00818 #endif
00819
00820 Ur_x2Face[j][i].d -= hdtodx1*(rsf*x1Flux[j ][i+1].d - lsf*x1Flux[j ][i].d );
00821 Ur_x2Face[j][i].Mx -= hdtodx1*(SQR(rsf)*x1Flux[j ][i+1].My - SQR(lsf)*x1Flux[j ][i].My);
00822 Ur_x2Face[j][i].My -= hdtodx1*(rsf*x1Flux[j ][i+1].Mz - lsf*x1Flux[j ][i].Mz);
00823 Ur_x2Face[j][i].Mz -= hdtodx1*(rsf*x1Flux[j ][i+1].Mx - lsf*x1Flux[j ][i].Mx);
00824 #ifndef BAROTROPIC
00825 Ur_x2Face[j][i].E -= hdtodx1*(rsf*x1Flux[j ][i+1].E - lsf*x1Flux[j ][i].E );
00826 #endif
00827 #ifdef MHD
00828 Ur_x2Face[j][i].By -= hdtodx1*(rsf*x1Flux[j ][i+1].Bz - lsf*x1Flux[j ][i].Bz);
00829 #endif
00830 #if (NSCALARS > 0)
00831 for (n=0; n<NSCALARS; n++) {
00832 Ul_x2Face[j][i].s[n]-=hdtodx1*(rsf*x1Flux[j-1][i+1].s[n]-lsf*x1Flux[j-1][i].s[n]);
00833 Ur_x2Face[j][i].s[n]-=hdtodx1*(rsf*x1Flux[j ][i+1].s[n]-lsf*x1Flux[j ][i].s[n]);
00834 }
00835 #endif
00836 }
00837 }
00838
00839
00840
00841
00842
00843
00844
00845 #ifdef MHD
00846 for (j=jl+1; j<=ju; j++) {
00847 for (i=il+1; i<=iu-1; i++) {
00848 #ifdef CYLINDRICAL
00849 hdtodx2 = hdt/(r[i]*pG->dx2);
00850 #endif
00851 dby = pG->B2i[ks][j][i] - pG->B2i[ks][j-1][i];
00852 B1 = pG->U[ks][j-1][i].B1c;
00853 B2 = pG->U[ks][j-1][i].B2c;
00854 B3 = pG->U[ks][j-1][i].B3c;
00855 V3 = pG->U[ks][j-1][i].M3/pG->U[ks][j-1][i].d;
00856
00857 Ul_x2Face[j][i].Mz += hdtodx2*B1*dby;
00858 Ul_x2Face[j][i].Mx += hdtodx2*B2*dby;
00859 Ul_x2Face[j][i].My += hdtodx2*B3*dby;
00860 Ul_x2Face[j][i].By += hdtodx2*V3*dby;
00861 #ifndef BAROTROPIC
00862 Ul_x2Face[j][i].E += hdtodx2*B3*V3*dby;
00863 #endif
00864
00865 dby = pG->B2i[ks][j+1][i] - pG->B2i[ks][j][i];
00866 B1 = pG->U[ks][j][i].B1c;
00867 B2 = pG->U[ks][j][i].B2c;
00868 B3 = pG->U[ks][j][i].B3c;
00869 V3 = pG->U[ks][j][i].M3/pG->U[ks][j][i].d;
00870
00871 Ur_x2Face[j][i].Mz += hdtodx2*B1*dby;
00872 Ur_x2Face[j][i].Mx += hdtodx2*B2*dby;
00873 Ur_x2Face[j][i].My += hdtodx2*B3*dby;
00874 Ur_x2Face[j][i].By += hdtodx2*V3*dby;
00875 #ifndef BAROTROPIC
00876 Ur_x2Face[j][i].E += hdtodx2*B3*V3*dby;
00877 #endif
00878 }
00879 }
00880 #endif
00881
00882
00883
00884
00885
00886
00887
00888
00889 if (StaticGravPot != NULL){
00890 for (j=jl+1; j<=ju; j++) {
00891 for (i=il+1; i<=iu-1; i++) {
00892 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00893 phic = (*StaticGravPot)((x1 ),x2,x3);
00894 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00895 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00896
00897 #ifdef CYLINDRICAL
00898 g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
00899 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
00900 Ur_x2Face[j][i].Mz -= hdt*pG->U[ks][j][i].d*g;
00901 #else
00902 Ur_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
00903 #endif
00904 #ifndef BAROTROPIC
00905 Ur_x2Face[j][i].E -= hdtodx1*(lsf*x1Flux[j ][i ].d*(phic - phil) +
00906 rsf*x1Flux[j ][i+1].d*(phir - phic));
00907 #endif
00908
00909 phic = (*StaticGravPot)((x1 ),(x2-pG->dx2),x3);
00910 phir = (*StaticGravPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
00911 phil = (*StaticGravPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
00912
00913 #ifdef CYLINDRICAL
00914 g = (*x1GravAcc)(x1vc(pG,i),(x2-pG->dx2),x3);
00915 Ul_x2Face[j][i].Mz -= hdt*pG->U[ks][j-1][i].d*g;
00916 #else
00917 Ul_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j-1][i].d;
00918 #endif
00919 #ifndef BAROTROPIC
00920 Ul_x2Face[j][i].E -= hdtodx1*(lsf*x1Flux[j-1][i ].d*(phic - phil) +
00921 rsf*x1Flux[j-1][i+1].d*(phir - phic));
00922 #endif
00923 }
00924 }
00925 }
00926
00927
00928
00929
00930
00931
00932 #ifdef SELF_GRAVITY
00933 for (j=jl+1; j<=ju; j++) {
00934 for (i=il+1; i<=iu-1; i++) {
00935 phic = pG->Phi[ks][j][i];
00936 phir = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j][i+1]);
00937 phil = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j][i-1]);
00938
00939 Ur_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
00940 #ifndef BAROTROPIC
00941 Ur_x2Face[j][i].E -= hdtodx1*(x1Flux[j ][i ].d*(phic - phil) +
00942 x1Flux[j ][i+1].d*(phir - phic));
00943 #endif
00944
00945 phic = pG->Phi[ks][j-1][i];
00946 phir = 0.5*(pG->Phi[ks][j-1][i] + pG->Phi[ks][j-1][i+1]);
00947 phil = 0.5*(pG->Phi[ks][j-1][i] + pG->Phi[ks][j-1][i-1]);
00948
00949 Ul_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j-1][i].d;
00950 #ifndef BAROTROPIC
00951 Ul_x2Face[j][i].E -= hdtodx1*(x1Flux[j-1][i ].d*(phic - phil) +
00952 x1Flux[j-1][i+1].d*(phir - phic));
00953 #endif
00954 }
00955 }
00956
00957 #endif
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968 #ifdef SHEARING_BOX
00969 if (ShearingBoxPot != NULL){
00970 for (j=jl+1; j<=ju; j++) {
00971 for (i=il+1; i<=iu-1; i++) {
00972 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00973 phic = (*ShearingBoxPot)((x1 ),x2,x3);
00974 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
00975 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
00976
00977 Ur_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
00978 #ifndef BAROTROPIC
00979 Ur_x2Face[j][i].E -= hdtodx1*(x1Flux[j ][i ].d*(phic - phil) +
00980 x1Flux[j ][i+1].d*(phir - phic));
00981 #endif
00982
00983 phic = (*ShearingBoxPot)((x1 ),(x2-pG->dx2),x3);
00984 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
00985 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
00986
00987 Ul_x2Face[j][i].Mz -= hdtodx1*(phir-phil)*pG->U[ks][j-1][i].d;
00988 #ifndef BAROTROPIC
00989 Ul_x2Face[j][i].E -= hdtodx1*(x1Flux[j-1][i ].d*(phic - phil) +
00990 x1Flux[j-1][i+1].d*(phir - phic));
00991 #endif
00992 }
00993 }
00994 }
00995
00996 if (ShBoxCoord == xz){
00997 for (j=jl+1; j<=ju; j++) {
00998 for (i=il+1; i<=iu-1; i++) {
00999 Ur_x2Face[j][i].Mz += pG->dt*Omega_0*pG->U[ks][j][i].M3;
01000 Ul_x2Face[j][i].Mz += pG->dt*Omega_0*pG->U[ks][j-1][i].M3;
01001 #ifdef FARGO
01002 Ur_x2Face[j][i].My += hdt*(qshear-2.)*Omega_0*pG->U[ks][j][i].M1;
01003 Ul_x2Face[j][i].My += hdt*(qshear-2.)*Omega_0*pG->U[ks][j-1][i].M1;
01004 #else
01005 Ur_x2Face[j][i].My -= pG->dt*Omega_0*pG->U[ks][j][i].M1;
01006 Ul_x2Face[j][i].My -= pG->dt*Omega_0*pG->U[ks][j-1][i].M1;
01007 #endif
01008 }
01009 }
01010 }
01011
01012 if (ShBoxCoord == xy){
01013 for (j=jl+1; j<=ju; j++) {
01014 for (i=il+1; i<=iu-1; i++) {
01015 Ur_x2Face[j][i].Mz += pG->dt*Omega_0*pG->U[ks][j][i].M2;
01016 Ul_x2Face[j][i].Mz += pG->dt*Omega_0*pG->U[ks][j-1][i].M2;
01017 #ifdef FARGO
01018 Ur_x2Face[j][i].Mx += hdt*(qshear-2.)*Omega_0*pG->U[ks][j][i].M1;
01019 Ul_x2Face[j][i].Mx += hdt*(qshear-2.)*Omega_0*pG->U[ks][j-1][i].M1;
01020 #else
01021 Ur_x2Face[j][i].Mx -= pG->dt*Omega_0*pG->U[ks][j][i].M1;
01022 Ul_x2Face[j][i].Mx -= pG->dt*Omega_0*pG->U[ks][j-1][i].M1;
01023 #endif
01024 }
01025 }
01026 }
01027 #endif
01028
01029
01030
01031
01032 #ifdef CYLINDRICAL
01033 for (j=js-1; j<=ju; j++) {
01034 for (i=is-1; i<=ie+1; i++) {
01035 Ur_x2Face[j][i].Mz += hdt*geom_src[j ][i];
01036 Ul_x2Face[j][i].Mz += hdt*geom_src[j-1][i];
01037 }
01038 }
01039 #endif
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049 #ifndef MHD
01050 #ifndef PARTICLES
01051 if ((StaticGravPot != NULL) || (CoolingFunc != NULL))
01052 #endif
01053 #endif
01054 {
01055 for (j=jl+1; j<=ju-1; j++) {
01056 for (i=il+1; i<=iu-1; i++) {
01057 #ifdef CYLINDRICAL
01058 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01059 hdtodx2 = hdt/(r[i]*pG->dx2);
01060 #endif
01061 dhalf[j][i] = pG->U[ks][j][i].d
01062 - hdtodx1*(rsf*x1Flux[j ][i+1].d - lsf*x1Flux[j][i].d)
01063 - hdtodx2*( x2Flux[j+1][i ].d - x2Flux[j][i].d);
01064 #ifdef PARTICLES
01065 pG->Coup[ks][j][i].grid_d = dhalf[j][i];
01066 #endif
01067 }
01068 }
01069 }
01070
01071
01072
01073
01074
01075 #ifndef MHD
01076 #ifndef PARTICLES
01077 if (CoolingFunc != NULL)
01078 #endif
01079 #endif
01080 {
01081 for (j=jl+1; j<=ju-1; j++) {
01082 for (i=il+1; i<=iu-1; i++) {
01083 #ifdef CYLINDRICAL
01084 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01085 hdtodx2 = hdt/(r[i]*pG->dx2);
01086 #endif
01087 M1h = pG->U[ks][j][i].M1
01088 - hdtodx1*(rsf*x1Flux[j][i+1].Mx - lsf*x1Flux[j][i].Mx)
01089 - hdtodx2*( x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
01090
01091 M2h = pG->U[ks][j][i].M2
01092 - hdtodx1*(SQR(rsf)*x1Flux[j][i+1].My - SQR(lsf)*x1Flux[j][i].My)
01093 - hdtodx2*( x2Flux[j+1][i].Mx - x2Flux[j][i].Mx);
01094
01095 M3h = pG->U[ks][j][i].M3
01096 - hdtodx1*(rsf*x1Flux[j][i+1].Mz - lsf*x1Flux[j][i].Mz)
01097 - hdtodx2*( x2Flux[j+1][i].My - x2Flux[j][i].My);
01098
01099 #ifndef BAROTROPIC
01100 Eh = pG->U[ks][j][i].E
01101 - hdtodx1*(rsf*x1Flux[j][i+1].E - lsf*x1Flux[j][i].E)
01102 - hdtodx2*( x2Flux[j+1][i].E - x2Flux[j][i].E);
01103 #endif
01104
01105
01106 if (StaticGravPot != NULL){
01107 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01108 #ifdef CYLINDRICAL
01109 g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
01110 M1h -= hdt*pG->U[ks][j][i].d*g;
01111 #else
01112 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
01113 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
01114 M1h -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
01115 #endif
01116 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01117 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01118 M2h -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
01119 }
01120
01121
01122 #ifdef SELF_GRAVITY
01123 phir = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j][i+1]);
01124 phil = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j][i-1]);
01125 M1h -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
01126
01127 phir = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j+1][i]);
01128 phil = 0.5*(pG->Phi[ks][j][i] + pG->Phi[ks][j-1][i]);
01129 M2h -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
01130 #endif
01131
01132
01133 #ifdef SHEARING_BOX
01134 if (ShearingBoxPot != NULL){
01135 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01136 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
01137 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
01138 M1h -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
01139
01140 phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
01141 phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
01142 M2h -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
01143 }
01144
01145 if (ShBoxCoord == xy) M1h += pG->dt*Omega_0*pG->U[ks][j][i].M2;
01146 if (ShBoxCoord == xz) M1h += pG->dt*Omega_0*pG->U[ks][j][i].M3;
01147 #ifdef FARGO
01148 if (ShBoxCoord == xy) M2h += hdt*(qshear-2.)*Omega_0*pG->U[ks][j][i].M1;
01149 if (ShBoxCoord == xz) M3h += hdt*(qshear-2.)*Omega_0*pG->U[ks][j][i].M1;
01150 #else
01151 if (ShBoxCoord == xy) M2h -= pG->dt*Omega_0*pG->U[ks][j][i].M1;
01152 if (ShBoxCoord == xz) M3h -= pG->dt*Omega_0*pG->U[ks][j][i].M1;
01153 #endif
01154 #endif
01155
01156
01157 #ifdef FEEDBACK
01158 M1h -= pG->Coup[ks][j][i].fb1;
01159 M2h -= pG->Coup[ks][j][i].fb2;
01160 M3h -= pG->Coup[ks][j][i].fb3;
01161 #endif
01162
01163
01164
01165 #ifdef CYLINDRICAL
01166 M1h += hdt*geom_src[j][i];
01167 #endif
01168
01169 #ifndef BAROTROPIC
01170 phalf[j][i] = Eh - 0.5*(M1h*M1h + M2h*M2h + M3h*M3h)/dhalf[j][i];
01171 #endif
01172
01173 #ifdef MHD
01174 B1ch = 0.5*(lsf*B1_x1Face[j][i] + rsf*B1_x1Face[j][i+1]);
01175 B2ch = 0.5*( B2_x2Face[j][i] + B2_x2Face[j+1][i]);
01176 B3ch = pG->U[ks][j][i].B3c
01177 - hdtodx1*(rsf*x1Flux[j][i+1].Bz - lsf*x1Flux[j][i].Bz)
01178 - hdtodx2*( x2Flux[j+1][i].By - x2Flux[j][i].By);
01179 emf3_cc[j][i] = (B1ch*M2h - B2ch*M1h)/dhalf[j][i];
01180
01181 #ifndef BAROTROPIC
01182 phalf[j][i] -= 0.5*(B1ch*B1ch + B2ch*B2ch + B3ch*B3ch);
01183 #endif
01184 #endif
01185
01186 #ifndef BAROTROPIC
01187 phalf[j][i] *= Gamma_1;
01188 #endif
01189
01190 #ifdef PARTICLES
01191 d1 = 1.0/dhalf[j][i];
01192 pG->Coup[ks][j][i].grid_v1 = M1h*d1;
01193 pG->Coup[ks][j][i].grid_v2 = M2h*d1;
01194 pG->Coup[ks][j][i].grid_v3 = M3h*d1;
01195 #ifndef BAROTROPIC
01196 pG->Coup[ks][j][i].grid_cs = sqrt(Gamma*phalf[j][i]*d1);
01197 #endif
01198 #endif
01199 }
01200 }
01201 }
01202
01203
01204
01205 #ifdef PARTICLES
01206 Integrate_Particles(pG,pD);
01207 #ifdef FEEDBACK
01208 exchange_feedback(pG, pD);
01209 #endif
01210 #endif
01211
01212
01213
01214
01215
01216
01217
01218
01219 #ifdef H_CORRECTION
01220 for (j=js-1; j<=je+1; j++) {
01221 for (i=is-1; i<=ie+2; i++) {
01222 #ifdef MHD
01223 Bx = B1_x1Face[j][i];
01224 #endif
01225 cfr = cfast(&(Ur_x1Face[j][i]),&Bx);
01226 cfl = cfast(&(Ul_x1Face[j][i]),&Bx);
01227 lambdar = Ur_x1Face[j][i].Mx/Ur_x1Face[j][i].d + cfr;
01228 lambdal = Ul_x1Face[j][i].Mx/Ul_x1Face[j][i].d - cfl;
01229 eta1[j][i] = 0.5*fabs(lambdar - lambdal);
01230 }
01231 }
01232
01233 for (j=js-1; j<=je+2; j++) {
01234 for (i=is-1; i<=ie+1; i++) {
01235 #ifdef MHD
01236 Bx = B2_x2Face[j][i];
01237 #endif
01238 cfr = cfast(&(Ur_x2Face[j][i]),&Bx);
01239 cfl = cfast(&(Ul_x2Face[j][i]),&Bx);
01240 lambdar = Ur_x2Face[j][i].Mx/Ur_x2Face[j][i].d + cfr;
01241 lambdal = Ul_x2Face[j][i].Mx/Ul_x2Face[j][i].d - cfl;
01242 eta2[j][i] = 0.5*fabs(lambdar - lambdal);
01243 }
01244 }
01245 #endif
01246
01247
01248
01249
01250
01251 for (j=js-1; j<=je+1; j++) {
01252 for (i=is; i<=ie+1; i++) {
01253 #ifdef H_CORRECTION
01254 etah = MAX(eta2[j][i-1],eta2[j][i]);
01255 etah = MAX(etah,eta2[j+1][i-1]);
01256 etah = MAX(etah,eta2[j+1][i ]);
01257 etah = MAX(etah,eta1[j ][i ]);
01258 #endif
01259 #ifdef MHD
01260 Bx = B1_x1Face[j][i];
01261 #endif
01262 Wl[i] = Cons1D_to_Prim1D(&Ul_x1Face[j][i],&Bx);
01263 Wr[i] = Cons1D_to_Prim1D(&Ur_x1Face[j][i],&Bx);
01264
01265 fluxes(Ul_x1Face[j][i],Ur_x1Face[j][i],Wl[i],Wr[i],Bx,&x1Flux[j][i]);
01266 }
01267 }
01268
01269
01270
01271
01272
01273 for (j=js; j<=je+1; j++) {
01274 for (i=is-1; i<=ie+1; i++) {
01275 #ifdef H_CORRECTION
01276 etah = MAX(eta1[j-1][i],eta1[j][i]);
01277 etah = MAX(etah,eta1[j-1][i+1]);
01278 etah = MAX(etah,eta1[j ][i+1]);
01279 etah = MAX(etah,eta2[j ][i ]);
01280 #endif
01281 #ifdef MHD
01282 Bx = B2_x2Face[j][i];
01283 #endif
01284 Wl[i] = Cons1D_to_Prim1D(&Ul_x2Face[j][i],&Bx);
01285 Wr[i] = Cons1D_to_Prim1D(&Ur_x2Face[j][i],&Bx);
01286
01287 fluxes(Ul_x2Face[j][i],Ur_x2Face[j][i],Wl[i],Wr[i],Bx,&x2Flux[j][i]);
01288 }
01289 }
01290
01291
01292
01293
01294
01295
01296
01297 #ifdef MHD
01298 integrate_emf3_corner(pG);
01299
01300
01301
01302
01303
01304 for (j=js; j<=je; j++) {
01305 for (i=is; i<=ie; i++) {
01306 #ifdef CYLINDRICAL
01307 dtodx2 = pG->dt/(ri[i]*pG->dx2);
01308 #endif
01309 pG->B1i[ks][j][i] -= dtodx2*(emf3[j+1][i ] - emf3[j][i]);
01310 pG->B2i[ks][j][i] += dtodx1*(emf3[j ][i+1] - emf3[j][i]);
01311 }
01312 #ifdef CYLINDRICAL
01313 dtodx2 = pG->dt/(ri[ie+1]*pG->dx2);
01314 #endif
01315 pG->B1i[ks][j][ie+1] -= dtodx2*(emf3[j+1][ie+1] - emf3[j][ie+1]);
01316 }
01317 for (i=is; i<=ie; i++) {
01318 pG->B2i[ks][je+1][i] += dtodx1*(emf3[je+1][i+1] - emf3[je+1][i]);
01319 }
01320 #endif
01321
01322
01323
01324
01325
01326
01327
01328 #ifdef CYLINDRICAL
01329 for (j=js; j<=je; j++) {
01330 for (i=is; i<=ie; i++) {
01331 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01332 hdtodx2 = hdt/(r[i]*pG->dx2);
01333
01334
01335 dhalf[j][i] = pG->U[ks][j][i].d
01336 - hdtodx1*(rsf*x1Flux[j ][i+1].d - lsf*x1Flux[j][i].d)
01337 - hdtodx2*( x2Flux[j+1][i ].d - x2Flux[j][i].d);
01338
01339
01340 M2h = pG->U[ks][j][i].M2
01341 - hdtodx1*(SQR(rsf)*x1Flux[j][i+1].My - SQR(lsf)*x1Flux[j][i].My)
01342 - hdtodx2*( x2Flux[j+1][i].Mx - x2Flux[j][i].Mx);
01343
01344
01345 if (StaticGravPot != NULL){
01346 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01347 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01348 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01349 M2h -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
01350 }
01351
01352
01353 geom_src[j][i] = SQR(M2h)/dhalf[j][i];
01354 #ifdef MHD
01355 B2ch = 0.5*(B2_x2Face[j][i] + B2_x2Face[j+1][i]);
01356 geom_src[j][i] -= SQR(B2ch);
01357 #endif
01358 #ifdef ISOTHERMAL
01359 geom_src[j][i] += Iso_csound2*dhalf[j][i];
01360 #ifdef MHD
01361 B1ch = 0.5*(rsf*B1_x1Face[j][i+1] + lsf*B1_x1Face[j][i]);
01362 B3ch = pG->U[ks][j][i].B3c
01363 - hdtodx1*(rsf*x1Flux[j ][i+1].Bz - lsf*x1Flux[j][i].Bz)
01364 - hdtodx2*( x2Flux[j+1][i ].By - x2Flux[j][i].By);
01365 geom_src[j][i] += 0.5*(SQR(B1ch)+SQR(B2ch)+SQR(B3ch));
01366 #endif
01367 #else
01368 Pavgh = 0.5*(lsf*x1Flux[j][i].Pflux + rsf*x1Flux[j][i+1].Pflux);
01369 geom_src[j][i] += Pavgh;
01370 #endif
01371 geom_src[j][i] /= x1vc(pG,i);
01372
01373
01374 pG->U[ks][j][i].M1 += pG->dt*geom_src[j][i];
01375 }
01376 }
01377 #endif
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388 #ifdef SHEARING_BOX
01389 fact = om_dt/(2. + (2.-qshear)*om_dt*om_dt);
01390 qom = qshear*Omega_0;
01391 for(j=js; j<=je; j++) {
01392 for(i=is; i<=ie; i++) {
01393 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01394
01395
01396 M1n = pG->U[ks][j][i].M1;
01397 #ifdef FARGO
01398 if (ShBoxCoord==xy) dM2n = pG->U[ks][j][i].M2;
01399 if (ShBoxCoord==xz) dM3n = pG->U[ks][j][i].M3;
01400 #else
01401 if (ShBoxCoord==xy) dM2n = pG->U[ks][j][i].M2 + qom*x1*pG->U[ks][j][i].d;
01402 if (ShBoxCoord==xz) dM3n = pG->U[ks][j][i].M3 + qom*x1*pG->U[ks][j][i].d;
01403 #endif
01404
01405
01406 if (ShBoxCoord==xy){
01407 frx1_dM2 = x1Flux[j][i+1].My;
01408 flx1_dM2 = x1Flux[j][i ].My;
01409 frx2_dM2 = x2Flux[j+1][i].Mx;
01410 flx2_dM2 = x2Flux[j ][i].Mx;
01411 }
01412 if (ShBoxCoord==xz){
01413 frx1_dM3 = x1Flux[j][i+1].Mz;
01414 flx1_dM3 = x1Flux[j][i ].Mz;
01415 frx2_dM3 = x2Flux[j+1][i].My;
01416 flx2_dM3 = x2Flux[j ][i].My;
01417 }
01418 #ifndef FARGO
01419 if (ShBoxCoord==xy){
01420 frx1_dM2 += qom*(x1+0.5*pG->dx1)*x1Flux[j][i+1].d;
01421 flx1_dM2 += qom*(x1-0.5*pG->dx1)*x1Flux[j][i ].d;
01422 frx2_dM2 += qom*(x1 )*x2Flux[j+1][i].d;
01423 flx2_dM2 += qom*(x1 )*x2Flux[j ][i].d;
01424 }
01425 if (ShBoxCoord==xz){
01426 frx1_dM3 += qom*(x1+0.5*pG->dx1)*x1Flux[j][i+1].d;
01427 flx1_dM3 += qom*(x1-0.5*pG->dx1)*x1Flux[j][i ].d;
01428 frx2_dM3 += qom*(x1 )*x2Flux[j+1][i].d;
01429 flx2_dM3 += qom*(x1 )*x2Flux[j ][i].d;
01430 }
01431 #endif
01432
01433
01434 M1e = M1n - hdtodx1*(x1Flux[j][i+1].Mx - x1Flux[j][i].Mx)
01435 - hdtodx2*(x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
01436
01437 if (ShBoxCoord==xy){
01438 dM2e = dM2n - hdtodx1*(frx1_dM2 - flx1_dM2)
01439 - hdtodx2*(frx2_dM2 - flx2_dM2);
01440 }
01441 if (ShBoxCoord==xz){
01442 dM3e = dM3n - hdtodx1*(frx1_dM3 - flx1_dM3)
01443 - hdtodx2*(frx2_dM3 - flx2_dM3);
01444 }
01445
01446 #ifdef FEEDBACK
01447 M1e -= 0.5*pG->Coup[ks][j][i].fb1;
01448 if (ShBoxCoord==xy) dM2e -= 0.5*pG->Coup[ks][j][i].fb2;
01449 if (ShBoxCoord==xz) dM3e -= 0.5*pG->Coup[ks][j][i].fb3;
01450 #endif
01451
01452
01453
01454
01455
01456 if (ShBoxCoord==xy){
01457 pG->U[ks][j][i].M1 += (4.0*dM2e + 2.0*(qshear-2.)*om_dt*M1e)*fact;
01458 pG->U[ks][j][i].M2 += 2.0*(qshear-2.)*(M1e + om_dt*dM2e)*fact;
01459 #ifndef FARGO
01460 pG->U[ks][j][i].M2 -=0.5*qshear*om_dt*(x1Flux[j][i].d+x1Flux[j][i+1].d);
01461 #endif
01462 }
01463 if (ShBoxCoord==xz){
01464 pG->U[ks][j][i].M1 += (4.0*dM3e + 2.0*(qshear-2.)*om_dt*M1e)*fact;
01465 pG->U[ks][j][i].M3 += 2.0*(qshear-2.)*(M1e + om_dt*dM3e)*fact;
01466 #ifndef FARGO
01467 pG->U[ks][j][i].M3 -=0.5*qshear*om_dt*(x1Flux[j][i].d+x1Flux[j][i+1].d);
01468 #endif
01469 }
01470
01471
01472
01473
01474 phic = (*ShearingBoxPot)((x1 ),x2,x3);
01475 phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
01476 phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
01477 #ifndef BAROTROPIC
01478 pG->U[ks][j][i].E -= dtodx1*(x1Flux[j][i ].d*(phic - phil) +
01479 x1Flux[j][i+1].d*(phir - phic));
01480 #endif
01481
01482 phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
01483 phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
01484 #ifndef BAROTROPIC
01485 pG->U[ks][j][i].E -= dtodx2*(x2Flux[j ][i].d*(phic - phil) +
01486 x2Flux[j+1][i].d*(phir - phic));
01487 #endif
01488 }
01489 }
01490
01491 #endif
01492
01493 if (StaticGravPot != NULL){
01494 for (j=js; j<=je; j++) {
01495 for (i=is; i<=ie; i++) {
01496 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
01497 phic = (*StaticGravPot)((x1 ),x2,x3);
01498 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
01499 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
01500
01501 #ifdef CYLINDRICAL
01502 g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
01503 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01504 dtodx2 = pG->dt/(r[i]*pG->dx2);
01505 pG->U[ks][j][i].M1 -= pG->dt*dhalf[j][i]*g;
01506 #else
01507 pG->U[ks][j][i].M1 -= dtodx1*dhalf[j][i]*(phir-phil);
01508 #endif
01509
01510 #ifndef BAROTROPIC
01511 pG->U[ks][j][i].E -= dtodx1*(lsf*x1Flux[j][i ].d*(phic - phil) +
01512 rsf*x1Flux[j][i+1].d*(phir - phic));
01513 #endif
01514 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01515 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01516
01517 pG->U[ks][j][i].M2 -= dtodx2*dhalf[j][i]*(phir-phil);
01518
01519 #ifndef BAROTROPIC
01520 pG->U[ks][j][i].E -= dtodx2*(x2Flux[j ][i].d*(phic - phil) +
01521 x2Flux[j+1][i].d*(phir - phic));
01522 #endif
01523 }
01524 }
01525 }
01526
01527
01528
01529
01530
01531
01532
01533
01534 #ifdef SELF_GRAVITY
01535
01536
01537 for (j=js; j<=je; j++){
01538 for (i=is; i<=ie; i++){
01539 phic = pG->Phi[ks][j][i];
01540 phil = 0.5*(pG->Phi[ks][j][i-1] + pG->Phi[ks][j][i ]);
01541 phir = 0.5*(pG->Phi[ks][j][i ] + pG->Phi[ks][j][i+1]);
01542
01543
01544 gxl = (pG->Phi[ks][j][i-1] - pG->Phi[ks][j][i ])/(pG->dx1);
01545 gxr = (pG->Phi[ks][j][i ] - pG->Phi[ks][j][i+1])/(pG->dx1);
01546
01547 gyl = 0.25*((pG->Phi[ks][j-1][i-1] - pG->Phi[ks][j+1][i-1]) +
01548 (pG->Phi[ks][j-1][i ] - pG->Phi[ks][j+1][i ]) )/(pG->dx2);
01549 gyr = 0.25*((pG->Phi[ks][j-1][i ] - pG->Phi[ks][j+1][i ]) +
01550 (pG->Phi[ks][j-1][i+1] - pG->Phi[ks][j+1][i+1]) )/(pG->dx2);
01551
01552
01553 flux_m1l = 0.5*(gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil;
01554 flux_m1r = 0.5*(gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir;
01555
01556 flux_m2l = gxl*gyl/four_pi_G;
01557 flux_m2r = gxr*gyr/four_pi_G;
01558
01559
01560 pG->U[ks][j][i].M1 -= dtodx1*(flux_m1r - flux_m1l);
01561 pG->U[ks][j][i].M2 -= dtodx1*(flux_m2r - flux_m2l);
01562 #ifndef BAROTROPIC
01563 pG->U[ks][j][i].E -= dtodx1*(x1Flux[j][i ].d*(phic - phil) +
01564 x1Flux[j][i+1].d*(phir - phic));
01565 #endif
01566 }
01567 }
01568
01569
01570
01571 for (j=js; j<=je; j++){
01572 for (i=is; i<=ie; i++){
01573 phic = pG->Phi[ks][j][i];
01574 phil = 0.5*(pG->Phi[ks][j-1][i] + pG->Phi[ks][j ][i]);
01575 phir = 0.5*(pG->Phi[ks][j ][i] + pG->Phi[ks][j+1][i]);
01576
01577
01578 gxl = 0.25*((pG->Phi[ks][j-1][i-1] - pG->Phi[ks][j-1][i+1]) +
01579 (pG->Phi[ks][j ][i-1] - pG->Phi[ks][j ][i+1]) )/(pG->dx1);
01580 gxr = 0.25*((pG->Phi[ks][j ][i-1] - pG->Phi[ks][j ][i+1]) +
01581 (pG->Phi[ks][j+1][i-1] - pG->Phi[ks][j+1][i+1]) )/(pG->dx1);
01582
01583 gyl = (pG->Phi[ks][j-1][i] - pG->Phi[ks][j ][i])/(pG->dx2);
01584 gyr = (pG->Phi[ks][j ][i] - pG->Phi[ks][j+1][i])/(pG->dx2);
01585
01586
01587 flux_m1l = gyl*gxl/four_pi_G;
01588 flux_m1r = gyr*gxr/four_pi_G;
01589
01590 flux_m2l = 0.5*(gyl*gyl-gxl*gxl)/four_pi_G + grav_mean_rho*phil;
01591 flux_m2r = 0.5*(gyr*gyr-gxr*gxr)/four_pi_G + grav_mean_rho*phir;
01592
01593
01594 pG->U[ks][j][i].M1 -= dtodx2*(flux_m1r - flux_m1l);
01595 pG->U[ks][j][i].M2 -= dtodx2*(flux_m2r - flux_m2l);
01596 #ifndef BAROTROPIC
01597 pG->U[ks][j][i].E -= dtodx2*(x2Flux[j ][i].d*(phic - phil) +
01598 x2Flux[j+1][i].d*(phir - phic));
01599 #endif
01600 }
01601 }
01602
01603
01604 for (j=js; j<=je+1; j++) {
01605 for (i=is; i<=ie+1; i++) {
01606 pG->x1MassFlux[ks][j][i] = x1Flux[j][i].d;
01607 pG->x2MassFlux[ks][j][i] = x2Flux[j][i].d;
01608 }
01609 }
01610 #endif
01611
01612
01613
01614
01615
01616 #ifndef BAROTROPIC
01617 if (CoolingFunc != NULL){
01618 for (j=js; j<=je; j++){
01619 for (i=is; i<=ie; i++){
01620 coolf = (*CoolingFunc)(dhalf[j][i],phalf[j][i],pG->dt);
01621 pG->U[ks][j][i].E -= pG->dt*coolf;
01622 }
01623 }
01624 }
01625 #endif
01626
01627
01628
01629
01630
01631 #ifdef FEEDBACK
01632 for (j=js; j<=je; j++)
01633 for (i=is; i<=ie; i++) {
01634 pG->U[ks][j][i].M1 -= pG->Coup[ks][j][i].fb1;
01635 pG->U[ks][j][i].M2 -= pG->Coup[ks][j][i].fb2;
01636 pG->U[ks][j][i].M3 -= pG->Coup[ks][j][i].fb3;
01637 #ifndef BAROTROPIC
01638 pG->U[ks][j][i].E += pG->Coup[ks][j][i].Eloss;
01639 pG->Coup[ks][j][i].Eloss *= dt1;
01640 pG->Eloss[ks][j][i] *= dt1;
01641 #endif
01642 }
01643 #endif
01644
01645
01646
01647
01648
01649
01650
01651 for (j=js; j<=je; j++) {
01652 for (i=is; i<=ie; i++) {
01653 #ifdef CYLINDRICAL
01654 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01655 #endif
01656 pG->U[ks][j][i].d -= dtodx1*(rsf*x1Flux[j][i+1].d - lsf*x1Flux[j][i].d );
01657 pG->U[ks][j][i].M1 -= dtodx1*(rsf*x1Flux[j][i+1].Mx - lsf*x1Flux[j][i].Mx);
01658 pG->U[ks][j][i].M2 -= dtodx1*(SQR(rsf)*x1Flux[j][i+1].My - SQR(lsf)*x1Flux[j][i].My);
01659 pG->U[ks][j][i].M3 -= dtodx1*(rsf*x1Flux[j][i+1].Mz - lsf*x1Flux[j][i].Mz);
01660 #ifndef BAROTROPIC
01661 pG->U[ks][j][i].E -= dtodx1*(rsf*x1Flux[j][i+1].E - lsf*x1Flux[j][i].E );
01662 #endif
01663 #ifdef MHD
01664 pG->U[ks][j][i].B2c -= dtodx1*(x1Flux[j][i+1].By - x1Flux[j][i].By);
01665 pG->U[ks][j][i].B3c -= dtodx1*(rsf*x1Flux[j][i+1].Bz - lsf*x1Flux[j][i].Bz);
01666 #endif
01667 #if (NSCALARS > 0)
01668 for (n=0; n<NSCALARS; n++)
01669 pG->U[ks][j][i].s[n] -= dtodx1*(rsf*x1Flux[j][i+1].s[n]
01670 - lsf*x1Flux[j][i ].s[n]);
01671 #endif
01672 }
01673 }
01674
01675
01676
01677
01678
01679 for (j=js; j<=je; j++) {
01680 for (i=is; i<=ie; i++) {
01681 #ifdef CYLINDRICAL
01682 dtodx2 = pG->dt/(r[i]*pG->dx2);
01683 #endif
01684 pG->U[ks][j][i].d -= dtodx2*(x2Flux[j+1][i].d - x2Flux[j][i].d );
01685 pG->U[ks][j][i].M1 -= dtodx2*(x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
01686 pG->U[ks][j][i].M2 -= dtodx2*(x2Flux[j+1][i].Mx - x2Flux[j][i].Mx);
01687 pG->U[ks][j][i].M3 -= dtodx2*(x2Flux[j+1][i].My - x2Flux[j][i].My);
01688 #ifndef BAROTROPIC
01689 pG->U[ks][j][i].E -= dtodx2*(x2Flux[j+1][i].E - x2Flux[j][i].E );
01690 #endif
01691 #ifdef MHD
01692 pG->U[ks][j][i].B3c -= dtodx2*(x2Flux[j+1][i].By - x2Flux[j][i].By);
01693 pG->U[ks][j][i].B1c -= dtodx2*(x2Flux[j+1][i].Bz - x2Flux[j][i].Bz);
01694 #endif
01695 #if (NSCALARS > 0)
01696 for (n=0; n<NSCALARS; n++)
01697 pG->U[ks][j][i].s[n] -= dtodx2*(x2Flux[j+1][i].s[n]
01698 - x2Flux[j ][i].s[n]);
01699 #endif
01700 }
01701 }
01702
01703
01704
01705
01706
01707
01708 #ifdef MHD
01709 for (j=js; j<=je; j++) {
01710 for (i=is; i<=ie; i++) {
01711 #ifdef CYLINDRICAL
01712 rsf = ri[i+1]/r[i]; lsf = ri[i]/r[i];
01713 #endif
01714 pG->U[ks][j][i].B1c =0.5*(lsf*pG->B1i[ks][j][i] + rsf*pG->B1i[ks][j][i+1]);
01715 pG->U[ks][j][i].B2c =0.5*( pG->B2i[ks][j][i] + pG->B2i[ks][j+1][i]);
01716
01717 pG->B3i[ks][j][i] = pG->U[ks][j][i].B3c;
01718 }
01719 }
01720 #endif
01721
01722 #ifdef STATIC_MESH_REFINEMENT
01723
01724
01725
01726
01727 for (ncg=0; ncg<pG->NCGrid; ncg++) {
01728
01729
01730
01731 for (dim=0; dim<2; dim++){
01732 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
01733
01734 if (dim==0) i = pG->CGrid[ncg].ijks[0];
01735 if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
01736 jcs = pG->CGrid[ncg].ijks[1];
01737 jce = pG->CGrid[ncg].ijke[1];
01738
01739 for (j=jcs, jj=0; j<=jce; j++, jj++){
01740 pG->CGrid[ncg].myFlx[dim][ks][jj].d = x1Flux[j][i].d;
01741 pG->CGrid[ncg].myFlx[dim][ks][jj].M1 = x1Flux[j][i].Mx;
01742 pG->CGrid[ncg].myFlx[dim][ks][jj].M2 = x1Flux[j][i].My;
01743 pG->CGrid[ncg].myFlx[dim][ks][jj].M3 = x1Flux[j][i].Mz;
01744 #ifndef BAROTROPIC
01745 pG->CGrid[ncg].myFlx[dim][ks][jj].E = x1Flux[j][i].E;
01746 #endif
01747 #ifdef MHD
01748 pG->CGrid[ncg].myFlx[dim][ks][jj].B1c = 0.0;
01749 pG->CGrid[ncg].myFlx[dim][ks][jj].B2c = x1Flux[j][i].By;
01750 pG->CGrid[ncg].myFlx[dim][ks][jj].B3c = x1Flux[j][i].Bz;
01751 #endif
01752 #if (NSCALARS > 0)
01753 for (n=0; n<NSCALARS; n++)
01754 pG->CGrid[ncg].myFlx[dim][ks][jj].s[n] = x1Flux[j][i].s[n];
01755 #endif
01756 }
01757 #ifdef MHD
01758 for (j=jcs, jj=0; j<=jce+1; j++, jj++){
01759 pG->CGrid[ncg].myEMF3[dim][ks][jj] = emf3[j][i];
01760 }
01761 #endif
01762 }
01763 }
01764
01765
01766
01767 for (dim=2; dim<4; dim++){
01768 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
01769
01770 ics = pG->CGrid[ncg].ijks[0];
01771 ice = pG->CGrid[ncg].ijke[0];
01772 if (dim==2) j = pG->CGrid[ncg].ijks[1];
01773 if (dim==3) j = pG->CGrid[ncg].ijke[1] + 1;
01774
01775 for (i=ics, ii=0; i<=ice; i++, ii++){
01776 pG->CGrid[ncg].myFlx[dim][ks][ii].d = x2Flux[j][i].d;
01777 pG->CGrid[ncg].myFlx[dim][ks][ii].M1 = x2Flux[j][i].Mz;
01778 pG->CGrid[ncg].myFlx[dim][ks][ii].M2 = x2Flux[j][i].Mx;
01779 pG->CGrid[ncg].myFlx[dim][ks][ii].M3 = x2Flux[j][i].My;
01780 #ifndef BAROTROPIC
01781 pG->CGrid[ncg].myFlx[dim][ks][ii].E = x2Flux[j][i].E;
01782 #endif
01783 #ifdef MHD
01784 pG->CGrid[ncg].myFlx[dim][ks][ii].B1c = x2Flux[j][i].Bz;
01785 pG->CGrid[ncg].myFlx[dim][ks][ii].B2c = 0.0;
01786 pG->CGrid[ncg].myFlx[dim][ks][ii].B3c = x2Flux[j][i].By;
01787 #endif
01788 #if (NSCALARS > 0)
01789 for (n=0; n<NSCALARS; n++)
01790 pG->CGrid[ncg].myFlx[dim][ks][ii].s[n] = x2Flux[j][i].s[n];
01791 #endif
01792 }
01793 #ifdef MHD
01794 for (i=ics, ii=0; i<=ice+1; i++, ii++){
01795 pG->CGrid[ncg].myEMF3[dim][ks][ii] = emf3[j][i];
01796 }
01797 #endif
01798 }
01799 }
01800 }
01801
01802 for (npg=0; npg<pG->NPGrid; npg++) {
01803
01804
01805
01806 for (dim=0; dim<2; dim++){
01807 if (pG->PGrid[npg].myFlx[dim] != NULL) {
01808
01809 if (dim==0) i = pG->PGrid[npg].ijks[0];
01810 if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
01811 jps = pG->PGrid[npg].ijks[1];
01812 jpe = pG->PGrid[npg].ijke[1];
01813
01814 for (j=jps, jj=0; j<=jpe; j++, jj++){
01815 pG->PGrid[npg].myFlx[dim][ks][jj].d = x1Flux[j][i].d;
01816 pG->PGrid[npg].myFlx[dim][ks][jj].M1 = x1Flux[j][i].Mx;
01817 pG->PGrid[npg].myFlx[dim][ks][jj].M2 = x1Flux[j][i].My;
01818 pG->PGrid[npg].myFlx[dim][ks][jj].M3 = x1Flux[j][i].Mz;
01819 #ifndef BAROTROPIC
01820 pG->PGrid[npg].myFlx[dim][ks][jj].E = x1Flux[j][i].E;
01821 #endif
01822 #ifdef MHD
01823 pG->PGrid[npg].myFlx[dim][ks][jj].B1c = 0.0;
01824 pG->PGrid[npg].myFlx[dim][ks][jj].B2c = x1Flux[j][i].By;
01825 pG->PGrid[npg].myFlx[dim][ks][jj].B3c = x1Flux[j][i].Bz;
01826 #endif
01827 #if (NSCALARS > 0)
01828 for (n=0; n<NSCALARS; n++)
01829 pG->PGrid[npg].myFlx[dim][ks][jj].s[n] = x1Flux[j][i].s[n];
01830 #endif
01831 }
01832 #ifdef MHD
01833 for (j=jps, jj=0; j<=jpe+1; j++, jj++){
01834 pG->PGrid[npg].myEMF3[dim][ks][jj] = emf3[j][i];
01835 }
01836 #endif
01837 }
01838 }
01839
01840
01841
01842 for (dim=2; dim<4; dim++){
01843 if (pG->PGrid[npg].myFlx[dim] != NULL) {
01844
01845 ips = pG->PGrid[npg].ijks[0];
01846 ipe = pG->PGrid[npg].ijke[0];
01847 if (dim==2) j = pG->PGrid[npg].ijks[1];
01848 if (dim==3) j = pG->PGrid[npg].ijke[1] + 1;
01849
01850 for (i=ips, ii=0; i<=ipe; i++, ii++){
01851 pG->PGrid[npg].myFlx[dim][ks][ii].d = x2Flux[j][i].d;
01852 pG->PGrid[npg].myFlx[dim][ks][ii].M1 = x2Flux[j][i].Mz;
01853 pG->PGrid[npg].myFlx[dim][ks][ii].M2 = x2Flux[j][i].Mx;
01854 pG->PGrid[npg].myFlx[dim][ks][ii].M3 = x2Flux[j][i].My;
01855 #ifndef BAROTROPIC
01856 pG->PGrid[npg].myFlx[dim][ks][ii].E = x2Flux[j][i].E;
01857 #endif
01858 #ifdef MHD
01859 pG->PGrid[npg].myFlx[dim][ks][ii].B1c = x2Flux[j][i].Bz;
01860 pG->PGrid[npg].myFlx[dim][ks][ii].B2c = 0.0;
01861 pG->PGrid[npg].myFlx[dim][ks][ii].B3c = x2Flux[j][i].By;
01862 #endif
01863 #if (NSCALARS > 0)
01864 for (n=0; n<NSCALARS; n++)
01865 pG->PGrid[npg].myFlx[dim][ks][ii].s[n] = x2Flux[j][i].s[n];
01866 #endif
01867 }
01868 #ifdef MHD
01869 for (i=ips, ii=0; i<=ipe+1; i++, ii++){
01870 pG->PGrid[npg].myEMF3[dim][ks][ii] = emf3[j][i];
01871 }
01872 #endif
01873 }
01874 }
01875 }
01876
01877 #endif
01878
01879 return;
01880 }
01881
01882
01883
01884
01885 void integrate_init_2d(MeshS *pM)
01886 {
01887 int nmax,size1=0,size2=0,nl,nd;
01888
01889
01890 for (nl=0; nl<(pM->NLevels); nl++){
01891 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01892 if (pM->Domain[nl][nd].Grid != NULL) {
01893 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
01894 size1 = pM->Domain[nl][nd].Grid->Nx[0];
01895 }
01896 if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
01897 size2 = pM->Domain[nl][nd].Grid->Nx[1];
01898 }
01899 }
01900 }
01901 }
01902
01903 size1 = size1 + 2*nghost;
01904 size2 = size2 + 2*nghost;
01905 nmax = MAX(size1,size2);
01906
01907 #ifdef MHD
01908 if ((emf3 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01909 goto on_error;
01910 if ((emf3_cc = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01911 goto on_error;
01912 #endif
01913 #ifdef H_CORRECTION
01914 if ((eta1 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01915 goto on_error;
01916 if ((eta2 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01917 goto on_error;
01918 #endif
01919
01920 if ((Bxc = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
01921 if ((Bxi = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
01922
01923 #ifdef MHD
01924 if ((B1_x1Face = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01925 goto on_error;
01926 if ((B2_x2Face = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01927 goto on_error;
01928 #endif
01929
01930 if ((U1d= (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01931 if ((Ul = (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01932 if ((Ur = (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01933 if ((W = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01934 if ((Wl = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01935 if ((Wr = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01936
01937 if ((Ul_x1Face=(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01938 goto on_error;
01939 if ((Ur_x1Face=(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01940 goto on_error;
01941 if ((Ul_x2Face=(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01942 goto on_error;
01943 if ((Ur_x2Face=(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01944 goto on_error;
01945
01946 if ((x1Flux =(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01947 goto on_error;
01948 if ((x2Flux =(Cons1DS**)calloc_2d_array(size2,size1,sizeof(Cons1DS)))==NULL)
01949 goto on_error;
01950
01951 #ifndef CYLINDRICAL
01952 #ifndef MHD
01953 #ifndef PARTICLES
01954 if((StaticGravPot != NULL) || (CoolingFunc != NULL))
01955 #endif
01956 #endif
01957 #endif
01958 {
01959 if ((dhalf = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01960 goto on_error;
01961 if ((phalf = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01962 goto on_error;
01963 }
01964
01965
01966 #ifdef CYLINDRICAL
01967 if ((geom_src = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01968 goto on_error;
01969 #endif
01970
01971 return;
01972
01973 on_error:
01974 integrate_destruct();
01975 ath_error("[integrate_init]: malloc returned a NULL pointer\n");
01976 }
01977
01978
01979
01980
01981 void integrate_destruct_2d(void)
01982 {
01983 #ifdef MHD
01984 if (emf3 != NULL) free_2d_array(emf3);
01985 if (emf3_cc != NULL) free_2d_array(emf3_cc);
01986 #endif
01987 #ifdef H_CORRECTION
01988 if (eta1 != NULL) free_2d_array(eta1);
01989 if (eta2 != NULL) free_2d_array(eta2);
01990 #endif
01991 if (Bxc != NULL) free(Bxc);
01992 if (Bxi != NULL) free(Bxi);
01993 #ifdef MHD
01994 if (B1_x1Face != NULL) free_2d_array(B1_x1Face);
01995 if (B2_x2Face != NULL) free_2d_array(B2_x2Face);
01996 #endif
01997
01998 if (U1d != NULL) free(U1d);
01999 if (Ul != NULL) free(Ul);
02000 if (Ur != NULL) free(Ur);
02001 if (W != NULL) free(W);
02002 if (Wl != NULL) free(Wl);
02003 if (Wr != NULL) free(Wr);
02004
02005 if (Ul_x1Face != NULL) free_2d_array(Ul_x1Face);
02006 if (Ur_x1Face != NULL) free_2d_array(Ur_x1Face);
02007 if (Ul_x2Face != NULL) free_2d_array(Ul_x2Face);
02008 if (Ur_x2Face != NULL) free_2d_array(Ur_x2Face);
02009 if (x1Flux != NULL) free_2d_array(x1Flux);
02010 if (x2Flux != NULL) free_2d_array(x2Flux);
02011 if (dhalf != NULL) free_2d_array(dhalf);
02012 if (phalf != NULL) free_2d_array(phalf);
02013
02014
02015 #ifdef CYLINDRICAL
02016 if (geom_src != NULL) free_2d_array(geom_src);
02017 #endif
02018
02019 return;
02020 }
02021
02022
02023
02024
02025
02026
02027 #ifdef MHD
02028 static void integrate_emf3_corner(GridS *pG)
02029 {
02030 int i,is,ie,j,js,je;
02031 Real emf_l1, emf_r1, emf_l2, emf_r2;
02032 Real rsf=1.0,lsf=1.0;
02033
02034 is = pG->is; ie = pG->ie;
02035 js = pG->js; je = pG->je;
02036
02037
02038 for (j=js-1; j<=je+2; j++) {
02039 for (i=is-1; i<=ie+2; i++) {
02040 #ifdef CYLINDRICAL
02041 rsf = pG->ri[i]/pG->r[i]; lsf = pG->ri[i]/pG->r[i-1];
02042 #endif
02043 if (x1Flux[j-1][i].d > 0.0) {
02044 emf_l2 = -x1Flux[j-1][i].By
02045 + (x2Flux[j][i-1].Bz - emf3_cc[j-1][i-1])*lsf;
02046 }
02047 else if (x1Flux[j-1][i].d < 0.0) {
02048 emf_l2 = -x1Flux[j-1][i].By
02049 + (x2Flux[j][i].Bz - emf3_cc[j-1][i])*rsf;
02050
02051 } else {
02052 emf_l2 = -x1Flux[j-1][i].By
02053 + 0.5*((x2Flux[j][i-1].Bz - emf3_cc[j-1][i-1])*lsf +
02054 (x2Flux[j][i ].Bz - emf3_cc[j-1][i ])*rsf );
02055 }
02056
02057 if (x1Flux[j][i].d > 0.0) {
02058 emf_r2 = -x1Flux[j][i].By
02059 + (x2Flux[j][i-1].Bz - emf3_cc[j][i-1])*lsf;
02060 }
02061 else if (x1Flux[j][i].d < 0.0) {
02062 emf_r2 = -x1Flux[j][i].By
02063 + (x2Flux[j][i].Bz - emf3_cc[j][i])*rsf;
02064
02065 } else {
02066 emf_r2 = -x1Flux[j][i].By
02067 + 0.5*((x2Flux[j][i-1].Bz - emf3_cc[j][i-1])*lsf +
02068 (x2Flux[j][i ].Bz - emf3_cc[j][i ])*rsf );
02069 }
02070
02071 if (x2Flux[j][i-1].d > 0.0) {
02072 emf_l1 = x2Flux[j][i-1].Bz
02073 + (-x1Flux[j-1][i].By - emf3_cc[j-1][i-1]);
02074 }
02075 else if (x2Flux[j][i-1].d < 0.0) {
02076 emf_l1 = x2Flux[j][i-1].Bz
02077 + (-x1Flux[j][i].By - emf3_cc[j][i-1]);
02078 } else {
02079 emf_l1 = x2Flux[j][i-1].Bz
02080 + 0.5*(-x1Flux[j-1][i].By - emf3_cc[j-1][i-1]
02081 -x1Flux[j ][i].By - emf3_cc[j ][i-1] );
02082 }
02083
02084 if (x2Flux[j][i].d > 0.0) {
02085 emf_r1 = x2Flux[j][i].Bz
02086 + (-x1Flux[j-1][i].By - emf3_cc[j-1][i]);
02087 }
02088 else if (x2Flux[j][i].d < 0.0) {
02089 emf_r1 = x2Flux[j][i].Bz
02090 + (-x1Flux[j][i].By - emf3_cc[j][i]);
02091 } else {
02092 emf_r1 = x2Flux[j][i].Bz
02093 + 0.5*(-x1Flux[j-1][i].By - emf3_cc[j-1][i]
02094 -x1Flux[j ][i].By - emf3_cc[j ][i] );
02095 }
02096
02097 emf3[j][i] = 0.25*(emf_l1 + emf_r1 + emf_l2 + emf_r2);
02098 }
02099 }
02100
02101 return;
02102 }
02103 #endif
02104
02105 #endif