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