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