Athena: prob/bubble.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file bubble.c
00004  *  \brief Problem generator for bubble in spherical isothermal atmosphere.
00005  */
00006 /*============================================================================*/
00007 
00008 #include <float.h>
00009 #include <math.h>
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012 #include "defs.h"
00013 #include "athena.h"
00014 #include "globals.h"
00015 #include "prototypes.h"
00016 
00017 /*==============================================================================
00018  * PRIVATE FUNCTION PROTOTYPES:
00019  * grav_pot() - gravitational potential for 2D problem (accn in Y)
00020  * outflow_ix1() - sets BCs on L-x1 (left edge) of grid used in 2D
00021  * outflow_ox1() - sets BCs on R-x1 (right edge) of grid used in 2D
00022  * outflow_ox2() - sets BCs on R-x2 (right edge) of grid used in 2D
00023  *============================================================================*/
00024 
00025 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00026 static void outflow_ix1(Grid *pGrid);
00027 static void outflow_ox1(Grid *pGrid);
00028 static void outflow_ox2(Grid *pGrid);
00029 static void outflow_ix3(Grid *pGrid);
00030 static void outflow_ox3(Grid *pGrid);
00031 
00032 /*=========================== PUBLIC FUNCTIONS ===============================*/
00033 /*----------------------------------------------------------------------------*/
00034 /* problem:  */
00035 
00036 void problem(Grid *pGrid, Domain *pDomain)
00037 {
00038   int i=0,j=0,k=0;
00039   int is,ie,js,je,ks,ke;
00040   int ifield;                                                                   //magentic configuration
00041   Real x1,x2,x3,r2,rad_bubble=0.25;;
00042 #ifdef MHD
00043   Real b0;
00044 #endif
00045   int Nx1, Nx2, Nx3;
00046   int ixs, jxs, kxs;
00047 
00048   is = pGrid->is;  ie = pGrid->ie;
00049   js = pGrid->js;  je = pGrid->je;
00050   ks = pGrid->ks;  ke = pGrid->ke;
00051 
00052   Nx1 = par_geti("grid","Nx1");
00053   Nx2 = par_geti("grid","Nx2");
00054   Nx3 = par_geti("grid","Nx3");
00055 
00056 /* Read magnetic field strength */
00057 #ifdef MHD
00058   b0 = par_getd("problem","b0");                                                //get ifield
00059   ifield = par_getd("problem","ifield");
00060 #endif
00061 
00062 /* Initialize atmosphere, and bubble centered at y=0.25 and radius=rad_bubble
00063  */
00064 
00065   for (k=ks; k<=ke; k++) {
00066     for (j=js; j<=je; j++) {
00067       for (i=is; i<=ie; i++) {
00068         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00069         r2 = x1*x1 + x2*x2 + x3*x3;
00070         pGrid->U[k][j][i].d = pow((1.0 + r2),-0.75);
00071         pGrid->U[k][j][i].M1 = 0.0;
00072         pGrid->U[k][j][i].M2 = 0.0;
00073         pGrid->U[k][j][i].M3 = 0.0;
00074         pGrid->U[k][j][i].E = pow((1.0 + r2),-0.75)/(Gamma*Gamma_1);
00075 
00076         r2 = x1*x1 + (x2-0.25)*(x2-0.25) + x3*x3;
00077         if (r2 < (rad_bubble*rad_bubble)) {
00078           pGrid->U[k][j][i].d = 0.01;
00079         }
00080 #ifdef MHD                                                                      //choose ifield
00081         if (ifield==1){
00082          pGrid->B1i[k][j][i] = b0;
00083          pGrid->U[k][j][i].B1c = b0;
00084          pGrid->U[k][j][i].E += 0.5*b0*b0;
00085          pGrid->B1i[k][j][ie+1] = b0;
00086         }
00087         if (ifield==2){
00088          pGrid->B2i[k][j][i] = b0;
00089          pGrid->U[k][j][i].B2c = b0;
00090          pGrid->U[k][j][i].E += 0.5*b0*b0;
00091          pGrid->B2i[k][je+1][i] = b0;
00092         }
00093         if (ifield==3){
00094          pGrid->B3i[k][j][i] = b0;
00095          pGrid->U[k][j][i].B3c = b0;
00096          pGrid->U[k][j][i].E += 0.5*b0*b0;
00097          pGrid->B3i[ke+1][je][i] = b0;
00098         }
00099 #endif
00100       }
00101     }
00102   }
00103 
00104 /* Enroll gravitational potential, and special BC fns */
00105 
00106   StaticGravPot = grav_pot;
00107   set_bvals_mhd_fun(left_x1,  outflow_ix1);
00108   set_bvals_mhd_fun(right_x1, outflow_ox1);
00109   set_bvals_mhd_fun(right_x2, outflow_ox2);
00110   set_bvals_mhd_fun(left_x3,  outflow_ix3);
00111   set_bvals_mhd_fun(right_x3, outflow_ox3);
00112 
00113 /* With viscosity and/or resistivity, read eta_Ohm and nu_V */
00114 
00115 #ifdef OHMIC
00116   eta_Ohm = par_getd("problem","eta");
00117 #endif
00118 #ifdef NAVIER_STOKES
00119   nu_V = par_getd("problem","nu");
00120 #endif
00121 #ifdef BRAGINSKII
00122   nu_V = par_getd("problem","nu");
00123 #endif
00124 
00125   return;
00126 }
00127 
00128 /*==============================================================================
00129  * PROBLEM USER FUNCTIONS:
00130  * problem_write_restart() - writes problem-specific user data to restart files
00131  * problem_read_restart()  - reads problem-specific user data from restart files
00132  * get_usr_expr()          - sets pointer to expression for special output data
00133  * get_usr_out_fun()       - returns a user defined output function pointer
00134  * get_usr_par_prop()      - returns a user defined particle selection function
00135  * Userwork_in_loop        - problem specific work IN     main loop
00136  * Userwork_after_loop     - problem specific work AFTER  main loop
00137  *----------------------------------------------------------------------------*/
00138 
00139 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00140 {
00141   return;
00142 }
00143 
00144 /*
00145  * 'problem_read_restart' must enroll special boundary value functions,
00146  *    and initialize gravity on restarts
00147  */
00148 
00149 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00150 {
00151   StaticGravPot = grav_pot;
00152   set_bvals_mhd_fun(left_x1,  outflow_ix1);
00153   set_bvals_mhd_fun(right_x1, outflow_ox1);
00154   set_bvals_mhd_fun(right_x2, outflow_ox2);
00155   set_bvals_mhd_fun(left_x3,  outflow_ix3);
00156   set_bvals_mhd_fun(right_x3, outflow_ox3);
00157 
00158   return;
00159 }
00160 
00161 Gasfun_t get_usr_expr(const char *expr)
00162 {
00163   return NULL;
00164 }
00165 
00166 VGFunout_t get_usr_out_fun(const char *name){
00167   return NULL;
00168 }
00169 
00170 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00171 {
00172 }
00173 
00174 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00175 {
00176 }
00177 
00178 /*=========================== PRIVATE FUNCTIONS ==============================*/
00179 /*! \fn static Real grav_pot(const Real x1, const Real x2, const Real x3)
00180  *  \brief Static gravitational potential */
00181 
00182 static Real grav_pot(const Real x1, const Real x2, const Real x3)
00183 {
00184   Real r2;
00185   r2 = x1*x1 + x2*x2 + x3*x3;
00186   return (0.75/Gamma)*log(1.0+r2);
00187 }
00188 
00189 /*----------------------------------------------------------------------------*/
00190 /*! \fn static void outflow_ix1(Grid *pGrid)
00191  *  \brief Special outflow boundary function
00192  */
00193 
00194 static void outflow_ix1(Grid *pGrid)
00195 {
00196   int is = pGrid->is;
00197   int js = pGrid->js, je = pGrid->je;
00198   int ks = pGrid->ks, ke = pGrid->ke;
00199   int i,j,k;
00200   Real r2,x1,x2,x3,d0,p0;
00201 #ifdef MHD
00202   int ibc_x1,ju,ku; /* j-upper, k-upper */
00203 #endif
00204 
00205   for (k=ks; k<=ke; k++) {
00206     for (j=js; j<=je; j++) {
00207       for (i=1; i<=nghost; i++) {
00208         cc_pos(pGrid,is,j,k,&x1,&x2,&x3);
00209         r2 = x1*x1 + x2*x2 + x3*x3;
00210         d0 = pGrid->U[k][j][is].d/pow((1.0 + r2),-0.75);
00211 
00212         p0 = pGrid->U[k][j][is].E - 0.5*(SQR(pGrid->U[k][j][is].M1)
00213           + SQR(pGrid->U[k][j][is].M2) + SQR(pGrid->U[k][j][is].M3))
00214              /pGrid->U[k][j][is].d;
00215         p0 /= pow((1.0 + r2),-0.75);
00216 
00217         cc_pos(pGrid,(is-i),j,k,&x1,&x2,&x3);
00218         r2 = x1*x1 + x2*x2 + x3*x3;
00219         pGrid->U[k][j][is-i].d = d0*pow((1.0 + r2),-0.75);
00220         pGrid->U[k][j][is-i].M1 = 0.0;
00221         pGrid->U[k][j][is-i].M2 = 0.0;
00222         pGrid->U[k][j][is-i].M3 = 0.0;
00223         pGrid->U[k][j][is-i].E = p0*pow((1.0 + r2),-0.75);
00224       }
00225     }
00226   }
00227 
00228 #ifdef MHD
00229   for (k=ks; k<=ke; k++) {
00230     for (j=js; j<=je; j++) {
00231       for (i=1; i<=nghost; i++) {
00232         pGrid->B1i[k][j][is-i]   = pGrid->B1i[k][j][is];
00233         pGrid->U[k][j][is-i].B1c = pGrid->U[k][j][is].B1c;
00234       }
00235     }
00236   }
00237 
00238   if (pGrid->Nx2 > 1) ju=je+1; else ju=je;
00239   for (k=ks; k<=ke; k++) {
00240     for (j=js; j<=ju; j++) {
00241       for (i=1; i<=nghost; i++) {
00242         pGrid->B2i[k][j][is-i]   = pGrid->B2i[k][j][is];
00243         pGrid->U[k][j][is-i].B2c = pGrid->U[k][j][is].B2c;
00244       }
00245     }
00246   }
00247 
00248   if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00249   for (k=ks; k<=ku; k++) {
00250     for (j=js; j<=je; j++) {
00251       for (i=1; i<=nghost; i++) {
00252         pGrid->B3i[k][j][is-i]   = pGrid->B3i[k][j][is];
00253         pGrid->U[k][j][is-i].B3c = pGrid->U[k][j][is].B3c;
00254       }
00255     }
00256   }
00257 #endif /* MHD */
00258 
00259   return;
00260 }
00261 
00262 /*----------------------------------------------------------------------------*/
00263 /*! \fn static void outflow_ox1(Grid *pGrid)
00264  *  \brief Special outflow boundary function
00265  */
00266 
00267 static void outflow_ox1(Grid *pGrid)
00268 {
00269   int ie = pGrid->ie;
00270   int js = pGrid->js, je = pGrid->je;
00271   int ks = pGrid->ks, ke = pGrid->ke;
00272   int i,j,k;
00273   Real r2,x1,x2,x3,d0,p0;
00274 #ifdef MHD
00275   int obc_x1,ju,ku; /* j-upper, k-upper */
00276 #endif
00277 
00278   for (k=ks; k<=ke; k++) {
00279     for (j=js; j<=je; j++) {
00280       for (i=1; i<=nghost; i++) {
00281         cc_pos(pGrid,ie,j,k,&x1,&x2,&x3);
00282         r2 = x1*x1 + x2*x2 + x3*x3;
00283         d0 = pGrid->U[k][j][ie].d/pow((1.0 + r2),-0.75);
00284 
00285         p0 = pGrid->U[k][j][ie].E - 0.5*(SQR(pGrid->U[k][j][ie].M1)
00286           + SQR(pGrid->U[k][j][ie].M2) + SQR(pGrid->U[k][j][ie].M3))
00287              /pGrid->U[k][j][ie].d;
00288         p0 /= pow((1.0 + r2),-0.75);
00289 
00290         cc_pos(pGrid,(ie+i),j,k,&x1,&x2,&x3);
00291         r2 = x1*x1 + x2*x2 + x3*x3;
00292         pGrid->U[k][j][ie+i].d = d0*pow((1.0 + r2),-0.75);
00293         pGrid->U[k][j][ie+i].M1 = 0.0;
00294         pGrid->U[k][j][ie+i].M2 = 0.0;
00295         pGrid->U[k][j][ie+i].M3 = 0.0;
00296         pGrid->U[k][j][ie+i].E = p0*pow((1.0 + r2),-0.75);
00297       }
00298     }
00299   }
00300 
00301 #ifdef MHD
00302 /* i=ie+1 is not set for the interface field B1i */
00303   for (k=ks; k<=ke; k++) {
00304     for (j=js; j<=je; j++) {
00305       for (i=1; i<=nghost; i++) {
00306         if (i>1) pGrid->B1i[k][j][ie+i]   = pGrid->B1i[k][j][ie+1];
00307         pGrid->U[k][j][ie+i].B1c = pGrid->U[k][j][ie].B1c;
00308       }
00309     }
00310   }
00311 
00312   if (pGrid->Nx2 > 1) ju=je+1; else ju=je;
00313   for (k=ks; k<=ke; k++) {
00314     for (j=js; j<=ju; j++) {
00315       for (i=1; i<=nghost; i++) {
00316         pGrid->B2i[k][j][ie+i]   = pGrid->B2i[k][j][ie];
00317         pGrid->U[k][j][ie+i].B2c = pGrid->U[k][j][ie].B2c;
00318       }
00319     }
00320   }
00321 
00322   if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00323   for (k=ks; k<=ku; k++) {
00324     for (j=js; j<=je; j++) {
00325       for (i=1; i<=nghost; i++) {
00326         pGrid->B3i[k][j][ie+i]   = pGrid->B3i[k][j][ie];
00327         pGrid->U[k][j][ie+i].B3c = pGrid->U[k][j][ie].B3c;
00328       }
00329     }
00330   }
00331 #endif /* MHD */
00332 
00333   return;
00334 }
00335 
00336 /*----------------------------------------------------------------------------*/
00337 /*! \fn static void outflow_ox2(Grid *pGrid)
00338  *  \brief Special outflow boundary function
00339  */
00340 
00341 static void outflow_ox2(Grid *pGrid)
00342 {
00343   int je = pGrid->je;
00344   int ks = pGrid->ks, ke = pGrid->ke;
00345   int i,j,k,il,iu; /* i-lower/upper */
00346   Real r2,x1,x2,x3,d0,p0;
00347 #ifdef MHD
00348   int obc_x2,ku; /* k-upper */
00349 #endif
00350 
00351   if (pGrid->Nx1 > 1){
00352     iu = pGrid->ie + nghost;
00353     il = pGrid->is - nghost;
00354   } else {
00355     iu = pGrid->ie;
00356     il = pGrid->is;
00357   }
00358 
00359   for (k=ks; k<=ke; k++) {
00360     for (j=1; j<=nghost; j++) {
00361       for (i=il; i<=iu; i++) {
00362         cc_pos(pGrid,i,je,k,&x1,&x2,&x3);
00363         r2 = x1*x1 + x2*x2 + x3*x3;
00364         d0 = pGrid->U[k][je][i].d/pow((1.0 + r2),-0.75);
00365 
00366         p0 = pGrid->U[k][je][i].E - 0.5*(SQR(pGrid->U[k][je][i].M1)
00367           + SQR(pGrid->U[k][je][i].M2) + SQR(pGrid->U[k][je][i].M3))
00368              /pGrid->U[k][je][i].d;
00369         p0 /= pow((1.0 + r2),-0.75);
00370 
00371         cc_pos(pGrid,i,(je+j),k,&x1,&x2,&x3);
00372         r2 = x1*x1 + x2*x2 + x3*x3;
00373         pGrid->U[k][je+j][i].d = d0*pow((1.0 + r2),-0.75);
00374         pGrid->U[k][je+j][i].M1 = 0.0;
00375         pGrid->U[k][je+j][i].M2 = 0.0;
00376         pGrid->U[k][je+j][i].M3 = 0.0;
00377         pGrid->U[k][je+j][i].E = p0*pow((1.0 + r2),-0.75);
00378       }
00379     }
00380   }
00381 
00382 #ifdef MHD
00383   for (k=ks; k<=ke; k++) {
00384     for (j=1; j<=nghost; j++) {
00385       for (i=il; i<=iu; i++) {
00386         pGrid->B1i[k][je+j][i]   = pGrid->B1i[k][je][i];
00387         pGrid->U[k][je+j][i].B1c = pGrid->U[k][je][i].B1c;
00388       }
00389     }
00390   }
00391 
00392 /* j=je+1 is not set for the interface field B2i */
00393   for (k=ks; k<=ke; k++) {
00394     for (j=1; j<=nghost; j++) {
00395       for (i=il; i<=iu; i++) {
00396          if (j>1) pGrid->B2i[k][je+j][i]   = pGrid->B2i[k][je+1][i];
00397         pGrid->U[k][je+j][i].B2c = pGrid->U[k][je][i].B2c;
00398       }
00399     }
00400   }
00401 
00402   if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00403   for (k=ks; k<=ku; k++) {
00404     for (j=1; j<=nghost; j++) {
00405       for (i=il; i<=iu; i++) {
00406         pGrid->B3i[k][je+j][i]   = pGrid->B3i[k][je][i];
00407         pGrid->U[k][je+j][i].B3c = pGrid->U[k][je][i].B3c;
00408       }
00409     }
00410   }
00411 #endif /* MHD */
00412 
00413   return;
00414 }
00415 
00416 /*----------------------------------------------------------------------------*/
00417 /*! \fn static void outflow_ix3(Grid *pGrid)
00418  *  \brief Special outflow boundary function
00419  */
00420 
00421 static void outflow_ix3(Grid *pGrid)
00422 {
00423   int ks = pGrid->ks;
00424   int i,j,k,il,iu,jl,ju; /* i-lower/upper;  j-lower/upper */
00425   Real r2,x1,x2,x3,d0,p0;
00426 
00427   if (pGrid->Nx1 > 1){
00428     iu = pGrid->ie + nghost;
00429     il = pGrid->is - nghost;
00430   } else {
00431     iu = pGrid->ie;
00432     il = pGrid->is;
00433   }
00434   if (pGrid->Nx2 > 1){
00435     ju = pGrid->je + nghost;
00436     jl = pGrid->js - nghost;
00437   } else {
00438     ju = pGrid->je;
00439     jl = pGrid->js;
00440   }
00441 
00442   for (k=1; k<=nghost; k++) {
00443     for (j=jl; j<=ju; j++) {
00444       for (i=il; i<=iu; i++) {
00445         cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00446         r2 = x1*x1 + x2*x2 + x3*x3;
00447         d0 = pGrid->U[ks][j][i].d/pow((1.0 + r2),-0.75);
00448 
00449         p0 = pGrid->U[ks][j][i].E - 0.5*(SQR(pGrid->U[ks][j][i].M1)
00450           + SQR(pGrid->U[ks][j][i].M2) + SQR(pGrid->U[ks][j][i].M3))
00451              /pGrid->U[ks][j][i].d;
00452         p0 /= pow((1.0 + r2),-0.75);
00453 
00454         cc_pos(pGrid,i,j,(ks-k),&x1,&x2,&x3);
00455         r2 = x1*x1 + x2*x2 + x3*x3;
00456         pGrid->U[ks-k][j][i].d = d0*pow((1.0 + r2),-0.75);
00457         pGrid->U[ks-k][j][i].M1 = 0.0;
00458         pGrid->U[ks-k][j][i].M2 = 0.0;
00459         pGrid->U[ks-k][j][i].M3 = 0.0;
00460         pGrid->U[ks-k][j][i].E = p0*pow((1.0 + r2),-0.75);
00461       }
00462     }
00463   }
00464 
00465 #ifdef MHD
00466   for (k=1; k<=nghost; k++) {
00467     for (j=jl; j<=ju; j++) {
00468       for (i=il; i<=iu; i++) {
00469         pGrid->B1i[ks-k][j][i]   = pGrid->B1i[ks][j][i];
00470         pGrid->U[ks-k][j][i].B1c = pGrid->U[ks][j][i].B1c;
00471       }
00472     }
00473   }
00474 
00475   for (k=1; k<=nghost; k++) {
00476     for (j=jl; j<=ju; j++) {
00477       for (i=il; i<=iu; i++) {
00478         pGrid->B2i[ks-k][j][i]   = pGrid->B2i[ks][j][i];
00479         pGrid->U[ks-k][j][i].B2c = pGrid->U[ks][j][i].B2c;
00480       }
00481     }
00482   }
00483 
00484   for (k=1; k<=nghost; k++) {
00485     for (j=jl; j<=ju; j++) {
00486       for (i=il; i<=iu; i++) {
00487         pGrid->B3i[ks-k][j][i]   = pGrid->B3i[ks][j][i];
00488         pGrid->U[ks-k][j][i].B3c = pGrid->U[ks][j][i].B3c;
00489       }
00490     }
00491   }
00492 #endif /* MHD */
00493 
00494   return;
00495 }
00496 
00497 /*----------------------------------------------------------------------------*/
00498 /*! \fn static void outflow_ox3(Grid *pGrid)
00499  *  \brief Special outflow boundary function
00500  */
00501 
00502 static void outflow_ox3(Grid *pGrid)
00503 {
00504   int ke = pGrid->ke;
00505   int i,j,k ,il,iu,jl,ju; /* i-lower/upper;  j-lower/upper */
00506   Real r2,x1,x2,x3,d0,p0;
00507 
00508   if (pGrid->Nx1 > 1){
00509     iu = pGrid->ie + nghost;
00510     il = pGrid->is - nghost;
00511   } else {
00512     iu = pGrid->ie;
00513     il = pGrid->is;
00514   }
00515   if (pGrid->Nx2 > 1){
00516     ju = pGrid->je + nghost;
00517     jl = pGrid->js - nghost;
00518   } else {
00519     ju = pGrid->je;
00520     jl = pGrid->js;
00521   }
00522 
00523   for (k=1; k<=nghost; k++) {
00524     for (j=jl; j<=ju; j++) {
00525       for (i=il; i<=iu; i++) {
00526         cc_pos(pGrid,i,j,ke,&x1,&x2,&x3);
00527         r2 = x1*x1 + x2*x2 + x3*x3;
00528         d0 = pGrid->U[ke][j][i].d/pow((1.0 + r2),-0.75);
00529 
00530         p0 = pGrid->U[ke][j][i].E - 0.5*(SQR(pGrid->U[ke][j][i].M1)
00531           + SQR(pGrid->U[ke][j][i].M2) + SQR(pGrid->U[ke][j][i].M3))
00532              /pGrid->U[ke][j][i].d;
00533         p0 /= pow((1.0 + r2),-0.75);
00534 
00535         cc_pos(pGrid,i,j,(ke+k),&x1,&x2,&x3);
00536         r2 = x1*x1 + x2*x2 + x3*x3;
00537         pGrid->U[ke+k][j][i].d = d0*pow((1.0 + r2),-0.75);
00538         pGrid->U[ke+k][j][i].M1 = 0.0;
00539         pGrid->U[ke+k][j][i].M2 = 0.0;
00540         pGrid->U[ke+k][j][i].M3 = 0.0;
00541         pGrid->U[ke+k][j][i].E = p0*pow((1.0 + r2),-0.75);
00542       }
00543     }
00544   }
00545 
00546 #ifdef MHD
00547   for (k=1; k<=nghost; k++) {
00548     for (j=jl; j<=ju; j++) {
00549       for (i=il; i<=iu; i++) {
00550         pGrid->B1i[ke+k][j][i]   = pGrid->B1i[ke][j][i];
00551         pGrid->U[ke+k][j][i].B1c = pGrid->U[ke][j][i].B1c;
00552       }
00553     }
00554   }
00555 
00556 
00557   for (k=1; k<=nghost; k++) {
00558     for (j=jl; j<=ju; j++) {
00559       for (i=il; i<=iu; i++) {
00560         pGrid->B2i[ke+k][j][i]   = pGrid->B2i[ke][j][i];
00561         pGrid->U[ke+k][j][i].B2c = pGrid->U[ke][j][i].B2c;
00562       }
00563     }
00564   }
00565 
00566 /* k=ke+1 is not set for the interface field B3i */
00567   for (k=1; k<=nghost; k++) {
00568     for (j=jl; j<=ju; j++) {
00569       for (i=il; i<=iu; i++) {
00570         if (k>1) pGrid->B3i[ke+k][j][i]   = pGrid->B3i[ke+1][j][i];
00571         pGrid->U[ke+k][j][i].B3c = pGrid->U[ke][j][i].B3c;
00572       }
00573     }
00574   }
00575 #endif /* MHD */
00576 
00577   return;
00578 }