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 #ifndef 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 *W=NULL, *Wl=NULL, *Wr=NULL;
00038 static Cons1DS *U1d=NULL, *Ul=NULL, *Ur=NULL;
00039
00040
00041 static ConsS *Uhalf=NULL;
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 void integrate_1d_vl(DomainS *pD)
00053 {
00054 GridS *pG=(pD->Grid);
00055 Real dtodx1=pG->dt/pG->dx1, hdtodx1=0.5*pG->dt/pG->dx1;
00056 int i, is = pG->is, ie = pG->ie;
00057 int js = pG->js;
00058 int ks = pG->ks;
00059 Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic;
00060 #if (NSCALARS > 0)
00061 int n;
00062 #endif
00063 #ifdef SELF_GRAVITY
00064 Real gxl,gxr,flx_m1l,flx_m1r;
00065 #endif
00066 #ifdef STATIC_MESH_REFINEMENT
00067 int ncg,npg,dim;
00068 #endif
00069
00070 int il=is-(nghost-1), iu=ie+(nghost-1);
00071
00072 for (i=is-nghost; i<=ie+nghost; i++) {
00073 Uhalf[i] = pG->U[ks][js][i];
00074 }
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 for (i=is-nghost; i<=ie+nghost; i++) {
00085 U1d[i].d = pG->U[ks][js][i].d;
00086 U1d[i].Mx = pG->U[ks][js][i].M1;
00087 U1d[i].My = pG->U[ks][js][i].M2;
00088 U1d[i].Mz = pG->U[ks][js][i].M3;
00089 #ifndef BAROTROPIC
00090 U1d[i].E = pG->U[ks][js][i].E;
00091 #endif
00092 #ifdef MHD
00093 U1d[i].By = pG->U[ks][js][i].B2c;
00094 U1d[i].Bz = pG->U[ks][js][i].B3c;
00095 Bxc[i] = pG->U[ks][js][i].B1c;
00096 Bxi[i] = pG->B1i[ks][js][i];
00097 #endif
00098 #if (NSCALARS > 0)
00099 for (n=0; n<NSCALARS; n++) U1d[i].s[n] = pG->U[ks][js][i].s[n];
00100 #endif
00101 }
00102
00103
00104
00105
00106 for (i=is-nghost; i<=ie+nghost; i++) {
00107 W[i] = Cons1D_to_Prim1D(&U1d[i],&Bxc[i]);
00108 }
00109
00110 for (i=il; i<=ie+nghost; i++) {
00111 Wl[i] = W[i-1];
00112 Wr[i] = W[i ];
00113
00114 Ul[i] = U1d[i-1];
00115 Ur[i] = U1d[i ];
00116 }
00117
00118
00119
00120
00121
00122
00123
00124 for (i=il; i<=ie+nghost; i++) {
00125 fluxes(Ul[i],Ur[i],Wl[i],Wr[i],Bxi[i],&x1Flux[i]);
00126 }
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 for (i=il; i<=iu; i++) {
00137 Uhalf[i].d -= hdtodx1*(x1Flux[i+1].d - x1Flux[i].d );
00138 Uhalf[i].M1 -= hdtodx1*(x1Flux[i+1].Mx - x1Flux[i].Mx);
00139 Uhalf[i].M2 -= hdtodx1*(x1Flux[i+1].My - x1Flux[i].My);
00140 Uhalf[i].M3 -= hdtodx1*(x1Flux[i+1].Mz - x1Flux[i].Mz);
00141 #ifndef BAROTROPIC
00142 Uhalf[i].E -= hdtodx1*(x1Flux[i+1].E - x1Flux[i].E );
00143 #endif
00144 #ifdef MHD
00145 Uhalf[i].B2c -= hdtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00146 Uhalf[i].B3c -= hdtodx1*(x1Flux[i+1].Bz - x1Flux[i].Bz);
00147 #endif
00148 #if (NSCALARS > 0)
00149 for (n=0; n<NSCALARS; n++)
00150 Uhalf[i].s[n] -= hdtodx1*(x1Flux[i+1].s[n] - x1Flux[i].s[n]);
00151 #endif
00152 }
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 if (StaticGravPot != NULL){
00164 for (i=il; i<=iu; i++) {
00165 cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00166 phic = (*StaticGravPot)((x1 ),x2,x3);
00167 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00168 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00169
00170 Uhalf[i].M1 -= hdtodx1*pG->U[ks][js][i].d*(phir-phil);
00171 #ifndef BAROTROPIC
00172 Uhalf[i].E -= hdtodx1*(x1Flux[i ].d*(phic - phil) +
00173 x1Flux[i+1].d*(phir - phic));
00174 #endif
00175 }
00176 }
00177
00178
00179
00180
00181
00182
00183 #ifdef SELF_GRAVITY
00184 for (i=il; i<=iu; i++) {
00185 phic = pG->Phi[ks][js][i];
00186 phir = 0.5*(pG->Phi[ks][js][i] + pG->Phi[ks][js][i+1]);
00187 phil = 0.5*(pG->Phi[ks][js][i] + pG->Phi[ks][js][i-1]);
00188
00189 Uhalf[i].M1 -= hdtodx1*pG->U[ks][js][i].d*(phir-phil);
00190 #ifndef BAROTROPIC
00191 Uhalf[i].E -= hdtodx1*(x1Flux[i ].d*(phic - phil) +
00192 x1Flux[i+1].d*(phir - phic));
00193 #endif
00194 }
00195 #endif
00196
00197
00198
00199
00200
00201
00202
00203
00204 for (i=il; i<=iu; i++) {
00205 U1d[i].d = Uhalf[i].d;
00206 U1d[i].Mx = Uhalf[i].M1;
00207 U1d[i].My = Uhalf[i].M2;
00208 U1d[i].Mz = Uhalf[i].M3;
00209 #ifndef BAROTROPIC
00210 U1d[i].E = Uhalf[i].E;
00211 #endif
00212 #ifdef MHD
00213 U1d[i].By = Uhalf[i].B2c;
00214 U1d[i].Bz = Uhalf[i].B3c;
00215 Bxc[i] = Uhalf[i].B1c;
00216 #endif
00217 #if (NSCALARS > 0)
00218 for (n=0; n<NSCALARS; n++) U1d[i].s[n] = Uhalf[i].s[n];
00219 #endif
00220 }
00221
00222
00223
00224
00225
00226 for (i=il; i<=iu; i++) {
00227 W[i] = Cons1D_to_Prim1D(&U1d[i],&Bxc[i]);
00228 }
00229
00230 lr_states(pG,W,Bxc,pG->dt,pG->dx1,is,ie,Wl,Wr,1);
00231
00232 for (i=is; i<=ie+1; i++) {
00233 Wl_x1Face[i] = Wl[i];
00234 Wr_x1Face[i] = Wr[i];
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 for (i=is; i<=ie+1; i++) {
00246 Ul[i] = Prim1D_to_Cons1D(&Wl_x1Face[i],&Bxi[i]);
00247 Ur[i] = Prim1D_to_Cons1D(&Wr_x1Face[i],&Bxi[i]);
00248
00249 fluxes(Ul[i],Ur[i],Wl_x1Face[i],Wr_x1Face[i],Bxi[i],&x1Flux[i]);
00250 }
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263 if (StaticGravPot != NULL){
00264 for (i=is; i<=ie; i++) {
00265 cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00266 phic = (*StaticGravPot)((x1 ),x2,x3);
00267 phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00268 phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00269
00270 pG->U[ks][js][i].M1 -= dtodx1*Uhalf[i].d*(phir-phil);
00271 #ifndef BAROTROPIC
00272 pG->U[ks][js][i].E -= dtodx1*(x1Flux[i ].d*(phic - phil) +
00273 x1Flux[i+1].d*(phir - phic));
00274 #endif
00275 }
00276 }
00277
00278
00279
00280
00281
00282
00283
00284 #ifdef SELF_GRAVITY
00285
00286
00287 for (i=is; i<=ie; i++){
00288 phic = pG->Phi[ks][js][i];
00289 phil = 0.5*(pG->Phi[ks][js][i-1] + pG->Phi[ks][js][i ]);
00290 phir = 0.5*(pG->Phi[ks][js][i ] + pG->Phi[ks][js][i+1]);
00291
00292
00293 gxl = (pG->Phi[ks][js][i-1] - pG->Phi[ks][js][i ])/(pG->dx1);
00294 gxr = (pG->Phi[ks][js][i ] - pG->Phi[ks][js][i+1])/(pG->dx1);
00295
00296
00297 flx_m1l = 0.5*(gxl*gxl)/four_pi_G + grav_mean_rho*phil;
00298 flx_m1r = 0.5*(gxr*gxr)/four_pi_G + grav_mean_rho*phir;
00299
00300
00301 pG->U[ks][js][i].M1 -= dtodx1*(flx_m1r - flx_m1l);
00302 #ifndef BAROTROPIC
00303 pG->U[ks][js][i].E -= dtodx1*(x1Flux[i ].d*(phic - phil) +
00304 x1Flux[i+1].d*(phir - phic));
00305 #endif
00306 }
00307
00308
00309
00310 for (i=is; i<=ie+1; i++) {
00311 pG->x1MassFlux[ks][js][i] = x1Flux[i].d;
00312 }
00313 #endif
00314
00315
00316
00317
00318
00319
00320
00321 for (i=is; i<=ie; i++) {
00322 pG->U[ks][js][i].d -= dtodx1*(x1Flux[i+1].d - x1Flux[i].d );
00323 pG->U[ks][js][i].M1 -= dtodx1*(x1Flux[i+1].Mx - x1Flux[i].Mx);
00324 pG->U[ks][js][i].M2 -= dtodx1*(x1Flux[i+1].My - x1Flux[i].My);
00325 pG->U[ks][js][i].M3 -= dtodx1*(x1Flux[i+1].Mz - x1Flux[i].Mz);
00326 #ifndef BAROTROPIC
00327 pG->U[ks][js][i].E -= dtodx1*(x1Flux[i+1].E - x1Flux[i].E );
00328 #endif
00329 #ifdef MHD
00330 pG->U[ks][js][i].B2c -= dtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00331 pG->U[ks][js][i].B3c -= dtodx1*(x1Flux[i+1].Bz - x1Flux[i].Bz);
00332
00333 pG->B2i[ks][js][i] = pG->U[ks][js][i].B2c;
00334 pG->B3i[ks][js][i] = pG->U[ks][js][i].B3c;
00335 #endif
00336 #if (NSCALARS > 0)
00337 for (n=0; n<NSCALARS; n++)
00338 pG->U[ks][js][i].s[n] -= dtodx1*(x1Flux[i+1].s[n] - x1Flux[i].s[n]);
00339 #endif
00340 }
00341
00342 #ifdef STATIC_MESH_REFINEMENT
00343
00344
00345
00346
00347
00348 for (ncg=0; ncg<pG->NCGrid; ncg++) {
00349 for (dim=0; dim<2; dim++){
00350 if (pG->CGrid[ncg].myFlx[dim] != NULL) {
00351
00352 if (dim==0) i = pG->CGrid[ncg].ijks[0];
00353 if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
00354
00355 pG->CGrid[ncg].myFlx[dim][ks][js].d = x1Flux[i].d;
00356 pG->CGrid[ncg].myFlx[dim][ks][js].M1 = x1Flux[i].Mx;
00357 pG->CGrid[ncg].myFlx[dim][ks][js].M2 = x1Flux[i].My;
00358 pG->CGrid[ncg].myFlx[dim][ks][js].M3 = x1Flux[i].Mz;
00359 #ifndef BAROTROPIC
00360 pG->CGrid[ncg].myFlx[dim][ks][js].E = x1Flux[i].E;
00361 #endif
00362 #ifdef MHD
00363 pG->CGrid[ncg].myFlx[dim][ks][js].B1c = 0.0;
00364 pG->CGrid[ncg].myFlx[dim][ks][js].B2c = x1Flux[i].By;
00365 pG->CGrid[ncg].myFlx[dim][ks][js].B3c = x1Flux[i].Bz;
00366 #endif
00367 #if (NSCALARS > 0)
00368 for (n=0; n<NSCALARS; n++)
00369 pG->CGrid[ncg].myFlx[dim][ks][js].s[n] = x1Flux[i].s[n];
00370 #endif
00371 }
00372 }
00373 }
00374
00375
00376 for (npg=0; npg<pG->NPGrid; npg++) {
00377 for (dim=0; dim<2; dim++){
00378 if (pG->PGrid[npg].myFlx[dim] != NULL) {
00379
00380 if (dim==0) i = pG->PGrid[npg].ijks[0];
00381 if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
00382
00383 pG->PGrid[npg].myFlx[dim][ks][js].d = x1Flux[i].d;
00384 pG->PGrid[npg].myFlx[dim][ks][js].M1 = x1Flux[i].Mx;
00385 pG->PGrid[npg].myFlx[dim][ks][js].M2 = x1Flux[i].My;
00386 pG->PGrid[npg].myFlx[dim][ks][js].M3 = x1Flux[i].Mz;
00387 #ifndef BAROTROPIC
00388 pG->PGrid[npg].myFlx[dim][ks][js].E = x1Flux[i].E;
00389 #endif
00390 #ifdef MHD
00391 pG->PGrid[npg].myFlx[dim][ks][js].B1c = 0.0;
00392 pG->PGrid[npg].myFlx[dim][ks][js].B2c = x1Flux[i].By;
00393 pG->PGrid[npg].myFlx[dim][ks][js].B3c = x1Flux[i].Bz;
00394 #endif
00395 #if (NSCALARS > 0)
00396 for (n=0; n<NSCALARS; n++)
00397 pG->PGrid[npg].myFlx[dim][ks][js].s[n] = x1Flux[i].s[n];
00398 #endif
00399 }
00400 }
00401 }
00402
00403 #endif
00404
00405 return;
00406 }
00407
00408
00409
00410
00411 void integrate_init_1d(MeshS *pM)
00412 {
00413 int size1=0,nl,nd;
00414
00415
00416 for (nl=0; nl<(pM->NLevels); nl++){
00417 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00418 if (pM->Domain[nl][nd].Grid != NULL) {
00419 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00420 size1 = pM->Domain[nl][nd].Grid->Nx[0];
00421 }
00422 }
00423 }
00424 }
00425
00426 size1 = size1 + 2*nghost;
00427
00428 if ((Wl_x1Face=(Prim1DS*)malloc(size1*sizeof(Prim1DS))) ==NULL) goto on_error;
00429 if ((Wr_x1Face=(Prim1DS*)malloc(size1*sizeof(Prim1DS))) ==NULL) goto on_error;
00430 if ((x1Flux =(Cons1DS*)malloc(size1*sizeof(Cons1DS))) ==NULL) goto on_error;
00431
00432 if ((Bxc = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00433 if ((Bxi = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00434
00435 if ((U1d= (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00436 if ((Ul = (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00437 if ((Ur = (Cons1DS*)malloc(size1*sizeof(Cons1DS))) == NULL) goto on_error;
00438 if ((W = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00439 if ((Wl = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00440 if ((Wr = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00441
00442 if ((Uhalf = (ConsS*)malloc(size1*sizeof(ConsS)))==NULL) goto on_error;
00443
00444 return;
00445
00446 on_error:
00447 integrate_destruct();
00448 ath_error("[integrate_init]: malloc returned a NULL pointer\n");
00449 }
00450
00451
00452
00453
00454 void integrate_destruct_1d(void)
00455 {
00456 if (Wl_x1Face != NULL) free(Wl_x1Face);
00457 if (Wr_x1Face != NULL) free(Wr_x1Face);
00458 if (x1Flux != NULL) free(x1Flux);
00459
00460 if (Bxc != NULL) free(Bxc);
00461 if (Bxi != NULL) free(Bxi);
00462
00463 if (U1d != NULL) free(U1d);
00464 if (Ul != NULL) free(Ul);
00465 if (Ur != NULL) free(Ur);
00466 if (W != NULL) free(W);
00467 if (Wl != NULL) free(Wl);
00468 if (Wr != NULL) free(Wr);
00469
00470 if (Uhalf != NULL) free(Uhalf);
00471
00472 return;
00473 }
00474 #endif
00475
00476 #endif