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 #include <math.h>
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 #include "../defs.h"
00033 #include "../athena.h"
00034 #include "../globals.h"
00035 #include "prototypes.h"
00036 #include "../prototypes.h"
00037
00038 #ifdef SPECIAL_RELATIVITY
00039
00040 #if defined(VL_INTEGRATOR) && defined(CARTESIAN)
00041
00042 #ifdef MHD
00043 #define USE_ENTROPY_FIX
00044 #endif
00045
00046
00047 static Prim1DS **Wl_x1Face=NULL, **Wr_x1Face=NULL;
00048 static Prim1DS **Wl_x2Face=NULL, **Wr_x2Face=NULL;
00049 static Cons1DS **x1Flux=NULL, **x2Flux=NULL;
00050 #ifdef FIRST_ORDER_FLUX_CORRECTION
00051 static Cons1DS **x1FluxP=NULL, **x2FluxP=NULL;
00052 static Real **emf3P=NULL;
00053 #endif
00054 #ifdef USE_ENTROPY_FIX
00055 static Real **S=NULL,**Shalf=NULL;
00056 static Real **x1FluxS=NULL,**x2FluxS=NULL;
00057 #ifdef FIRST_ORDER_FLUX_CORRECTION
00058 static Real **x1FluxSP=NULL,**x2FluxSP=NULL;
00059 #endif
00060 #endif
00061
00062
00063 #ifdef MHD
00064 static Real **B1_x1Face=NULL, **B2_x2Face=NULL;
00065 static Real **emf3=NULL, **emf3_cc=NULL;
00066 #endif
00067
00068
00069 static Real *Bxc=NULL, *Bxi=NULL;
00070 static Prim1DS *W1d=NULL, *Wl=NULL, *Wr=NULL;
00071 static Cons1DS *U1d=NULL, *Ul=NULL, *Ur=NULL;
00072
00073
00074 static PrimS **W=NULL;
00075
00076
00077 static ConsS **Uhalf=NULL;
00078 static PrimS **Whalf=NULL;
00079
00080
00081 extern Real etah;
00082 #ifdef H_CORRECTION
00083 static Real **eta1=NULL, **eta2=NULL;
00084 #endif
00085
00086
00087
00088
00089
00090
00091 #ifdef MHD
00092 static void integrate_emf3_corner(GridS *pG);
00093 #endif
00094 #ifdef FIRST_ORDER_FLUX_CORRECTION
00095 static void FixCell(GridS *pG, Int3Vect);
00096 #endif
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 void integrate_2d_vl(DomainS *pD)
00107 {
00108 GridS *pG=(pD->Grid);
00109 Real dtodx1=pG->dt/pG->dx1, dtodx2=pG->dt/pG->dx2;
00110 Real hdtodx1 = 0.5*dtodx1, hdtodx2 = 0.5*dtodx2;
00111 int i, is = pG->is, ie = pG->ie;
00112 int j, js = pG->js, je = pG->je;
00113 int ks = pG->ks;
00114 int cart_x1 = 1, cart_x2 = 2, cart_x3 = 3;
00115 Real pmhalf,pmnew;
00116 Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic,Bx;
00117 #if (NSCALARS > 0)
00118 int n;
00119 #endif
00120 #ifdef SELF_GRAVITY
00121 Real gxl,gxr,gyl,gyr,flx_m1l,flx_m1r,flx_m2l,flx_m2r;
00122 #endif
00123 #ifdef H_CORRECTION
00124 Real cfr,cfl,lambdar,lambdal;
00125 #endif
00126 #ifdef STATIC_MESH_REFINEMENT
00127 int ncg,npg,dim;
00128 int ii,ics,ice,jj,jcs,jce,ips,ipe,jps,jpe;
00129 #endif
00130 #ifdef FIRST_ORDER_FLUX_CORRECTION
00131 ConsS Ucheck;
00132 PrimS Wcheck;
00133 int flag_cell=0,negd=0,negP=0,superl=0,NaNFlux=0;
00134 int entropy,final,fail;
00135 Real Vsq;
00136 Int3Vect BadCell;
00137 #endif
00138 int il=is-(nghost-1), iu=ie+(nghost-1);
00139 int jl=js-(nghost-1), ju=je+(nghost-1);
00140
00141
00142 etah = 0.0;
00143
00144 for (j=js-nghost; j<=je+nghost; j++) {
00145 for (i=is-nghost; i<=ie+nghost; i++) {
00146 Uhalf[j][i] = pG->U[ks][j][i];
00147 W[j][i] = Cons_to_Prim(&(pG->U[ks][j][i]));
00148 #ifdef MHD
00149 B1_x1Face[j][i] = pG->B1i[ks][j][i];
00150 B2_x2Face[j][i] = pG->B2i[ks][j][i];
00151 #endif
00152 #ifdef USE_ENTROPY_FIX
00153 S[j][i] = W[j][i].P * pow(W[j][i].d,1.0-Gamma);
00154 S[j][i]*= pG->U[ks][j][i].d / W[j][i].d;
00155 Shalf[j][i] = S[j][i];
00156 #endif
00157 }
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168 for (j=js-nghost; j<=je+nghost; j++) {
00169 for (i=is-nghost; i<=ie+nghost; i++) {
00170 W1d[i].d = W[j][i].d;
00171 W1d[i].Vx = W[j][i].V1;
00172 W1d[i].Vy = W[j][i].V2;
00173 W1d[i].Vz = W[j][i].V3;
00174 W1d[i].P = W[j][i].P;
00175 #ifdef MHD
00176 W1d[i].By = W[j][i].B2c;
00177 W1d[i].Bz = W[j][i].B3c;
00178 Bxc[i] = W[j][i].B1c;
00179 Bxi[i] = pG->B1i[ks][j][i];
00180 #endif
00181 #if (NSCALARS > 0)
00182 for (n=0; n<NSCALARS; n++) W1d[i].s[n] = W[j][i].s[n];
00183 #endif
00184 }
00185
00186
00187
00188 for (i=il; i<=ie+nghost; i++) {
00189 Wl[i] = W1d[i-1];
00190 Wr[i] = W1d[i ];
00191
00192
00193 Ul[i] = Prim1D_to_Cons1D(&Wl[i], &Bxc[i-1]);
00194 Ur[i] = Prim1D_to_Cons1D(&Wr[i], &Bxc[i ]);
00195 }
00196
00197
00198
00199
00200
00201
00202
00203 for (i=il; i<=ie+nghost; i++) {
00204 fluxes(Ul[i],Ur[i],Wl[i],Wr[i],Bxi[i],&x1Flux[j][i]);
00205 #ifdef USE_ENTROPY_FIX
00206 entropy_flux(Ul[i],Ur[i],Wl[i],Wr[i],Bxi[i],&x1FluxS[j][i]);
00207 #endif
00208 }
00209 }
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219 for (i=is-nghost; i<=ie+nghost; i++) {
00220 for (j=js-nghost; j<=je+nghost; j++) {
00221 W1d[j].d = W[j][i].d;
00222 W1d[j].Vx = W[j][i].V2;
00223 W1d[j].Vy = W[j][i].V3;
00224 W1d[j].Vz = W[j][i].V1;
00225 W1d[j].P = W[j][i].P;
00226 #ifdef MHD
00227 W1d[j].By = W[j][i].B3c;
00228 W1d[j].Bz = W[j][i].B1c;
00229 Bxc[j] = W[j][i].B2c;
00230 Bxi[j] = pG->B2i[ks][j][i];
00231 #endif
00232 #if (NSCALARS > 0)
00233 for (n=0; n<NSCALARS; n++) W1d[j].s[n] = W[j][i].s[n];
00234 #endif
00235 }
00236
00237
00238
00239
00240 for (j=jl; j<=je+nghost; j++) {
00241 Wl[j] = W1d[j-1];
00242 Wr[j] = W1d[j ];
00243
00244
00245 Ul[j] = Prim1D_to_Cons1D(&Wl[j], &Bxc[j-1]);
00246 Ur[j] = Prim1D_to_Cons1D(&Wr[j], &Bxc[j ]);
00247 }
00248
00249
00250
00251
00252
00253
00254
00255 for (j=jl; j<=je+nghost; j++) {
00256 fluxes(Ul[j],Ur[j],Wl[j],Wr[j],Bxi[j],&x2Flux[j][i]);
00257 #ifdef USE_ENTROPY_FIX
00258 entropy_flux(Ul[j],Ur[j],Wl[j],Wr[j],Bxi[j],&x2FluxS[j][i]);
00259 #endif
00260 }
00261 }
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272 #ifdef MHD
00273 for (j=js-nghost; j<=je+nghost; j++) {
00274 for (i=is-nghost; i<=ie+nghost; i++) {
00275 emf3_cc[j][i] = (W[j][i].B1c*W[j][i].V2 - W[j][i].B2c*W[j][i].V1);
00276 }
00277 }
00278 integrate_emf3_corner(pG);
00279
00280
00281
00282
00283
00284 for (j=jl; j<=ju; j++) {
00285 for (i=il; i<=iu; i++) {
00286 B1_x1Face[j][i] -= hdtodx2*(emf3[j+1][i ] - emf3[j][i]);
00287 B2_x2Face[j][i] += hdtodx1*(emf3[j ][i+1] - emf3[j][i]);
00288 }
00289 B1_x1Face[j][iu+1] -= hdtodx2*(emf3[j+1][iu+1]-emf3[j][iu+1]);
00290 }
00291 for (i=il; i<=iu; i++) {
00292 B2_x2Face[ju+1][i] += hdtodx1*(emf3[ju+1][i+1]-emf3[ju+1][i]);
00293 }
00294
00295
00296
00297
00298
00299
00300 for (j=jl; j<=ju; j++) {
00301 for (i=il; i<=iu; i++) {
00302 Uhalf[j][i].B1c = 0.5*(B1_x1Face[j][i] + B1_x1Face[j][i+1]);
00303 Uhalf[j][i].B2c = 0.5*(B2_x2Face[j][i] + B2_x2Face[j+1][i]);
00304 }
00305 }
00306 #endif
00307
00308
00309
00310
00311
00312
00313
00314 for (j=jl; j<=ju; j++) {
00315 for (i=il; i<=iu; i++) {
00316 Uhalf[j][i].d -= hdtodx1*(x1Flux[j][i+1].d - x1Flux[j][i].d );
00317 Uhalf[j][i].M1 -= hdtodx1*(x1Flux[j][i+1].Mx - x1Flux[j][i].Mx);
00318 Uhalf[j][i].M2 -= hdtodx1*(x1Flux[j][i+1].My - x1Flux[j][i].My);
00319 Uhalf[j][i].M3 -= hdtodx1*(x1Flux[j][i+1].Mz - x1Flux[j][i].Mz);
00320 #ifndef BAROTROPIC
00321 Uhalf[j][i].E -= hdtodx1*(x1Flux[j][i+1].E - x1Flux[j][i].E );
00322 #endif
00323 #ifdef MHD
00324 Uhalf[j][i].B3c -= hdtodx1*(x1Flux[j][i+1].Bz - x1Flux[j][i].Bz);
00325 #endif
00326 #if (NSCALARS > 0)
00327 for (n=0; n<NSCALARS; n++)
00328 Uhalf[j][i].s[n] -= hdtodx1*(x1Flux[j][i+1].s[n] - x1Flux[j][i].s[n]);
00329 #endif
00330 #ifdef USE_ENTROPY_FIX
00331 Shalf[j][i] -= hdtodx1*(x1FluxS[j][i+1] - x1FluxS[j][i]);
00332 #endif
00333 }
00334 }
00335
00336
00337
00338
00339
00340 for (j=jl; j<=ju; j++) {
00341 for (i=il; i<=iu; i++) {
00342 Uhalf[j][i].d -= hdtodx2*(x2Flux[j+1][i].d - x2Flux[j][i].d );
00343 Uhalf[j][i].M1 -= hdtodx2*(x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
00344 Uhalf[j][i].M2 -= hdtodx2*(x2Flux[j+1][i].Mx - x2Flux[j][i].Mx);
00345 Uhalf[j][i].M3 -= hdtodx2*(x2Flux[j+1][i].My - x2Flux[j][i].My);
00346 #ifndef BAROTROPIC
00347 Uhalf[j][i].E -= hdtodx2*(x2Flux[j+1][i].E - x2Flux[j][i].E );
00348 #endif
00349 #ifdef MHD
00350 Uhalf[j][i].B3c -= hdtodx2*(x2Flux[j+1][i].By - x2Flux[j][i].By);
00351 #endif
00352 #if (NSCALARS > 0)
00353 for (n=0; n<NSCALARS; n++)
00354 Uhalf[j][i].s[n] -= hdtodx2*(x2Flux[j+1][i].s[n] - x2Flux[j][i].s[n]);
00355 #endif
00356 #ifdef USE_ENTROPY_FIX
00357 Shalf[j][i] -= hdtodx2*(x2FluxS[j+1][i] - x2FluxS[j][i] );
00358 #endif
00359 }
00360 }
00361
00362 #ifdef FIRST_ORDER_FLUX_CORRECTION
00363
00364
00365
00366
00367
00368 NaNFlux = 0;
00369 for (j=js; j<=je+1; j++) {
00370 for (i=is; i<=ie+1; i++) {
00371 if ((x1Flux[j][i].d != x1Flux[j][i].d) ||
00372 #ifndef BAROTROPIC
00373 (x1Flux[j][i].E != x1Flux[j][i].E) ||
00374 #endif
00375 #ifdef MHD
00376 (x1Flux[j][i].By != x1Flux[j][i].By) ||
00377 (x1Flux[j][i].Bz != x1Flux[j][i].Bz) ||
00378 #endif
00379 (x1Flux[j][i].Mx != x1Flux[j][i].Mx) ||
00380 (x1Flux[j][i].My != x1Flux[j][i].My) ||
00381 (x1Flux[j][i].Mz != x1Flux[j][i].Mz)) {
00382 NaNFlux++;
00383 }
00384 if ((x2Flux[j][i].d != x2Flux[j][i].d) ||
00385 #ifndef BAROTROPIC
00386 (x2Flux[j][i].E != x2Flux[j][i].E) ||
00387 #endif
00388 #ifdef MHD
00389 (x2Flux[j][i].By != x2Flux[j][i].By) ||
00390 (x2Flux[j][i].Bz != x2Flux[j][i].Bz) ||
00391 #endif
00392 (x2Flux[j][i].Mx != x2Flux[j][i].Mx) ||
00393 (x2Flux[j][i].My != x2Flux[j][i].My) ||
00394 (x2Flux[j][i].Mz != x2Flux[j][i].Mz)) {
00395 NaNFlux++;
00396 }
00397 x1FluxP[j][i] = x1Flux[j][i];
00398 x2FluxP[j][i] = x2Flux[j][i];
00399 #ifdef MHD
00400 emf3P[j][i] = emf3[j][i];
00401 #endif
00402 #ifdef USE_ENTROPY_FIX
00403 x1FluxSP[j][i] = x1FluxS[j][i];
00404 x2FluxSP[j][i] = x2FluxS[j][i];
00405 #endif
00406 }
00407 }
00408 if (NaNFlux != 0) {
00409 ath_error("[Step5d] %i first-order fluxes are NaN!\n",NaNFlux);
00410 }
00411 #endif
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422 if (StaticGravPot != NULL){
00423 for (j=jl; j<=ju; j++) {
00424 for (i=il; i<=iu; i++) {
00425 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00426 phic = (*StaticGravPot)( x1, x2,x3);
00427 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00428 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00429
00430 Uhalf[j][i].M1 -= hdtodx1*(phir-phil)*pG->U[ks][j][i].d;
00431 #ifndef BAROTROPIC
00432 Uhalf[j][i].E -= hdtodx1*(x1Flux[j][i ].d*(phic - phil)
00433 + x1Flux[j][i+1].d*(phir - phic));
00434 #endif
00435 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
00436 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00437
00438 Uhalf[j][i].M2 -= hdtodx2*(phir-phil)*pG->U[ks][j][i].d;
00439 #ifndef BAROTROPIC
00440 Uhalf[j][i].E -= hdtodx2*(x2Flux[j ][i].d*(phic - phil)
00441 + x2Flux[j+1][i].d*(phir - phic));
00442 #endif
00443 }
00444 }
00445 }
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457 #ifdef FIRST_ORDER_FLUX_CORRECTION
00458 negd = 0;
00459 negP = 0;
00460 superl = 0;
00461 flag_cell = 0;
00462 #endif
00463 for (j=js-nghost; j<=je+nghost; j++) {
00464 for (i=is-nghost; i<=ie+nghost; i++) {
00465 Whalf[j][i] = check_Prim(&(Uhalf[j][i]));
00466 #ifdef FIRST_ORDER_FLUX_CORRECTION
00467 if (Whalf[j][i].d < 0.0) {
00468 flag_cell = 1;
00469 negd++;
00470 }
00471 if (Whalf[j][i].P < 0.0) {
00472 flag_cell = 1;
00473 negP++;
00474 }
00475 Vsq = SQR(Whalf[j][i].V1) + SQR(Whalf[j][i].V2) + SQR(Whalf[j][i].V3);
00476 if (Vsq > 1.0) {
00477 flag_cell = 1;
00478 superl++;
00479 }
00480 if (flag_cell != 0) {
00481 #ifdef USE_ENTROPY_FIX
00482 Wcheck = entropy_fix (&(Uhalf[j][i]),&(Shalf[j][i]));
00483 Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00484 if (Wcheck.d > 0.0 && Wcheck.P > 0.0 && Vsq < 1.0){
00485 entropy++;
00486 Whalf[j][i].d = Wcheck.d;
00487 Whalf[j][i].P = Wcheck.P;
00488 Whalf[j][i].V1 = Wcheck.V1;
00489 Whalf[j][i].V2 = Wcheck.V2;
00490 Whalf[j][i].V3 = Wcheck.V3;
00491 flag_cell=0;
00492 } else {
00493 #endif
00494 Whalf[j][i].d = W[j][i].d;
00495 Whalf[j][i].V1 = W[j][i].V1;
00496 Whalf[j][i].V2 = W[j][i].V2;
00497 Whalf[j][i].V3 = W[j][i].V3;
00498 Whalf[j][i].P = W[j][i].P;
00499 flag_cell=0;
00500 #ifdef USE_ENTROPY_FIX
00501 }
00502 #endif
00503 }
00504 #endif
00505 }
00506 }
00507
00508 #ifdef FIRST_ORDER_FLUX_CORRECTION
00509 if (negd > 0 || negP > 0 || superl > 0){
00510 printf("[Step7]: %i cells had d<0; %i cells had P<0\n",negd,negP);
00511 printf("[Step7]: %i cells had v>1 at t_half\n",superl);
00512 #ifdef USE_ENTROPY_FIX
00513 printf("[Step7]: %i cells fixed using entropy at t_half\n",entropy);
00514 #endif
00515 }
00516 #endif
00517
00518
00519
00520
00521
00522
00523
00524
00525 for (j=js-1; j<=je+1; j++) {
00526 for (i=il; i<=iu; i++) {
00527 W1d[i].d = Whalf[j][i].d;
00528 W1d[i].Vx = Whalf[j][i].V1;
00529 W1d[i].Vy = Whalf[j][i].V2;
00530 W1d[i].Vz = Whalf[j][i].V3;
00531 #ifndef BAROTROPIC
00532 W1d[i].P = Whalf[j][i].P;
00533 #endif
00534 #ifdef MHD
00535 W1d[i].By = Whalf[j][i].B2c;
00536 W1d[i].Bz = Whalf[j][i].B3c;
00537 Bxc[i] = Whalf[j][i].B1c;
00538 #endif
00539 #if (NSCALARS > 0)
00540 for (n=0; n<NSCALARS; n++) W1d[i].s[n] = Whalf[j][i].s[n];
00541 #endif
00542 }
00543
00544
00545
00546
00547
00548 lr_states(pG,W1d,Bxc,pG->dt,pG->dx1,is,ie,Wl,Wr,cart_x1);
00549
00550 #ifdef FIRST_ORDER_FLUX_CORRECTION
00551 for (i=il; i<=iu; i++) {
00552 Vsq = SQR(Wl[i].Vx) + SQR(Wl[i].Vy) + SQR(Wl[i].Vz);
00553 if (Vsq > 1.0){
00554 Wl[i] = W1d[i];
00555 Wr[i] = W1d[i];
00556 }
00557 Vsq = SQR(Wr[i].Vx) + SQR(Wr[i].Vy) + SQR(Wr[i].Vz);
00558 if (Vsq > 1.0){
00559 Wl[i] = W1d[i];
00560 Wr[i] = W1d[i];
00561 }
00562 }
00563 #endif
00564
00565 for (i=is; i<=ie+1; i++) {
00566 Wl_x1Face[j][i] = Wl[i];
00567 Wr_x1Face[j][i] = Wr[i];
00568 }
00569 }
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580 for (i=is-1; i<=ie+1; i++) {
00581 for (j=jl; j<=ju; j++) {
00582 W1d[j].d = Whalf[j][i].d;
00583 W1d[j].Vx = Whalf[j][i].V2;
00584 W1d[j].Vy = Whalf[j][i].V3;
00585 W1d[j].Vz = Whalf[j][i].V1;
00586 #ifndef BAROTROPIC
00587 W1d[j].P = Whalf[j][i].P;
00588 #endif
00589 #ifdef MHD
00590 W1d[j].By = Whalf[j][i].B3c;
00591 W1d[j].Bz = Whalf[j][i].B1c;
00592 Bxc[j] = Whalf[j][i].B2c;
00593 #endif
00594 #if (NSCALARS > 0)
00595 for (n=0; n<NSCALARS; n++) W1d[j].s[n] = Whalf[j][i].s[n];
00596 #endif
00597 }
00598
00599
00600
00601
00602
00603 lr_states(pG,W1d,Bxc,pG->dt,pG->dx2,js,je,Wl,Wr,cart_x2);
00604
00605 #ifdef FIRST_ORDER_FLUX_CORRECTION
00606 for (j=jl; j<=ju; j++) {
00607 Vsq = SQR(Wl[j].Vx) + SQR(Wl[j].Vy) + SQR(Wl[j].Vz);
00608 if (Vsq > 1.0){
00609 Wl[j] = W1d[j];
00610 Wr[j] = W1d[j];
00611 }
00612 Vsq = SQR(Wr[j].Vx) + SQR(Wr[j].Vy) + SQR(Wr[j].Vz);
00613 if (Vsq > 1.0){
00614 Wl[j] = W1d[j];
00615 Wr[j] = W1d[j];
00616 }
00617 }
00618 #endif
00619
00620 for (j=js; j<=je+1; j++) {
00621 Wl_x2Face[j][i] = Wl[j];
00622 Wr_x2Face[j][i] = Wr[j];
00623 }
00624 }
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635 #ifdef H_CORRECTION
00636 for (j=js-1; j<=je+1; j++) {
00637 for (i=is-1; i<=iu; i++) {
00638 #ifdef MHD
00639 Bx = B1_x1Face[j][i];
00640 #endif
00641 cfr = cfast(&(Ur_x1Face[j][i]),&Bx);
00642 cfl = cfast(&(Ul_x1Face[j][i]),&Bx);
00643 lambdar = Ur_x1Face[j][i].Mx/Ur_x1Face[j][i].d + cfr;
00644 lambdal = Ul_x1Face[j][i].Mx/Ul_x1Face[j][i].d - cfl;
00645 eta1[j][i] = 0.5*fabs(lambdar - lambdal);
00646 }
00647 }
00648
00649 for (j=js-1; j<=ju; j++) {
00650 for (i=is-1; i<=ie+1; i++) {
00651 #ifdef MHD
00652 Bx = B2_x2Face[j][i];
00653 #endif
00654 cfr = cfast(&(Ur_x2Face[j][i]),&Bx);
00655 cfl = cfast(&(Ul_x2Face[j][i]),&Bx);
00656 lambdar = Ur_x2Face[j][i].Mx/Ur_x2Face[j][i].d + cfr;
00657 lambdal = Ul_x2Face[j][i].Mx/Ul_x2Face[j][i].d - cfl;
00658 eta2[j][i] = 0.5*fabs(lambdar - lambdal);
00659 }
00660 }
00661 #endif
00662
00663
00664
00665
00666
00667 for (j=js-1; j<=je+1; j++) {
00668 for (i=is; i<=ie+1; i++) {
00669 #ifdef H_CORRECTION
00670 etah = MAX(eta2[j][i-1],eta2[j][i]);
00671 etah = MAX(etah,eta2[j+1][i-1]);
00672 etah = MAX(etah,eta2[j+1][i ]);
00673 etah = MAX(etah,eta1[j ][i ]);
00674 #endif
00675 #ifdef MHD
00676 Bx = B1_x1Face[j][i];
00677 #endif
00678 Ul[i] = Prim1D_to_Cons1D(&Wl_x1Face[j][i],&Bx);
00679 Ur[i] = Prim1D_to_Cons1D(&Wr_x1Face[j][i],&Bx);
00680
00681 fluxes(Ul[i],Ur[i],Wl_x1Face[j][i],Wr_x1Face[j][i],Bx,&x1Flux[j][i]);
00682
00683 #ifdef USE_ENTROPY_FIX
00684 entropy_flux(Ul[i], Ur[i],
00685 Wl_x1Face[j][i],Wr_x1Face[j][i],
00686 Bx, &x1FluxS[j][i]);
00687 #endif
00688
00689 #ifdef FIRST_ORDER_FLUX_CORRECTION
00690
00691 if ((x1Flux[j][i].d != x1Flux[j][i].d) ||
00692 #ifndef BAROTROPIC
00693 (x1Flux[j][i].E != x1Flux[j][i].E) ||
00694 #endif
00695 #ifdef MHD
00696 (x1Flux[j][i].By != x1Flux[j][i].By) ||
00697 (x1Flux[j][i].Bz != x1Flux[j][i].Bz) ||
00698 #endif
00699 (x1Flux[j][i].Mx != x1Flux[j][i].Mx) ||
00700 (x1Flux[j][i].My != x1Flux[j][i].My) ||
00701 (x1Flux[j][i].Mz != x1Flux[j][i].Mz)) {
00702 x1Flux[j][i] = x1FluxP[j][i];
00703 NaNFlux++;
00704 }
00705 #endif
00706
00707 }
00708 }
00709
00710
00711
00712
00713
00714 for (j=js; j<=je+1; j++) {
00715 for (i=is-1; i<=ie+1; i++) {
00716 #ifdef H_CORRECTION
00717 etah = MAX(eta1[j-1][i],eta1[j][i]);
00718 etah = MAX(etah,eta1[j-1][i+1]);
00719 etah = MAX(etah,eta1[j ][i+1]);
00720 etah = MAX(etah,eta2[j ][i]);
00721 #endif
00722 #ifdef MHD
00723 Bx = B2_x2Face[j][i];
00724 #endif
00725 Ul[i] = Prim1D_to_Cons1D(&Wl_x2Face[j][i],&Bx);
00726 Ur[i] = Prim1D_to_Cons1D(&Wr_x2Face[j][i],&Bx);
00727
00728 fluxes(Ul[i],Ur[i],Wl_x2Face[j][i],Wr_x2Face[j][i],Bx,&x2Flux[j][i]);
00729
00730 #ifdef USE_ENTROPY_FIX
00731 entropy_flux(Ul[i], Ur[i],
00732 Wl_x2Face[j][i],Wr_x2Face[j][i],
00733 Bx, &x2FluxS[j][i]);
00734 #endif
00735
00736 #ifdef FIRST_ORDER_FLUX_CORRECTION
00737
00738 if ((x2Flux[j][i].d != x2Flux[j][i].d) ||
00739 #ifndef BAROTROPIC
00740 (x2Flux[j][i].E != x2Flux[j][i].E) ||
00741 #endif
00742 #ifdef MHD
00743 (x2Flux[j][i].By != x2Flux[j][i].By) ||
00744 (x2Flux[j][i].Bz != x2Flux[j][i].Bz) ||
00745 #endif
00746 (x2Flux[j][i].Mx != x2Flux[j][i].Mx) ||
00747 (x2Flux[j][i].My != x2Flux[j][i].My) ||
00748 (x2Flux[j][i].Mz != x2Flux[j][i].Mz)) {
00749 x2Flux[j][i] = x2FluxP[j][i];
00750 NaNFlux++;
00751 }
00752 #endif
00753
00754 }
00755 }
00756
00757 #ifdef FIRST_ORDER_FLUX_CORRECTION
00758 if (NaNFlux != 0) {
00759 printf("[Step11] %i second-order fluxes replaced\n",NaNFlux);
00760 NaNFlux=0;
00761 }
00762 #endif
00763
00764
00765
00766
00767
00768
00769
00770
00771 #ifdef MHD
00772 for (j=js-1; j<=je+1; j++) {
00773 for (i=is-1; i<=ie+1; i++) {
00774 emf3_cc[j][i] = (Whalf[j][i].B1c*Whalf[j][i].V2 -
00775 Whalf[j][i].B2c*Whalf[j][i].V1);
00776 }
00777 }
00778
00779
00780
00781
00782
00783 integrate_emf3_corner(pG);
00784
00785
00786
00787
00788
00789 for (j=js; j<=je; j++) {
00790 for (i=is; i<=ie; i++) {
00791 pG->B1i[ks][j][i] -= dtodx2*(emf3[j+1][i ] - emf3[j][i]);
00792 pG->B2i[ks][j][i] += dtodx1*(emf3[j ][i+1] - emf3[j][i]);
00793 }
00794 pG->B1i[ks][j][ie+1] -= dtodx2*(emf3[j+1][ie+1] - emf3[j][ie+1]);
00795 }
00796 for (i=is; i<=ie; i++) {
00797 pG->B2i[ks][je+1][i] += dtodx1*(emf3[je+1][i+1] - emf3[je+1][i]);
00798 }
00799
00800
00801
00802
00803
00804 for (j=js; j<=je; j++) {
00805 for (i=is; i<=ie; i++) {
00806 pG->U[ks][j][i].B1c = 0.5*(pG->B1i[ks][j][i]+pG->B1i[ks][j][i+1]);
00807 pG->U[ks][j][i].B2c = 0.5*(pG->B2i[ks][j][i]+pG->B2i[ks][j+1][i]);
00808 }
00809 }
00810 #endif
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821 if (StaticGravPot != NULL){
00822 for (j=js; j<=je; j++) {
00823 for (i=is; i<=ie; i++) {
00824 cc_pos(pG,i,j,ks,&x1,&x2,&x3);
00825 phic = (*StaticGravPot)( x1, x2,x3);
00826 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00827 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00828
00829 pG->U[ks][j][i].M1 -= dtodx1*(phir-phil)*Uhalf[j][i].d;
00830 #ifndef BAROTROPIC
00831 pG->U[ks][j][i].E -= dtodx1*(x1Flux[j][i ].d*(phic - phil)
00832 + x1Flux[j][i+1].d*(phir - phic));
00833 #endif
00834 phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
00835 phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00836
00837 pG->U[ks][j][i].M2 -= dtodx2*(phir-phil)*Uhalf[j][i].d;
00838 #ifndef BAROTROPIC
00839 pG->U[ks][j][i].E -= dtodx2*(x2Flux[j ][i].d*(phic - phil)
00840 + x2Flux[j+1][i].d*(phir - phic));
00841 #endif
00842 }
00843 }
00844 }
00845
00846
00847
00848
00849
00850
00851
00852 for (j=js; j<=je; j++) {
00853 for (i=is; i<=ie; i++) {
00854 pG->U[ks][j][i].d -= dtodx1*(x1Flux[j][i+1].d - x1Flux[j][i].d );
00855 pG->U[ks][j][i].M1 -= dtodx1*(x1Flux[j][i+1].Mx - x1Flux[j][i].Mx);
00856 pG->U[ks][j][i].M2 -= dtodx1*(x1Flux[j][i+1].My - x1Flux[j][i].My);
00857 pG->U[ks][j][i].M3 -= dtodx1*(x1Flux[j][i+1].Mz - x1Flux[j][i].Mz);
00858 #ifndef BAROTROPIC
00859 pG->U[ks][j][i].E -= dtodx1*(x1Flux[j][i+1].E - x1Flux[j][i].E );
00860 #endif
00861 #ifdef MHD
00862 pG->U[ks][j][i].B3c -= dtodx1*(x1Flux[j][i+1].Bz - x1Flux[j][i].Bz);
00863 #endif
00864 #if (NSCALARS > 0)
00865 for (n=0; n<NSCALARS; n++)
00866 pG->U[ks][j][i].s[n] -= dtodx1*(x1Flux[j][i+1].s[n]
00867 - x1Flux[j][i ].s[n]);
00868 #endif
00869 #ifdef USE_ENTROPY_FIX
00870 S[j][i] -= dtodx1*(x1FluxS[j][i+1] - x1FluxS[j][i]);
00871 #endif
00872 }
00873 }
00874
00875
00876
00877
00878
00879 for (j=js; j<=je; j++) {
00880 for (i=is; i<=ie; i++) {
00881 pG->U[ks][j][i].d -= dtodx2*(x2Flux[j+1][i].d - x2Flux[j][i].d );
00882 pG->U[ks][j][i].M1 -= dtodx2*(x2Flux[j+1][i].Mz - x2Flux[j][i].Mz);
00883 pG->U[ks][j][i].M2 -= dtodx2*(x2Flux[j+1][i].Mx - x2Flux[j][i].Mx);
00884 pG->U[ks][j][i].M3 -= dtodx2*(x2Flux[j+1][i].My - x2Flux[j][i].My);
00885 #ifndef BAROTROPIC
00886 pG->U[ks][j][i].E -= dtodx2*(x2Flux[j+1][i].E - x2Flux[j][i].E );
00887 #endif
00888 #ifdef MHD
00889 pG->U[ks][j][i].B3c -= dtodx2*(x2Flux[j+1][i].By - x2Flux[j][i].By);
00890 #endif
00891 #if (NSCALARS > 0)
00892 for (n=0; n<NSCALARS; n++)
00893 pG->U[ks][j][i].s[n] -= dtodx2*(x2Flux[j+1][i].s[n]
00894 - x2Flux[j ][i].s[n]);
00895 #endif
00896 #ifdef USE_ENTROPY_FIX
00897 S[j][i] -= dtodx2*(x2FluxS[j+1][i] - x2FluxS[j][i]);
00898 #endif
00899 }
00900 }
00901
00902 #ifdef FIRST_ORDER_FLUX_CORRECTION
00903
00904
00905
00906
00907
00908
00909
00910 for (j=js; j<=je; j++) {
00911 for (i=is; i<=ie; i++) {
00912 Wcheck = check_Prim(&(pG->U[ks][j][i]));
00913 if (Wcheck.d < 0.0) {
00914 flag_cell = 1;
00915 BadCell.i = i;
00916 BadCell.j = j;
00917 BadCell.k = ks;
00918 negd++;
00919 }
00920 if (Wcheck.P < 0.0) {
00921 flag_cell = 1;
00922 BadCell.i = i;
00923 BadCell.j = j;
00924 BadCell.k = ks;
00925 negP++;
00926 }
00927 Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00928 if (Vsq > 1.0) {
00929 flag_cell = 1;
00930 BadCell.i = i;
00931 BadCell.j = j;
00932 BadCell.k = ks;
00933 superl++;
00934 }
00935 if (flag_cell != 0) {
00936 FixCell(pG, BadCell);
00937 flag_cell=0;
00938 }
00939 }
00940 }
00941
00942 if (negd > 0 || negP > 0 || superl > 0){
00943 printf("[Step15a]: %i cells had d<0; %i cells had P<0;\n",negd,negP);
00944 printf("[Step15a]: %i cells had v>1 at 1st order correction\n",superl);
00945 }
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958 #ifdef MHD
00959 negd = 0;
00960 negP = 0;
00961 superl = 0;
00962 entropy = 0;
00963 final = 0;
00964 fail = 0;
00965 for (j=js; j<=je; j++) {
00966 for (i=is; i<=ie; i++) {
00967 flag_cell=0;
00968 Wcheck = check_Prim(&(pG->U[ks][j][i]));
00969 if (Wcheck.d < 0.0) {
00970 flag_cell = 1;
00971 negd++;
00972 }
00973 if (Wcheck.P < 0.0) {
00974 flag_cell = 1;
00975 negP++;
00976 }
00977 Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00978 if (Vsq > 1.0) {
00979 flag_cell = 1;
00980 superl++;
00981 }
00982 #ifdef USE_ENTROPY_FIX
00983 if (flag_cell != 0) {
00984 entropy++;
00985 Wcheck = entropy_fix (&(pG->U[ks][j][i]),&(S[j][i]));
00986 Ucheck = Prim_to_Cons(&Wcheck);
00987 Wcheck = check_Prim(&Ucheck);
00988 Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00989 if (Wcheck.d > 0.0 && Wcheck.P > 0.0 && Vsq < 1.0){
00990 pG->U[ks][j][i].d = Ucheck.d;
00991 pG->U[ks][j][i].M1 = Ucheck.M1;
00992 pG->U[ks][j][i].M2 = Ucheck.M2;
00993 pG->U[ks][j][i].M3 = Ucheck.M3;
00994 pG->U[ks][j][i].E = Ucheck.E;
00995 flag_cell = 0;
00996 }
00997 }
00998 #endif
00999 if (flag_cell != 0) {
01000 final++;
01001 Wcheck = fix_vsq (&(pG->U[ks][j][i]));
01002 Ucheck = Prim_to_Cons(&Wcheck);
01003 pG->U[ks][j][i].d = Ucheck.d;
01004 pG->U[ks][j][i].M1 = Ucheck.M1;
01005 pG->U[ks][j][i].M2 = Ucheck.M2;
01006 pG->U[ks][j][i].M3 = Ucheck.M3;
01007 pG->U[ks][j][i].E = Ucheck.E;
01008 Wcheck = check_Prim(&(pG->U[ks][j][i]));
01009 Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
01010 if (Wcheck.d < 0.0 || Wcheck.P < 0.0 || Vsq > 1.0){
01011 fail++;
01012 }
01013 }
01014 }
01015 }
01016
01017 if (negd > 0 || negP > 0 || superl > 0) {
01018 printf("[Step15b]: %i cells had d<0; %i cells had P<0;\n",negd,negP);
01019 printf("[Step15b]: %i cells had v>1 after 1st order correction\n",superl);
01020 printf("[Step15b]: %i cells required an entropy fix\n",entropy);
01021 printf("[Step15b]: %i cells required a final fix\n",final);
01022 printf("[Step15b]: %i cells had an unphysical state\n",fail);
01023 }
01024 #endif
01025 #endif
01026
01027 #ifdef STATIC_MESH_REFINEMENT
01028
01029
01030
01031
01032
01033
01034 for (ncg=0; ncg<pG->NCGrid; ncg++) {
01035
01036
01037
01038 for (dim=0; dim<2; dim++){
01039 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
01040
01041 if (dim==0) i = pG->CGrid[ncg].ijks[0];
01042 if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
01043 jcs = pG->CGrid[ncg].ijks[1];
01044 jce = pG->CGrid[ncg].ijke[1];
01045
01046 for (j=jcs, jj=0; j<=jce; j++, jj++){
01047 pG->CGrid[ncg].myFlx[dim][ks][jj].d = x1Flux[j][i].d;
01048 pG->CGrid[ncg].myFlx[dim][ks][jj].M1 = x1Flux[j][i].Mx;
01049 pG->CGrid[ncg].myFlx[dim][ks][jj].M2 = x1Flux[j][i].My;
01050 pG->CGrid[ncg].myFlx[dim][ks][jj].M3 = x1Flux[j][i].Mz;
01051 #ifndef BAROTROPIC
01052 pG->CGrid[ncg].myFlx[dim][ks][jj].E = x1Flux[j][i].E;
01053 #endif
01054 #ifdef MHD
01055 pG->CGrid[ncg].myFlx[dim][ks][jj].B1c = 0.0;
01056 pG->CGrid[ncg].myFlx[dim][ks][jj].B2c = x1Flux[j][i].By;
01057 pG->CGrid[ncg].myFlx[dim][ks][jj].B3c = x1Flux[j][i].Bz;
01058 #endif
01059 #if (NSCALARS > 0)
01060 for (n=0; n<NSCALARS; n++)
01061 pG->CGrid[ncg].myFlx[dim][ks][jj].s[n] = x1Flux[j][i].s[n];
01062 #endif
01063 }
01064 #ifdef MHD
01065 for (j=jcs, jj=0; j<=jce+1; j++, jj++){
01066 pG->CGrid[ncg].myEMF3[dim][ks][jj] = emf3[j][i];
01067 }
01068 #endif
01069 }
01070 }
01071
01072
01073
01074 for (dim=2; dim<4; dim++){
01075 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
01076
01077 ics = pG->CGrid[ncg].ijks[0];
01078 ice = pG->CGrid[ncg].ijke[0];
01079 if (dim==2) j = pG->CGrid[ncg].ijks[1];
01080 if (dim==3) j = pG->CGrid[ncg].ijke[1] + 1;
01081
01082 for (i=ics, ii=0; i<=ice; i++, ii++){
01083 pG->CGrid[ncg].myFlx[dim][ks][ii].d = x2Flux[j][i].d;
01084 pG->CGrid[ncg].myFlx[dim][ks][ii].M1 = x2Flux[j][i].Mz;
01085 pG->CGrid[ncg].myFlx[dim][ks][ii].M2 = x2Flux[j][i].Mx;
01086 pG->CGrid[ncg].myFlx[dim][ks][ii].M3 = x2Flux[j][i].My;
01087 #ifndef BAROTROPIC
01088 pG->CGrid[ncg].myFlx[dim][ks][ii].E = x2Flux[j][i].E;
01089 #endif
01090 #ifdef MHD
01091 pG->CGrid[ncg].myFlx[dim][ks][ii].B1c = x2Flux[j][i].Bz;
01092 pG->CGrid[ncg].myFlx[dim][ks][ii].B2c = 0.0;
01093 pG->CGrid[ncg].myFlx[dim][ks][ii].B3c = x2Flux[j][i].By;
01094 #endif
01095 #if (NSCALARS > 0)
01096 for (n=0; n<NSCALARS; n++)
01097 pG->CGrid[ncg].myFlx[dim][ks][ii].s[n] = x2Flux[j][i].s[n];
01098 #endif
01099 }
01100 #ifdef MHD
01101 for (i=ics, ii=0; i<=ice+1; i++, ii++){
01102 pG->CGrid[ncg].myEMF3[dim][ks][ii] = emf3[j][i];
01103 }
01104 #endif
01105 }
01106 }
01107 }
01108
01109
01110
01111
01112
01113 for (npg=0; npg<pG->NPGrid; npg++) {
01114
01115
01116
01117 for (dim=0; dim<2; dim++){
01118 if (pG->PGrid[npg].myFlx[dim] != NULL) {
01119
01120 if (dim==0) i = pG->PGrid[npg].ijks[0];
01121 if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
01122 jps = pG->PGrid[npg].ijks[1];
01123 jpe = pG->PGrid[npg].ijke[1];
01124
01125 for (j=jps, jj=0; j<=jpe; j++, jj++){
01126 pG->PGrid[npg].myFlx[dim][ks][jj].d = x1Flux[j][i].d;
01127 pG->PGrid[npg].myFlx[dim][ks][jj].M1 = x1Flux[j][i].Mx;
01128 pG->PGrid[npg].myFlx[dim][ks][jj].M2 = x1Flux[j][i].My;
01129 pG->PGrid[npg].myFlx[dim][ks][jj].M3 = x1Flux[j][i].Mz;
01130 #ifndef BAROTROPIC
01131 pG->PGrid[npg].myFlx[dim][ks][jj].E = x1Flux[j][i].E;
01132 #endif
01133 #ifdef MHD
01134 pG->PGrid[npg].myFlx[dim][ks][jj].B1c = 0.0;
01135 pG->PGrid[npg].myFlx[dim][ks][jj].B2c = x1Flux[j][i].By;
01136 pG->PGrid[npg].myFlx[dim][ks][jj].B3c = x1Flux[j][i].Bz;
01137 #endif
01138 #if (NSCALARS > 0)
01139 for (n=0; n<NSCALARS; n++)
01140 pG->PGrid[npg].myFlx[dim][ks][jj].s[n] = x1Flux[j][i].s[n];
01141 #endif
01142 }
01143 #ifdef MHD
01144 for (j=jps, jj=0; j<=jpe+1; j++, jj++){
01145 pG->PGrid[npg].myEMF3[dim][ks][jj] = emf3[j][i];
01146 }
01147 #endif
01148 }
01149 }
01150
01151
01152
01153 for (dim=2; dim<4; dim++){
01154 if (pG->PGrid[npg].myFlx[dim] != NULL) {
01155
01156 ips = pG->PGrid[npg].ijks[0];
01157 ipe = pG->PGrid[npg].ijke[0];
01158 if (dim==2) j = pG->PGrid[npg].ijks[1];
01159 if (dim==3) j = pG->PGrid[npg].ijke[1] + 1;
01160
01161 for (i=ips, ii=0; i<=ipe; i++, ii++){
01162 pG->PGrid[npg].myFlx[dim][ks][ii].d = x2Flux[j][i].d;
01163 pG->PGrid[npg].myFlx[dim][ks][ii].M1 = x2Flux[j][i].Mz;
01164 pG->PGrid[npg].myFlx[dim][ks][ii].M2 = x2Flux[j][i].Mx;
01165 pG->PGrid[npg].myFlx[dim][ks][ii].M3 = x2Flux[j][i].My;
01166 #ifndef BAROTROPIC
01167 pG->PGrid[npg].myFlx[dim][ks][ii].E = x2Flux[j][i].E;
01168 #endif
01169 #ifdef MHD
01170 pG->PGrid[npg].myFlx[dim][ks][ii].B1c = x2Flux[j][i].Bz;
01171 pG->PGrid[npg].myFlx[dim][ks][ii].B2c = 0.0;
01172 pG->PGrid[npg].myFlx[dim][ks][ii].B3c = x2Flux[j][i].By;
01173 #endif
01174 #if (NSCALARS > 0)
01175 for (n=0; n<NSCALARS; n++)
01176 pG->PGrid[npg].myFlx[dim][ks][ii].s[n] = x2Flux[j][i].s[n];
01177 #endif
01178 }
01179 #ifdef MHD
01180 for (i=ips, ii=0; i<=ipe+1; i++, ii++){
01181 pG->PGrid[npg].myEMF3[dim][ks][ii] = emf3[j][i];
01182 }
01183 #endif
01184 }
01185 }
01186 }
01187
01188 #endif
01189
01190 return;
01191 }
01192
01193
01194
01195
01196 void integrate_init_2d(MeshS *pM)
01197 {
01198 int nmax,size1=0,size2=0,nl,nd;
01199
01200
01201 for (nl=0; nl<(pM->NLevels); nl++){
01202 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01203 if (pM->Domain[nl][nd].Grid != NULL) {
01204 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
01205 size1 = pM->Domain[nl][nd].Grid->Nx[0];
01206 }
01207 if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
01208 size2 = pM->Domain[nl][nd].Grid->Nx[1];
01209 }
01210 }
01211 }
01212 }
01213
01214 size1 = size1 + 2*nghost;
01215 size2 = size2 + 2*nghost;
01216 nmax = MAX(size1,size2);
01217
01218 #ifdef MHD
01219 if ((emf3 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01220 goto on_error;
01221 if ((emf3_cc = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01222 goto on_error;
01223 #endif
01224 #ifdef H_CORRECTION
01225 if ((eta1 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01226 goto on_error;
01227 if ((eta2 = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01228 goto on_error;
01229 #endif
01230 if ((Bxc = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
01231 if ((Bxi = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
01232 #ifdef MHD
01233 if ((B1_x1Face = (Real**)calloc_2d_array(size2,size1, sizeof(Real))) == NULL)
01234 goto on_error;
01235 if ((B2_x2Face = (Real**)calloc_2d_array(size2,size1, sizeof(Real))) == NULL)
01236 goto on_error;
01237 #endif
01238
01239 if ((U1d= (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01240 if ((Ul = (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01241 if ((Ur = (Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
01242 if ((W1d= (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01243 if ((Wl = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01244 if ((Wr = (Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
01245
01246 if ((Wl_x1Face = (Prim1DS**)calloc_2d_array(size2,size1, sizeof(Prim1DS)))
01247 == NULL) goto on_error;
01248 if ((Wr_x1Face = (Prim1DS**)calloc_2d_array(size2,size1, sizeof(Prim1DS)))
01249 == NULL) goto on_error;
01250 if ((Wl_x2Face = (Prim1DS**)calloc_2d_array(size2,size1, sizeof(Prim1DS)))
01251 == NULL) goto on_error;
01252 if ((Wr_x2Face = (Prim1DS**)calloc_2d_array(size2,size1, sizeof(Prim1DS)))
01253 == NULL) goto on_error;
01254 if ((x1Flux = (Cons1DS**)calloc_2d_array(size2,size1, sizeof(Cons1DS)))
01255 == NULL) goto on_error;
01256 if ((x2Flux = (Cons1DS**)calloc_2d_array(size2,size1, sizeof(Cons1DS)))
01257 == NULL) goto on_error;
01258 #ifdef FIRST_ORDER_FLUX_CORRECTION
01259 if ((x1FluxP = (Cons1DS**)calloc_2d_array(size2,size1, sizeof(Cons1DS)))
01260 == NULL) goto on_error;
01261 if ((x2FluxP = (Cons1DS**)calloc_2d_array(size2,size1, sizeof(Cons1DS)))
01262 == NULL) goto on_error;
01263 #ifdef MHD
01264 if ((emf3P = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01265 goto on_error;
01266 #endif
01267 #endif
01268 #ifdef USE_ENTROPY_FIX
01269 if ((S = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01270 goto on_error;
01271 if ((Shalf = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01272 goto on_error;
01273 if ((x1FluxS = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01274 goto on_error;
01275 if ((x2FluxS = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01276 goto on_error;
01277 #ifdef FIRST_ORDER_FLUX_CORRECTION
01278 if ((x1FluxSP = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01279 goto on_error;
01280 if ((x2FluxSP = (Real**)calloc_2d_array(size2, size1, sizeof(Real))) == NULL)
01281 goto on_error;
01282 #endif
01283 #endif
01284
01285 if ((Uhalf = (ConsS**)calloc_2d_array(size2,size1,sizeof(ConsS)))==NULL)
01286 goto on_error;
01287 if ((Whalf = (PrimS**)calloc_2d_array(size2,size1,sizeof(PrimS)))==NULL)
01288 goto on_error;
01289 if ((W = (PrimS**)calloc_2d_array(size2,size1,sizeof(PrimS)))==NULL)
01290 goto on_error;
01291
01292 return;
01293
01294 on_error:
01295 integrate_destruct();
01296 ath_error("[integrate_init]: malloc returned a NULL pointer\n");
01297 }
01298
01299
01300
01301
01302 void integrate_destruct_2d(void)
01303 {
01304 #ifdef MHD
01305 if (emf3 != NULL) free_2d_array(emf3);
01306 if (emf3_cc != NULL) free_2d_array(emf3_cc);
01307 #endif
01308 #ifdef H_CORRECTION
01309 if (eta1 != NULL) free_2d_array(eta1);
01310 if (eta2 != NULL) free_2d_array(eta2);
01311 #endif
01312 if (Bxc != NULL) free(Bxc);
01313 if (Bxi != NULL) free(Bxi);
01314 #ifdef MHD
01315 if (B1_x1Face != NULL) free_2d_array(B1_x1Face);
01316 if (B2_x2Face != NULL) free_2d_array(B2_x2Face);
01317 #endif
01318
01319 if (U1d != NULL) free(U1d);
01320 if (Ul != NULL) free(Ul);
01321 if (Ur != NULL) free(Ur);
01322 if (W1d != NULL) free(W1d);
01323 if (Wl != NULL) free(Wl);
01324 if (Wr != NULL) free(Wr);
01325
01326 if (Wl_x1Face != NULL) free_2d_array(Wl_x1Face);
01327 if (Wr_x1Face != NULL) free_2d_array(Wr_x1Face);
01328 if (Wl_x2Face != NULL) free_2d_array(Wl_x2Face);
01329 if (Wr_x2Face != NULL) free_2d_array(Wr_x2Face);
01330 if (x1Flux != NULL) free_2d_array(x1Flux);
01331 if (x2Flux != NULL) free_2d_array(x2Flux);
01332 #ifdef FIRST_ORDER_FLUX_CORRECTION
01333 if (x1FluxP != NULL) free_2d_array(x1FluxP);
01334 if (x2FluxP != NULL) free_2d_array(x2FluxP);
01335 #ifdef MHD
01336 if (emf3P != NULL) free_2d_array(emf3P);
01337 #endif
01338 #endif
01339 #ifdef USE_ENTROPY_FIX
01340 if (S != NULL) free_2d_array(S);
01341 if (Shalf != NULL) free_2d_array(Shalf);
01342 if (x1FluxS != NULL) free_2d_array(x1FluxS);
01343 if (x2FluxS != NULL) free_2d_array(x2FluxS);
01344 #ifdef FIRST_ORDER_FLUX_CORRECTION
01345 if (x1FluxSP != NULL) free_2d_array(x1FluxSP);
01346 if (x2FluxSP != NULL) free_2d_array(x2FluxSP);
01347 #endif
01348 #endif
01349
01350 if (Uhalf != NULL) free_2d_array(Uhalf);
01351 if (Whalf != NULL) free_2d_array(Whalf);
01352 if (W != NULL) free_2d_array(W);
01353
01354 return;
01355 }
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370 #ifdef MHD
01371 static void integrate_emf3_corner(GridS *pG)
01372 {
01373 int i,il,iu,j,jl,ju;
01374 Real emf_l1, emf_r1, emf_l2, emf_r2;
01375
01376 il = pG->is-(nghost-1); iu = pG->ie+(nghost-1);
01377 jl = pG->js-(nghost-1); ju = pG->je+(nghost-1);
01378
01379
01380 for (j=jl; j<=ju+1; j++) {
01381 for (i=il; i<=iu+1; i++) {
01382 if (x1Flux[j-1][i].d > 0.0) {
01383 emf_l2 = -x1Flux[j-1][i].By
01384 + (x2Flux[j][i-1].Bz - emf3_cc[j-1][i-1]);
01385 }
01386 else if (x1Flux[j-1][i].d < 0.0) {
01387 emf_l2 = -x1Flux[j-1][i].By
01388 + (x2Flux[j][i].Bz - emf3_cc[j-1][i]);
01389
01390 } else {
01391 emf_l2 = -x1Flux[j-1][i].By
01392 + 0.5*(x2Flux[j][i-1].Bz - emf3_cc[j-1][i-1] +
01393 x2Flux[j][i ].Bz - emf3_cc[j-1][i ] );
01394 }
01395
01396 if (x1Flux[j][i].d > 0.0) {
01397 emf_r2 = -x1Flux[j][i].By
01398 + (x2Flux[j][i-1].Bz - emf3_cc[j][i-1]);
01399 }
01400 else if (x1Flux[j][i].d < 0.0) {
01401 emf_r2 = -x1Flux[j][i].By
01402 + (x2Flux[j][i].Bz - emf3_cc[j][i]);
01403
01404 } else {
01405 emf_r2 = -x1Flux[j][i].By
01406 + 0.5*(x2Flux[j][i-1].Bz - emf3_cc[j][i-1] +
01407 x2Flux[j][i ].Bz - emf3_cc[j][i ] );
01408 }
01409
01410 if (x2Flux[j][i-1].d > 0.0) {
01411 emf_l1 = x2Flux[j][i-1].Bz
01412 + (-x1Flux[j-1][i].By - emf3_cc[j-1][i-1]);
01413 }
01414 else if (x2Flux[j][i-1].d < 0.0) {
01415 emf_l1 = x2Flux[j][i-1].Bz
01416 + (-x1Flux[j][i].By - emf3_cc[j][i-1]);
01417 } else {
01418 emf_l1 = x2Flux[j][i-1].Bz
01419 + 0.5*(-x1Flux[j-1][i].By - emf3_cc[j-1][i-1]
01420 -x1Flux[j ][i].By - emf3_cc[j ][i-1] );
01421 }
01422
01423 if (x2Flux[j][i].d > 0.0) {
01424 emf_r1 = x2Flux[j][i].Bz
01425 + (-x1Flux[j-1][i].By - emf3_cc[j-1][i]);
01426 }
01427 else if (x2Flux[j][i].d < 0.0) {
01428 emf_r1 = x2Flux[j][i].Bz
01429 + (-x1Flux[j][i].By - emf3_cc[j][i]);
01430 } else {
01431 emf_r1 = x2Flux[j][i].Bz
01432 + 0.5*(-x1Flux[j-1][i].By - emf3_cc[j-1][i]
01433 -x1Flux[j ][i].By - emf3_cc[j ][i] );
01434 }
01435
01436 emf3[j][i] = 0.25*(emf_l1 + emf_r1 + emf_l2 + emf_r2);
01437 }
01438 }
01439
01440 return;
01441 }
01442 #endif
01443
01444 #ifdef FIRST_ORDER_FLUX_CORRECTION
01445
01446
01447
01448
01449 static void FixCell(GridS *pG, Int3Vect ix)
01450 {
01451 int ks=pG->ks;
01452 Real dtodx1=pG->dt/pG->dx1, dtodx2=pG->dt/pG->dx2;
01453 Cons1DS x1FD_i, x1FD_ip1, x2FD_j, x2FD_jp1;
01454 #ifdef MHD
01455 int i,j;
01456 Real emf3D_ji, emf3D_jip1, emf3D_jp1i, emf3D_jp1ip1;
01457 #endif
01458 #ifdef USE_ENTROPY_FIX
01459 Real x1SFD_i, x1SFD_ip1, x2SFD_j, x2SFD_jp1;
01460 #endif
01461
01462
01463
01464
01465 x1FD_i.d = x1Flux[ix.j][ix.i].d - x1FluxP[ix.j][ix.i].d;
01466 x2FD_j.d = x2Flux[ix.j][ix.i].d - x2FluxP[ix.j][ix.i].d;
01467
01468 x1FD_ip1.d = x1Flux[ix.j][ix.i+1].d - x1FluxP[ix.j][ix.i+1].d;
01469 x2FD_jp1.d = x2Flux[ix.j+1][ix.i].d - x2FluxP[ix.j+1][ix.i].d;
01470
01471 x1FD_i.Mx = x1Flux[ix.j][ix.i].Mx - x1FluxP[ix.j][ix.i].Mx;
01472 x2FD_j.Mx = x2Flux[ix.j][ix.i].Mx - x2FluxP[ix.j][ix.i].Mx;
01473
01474 x1FD_ip1.Mx = x1Flux[ix.j][ix.i+1].Mx - x1FluxP[ix.j][ix.i+1].Mx;
01475 x2FD_jp1.Mx = x2Flux[ix.j+1][ix.i].Mx - x2FluxP[ix.j+1][ix.i].Mx;
01476
01477 x1FD_i.My = x1Flux[ix.j][ix.i].My - x1FluxP[ix.j][ix.i].My;
01478 x2FD_j.My = x2Flux[ix.j][ix.i].My - x2FluxP[ix.j][ix.i].My;
01479
01480 x1FD_ip1.My = x1Flux[ix.j][ix.i+1].My - x1FluxP[ix.j][ix.i+1].My;
01481 x2FD_jp1.My = x2Flux[ix.j+1][ix.i].My - x2FluxP[ix.j+1][ix.i].My;
01482
01483 x1FD_i.Mz = x1Flux[ix.j][ix.i].Mz - x1FluxP[ix.j][ix.i].Mz;
01484 x2FD_j.Mz = x2Flux[ix.j][ix.i].Mz - x2FluxP[ix.j][ix.i].Mz;
01485
01486 x1FD_ip1.Mz = x1Flux[ix.j][ix.i+1].Mz - x1FluxP[ix.j][ix.i+1].Mz;
01487 x2FD_jp1.Mz = x2Flux[ix.j+1][ix.i].Mz - x2FluxP[ix.j+1][ix.i].Mz;
01488
01489 #ifndef BAROTROPIC
01490 x1FD_i.E = x1Flux[ix.j][ix.i].E - x1FluxP[ix.j][ix.i].E;
01491 x2FD_j.E = x2Flux[ix.j][ix.i].E - x2FluxP[ix.j][ix.i].E;
01492
01493 x1FD_ip1.E = x1Flux[ix.j][ix.i+1].E - x1FluxP[ix.j][ix.i+1].E;
01494 x2FD_jp1.E = x2Flux[ix.j+1][ix.i].E - x2FluxP[ix.j+1][ix.i].E;
01495 #endif
01496
01497 #ifdef MHD
01498 x1FD_i.Bz = x1Flux[ix.j][ix.i].Bz - x1FluxP[ix.j][ix.i].Bz;
01499 x2FD_j.By = x2Flux[ix.j][ix.i].By - x2FluxP[ix.j][ix.i].By;
01500
01501 x1FD_ip1.Bz = x1Flux[ix.j][ix.i+1].Bz - x1FluxP[ix.j][ix.i+1].Bz;
01502 x2FD_jp1.By = x2Flux[ix.j+1][ix.i].By - x2FluxP[ix.j+1][ix.i].By;
01503
01504 emf3D_ji = emf3[ix.j ][ix.i ] - emf3P[ix.j ][ix.i ];
01505 emf3D_jip1 = emf3[ix.j ][ix.i+1] - emf3P[ix.j ][ix.i+1];
01506 emf3D_jp1i = emf3[ix.j+1][ix.i ] - emf3P[ix.j+1][ix.i ];
01507 emf3D_jp1ip1 = emf3[ix.j+1][ix.i+1] - emf3P[ix.j+1][ix.i+1];
01508 #endif
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520 pG->U[ks][ix.j][ix.i].d += dtodx1*(x1FD_ip1.d - x1FD_i.d );
01521 pG->U[ks][ix.j][ix.i].M1 += dtodx1*(x1FD_ip1.Mx - x1FD_i.Mx);
01522 pG->U[ks][ix.j][ix.i].M2 += dtodx1*(x1FD_ip1.My - x1FD_i.My);
01523 pG->U[ks][ix.j][ix.i].M3 += dtodx1*(x1FD_ip1.Mz - x1FD_i.Mz);
01524 #ifndef BAROTROPIC
01525 pG->U[ks][ix.j][ix.i].E += dtodx1*(x1FD_ip1.E - x1FD_i.E );
01526 #endif
01527 #ifdef MHD
01528 pG->U[ks][ix.j][ix.i].B3c += dtodx1*(x1FD_ip1.Bz - x1FD_i.Bz);
01529 #endif
01530
01531
01532
01533
01534 pG->U[ks][ix.j][ix.i].d += dtodx2*(x2FD_jp1.d - x2FD_j.d );
01535 pG->U[ks][ix.j][ix.i].M1 += dtodx2*(x2FD_jp1.Mz - x2FD_j.Mz);
01536 pG->U[ks][ix.j][ix.i].M2 += dtodx2*(x2FD_jp1.Mx - x2FD_j.Mx);
01537 pG->U[ks][ix.j][ix.i].M3 += dtodx2*(x2FD_jp1.My - x2FD_j.My);
01538 #ifndef BAROTROPIC
01539 pG->U[ks][ix.j][ix.i].E += dtodx2*(x2FD_jp1.E - x2FD_j.E );
01540 #endif
01541 #ifdef MHD
01542 pG->U[ks][ix.j][ix.i].B3c += dtodx2*(x2FD_jp1.By - x2FD_j.By);
01543 #endif
01544
01545
01546
01547
01548
01549
01550 if (ix.i > pG->is) {
01551 pG->U[ks][ix.j][ix.i-1].d += dtodx1*(x1FD_i.d );
01552 pG->U[ks][ix.j][ix.i-1].M1 += dtodx1*(x1FD_i.Mx);
01553 pG->U[ks][ix.j][ix.i-1].M2 += dtodx1*(x1FD_i.My);
01554 pG->U[ks][ix.j][ix.i-1].M3 += dtodx1*(x1FD_i.Mz);
01555 #ifndef BAROTROPIC
01556 pG->U[ks][ix.j][ix.i-1].E += dtodx1*(x1FD_i.E );
01557 #endif
01558 #ifdef MHD
01559 pG->U[ks][ix.j][ix.i-1].B3c += dtodx1*(x1FD_i.Bz);
01560 #endif
01561
01562
01563
01564 }
01565
01566 if (ix.i < pG->ie) {
01567 pG->U[ks][ix.j][ix.i+1].d -= dtodx1*(x1FD_ip1.d );
01568 pG->U[ks][ix.j][ix.i+1].M1 -= dtodx1*(x1FD_ip1.Mx);
01569 pG->U[ks][ix.j][ix.i+1].M2 -= dtodx1*(x1FD_ip1.My);
01570 pG->U[ks][ix.j][ix.i+1].M3 -= dtodx1*(x1FD_ip1.Mz);
01571 #ifndef BAROTROPIC
01572 pG->U[ks][ix.j][ix.i+1].E -= dtodx1*(x1FD_ip1.E );
01573 #endif
01574 #ifdef MHD
01575 pG->U[ks][ix.j][ix.i+1].B3c -= dtodx1*(x1FD_ip1.Bz);
01576 #endif
01577
01578
01579
01580 }
01581
01582
01583
01584 if (ix.j > pG->js) {
01585 pG->U[ks][ix.j-1][ix.i].d += dtodx2*(x2FD_j.d );
01586 pG->U[ks][ix.j-1][ix.i].M1 += dtodx2*(x2FD_j.Mz);
01587 pG->U[ks][ix.j-1][ix.i].M2 += dtodx2*(x2FD_j.Mx);
01588 pG->U[ks][ix.j-1][ix.i].M3 += dtodx2*(x2FD_j.My);
01589 #ifndef BAROTROPIC
01590 pG->U[ks][ix.j-1][ix.i].E += dtodx2*(x2FD_j.E );
01591 #endif
01592 #ifdef MHD
01593 pG->U[ks][ix.j-1][ix.i].B3c += dtodx2*(x2FD_j.By);
01594 #endif
01595
01596
01597
01598 }
01599
01600 if (ix.j < pG->je) {
01601 pG->U[ks][ix.j+1][ix.i].d -= dtodx2*(x2FD_jp1.d );
01602 pG->U[ks][ix.j+1][ix.i].M1 -= dtodx2*(x2FD_jp1.Mz);
01603 pG->U[ks][ix.j+1][ix.i].M2 -= dtodx2*(x2FD_jp1.Mx);
01604 pG->U[ks][ix.j+1][ix.i].M3 -= dtodx2*(x2FD_jp1.My);
01605 #ifndef BAROTROPIC
01606 pG->U[ks][ix.j+1][ix.i].E -= dtodx2*(x2FD_jp1.E );
01607 #endif
01608 #ifdef MHD
01609 pG->U[ks][ix.j+1][ix.i].B3c -= dtodx2*(x2FD_jp1.By);
01610 #endif
01611
01612
01613
01614 }
01615
01616
01617 #ifdef MHD
01618 pG->B1i[ks][ix.j][ix.i ] += dtodx2*(emf3D_jp1i - emf3D_ji);
01619 pG->B1i[ks][ix.j][ix.i+1] += dtodx2*(emf3D_jp1ip1 - emf3D_jip1);
01620 pG->B2i[ks][ix.j ][ix.i] -= dtodx1*(emf3D_jip1 - emf3D_ji);
01621 pG->B2i[ks][ix.j+1][ix.i] -= dtodx1*(emf3D_jp1ip1 - emf3D_jp1i);
01622
01623
01624 if (ix.j > pG->js) {
01625 pG->B1i[ks][ix.j-1][ix.i ] -= dtodx2*(emf3D_ji);
01626 pG->B1i[ks][ix.j-1][ix.i+1] -= dtodx2*(emf3D_jip1);
01627 }
01628 if (ix.j < pG->je) {
01629 pG->B1i[ks][ix.j+1][ix.i ] += dtodx2*(emf3D_jp1i);
01630 pG->B1i[ks][ix.j+1][ix.i+1] += dtodx2*(emf3D_jp1ip1);
01631 }
01632
01633 if (ix.i > pG->is) {
01634 pG->B2i[ks][ix.j ][ix.i-1] += dtodx1*(emf3D_ji);
01635 pG->B2i[ks][ix.j+1][ix.i-1] += dtodx1*(emf3D_jp1i);
01636 }
01637 if (ix.i < pG->ie) {
01638 pG->B2i[ks][ix.j ][ix.i+1] -= dtodx1*(emf3D_jip1);
01639 pG->B2i[ks][ix.j+1][ix.i+1] -= dtodx1*(emf3D_jp1ip1);
01640 }
01641
01642
01643 for (j=(ix.j-1); j<=(ix.j+1); j++) {
01644 for (i=(ix.i-1); i<=(ix.i+1); i++) {
01645 pG->U[ks][j][i].B1c = 0.5*(pG->B1i[ks][j][i] + pG->B1i[ks][j][i+1]);
01646 pG->U[ks][j][i].B2c = 0.5*(pG->B2i[ks][j][i] + pG->B2i[ks][j+1][i]);
01647 }}
01648 #endif
01649
01650 #ifdef STATIC_MESH_REFINEMENT
01651
01652
01653
01654 x1Flux[ix.j][ix.i] = x1FluxP[ix.j][ix.i];
01655 x2Flux[ix.j][ix.i] = x2FluxP[ix.j][ix.i];
01656 #endif
01657
01658 }
01659 #endif
01660
01661 #endif
01662
01663 #endif