Athena: integrators/integrate_1d

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file integrate_1d_vl.c
00004  *  \brief Integrate MHD equations using 1D version of MUSCL-Hancock (VL)
00005  *   integrator.
00006  *
00007  * PURPOSE: Integrate MHD equations using 1D version of MUSCL-Hancock (VL)
00008  *   integrator.  Updates U.[d,M1,M2,M3,E,B2c,B3c,s] in Grid structure.
00009  *   Adds gravitational source terms, self-gravity.
00010  *        
00011  * CONTAINS PUBLIC FUNCTIONS: 
00012  * - integrate_1d_vl()
00013  * - integrate_init_1d()
00014  * - integrate_destruct_1d() */
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 /* The L/R states of primitive variables and fluxes at each cell face */
00032 static Prim1DS *Wl_x1Face=NULL, *Wr_x1Face=NULL;
00033 static Cons1DS *x1Flux=NULL;
00034 
00035 /* 1D scratch vectors used by lr_states and flux functions */
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 /* conserved variables at t^{n+1/2} computed in predict step */
00041 static ConsS *Uhalf=NULL;
00042 
00043 /*=========================== PUBLIC FUNCTIONS ===============================*/
00044 /*----------------------------------------------------------------------------*/
00045 /*! \fn void integrate_1d_vl(DomainS *pD)
00046  *  \brief 1D version of van Leer unsplit integrator for MHD.  
00047  *
00048  *   The numbering of steps follows the numbering in the 3D version.
00049  *   NOT ALL STEPS ARE NEEDED IN 1D.
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 /*=== STEP 1: Compute first-order fluxes at t^{n} in x1-direction ============*/
00077 /* No source terms are needed since there is no temporal evolution */
00078 
00079 /*--- Step 1a ------------------------------------------------------------------
00080  * Load 1D vector of conserved variables;  
00081  * U1d = (d, M1, M2, M3, E, B2c, B3c, s[n])
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 /* BAROTROPIC */
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 /* MHD */
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 /*--- Step 1b ------------------------------------------------------------------
00104  * Compute first-order L/R states */
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 /*--- Step 1c ------------------------------------------------------------------
00119  * No source terms needed */
00120 
00121 /*--- Step 1d ------------------------------------------------------------------
00122  * Compute flux in x1-direction */
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 /*=== STEPS 2-4: Not needed in 1D ===*/
00129 
00130 /*=== STEP 5: Update cell-centered variables to half-timestep ================*/
00131 
00132 /*--- Step 5a ------------------------------------------------------------------
00133  * Update cell-centered variables (including B2c and B3c) to half-timestep
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 /* BAROTROPIC */
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 /* MHD */
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 /*=== STEP 6: Add source terms to predict values at half-timestep ============*/
00155 
00156 /*--- Step 6a ------------------------------------------------------------------
00157  * Add source terms from a static gravitational potential for 0.5*dt to predict
00158  * step.  To improve conservation of total energy, we average the energy
00159  * source term computed at cell faces.
00160  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
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 /*--- Step 6b ------------------------------------------------------------------
00179  * Add source terms for self gravity for 0.5*dt to predict step.
00180  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
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 /* SELF_GRAVITY */
00196 
00197 /*=== STEP 7: Compute second-order L/R x1-interface states ===================*/
00198 
00199 /*--- Step 7a ------------------------------------------------------------------
00200  * Load 1D vector of conserved variables;
00201  * U = (d, M1, M2, M3, E, B2c, B3c, s[n])
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 /* BAROTROPIC */
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 /* MHD */
00217 #if (NSCALARS > 0)
00218     for (n=0; n<NSCALARS; n++) U1d[i].s[n] = Uhalf[i].s[n];
00219 #endif /* NSCALARS */
00220   }
00221 
00222 /*--- Step 7b ------------------------------------------------------------------
00223  * Compute L/R states on x1-interfaces, store into arrays
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 /*=== STEPS 8-9: Not needed in 1D ===*/
00238 
00239 /*=== STEP 10: Compute x1-Flux ===============================================*/
00240 
00241 /*--- Step 10b -----------------------------------------------------------------
00242  * Compute second-order fluxes in x1-direction
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 /*=== STEP 11: Not needed in 1D ===*/
00253         
00254 /*=== STEP 12: Add source terms for a full timestep using n+1/2 states =======*/
00255        
00256 /*--- Step 12a -----------------------------------------------------------------
00257  * Add gravitational source terms due to a Static Potential
00258  * To improve conservation of total energy, we average the energy
00259  * source term computed at cell faces.
00260  *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
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 /*--- Step 12b -----------------------------------------------------------------
00279  * Add gravitational source terms for self-gravity.
00280  * A flux correction using Phi^{n+1} in the main loop is required to make
00281  * the source terms 2nd order: see selfg_flux_correction().
00282  */
00283 
00284 #ifdef SELF_GRAVITY
00285 /* Add fluxes and source terms due to (d/dx1) terms  */
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 /* gx, gy and gz centered at L and R x1-faces */
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 /* momentum fluxes in x1.  2nd term is needed only if Jean's swindle used */
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 /* Update momenta and energy with d/dx1 terms  */
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 /* BAROTROPIC */
00306   }
00307 
00308 /* Save mass fluxes in Grid structure for source term correction in main loop */
00309 
00310   for (i=is; i<=ie+1; i++) {
00311     pG->x1MassFlux[ks][js][i] = x1Flux[i].d;
00312   }
00313 #endif /* SELF_GRAVITY */
00314 
00315 /*=== STEP 13: Update cell-centered values for a full timestep ===============*/
00316 
00317 /*--- Step 13a -----------------------------------------------------------------
00318  * Update cell-centered variables in pG (including B2c and B3c) using x1-Fluxes
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 /* BAROTROPIC */
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 /* For consistency, set B2i and B3i to cell-centered values.  */
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 /* MHD */
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 /*--- Step 13d -----------------------------------------------------------------
00344  * With SMR, store fluxes at boundaries of child and parent grids.
00345  */
00346 
00347 /* x1-boundaries of child Grids (interior to THIS Grid) */
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 /* BAROTROPIC */
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 /* MHD */
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 /* x1-boundaries of parent Grids (at boundaries of THIS Grid)  */
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 /* BAROTROPIC */
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 /* MHD */
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 /* STATIC_MESH_REFINEMENT */
00404 
00405   return;
00406 }
00407 
00408 /*----------------------------------------------------------------------------*/
00409 /*! \fn void integrate_init_1d(MeshS *pM)
00410  *  \brief Allocate temporary integration arrays */
00411 void integrate_init_1d(MeshS *pM)
00412 {
00413   int size1=0,nl,nd;
00414 
00415 /* Cycle over all Grids on this processor to find maximum Nx1 */
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 /*! \fn void integrate_destruct_1d(void)
00453  *  \brief Free temporary integration arrays */
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 /* VL_INTEGRATOR */
00475 
00476 #endif /* SPECIAL_RELATIVITY */
integrators/integrate_1d_vl.c
