00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include <string.h>
00021 #include "../defs.h"
00022 #include "../athena.h"
00023 #include "../globals.h"
00024 #include "prototypes.h"
00025 #include "../prototypes.h"
00026
00027 #ifdef SPECIAL_RELATIVITY
00028
00029 #if defined(VL_INTEGRATOR) && defined(CARTESIAN)
00030
00031
00032 static Prim1DS *Wl_x1Face=NULL, *Wr_x1Face=NULL;
00033 static Cons1DS *x1Flux=NULL;
00034
00035
00036 static Real *Bxc=NULL, *Bxi=NULL;
00037 static Prim1DS *W1d=NULL, *Wl=NULL, *Wr=NULL;
00038 static Cons1DS *U1d=NULL, *Ul=NULL, *Ur=NULL;
00039 #ifdef FIRST_ORDER_FLUX_CORRECTION
00040 static Cons1DS *x1FluxP=NULL;
00041 #endif
00042
00043
00044 static PrimS *W=NULL;
00045
00046
00047 static ConsS *Uhalf=NULL;
00048 static PrimS *Whalf=NULL;
00049
00050 #ifdef FIRST_ORDER_FLUX_CORRECTION
00051 static void FixCell(GridS *pG, Int3Vect);
00052 #endif
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 void integrate_1d_vl(DomainS *pD)
00064 {
00065 GridS *pG=(pD->Grid);
00066 ConsS U;
00067 Real dtodx1=pG->dt/pG->dx1, hdtodx1=0.5*pG->dt/pG->dx1;
00068 int i, is = pG->is, ie = pG->ie;
00069 int js = pG->js;
00070 int ks = pG->ks;
00071 int cart_x1 = 1, cart_x2 = 2, cart_x3 = 3;
00072 Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic;
00073 #if (NSCALARS > 0)
00074 int n;
00075 #endif
00076 #ifdef SELF_GRAVITY
00077 Real gxl,gxr,flx_m1l,flx_m1r;
00078 #endif
00079 #ifdef STATIC_MESH_REFINEMENT
00080 int ncg,npg,dim;
00081 #endif
00082 #ifdef FIRST_ORDER_FLUX_CORRECTION
00083 int flag_cell=0,negd=0,negP=0,superl=0,NaNFlux=0;
00084 int fail=0,final=0;
00085 Real Vsq,Bx;
00086 PrimS Wcheck;
00087 Int3Vect BadCell;
00088 #endif
00089
00090
00091 int il=is-(nghost-1), iu=ie+(nghost-1);
00092
00093 for (i=is-nghost; i<=ie+nghost; i++) {
00094 Uhalf[i] = pG->U[ks][js][i];
00095 W[i] = Cons_to_Prim(&(pG->U[ks][js][i]));
00096 }
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 for (i=is-nghost; i<=ie+nghost; i++) {
00107 W1d[i].d = W[i].d;
00108 W1d[i].Vx = W[i].V1;
00109 W1d[i].Vy = W[i].V2;
00110 W1d[i].Vz = W[i].V3;
00111 W1d[i].P = W[i].P;
00112 #ifdef MHD
00113 W1d[i].By = W[i].B2c;
00114 W1d[i].Bz = W[i].B3c;
00115 Bxc[i] = W[i].B1c;
00116 Bxi[i] = pG->B1i[ks][js][i];
00117 #endif
00118 #if (NSCALARS > 0)
00119 for (n=0; n<NSCALARS; n++) W1d[i].s[n] = W[i].s[n];
00120 #endif
00121 }
00122
00123
00124
00125
00126
00127 for (i=is-nghost; i<=ie+nghost; i++) {
00128 U1d[i] = Prim1D_to_Cons1D(&W1d[i],&Bxc[i]);
00129 }
00130
00131 for (i=il; i<=ie+nghost; i++) {
00132 Wl[i] = W1d[i-1];
00133 Wr[i] = W1d[i ];
00134
00135 Ul[i] = U1d[i-1];
00136 Ur[i] = U1d[i ];
00137 }
00138
00139
00140
00141
00142
00143
00144
00145 for (i=il; i<=ie+nghost; i++) {
00146 fluxes(Ul[i],Ur[i],Wl[i],Wr[i],Bxi[i],&x1Flux[i]);
00147 }
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 for (i=il; i<=iu; i++) {
00158 Uhalf[i].d -= hdtodx1*(x1Flux[i+1].d - x1Flux[i].d );
00159 Uhalf[i].M1 -= hdtodx1*(x1Flux[i+1].Mx - x1Flux[i].Mx);
00160 Uhalf[i].M2 -= hdtodx1*(x1Flux[i+1].My - x1Flux[i].My);
00161 Uhalf[i].M3 -= hdtodx1*(x1Flux[i+1].Mz - x1Flux[i].Mz);
00162 Uhalf[i].E -= hdtodx1*(x1Flux[i+1].E - x1Flux[i].E );
00163 #ifdef MHD
00164 Uhalf[i].B2c -= hdtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00165 Uhalf[i].B3c -= hdtodx1*(x1Flux[i+1].Bz - x1Flux[i].Bz);
00166 #endif
00167 #if (NSCALARS > 0)
00168 for (n=0; n<NSCALARS; n++)
00169 Uhalf[i].s[n] -= hdtodx1*(x1Flux[i+1].s[n] - x1Flux[i].s[n]);
00170 #endif
00171 #ifdef FIRST_ORDER_FLUX_CORRECTION
00172 x1FluxP[i] = x1Flux[i];
00173 #endif
00174 }
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185 if (StaticGravPot != NULL){
00186 for (i=il; i<=iu; i++) {
00187 cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00188 phic = (*StaticGravPot)((x1 ),x2,x3);
00189 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00190 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00191
00192 Uhalf[i].M1 -= hdtodx1*pG->U[ks][js][i].d*(phir-phil);
00193 Uhalf[i].E -= hdtodx1*(x1Flux[i ].d*(phic - phil) +
00194 x1Flux[i+1].d*(phir - phic));
00195 }
00196 }
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206 #ifdef FIRST_ORDER_FLUX_CORRECTION
00207 negd = 0;
00208 negP = 0;
00209 superl = 0;
00210 flag_cell = 0;
00211 #endif
00212 for (i=il; i<=iu; i++) {
00213 Whalf[i] = Cons_to_Prim(&Uhalf[i]);
00214 #ifdef FIRST_ORDER_FLUX_CORRECTION
00215 if (Whalf[i].d < 0.0) {
00216 flag_cell = 1;
00217 negd++;
00218 }
00219 if (Whalf[i].P < 0.0) {
00220 flag_cell = 1;
00221 negP++;
00222 }
00223 Vsq = SQR(Whalf[i].V1) + SQR(Whalf[i].V2) + SQR(Whalf[i].V3);
00224 if (Vsq > 1.0) {
00225 flag_cell = 1;
00226 superl++;
00227 }
00228 if (flag_cell != 0) {
00229 Whalf[i].d = W[i].d;
00230 Whalf[i].V1 = W[i].V1;
00231 Whalf[i].V2 = W[i].V2;
00232 Whalf[i].V3 = W[i].V3;
00233 Whalf[i].P = W[i].P;
00234 flag_cell=0;
00235 }
00236 #endif
00237 }
00238
00239 #ifdef FIRST_ORDER_FLUX_CORRECTION
00240 if (negd > 0 || negP > 0 || superl > 0)
00241 printf("[Step7]: %i cells had d<0; %i cells had P<0; %i cells had v>1\n"
00242 ,negd,negP,superl);
00243 #endif
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253 for (i=il; i<=iu; i++) {
00254 W1d[i].d = Whalf[i].d;
00255 W1d[i].Vx = Whalf[i].V1;
00256 W1d[i].Vy = Whalf[i].V2;
00257 W1d[i].Vz = Whalf[i].V3;
00258 W1d[i].P = Whalf[i].P;
00259 #ifdef MHD
00260 W1d[i].By = Whalf[i].B2c;
00261 W1d[i].Bz = Whalf[i].B3c;
00262 Bxc[i] = Whalf[i].B1c;
00263 #endif
00264 #if (NSCALARS > 0)
00265 for (n=0; n<NSCALARS; n++) W1d[i].s[n] = Whalf[i].s[n];
00266 #endif
00267 }
00268
00269
00270
00271
00272
00273 lr_states(pG,W1d,Bxc,pG->dt,pG->dx1,is,ie,Wl,Wr,cart_x1);
00274
00275 #ifdef FIRST_ORDER_FLUX_CORRECTION
00276 for (i=il; i<=iu; i++) {
00277 Vsq = SQR(Wl[i].Vx) + SQR(Wl[i].Vy) + SQR(Wl[i].Vz);
00278 if (Vsq > 1.0){
00279 Wl[i] = W1d[i];
00280 Wr[i] = W1d[i];
00281 }
00282 Vsq = SQR(Wr[i].Vx) + SQR(Wr[i].Vy) + SQR(Wr[i].Vz);
00283 if (Vsq > 1.0){
00284 Wl[i] = W1d[i];
00285 Wr[i] = W1d[i];
00286 }
00287 }
00288 #endif
00289
00290 for (i=is; i<=ie+1; i++) {
00291 Wl_x1Face[i] = Wl[i];
00292 Wr_x1Face[i] = Wr[i];
00293 }
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 for (i=is; i<=ie+1; i++) {
00304 Ul[i] = Prim1D_to_Cons1D(&Wl_x1Face[i],&Bxi[i]);
00305 Ur[i] = Prim1D_to_Cons1D(&Wr_x1Face[i],&Bxi[i]);
00306
00307 fluxes(Ul[i],Ur[i],Wl_x1Face[i],Wr_x1Face[i],Bxi[i],&x1Flux[i]);
00308 }
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321 if (StaticGravPot != NULL){
00322 for (i=is; i<=ie; i++) {
00323 cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00324 phic = (*StaticGravPot)((x1 ),x2,x3);
00325 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00326 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00327
00328 pG->U[ks][js][i].M1 -= dtodx1*Uhalf[i].d*(phir-phil);
00329 #ifndef BAROTROPIC
00330 pG->U[ks][js][i].E -= dtodx1*(x1Flux[i ].d*(phic - phil) +
00331 x1Flux[i+1].d*(phir - phic));
00332 #endif
00333 }
00334 }
00335
00336
00337
00338
00339
00340
00341
00342 for (i=is; i<=ie; i++) {
00343 pG->U[ks][js][i].d -= dtodx1*(x1Flux[i+1].d - x1Flux[i].d );
00344 pG->U[ks][js][i].M1 -= dtodx1*(x1Flux[i+1].Mx - x1Flux[i].Mx);
00345 pG->U[ks][js][i].M2 -= dtodx1*(x1Flux[i+1].My - x1Flux[i].My);
00346 pG->U[ks][js][i].M3 -= dtodx1*(x1Flux[i+1].Mz - x1Flux[i].Mz);
00347 #ifndef BAROTROPIC
00348 pG->U[ks][js][i].E -= dtodx1*(x1Flux[i+1].E - x1Flux[i].E );
00349 #endif
00350 #ifdef MHD
00351 pG->U[ks][js][i].B2c -= dtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00352 pG->U[ks][js][i].B3c -= dtodx1*(x1Flux[i+1].Bz - x1Flux[i].Bz);
00353
00354 pG->B2i[ks][js][i] = pG->U[ks][js][i].B2c;
00355 pG->B3i[ks][js][i] = pG->U[ks][js][i].B3c;
00356 #endif
00357 #if (NSCALARS > 0)
00358 for (n=0; n<NSCALARS; n++)
00359 pG->U[ks][js][i].s[n] -= dtodx1*(x1Flux[i+1].s[n] - x1Flux[i].s[n]);
00360 #endif
00361 }
00362
00363 #ifdef FIRST_ORDER_FLUX_CORRECTION
00364
00365
00366
00367
00368
00369
00370 for (i=is; i<=ie; i++) {
00371 Wcheck = check_Prim(&(pG->U[ks][js][i]));
00372 if (Wcheck.d < 0.0) {
00373 flag_cell = 1;
00374 BadCell.i = i;
00375 BadCell.j = js;
00376 BadCell.k = ks;
00377 negd++;
00378 }
00379 if (Wcheck.P < 0.0) {
00380 flag_cell = 1;
00381 BadCell.i = i;
00382 BadCell.j = js;
00383 BadCell.k = ks;
00384 negP++;
00385 }
00386 Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00387 if (Vsq > 1.0) {
00388 flag_cell = 1;
00389 BadCell.i = i;
00390 BadCell.j = js;
00391 BadCell.k = ks;
00392 superl++;
00393 }
00394 if (flag_cell != 0) {
00395 FixCell(pG, BadCell);
00396 flag_cell=0;
00397 }
00398
00399 }
00400
00401 if (negd > 0 || negP > 0 || superl > 0){
00402 printf("[Step15a]: %i cells had d<0; %i cells had P<0;\n",negd,negP);
00403 printf("[Step15a]: %i cells had v>1 at 1st correction\n",superl);
00404 }
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417 #ifdef MHD
00418 fail = 0;
00419 negd = 0;
00420 negP = 0;
00421 final = 0;
00422 superl = 0;
00423 for (i=is; i<=ie; i++) {
00424 flag_cell=0;
00425 Wcheck = check_Prim(&(pG->U[ks][js][i]));
00426 if (Wcheck.d < 0.0) {
00427 flag_cell = 1;
00428 negd++;
00429 }
00430 if (Wcheck.P < 0.0) {
00431 flag_cell = 1;
00432 negP++;
00433 }
00434 Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00435 if (Vsq > 1.0) {
00436 flag_cell = 1;
00437 superl++;
00438 }
00439 if (flag_cell != 0) {
00440 final++;
00441 Wcheck = fix_vsq (&(pG->U[ks][js][i]));
00442 U = Prim_to_Cons(&Wcheck);
00443 pG->U[ks][js][i].d = U.d;
00444 pG->U[ks][js][i].M1 = U.M1;
00445 pG->U[ks][js][i].M2 = U.M2;
00446 pG->U[ks][js][i].M3 = U.M3;
00447 pG->U[ks][js][i].E = U.E;
00448 Wcheck = check_Prim(&(pG->U[ks][js][i]));
00449 Vsq = SQR(Wcheck.V1) + SQR(Wcheck.V2) + SQR(Wcheck.V3);
00450 if (Wcheck.d < 0.0 || Wcheck.P < 0.0 || Vsq > 1.0){
00451 fail++;
00452 }
00453 }
00454 }
00455
00456 if (negd > 0 || negP > 0 || superl > 0) {
00457 printf("[Step15b]: %i cells had d<0; %i cells had P<0;\n",negd,negP);
00458 printf("[Step15b]: %i cells had v>1 after 1st order correction\n",superl);
00459 printf("[Step15b]: %i cells required a final fix\n",final);
00460 printf("[Step15b]: %i cells had an unphysical state\n",fail);
00461 }
00462 #endif
00463 #endif
00464
00465
00466
00467 #ifdef STATIC_MESH_REFINEMENT
00468
00469
00470
00471
00472
00473 for (ncg=0; ncg<pG->NCGrid; ncg++) {
00474 for (dim=0; dim<2; dim++){
00475 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
00476
00477 if (dim==0) i = pG->CGrid[ncg].ijks[0];
00478 if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
00479
00480 pG->CGrid[ncg].myFlx[dim][ks][js].d = x1Flux[i].d;
00481 pG->CGrid[ncg].myFlx[dim][ks][js].M1 = x1Flux[i].Mx;
00482 pG->CGrid[ncg].myFlx[dim][ks][js].M2 = x1Flux[i].My;
00483 pG->CGrid[ncg].myFlx[dim][ks][js].M3 = x1Flux[i].Mz;
00484 #ifndef BAROTROPIC
00485 pG->CGrid[ncg].myFlx[dim][ks][js].E = x1Flux[i].E;
00486 #endif
00487 #ifdef MHD
00488 pG->CGrid[ncg].myFlx[dim][ks][js].B1c = 0.0;
00489 pG->CGrid[ncg].myFlx[dim][ks][js].B2c = x1Flux[i].By;
00490 pG->CGrid[ncg].myFlx[dim][ks][js].B3c = x1Flux[i].Bz;
00491 #endif
00492 #if (NSCALARS > 0)
00493 for (n=0; n<NSCALARS; n++)
00494 pG->CGrid[ncg].myFlx[dim][ks][js].s[n] = x1Flux[i].s[n];
00495 #endif
00496 }
00497 }
00498 }
00499
00500
00501 for (npg=0; npg<pG->NPGrid; npg++) {
00502 for (dim=0; dim<2; dim++){
00503 if (pG->PGrid[npg].myFlx[dim] != NULL) {
00504
00505 if (dim==0) i = pG->PGrid[npg].ijks[0];
00506 if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
00507
00508 pG->PGrid[npg].myFlx[dim][ks][js].d = x1Flux[i].d;
00509 pG->PGrid[npg].myFlx[dim][ks][js].M1 = x1Flux[i].Mx;
00510 pG->PGrid[npg].myFlx[dim][ks][js].M2 = x1Flux[i].My;
00511 pG->PGrid[npg].myFlx[dim][ks][js].M3 = x1Flux[i].Mz;
00512 #ifndef BAROTROPIC
00513 pG->PGrid[npg].myFlx[dim][ks][js].E = x1Flux[i].E;
00514 #endif
00515 #ifdef MHD
00516 pG->PGrid[npg].myFlx[dim][ks][js].B1c = 0.0;
00517 pG->PGrid[npg].myFlx[dim][ks][js].B2c = x1Flux[i].By;
00518 pG->PGrid[npg].myFlx[dim][ks][js].B3c = x1Flux[i].Bz;
00519 #endif
00520 #if (NSCALARS > 0)
00521 for (n=0; n<NSCALARS; n++)
00522 pG->PGrid[npg].myFlx[dim][ks][js].s[n] = x1Flux[i].s[n];
00523 #endif
00524 }
00525 }
00526 }
00527
00528 #endif
00529
00530 return;
00531 }
00532
00533
00534
00535
00536 void integrate_init_1d(MeshS *pM)
00537 {
00538 int size1=0,nl,nd;
00539
00540
00541 for (nl=0; nl<(pM->NLevels); nl++){
00542 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00543 if (pM->Domain[nl][nd].Grid != NULL) {
00544 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00545 size1 = pM->Domain[nl][nd].Grid->Nx[0];
00546 }
00547 }
00548 }
00549 }
00550
00551 size1 = size1 + 2*nghost;
00552
00553 if ((Wl_x1Face=(Prim1DS*)malloc(size1*sizeof(Prim1DS))) ==NULL) goto on_error;
00554 if ((Wr_x1Face=(Prim1DS*)malloc(size1*sizeof(Prim1DS))) ==NULL) goto on_error;
00555 if ((x1Flux =(Cons1DS*)malloc(size1*sizeof(Cons1DS))) ==NULL) goto on_error;
00556 #ifdef FIRST_ORDER_FLUX_CORRECTION
00557 if ((x1FluxP =(Cons1DS*)malloc(size1*sizeof(Cons1DS))) ==NULL) goto on_error;
00558 #endif
00559
00560 if ((Bxc = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00561 if ((Bxi = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00562
00563 if ((U1d= (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00564 if ((Ul = (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00565 if ((Ur = (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00566 if ((W1d= (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00567 if ((Wl = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00568 if ((Wr = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00569
00570 if ((Uhalf = (ConsS*)malloc(size1*sizeof(ConsS)))==NULL) goto on_error;
00571 if ((Whalf = (PrimS*)malloc(size1*sizeof(PrimS)))==NULL) goto on_error;
00572 if ((W = (PrimS*)malloc(size1*sizeof(PrimS)))==NULL) goto on_error;
00573
00574 return;
00575
00576 on_error:
00577 integrate_destruct();
00578 ath_error("[integrate_init]: malloc returned a NULL pointer\n");
00579 }
00580
00581
00582
00583
00584 void integrate_destruct_1d(void)
00585 {
00586 if (Wl_x1Face != NULL) free(Wl_x1Face);
00587 if (Wr_x1Face != NULL) free(Wr_x1Face);
00588 if (x1Flux != NULL) free(x1Flux);
00589 #ifdef FIRST_ORDER_FLUX_CORRECTION
00590 if (x1FluxP != NULL) free(x1FluxP);
00591 #endif
00592
00593
00594 if (Bxc != NULL) free(Bxc);
00595 if (Bxi != NULL) free(Bxi);
00596
00597 if (U1d != NULL) free(U1d);
00598 if (Ul != NULL) free(Ul);
00599 if (Ur != NULL) free(Ur);
00600 if (W1d != NULL) free(W1d);
00601 if (Wl != NULL) free(Wl);
00602 if (Wr != NULL) free(Wr);
00603
00604 if (Uhalf != NULL) free(Uhalf);
00605 if (Whalf != NULL) free(Whalf);
00606 if (W != NULL) free(W);
00607
00608 return;
00609 }
00610
00611 #ifdef FIRST_ORDER_FLUX_CORRECTION
00612
00613
00614
00615
00616 static void FixCell(GridS *pG, Int3Vect indx)
00617 {
00618 int ks=pG->ks,js=pG->js;
00619 Real dtodx1=pG->dt/pG->dx1;
00620 Cons1DS x1FD_i, x1FD_ip1, x2FD_j, x2FD_jp1;
00621
00622
00623
00624 x1FD_i.d = x1Flux[indx.i ].d - x1FluxP[indx.i ].d;
00625 x1FD_ip1.d = x1Flux[indx.i+1].d - x1FluxP[indx.i+1].d;
00626
00627 x1FD_i.Mx = x1Flux[indx.i ].Mx - x1FluxP[indx.i ].Mx;
00628 x1FD_ip1.Mx = x1Flux[indx.i+1].Mx - x1FluxP[indx.i+1].Mx;
00629
00630 x1FD_i.My = x1Flux[indx.i ].My - x1FluxP[indx.i ].My;
00631 x1FD_ip1.My = x1Flux[indx.i+1].My - x1FluxP[indx.i+1].My;
00632
00633 x1FD_i.Mz = x1Flux[indx.i ].Mz - x1FluxP[indx.i ].Mz;
00634 x1FD_ip1.Mz = x1Flux[indx.i+1].Mz - x1FluxP[indx.i+1].Mz;
00635
00636 #ifdef MHD
00637 x1FD_i.By = x1Flux[indx.i ].By - x1FluxP[indx.i ].By;
00638 x1FD_ip1.By = x1Flux[indx.i+1].By - x1FluxP[indx.i+1].By;
00639
00640 x1FD_i.Bz = x1Flux[indx.i ].Bz - x1FluxP[indx.i ].Bz;
00641 x1FD_ip1.Bz = x1Flux[indx.i+1].Bz - x1FluxP[indx.i+1].Bz;
00642 #endif
00643
00644 #ifndef BAROTROPIC
00645 x1FD_i.E = x1Flux[indx.i ].E - x1FluxP[indx.i ].E;
00646 x1FD_ip1.E = x1Flux[indx.i+1].E - x1FluxP[indx.i+1].E;
00647 #endif
00648
00649
00650 pG->U[ks][js][indx.i].d += dtodx1*(x1FD_ip1.d - x1FD_i.d );
00651 pG->U[ks][js][indx.i].M1 += dtodx1*(x1FD_ip1.Mx - x1FD_i.Mx);
00652 pG->U[ks][js][indx.i].M2 += dtodx1*(x1FD_ip1.My - x1FD_i.My);
00653 pG->U[ks][js][indx.i].M3 += dtodx1*(x1FD_ip1.Mz - x1FD_i.Mz);
00654 #ifdef MHD
00655 pG->U[ks][js][indx.i].B2c += dtodx1*(x1FD_ip1.By - x1FD_i.By);
00656 pG->U[ks][js][indx.i].B3c += dtodx1*(x1FD_ip1.Bz - x1FD_i.Bz);
00657
00658 pG->B2i[ks][js][indx.i] = pG->U[ks][js][indx.i].B2c;
00659 pG->B3i[ks][js][indx.i] = pG->U[ks][js][indx.i].B3c;
00660 #endif
00661 #ifndef BAROTROPIC
00662 pG->U[ks][js][indx.i].E += dtodx1*(x1FD_ip1.E - x1FD_i.E );
00663 #endif
00664
00665
00666
00667 if (indx.i > pG->is) {
00668 pG->U[ks][js][indx.i-1].d += dtodx1*(x1FD_i.d );
00669 pG->U[ks][js][indx.i-1].M1 += dtodx1*(x1FD_i.Mx);
00670 pG->U[ks][js][indx.i-1].M2 += dtodx1*(x1FD_i.My);
00671 pG->U[ks][js][indx.i-1].M3 += dtodx1*(x1FD_i.Mz);
00672 #ifdef MHD
00673 pG->U[ks][js][indx.i-1].B2c += dtodx1*(x1FD_i.By);
00674 pG->U[ks][js][indx.i-1].B3c += dtodx1*(x1FD_i.Bz);
00675
00676 pG->B2i[ks][js][indx.i-1] = pG->U[ks][js][indx.i-1].B2c;
00677 pG->B3i[ks][js][indx.i-1] = pG->U[ks][js][indx.i-1].B3c;
00678 #endif
00679 #ifndef BAROTROPIC
00680 pG->U[ks][js][indx.i-1].E += dtodx1*(x1FD_i.E );
00681 #endif
00682 }
00683
00684 if (indx.i < pG->ie) {
00685 pG->U[ks][js][indx.i+1].d -= dtodx1*(x1FD_ip1.d );
00686 pG->U[ks][js][indx.i+1].M1 -= dtodx1*(x1FD_ip1.Mx);
00687 pG->U[ks][js][indx.i+1].M2 -= dtodx1*(x1FD_ip1.My);
00688 pG->U[ks][js][indx.i+1].M3 -= dtodx1*(x1FD_ip1.Mz);
00689 #ifdef MHD
00690 pG->U[ks][js][indx.i+1].B2c -= dtodx1*(x1FD_ip1.By);
00691 pG->U[ks][js][indx.i+1].B3c -= dtodx1*(x1FD_ip1.Bz);
00692
00693 pG->B2i[ks][js][indx.i+1] = pG->U[ks][js][indx.i+1].B2c;
00694 pG->B3i[ks][js][indx.i+1] = pG->U[ks][js][indx.i+1].B3c;
00695 #endif MHD
00696 #ifndef BAROTROPIC
00697 pG->U[ks][js][indx.i+1].E -= dtodx1*(x1FD_ip1.E );
00698 #endif
00699 }
00700
00701 }
00702 #endif
00703
00704 #endif
00705
00706 #endif