Athena: prob/shkset3d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file shkset3d.c
00004  *  \brief Sets up shock at angle to grid to test 3D algorithms.
00005  *
00006  * PURPOSE: Sets up shock at angle to grid to test 3D algorithms.  Setting up
00007  *   the initial conditions to minimize grid noise is very complex, so this
00008  *   special function has been written to handle only 3D problems.
00009  *
00010  * This problem cannot be run with Static Mesh Refinement.
00011  *
00012  * This code is most easily understood in terms of a one dimensional
00013  * problem in the coordinate system (x1,x2,x3).  Two coordinate rotations are
00014  * applied to obtain a new wave vector in a 3D space in the (x,y,z)
00015  * coordinate system.
00016  *
00017  *   First rotate about the x2 axis:
00018  *  -  x1' = x1*cos(ang_2) - x3*sin(ang_2)
00019  *  -  x2' = x2
00020  *  -  x3' = x1*sin(ang_2) + x3*cos(ang_2)
00021  *
00022  *   Next rotate about the x3' axis:
00023  *   - x = x1'*cos(ang_3) - x2'*sin(ang_3)
00024  *   - y = x1'*sin(ang_3) + x2'*cos(ang_3)
00025  *   - z = x3'
00026  *
00027  *   Expanding this out we get:
00028  *   - x = x1*cos(ang_2)*cos(ang_3) - x2*sin(ang_3) - x3*sin(ang_2)*cos(ang_3)
00029  *   - y = x1*cos(ang_2)*sin(ang_3) + x2*cos(ang_3) - x3*sin(ang_2)*sin(ang_3)
00030  *   - z = x1*sin(ang_2)                            + x3*cos(ang_2)
00031  *
00032  *   This inverts to:
00033  *   - x1 =  x*cos(ang_2)*cos(ang_3) + y*cos(ang_2)*sin(ang_3) + z*sin(ang_2)
00034  *   - x2 = -x*sin(ang_3)            + y*cos(ang_3)
00035  *   - x3 = -x*sin(ang_2)*cos(ang_3) - y*sin(ang_2)*sin(ang_3) + z*cos(ang_2)
00036  * Note these transformations are the same used in linear_wave3d and cpaw3d
00037  *
00038  * The initial conditions must be translation invariant, i.e. 
00039  * q(x,y,z) = q(x+tx, y+ty, z+tz) for (tx,ty,tz) such that 
00040  * x1(x,y,z) = x1(x+tx, y+ty, z+ty).  Inserting this last expression for
00041  * x1(x,y,z) into the transformation above leads to
00042  *
00043  * - tx*cos(ang_2)*cos(ang_3) + ty*cos(ang_2)*sin(ang_3) + tz*sin(ang_2) = 0.
00044  *
00045  * Now, restrict the translation symmetry to be discrete by inserting
00046  * (tx, ty, tz) = (nx*dx, ny*dy, nz*dz) where (nx, ny, nz) are integers (not 
00047  * the number of grid cells in each direction) and (dx, dy, dz) are the grid
00048  * cell sizes.  With some simplification the translation symmetry becomes
00049  *
00050  * - nx + ny*tan(ang_3)*(dy/dx) + nz*(tan(ang_2)/cos(ang_3))*(dz/dx) = 0.
00051  *
00052  * Now, choose an integer unit cell size (rx, ry, rz) where
00053  *
00054  * - tan(ang_3)*(dy/dx) = (rx/ry)
00055  *
00056  * and
00057  *
00058  * - (tan(ang_2)/cos(ang_3))*(dz/dx) = (rx/rz).
00059  *
00060  * With this choice, or equation for discrete translation symmetry becomes
00061  *
00062  * - (nx/rx) + (ny/ry) + (nz/rz) = 0.
00063  *
00064  * PRIVATE FUNCTION PROTOTYPES:
00065  * - lx_bc() - use periodicity of "unit cell" to apply BCs to left-x edge
00066  * - rx_bc() - use periodicity of "unit cell" to apply BCs to right-x edge
00067  * - ly_bc() - use periodicity of "unit cell" to apply BCs to left-y edge
00068  * - ry_bc() - use periodicity of "unit cell" to apply BCs to right-y edge
00069  * - lz_bc() - use periodicity of "unit cell" to apply BCs to left-z edge
00070  * - rz_bc() - use periodicity of "unit cell" to apply BCs to right-z edge
00071  * - Ax() - x-component of vector potential for initial conditions
00072  * - Ay() - y-component of vector potential for initial conditions
00073  * - Az() - z-component of vector potential for initial conditions
00074  *
00075  * This is the equation for a set of discrete points which lie in a
00076  * plane and is a key relation for setting ghost cells in the boundary
00077  * condition routines below.  --  T. A. Gardiner  --  7/21/2006               */
00078 /*============================================================================*/
00079 
00080 #include <math.h>
00081 #include <stdio.h>
00082 #include <stdlib.h>
00083 #include "defs.h"
00084 #include "athena.h"
00085 #include "globals.h"
00086 #include "prototypes.h"
00087 
00088 #ifdef STATIC_MESH_REFINEMENT
00089 #error The shkset3d test does not work with SMR.
00090 #endif
00091 
00092 /* Parameters which define initial solution -- made global so that they can be
00093  * shared with all private functions in this file */
00094 
00095 static int rx, ry, rz;
00096 static Real ang_2, ang_3; /* Rotation angles about the 2 and 3' axis */
00097 static Real sin_a2, cos_a2, sin_a3, cos_a3;
00098 /* The Riemann Problem Left and Right states */
00099 static  Real dl, vxl, vyl, vzl;
00100 static  Real dr, vxr, vyr, vzr;
00101 #ifdef MHD
00102 static  Real Bxl, Byl, Bzl, Bxr, Byr, Bzr;
00103 #endif /* MHD */
00104 #ifdef ADIABATIC
00105 static  Real Pl, Pr;
00106 #endif /* ADIABATIC */
00107 
00108 /*==============================================================================
00109  * PRIVATE FUNCTION PROTOTYPES:
00110  * lx_bc() - use periodicity of "unit cell" to apply BCs to left-x edge
00111  * rx_bc() - use periodicity of "unit cell" to apply BCs to right-x edge
00112  * ly_bc() - use periodicity of "unit cell" to apply BCs to left-y edge
00113  * ry_bc() - use periodicity of "unit cell" to apply BCs to right-y edge
00114  * lz_bc() - use periodicity of "unit cell" to apply BCs to left-z edge
00115  * rz_bc() - use periodicity of "unit cell" to apply BCs to right-z edge
00116  * Ax() - x-component of vector potential for initial conditions
00117  * Ay() - y-component of vector potential for initial conditions
00118  * Az() - z-component of vector potential for initial conditions
00119  *============================================================================*/
00120 
00121 static void lx_bc(GridS *pG);
00122 static void rx_bc(GridS *pG);
00123 static void ly_bc(GridS *pG);
00124 static void ry_bc(GridS *pG);
00125 static void lz_bc(GridS *pG);
00126 static void rz_bc(GridS *pG);
00127 static Real Ax(const Real x1, const Real x2, const Real x3);
00128 static Real Ay(const Real x1, const Real x2, const Real x3);
00129 static Real Az(const Real x1, const Real x2, const Real x3);
00130 
00131 /*=========================== PUBLIC FUNCTIONS ===============================*/
00132 /*----------------------------------------------------------------------------*/
00133 /* problem:  */
00134 
00135 void problem(DomainS *pDomain)
00136 {
00137   GridS *pGrid = pDomain->Grid;
00138   int n;
00139   int i, j, k, ix, jx, kx, mix, mjx, mkx;
00140 
00141   int Nx1, Nx2, Nx3; /* Complete Grid dimensions */
00142   int sx, sy, sz; /* Number of unit cells in each direction */
00143   div_t id;
00144   double d_ix;
00145   Real xl, x, xr, dt_usr;
00146   Real lx1,lx2,lx3,cx1,cx2,cx3,rx1,rx2,rx3;
00147   Real sp0_x1, sp0_x2, sp0_x3;
00148 
00149 /* Variables assocaiated with a cubic array on which the Riemann problem
00150  * interface normal is aligned along the (1,1,1) direction. */
00151   int N, isqa, jsqa, ksqa;
00152   int scx, scy, scz; /* scx = N/rx, scy = N/ry, scz = N/rz */
00153   Real sdx, sdy, sdz;
00154   ConsS *pq, ***sqa=NULL; /* allocated to dimensions sqa[N][N][2N] */
00155 #ifdef MHD
00156   Real ***sB1i=NULL, ***sB2i=NULL, ***sB3i=NULL;
00157 #endif /* MHD */
00158 
00159 /* The Riemann Problem Left and Right states in the Grid coordinate system */
00160   ConsS ql, qr;
00161 /* The unit cell enclosing the Riemann problem interface. */
00162   ConsS ***qa=NULL; /* allocated to dimensions qa[rz][ry][rx] */
00163 #ifdef MHD
00164   Real ***aB1i=NULL, ***aB2i=NULL, ***aB3i=NULL;
00165 #endif /* MHD */
00166   int qa_min_ix, qa_min_jx, qa_min_kx;
00167   int qa_max_ix, qa_max_jx, qa_max_kx;
00168 
00169 /*--- Step 1 -------------------------------------------------------------------
00170  * Read initial conditions, check that inputs are within valid range, allocate
00171  * memory for arrays.  We know this Domain must be the root Grid, since this
00172  * problem generator does not work with SMR.
00173  */
00174 
00175   Nx1 = pDomain->Nx[0];
00176   Nx2 = pDomain->Nx[1];
00177   Nx3 = pDomain->Nx[2];
00178   if(pDomain->Nx[2] <= 0)
00179     ath_error("[get_set_input]: This problem assumes a 3D domain\n");
00180 
00181 /* Get dimensions of unit cell, check that total grid size is an integer number
00182  * of cells.  Reminder: the expression x % y in C returns the remainder when
00183  * x is divided by y */
00184 
00185   rx = par_geti("problem","rx");
00186   ry = par_geti("problem","ry");
00187   rz = par_geti("problem","rz");
00188 
00189   if(rx <= 0 || ry <= 0 || rz <= 0)
00190     ath_error("[get_set_input]: Only rx > 0, ry > 0 and rz > 0 is allowed\n");
00191 
00192   id = div(Nx1,rx);
00193   sx = id.quot;
00194   if(id.rem != 0)
00195     ath_error("[get_set_input]: Choose Nx1 % rx = 0 instead of %d\n",Nx1%rx);
00196 
00197   id = div(Nx2,ry);
00198   sy = id.quot;
00199   if(id.rem != 0)
00200     ath_error("[get_set_input]: Choose Nx2 % ry = 0 instead of %d\n",Nx2%ry);
00201 
00202   id = div(Nx3,rz);
00203   sz = id.quot;
00204   if(id.rem != 0)
00205     ath_error("[get_set_input]: Choose Nx3 % rz = 0 instead of %d\n",Nx3%rz);
00206 
00207   if(sy > sx || sz > sx)
00208     ath_error("Nx1/rx is assumed to be larger than the y- and z-dir. ratios\n");
00209 
00210 /* Imposing translation symmetry over the domain (rx,ry,rz) gives */
00211 
00212   ang_3 = atan((rx*pGrid->dx1)/(ry*pGrid->dx2));
00213   sin_a3 = sin(ang_3);
00214   cos_a3 = cos(ang_3);
00215 
00216   ang_2 = atan((rx*pGrid->dx1*cos_a3)/(rz*pGrid->dx3));
00217   sin_a2 = sin(ang_2);
00218   cos_a2 = cos(ang_2);
00219 
00220 /* Read in the left and right states in the Riemann problem */
00221 
00222   dl  = par_getd("problem","dl");
00223   vxl = par_getd("problem","vxl");
00224   vyl = par_getd("problem","vyl");
00225   vzl = par_getd("problem","vzl");
00226 #ifdef MHD
00227   Bxl = par_getd("problem","Bxl");
00228   Byl = par_getd("problem","Byl");
00229   Bzl = par_getd("problem","Bzl");
00230 #endif /* MHD */
00231 #ifdef ADIABATIC
00232   Pl  = par_getd("problem","pl");
00233 #endif /* ADIABATIC */
00234 
00235   dr  = par_getd("problem","dr");
00236   vxr = par_getd("problem","vxr");
00237   vyr = par_getd("problem","vyr");
00238   vzr = par_getd("problem","vzr");
00239 #ifdef MHD
00240   Bxr = par_getd("problem","Bxr");
00241   Byr = par_getd("problem","Byr");
00242   Bzr = par_getd("problem","Bzr");
00243 #endif /* MHD */
00244 #ifdef ADIABATIC
00245   Pr  = par_getd("problem","pr");
00246 #endif /* ADIABATIC */
00247 
00248 /* Calculate the conservative (ConsS) L/R states in the Grid's coordinates */
00249 
00250 /* This is the left state */
00251   ql.d  = dl;
00252   ql.M1 = dl*(vxl*cos_a2*cos_a3 - vyl*sin_a3 - vzl*sin_a2*cos_a3);
00253   ql.M2 = dl*(vxl*cos_a2*sin_a3 + vyl*cos_a3 - vzl*sin_a2*sin_a3);
00254   ql.M3 = dl*(vxl*sin_a2                     + vzl*cos_a2);
00255 
00256 #ifdef MHD
00257   ql.B1c = (Bxl*cos_a2*cos_a3 - Byl*sin_a3 - Bzl*sin_a2*cos_a3);
00258   ql.B2c = (Bxl*cos_a2*sin_a3 + Byl*cos_a3 - Bzl*sin_a2*sin_a3);
00259   ql.B3c = (Bxl*sin_a2                     + Bzl*cos_a2);
00260 #endif /* MHD */
00261 
00262 #ifdef ADIABATIC
00263   ql.E  = Pl/Gamma_1
00264 #ifdef MHD
00265     + 0.5*(Bxl*Bxl + Byl*Byl + Bzl*Bzl)
00266 #endif /* MHD */
00267     + 0.5*dl*(vxl*vxl + vyl*vyl + vzl*vzl);
00268 #endif
00269 
00270 /* This is the right state */
00271   qr.d  = dr;
00272   qr.M1 = dr*(vxr*cos_a2*cos_a3 - vyr*sin_a3 - vzr*sin_a2*cos_a3);
00273   qr.M2 = dr*(vxr*cos_a2*sin_a3 + vyr*cos_a3 - vzr*sin_a2*sin_a3);
00274   qr.M3 = dr*(vxr*sin_a2                     + vzr*cos_a2);
00275 
00276 #ifdef MHD
00277   qr.B1c = (Bxr*cos_a2*cos_a3 - Byr*sin_a3 - Bzr*sin_a2*cos_a3);
00278   qr.B2c = (Bxr*cos_a2*sin_a3 + Byr*cos_a3 - Bzr*sin_a2*sin_a3);
00279   qr.B3c = (Bxr*sin_a2                     + Bzr*cos_a2);
00280 #endif /* MHD */
00281 
00282 #ifdef ADIABATIC
00283   qr.E  = Pr/Gamma_1
00284 #ifdef MHD
00285     + 0.5*(Bxr*Bxr + Byr*Byr + Bzr*Bzr)
00286 #endif /* MHD */
00287     + 0.5*dr*(vxr*vxr + vyr*vyr + vzr*vzr);
00288 #endif
00289 
00290 /* Allocate 3D arrays with size of UNIT CELL */
00291 
00292   if((qa = (ConsS***)calloc_3d_array(rz, ry, rx+rx, sizeof(ConsS)))==NULL)
00293     ath_error("[shkset3d]: Error allocating qa[%d][%d][%d]\n",rz,ry,rx+rx);
00294 
00295 #ifdef MHD
00296   if((aB1i = (Real***)calloc_3d_array(rz, ry, rx+rx+1, sizeof(Real))) == NULL)
00297     ath_error("[shkset3d]: Error allocating aB1i[%d][%d][%d]\n",rz,ry,rx+rx+1);
00298 
00299   if((aB2i = (Real***)calloc_3d_array(rz, ry+1, rx+rx, sizeof(Real))) == NULL)
00300     ath_error("[shkset3d]: Error allocating aB2i[%d][%d][%d]\n",rz,ry+1,rx+rx);
00301 
00302   if((aB3i = (Real***)calloc_3d_array(rz+1, ry, rx+rx, sizeof(Real))) == NULL)
00303     ath_error("[shkset3d]: Error allocating aB3i[%d][%d][%d]\n",rz+1,ry,rx+rx);
00304 #endif /* MHD */
00305 
00306 /* Initialize qa_min_co and qa_max_co */
00307   qa_min_jx = 0;  qa_max_jx = ry;
00308   qa_min_kx = 0;  qa_max_kx = rz;
00309 
00310 /* qa_min_ix comes from calculating x=0 at jx=ry, kx=rz */
00311 /* qa_max_ix comes from calculating x=0 at jx=0, kx=0 */
00312 /* qa_max_ix - qa_min_ix = 2*rx! */
00313 
00314   d_ix = -pGrid->MinX[0]/pGrid->dx1 - rx*(pGrid->MinX[1]/(ry*pGrid->dx2) 
00315     + pGrid->MinX[2]/(rz*pGrid->dx3));
00316   qa_max_ix = ceil(d_ix);
00317   qa_min_ix = qa_max_ix - rx - rx;
00318 
00319 /* Note that this code has the undesireable "feature" that when I wrote it I
00320  * assumed that d_ix would be equal to an integer.  This is not always the 
00321  * case..., so if it fails we should quit and give the user advice
00322  * T. A. Gardiner  --  7/21/2006  */
00323 
00324   if(qa_max_ix - d_ix > 1.0e-12){
00325     ath_perr(-1,"[shkset3d]: qa_max_ix - d_ix = %g\n",
00326             qa_max_ix - d_ix);
00327     ath_error("[shkset3d]: Try setting the x2min and x3min = 0.0\n");
00328   }
00329 
00330 /*--- Step 2 -------------------------------------------------------------------
00331  * Calculate the size of a square domain which can be restriced in an integral
00332  * fashion to the relevant grid qa[][][].  Note, one can multiply N, scx, scy,
00333  * and scz by some common integer to increase the resolution to the
00334  * approximation of the volume and interface integral averages.  Tests up to
00335  * now have shown no difference in the initial conditions when multiplying by
00336  * 3 or 4.*/
00337 
00338   N = rx*ry*rz;
00339   scx = ry*rz;
00340   scy = rx*rz;
00341   scz = rx*ry;
00342 
00343   sdx = pGrid->dx1/scx;
00344   sdy = pGrid->dx2/scy;
00345   sdz = pGrid->dx3/scz;
00346 
00347   if((sqa = (ConsS***)calloc_3d_array(N, N, N+N, sizeof(ConsS))) == NULL)
00348     ath_error("[shkset3d]: Error allocating sqa[%d][%d][%d]\n",N,N,N+N);
00349 
00350 #ifdef MHD
00351   if((sB1i = (Real***)calloc_3d_array(N, N, N+N+1, sizeof(Real))) == NULL)
00352     ath_error("[shkset3d]: Error allocating sB1i[%d][%d][%d]\n",N,N,N+N+1);
00353 
00354   if((sB2i = (Real***)calloc_3d_array(N, N+1, N+N, sizeof(Real))) == NULL)
00355     ath_error("[shkset3d]: Error allocating sB2i[%d][%d][%d]\n",N,N+1,N+N);
00356 
00357   if((sB3i = (Real***)calloc_3d_array(N+1, N, N+N, sizeof(Real))) == NULL)
00358     ath_error("[shkset3d]: Error allocating sB3i[%d][%d][%d]\n",N+1,N,N+N);
00359 #endif /* MHD */
00360 
00361 /* Calculate the Position of the left-most corner of the sqa and qa grids. */
00362   sp0_x1 = pGrid->MinX[0] + qa_min_ix*pGrid->dx1;
00363   sp0_x2 = pGrid->MinX[1] + qa_min_jx*pGrid->dx2;
00364   sp0_x3 = pGrid->MinX[2] + qa_min_kx*pGrid->dx3;
00365 
00366 /*--- Step 3 -------------------------------------------------------------------
00367  * First, initialize the interface magnetic fields in sqa
00368  */
00369 
00370 #ifdef MHD
00371   for(k=0; k<=N; k++){
00372     for(j=0; j<=N; j++){
00373       for(i=0; i<=(N+N); i++){
00374 /* Calculate the right corner position in sqa[k][j][i] */
00375         rx1 = sp0_x1 + (i+1)*sdx;
00376         rx2 = sp0_x2 + (j+1)*sdy;
00377         rx3 = sp0_x3 + (k+1)*sdz;
00378 
00379 /* Calculate the left corner position in sqa[k][j][i] */
00380         lx1 = sp0_x1 + i*sdx;
00381         lx2 = sp0_x2 + j*sdy;
00382         lx3 = sp0_x3 + k*sdz;
00383 
00384 /* Calculate the cell-center position and the x-position */
00385         cx1 = lx1 + 0.5*sdx;
00386         cx2 = lx2 + 0.5*sdy;
00387         cx3 = lx3 + 0.5*sdz;
00388 
00389 /* Initialize B1i */
00390         if(j<N && k<N){ 
00391 /* Calculate the left- and right-most x-position on this x1-face */
00392           xl = lx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + lx3*sin_a2;
00393           xr = lx1*cos_a2*cos_a3 + rx2*cos_a2*sin_a3 + rx3*sin_a2;
00394 
00395           if(xl >= 0.0){ /* This cell is in the right state */
00396             sB1i[k][j][i] = qr.B1c;
00397           }
00398           else if(xr <= 0.0){ /* This cell is in the left state */
00399             sB1i[k][j][i] = ql.B1c;
00400           }
00401           else{
00402             sB1i[k][j][i] = Bxl*cos_a2*cos_a3 +
00403               (Az(lx1, rx2, cx3) - Az(lx1, lx2, cx3) )/sdy -
00404               (Ay(lx1, cx2, rx3) - Ay(lx1, cx2, lx3) )/sdz;
00405           }
00406         }
00407 
00408 /* Initialize B2i */
00409         if(i<(N+N) && k<N){ 
00410 /* Calculate the left- and right-most x-position on this x2-face */
00411           xl = lx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + lx3*sin_a2;
00412           xr = rx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + rx3*sin_a2;
00413 
00414           if(xl >= 0.0){ /* This cell is in the right state */
00415             sB2i[k][j][i] = qr.B2c;
00416           }
00417           else if(xr <= 0.0){ /* This cell is in the left state */
00418             sB2i[k][j][i] = ql.B2c;
00419           }
00420           else{
00421             sB2i[k][j][i] = Bxl*cos_a2*sin_a3 +
00422               (Ax(cx1, lx2, rx3) - Ax(cx1, lx2, lx3) )/sdz -
00423               (Az(rx1, lx2, cx3) - Az(lx1, lx2, cx3) )/sdx;
00424           }
00425         }
00426 
00427 /* Initialize B3i */
00428         if(i<(N+N) && j<N){ 
00429 /* Calculate the left- and right-most x-position on this x3-face */
00430           xl = lx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + lx3*sin_a2;
00431           xr = rx1*cos_a2*cos_a3 + rx2*cos_a2*sin_a3 + lx3*sin_a2;
00432 
00433           if(xl >= 0.0){ /* This cell is in the right state */
00434             sB3i[k][j][i] = qr.B3c;
00435           }
00436           else if(xr <= 0.0){ /* This cell is in the left state */
00437             sB3i[k][j][i] = ql.B3c;
00438           }
00439           else{
00440             sB3i[k][j][i] = Bxl*sin_a2 +
00441               (Ay(rx1, cx2, lx3) - Ay(lx1, cx2, lx3) )/sdx -
00442               (Ax(cx1, rx2, lx3) - Ax(cx1, lx2, lx3) )/sdy;
00443           }
00444         }
00445       }
00446     }
00447   }
00448 #endif /* MHD */
00449 
00450 /*--- Step 4 -------------------------------------------------------------------
00451  * Next, initialize the cell-centered quantities (density, momenta, total
00452  * energy, and cell-centered B) in sqa
00453  */
00454 
00455   for(k=0; k<N; k++){
00456     for(j=0; j<N; j++){
00457       for(i=0; i<(N+N); i++){
00458 /* Calculate the right corner position in sqa[k][j][i] */
00459         rx1 = sp0_x1 + (i+1)*sdx;
00460         rx2 = sp0_x2 + (j+1)*sdy;
00461         rx3 = sp0_x3 + (k+1)*sdz;
00462 
00463         xr = rx1*cos_a2*cos_a3 + rx2*cos_a2*sin_a3 + rx3*sin_a2;
00464 
00465 /* Calculate the left corner position in sqa[k][j][i] */
00466         lx1 = sp0_x1 + i*sdx;
00467         lx2 = sp0_x2 + j*sdy;
00468         lx3 = sp0_x3 + k*sdz;
00469 
00470         xl = lx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + lx3*sin_a2;
00471 
00472 /* Calculate the cell-center position and the x-position */
00473         cx1 = lx1 + 0.5*sdx;
00474         cx2 = lx2 + 0.5*sdy;
00475         cx3 = lx3 + 0.5*sdz;
00476         x = cx1*cos_a2*cos_a3 + cx2*cos_a2*sin_a3 + cx3*sin_a2;
00477 
00478         if(xr <= 0.0){ /* This cell is in the left state */
00479           sqa[k][j][i] = ql;
00480           continue;
00481         }
00482 
00483         if(xl >= 0.0){ /* This cell is in the right state */
00484           sqa[k][j][i] = qr;
00485           continue;
00486         }
00487 
00488         pq = &(sqa[k][j][i]); /* get a pointer to the fluid element */
00489 
00490 /* Initialize the density, the momenta and the total energy */
00491         if(x < 0.0){ /* Left state */
00492           pq->d  = dl;
00493           pq->M1 = dl*(vxl*cos_a2*cos_a3 - vyl*sin_a3 - vzl*sin_a2*cos_a3);
00494           pq->M2 = dl*(vxl*cos_a2*sin_a3 + vyl*cos_a3 - vzl*sin_a2*sin_a3);
00495           pq->M3 = dl*(vxl*sin_a2                     + vzl*cos_a2);
00496 #ifdef ADIABATIC
00497           pq->E  = Pl/Gamma_1 + 0.5*dl*(vxl*vxl + vyl*vyl + vzl*vzl);
00498 #endif
00499         }
00500         else{ /* Right state */
00501           pq->d  = dr;
00502           pq->M1 = dr*(vxr*cos_a2*cos_a3 - vyr*sin_a3 - vzr*sin_a2*cos_a3);
00503           pq->M2 = dr*(vxr*cos_a2*sin_a3 + vyr*cos_a3 - vzr*sin_a2*sin_a3);
00504           pq->M3 = dr*(vxr*sin_a2                     + vzr*cos_a2);
00505 #ifdef ADIABATIC
00506           pq->E  = Pr/Gamma_1 + 0.5*dr*(vxr*vxr + vyr*vyr + vzr*vzr);
00507 #endif
00508         }
00509 
00510 /* Initialize the cell-center magnetic field components */
00511 #ifdef MHD
00512         pq->B1c = 0.5*(sB1i[k][j][i] + sB1i[k  ][j  ][i+1]);
00513         pq->B2c = 0.5*(sB2i[k][j][i] + sB2i[k  ][j+1][i  ]);
00514         pq->B3c = 0.5*(sB3i[k][j][i] + sB3i[k+1][j  ][i  ]);
00515 
00516 #ifdef ADIABATIC
00517         pq->E += 0.5*(pq->B1c*pq->B1c + pq->B2c*pq->B2c + pq->B3c*pq->B3c);
00518 #endif /* ADIABATIC */
00519 #endif /* MHD */
00520       }
00521     }
00522   }
00523 
00524 /*--- Step 5 -------------------------------------------------------------------
00525  *  Now conservatively restrict sqa[][][] into qa[][][] 
00526  */
00527 
00528 /* First, the interface magnetic fields */
00529 
00530 #ifdef MHD
00531   for(k=0; k<=rz; k++){
00532     for(j=0; j<=ry; j++){
00533       for(i=0; i<=(rx+rx); i++){
00534 
00535 /* Initialize B1i */
00536         if(j<ry && k<rz){ 
00537           aB1i[k][j][i] = 0.0;
00538           isqa = i*scx; /* The left-x1 interface field */
00539           for(ksqa=k*scz; ksqa<(k+1)*scz; ksqa++){
00540             for(jsqa=j*scy; jsqa<(j+1)*scy; jsqa++){
00541               aB1i[k][j][i] += sB1i[ksqa][jsqa][isqa];
00542             }
00543           }
00544           aB1i[k][j][i] /= (Real)(scy*scz);
00545         }
00546 
00547 /* Initialize B2i */
00548         if(i<(rx+rx) && k<rz){ 
00549           aB2i[k][j][i] = 0.0;
00550           jsqa = j*scy; /* The left-x2 interface field */
00551           for(ksqa=k*scz; ksqa<(k+1)*scz; ksqa++){
00552             for(isqa=i*scx; isqa<(i+1)*scx; isqa++){
00553               aB2i[k][j][i] += sB2i[ksqa][jsqa][isqa];
00554             }
00555           }
00556           aB2i[k][j][i] /= (Real)(scx*scz);
00557         }
00558 
00559 /* Initialize B3i */
00560         if(i<(rx+rx) && j<ry){ 
00561           aB3i[k][j][i] = 0.0;
00562           ksqa = k*scz; /* The left-x3 interface field */
00563           for(jsqa=j*scy; jsqa<(j+1)*scy; jsqa++){
00564             for(isqa=i*scx; isqa<(i+1)*scx; isqa++){
00565               aB3i[k][j][i] += sB3i[ksqa][jsqa][isqa];
00566             }
00567           }
00568           aB3i[k][j][i] /= (Real)(scx*scy);
00569         }
00570       }
00571     }
00572   }
00573 #endif /* MHD */
00574 
00575 /* Next, the volume averaged quantities */
00576   for(k=0; k<rz; k++){
00577     for(j=0; j<ry; j++){
00578       for(i=0; i<(rx+rx); i++){
00579         qa[k][j][i].d = 0.0;
00580         qa[k][j][i].M1 = qa[k][j][i].M2 = qa[k][j][i].M3 = 0.0;
00581 #ifdef ADIABATIC
00582         qa[k][j][i].E = 0.0;
00583 #endif
00584 
00585 /* loop over the overlapping elements in the sqa array and sum up the
00586  * conservative variables. */
00587 
00588         for(ksqa=k*scz; ksqa<(k+1)*scz; ksqa++){
00589           for(jsqa=j*scy; jsqa<(j+1)*scy; jsqa++){
00590             for(isqa=i*scx; isqa<(i+1)*scx; isqa++){
00591               qa[k][j][i].d  += sqa[ksqa][jsqa][isqa].d;
00592               qa[k][j][i].M1 += sqa[ksqa][jsqa][isqa].M1;
00593               qa[k][j][i].M2 += sqa[ksqa][jsqa][isqa].M2;
00594               qa[k][j][i].M3 += sqa[ksqa][jsqa][isqa].M3;
00595 #ifdef ADIABATIC
00596               qa[k][j][i].E  += sqa[ksqa][jsqa][isqa].E;
00597 #endif
00598             }
00599           }
00600         }
00601 
00602 /* Normalize them - Note, scx*scy*scz = N*N; */
00603 
00604         qa[k][j][i].d  /= (Real)(N*N);
00605         qa[k][j][i].M1 /= (Real)(N*N);
00606         qa[k][j][i].M2 /= (Real)(N*N);
00607         qa[k][j][i].M3 /= (Real)(N*N);
00608 #ifdef ADIABATIC
00609         qa[k][j][i].E  /= (Real)(N*N);
00610 #endif
00611 
00612 /* Finally initialize the cell center field components. */
00613 #ifdef MHD
00614         qa[k][j][i].B1c = 0.5*(aB1i[k][j][i] + aB1i[k  ][j  ][i+1]);
00615         qa[k][j][i].B2c = 0.5*(aB2i[k][j][i] + aB2i[k  ][j+1][i  ]);
00616         qa[k][j][i].B3c = 0.5*(aB3i[k][j][i] + aB3i[k+1][j  ][i  ]);
00617 #endif /* MHD */
00618       }
00619     }
00620   }
00621 
00622 /*--- Step 6 -------------------------------------------------------------------
00623  * Now initialize variables over the whole grid
00624  */
00625 
00626   for (k=0; k<=(pGrid->ke)+nghost; k++) {
00627     for (j=0; j<=(pGrid->je)+nghost; j++) {
00628       for (i=0; i<=(pGrid->ie)+nghost; i++) {
00629 
00630 /* Calculate the coordinate of this fluid element and the mapped coordinate to
00631  * the equivalent fluid element on the interval (0 < iy <= ry) and
00632  * (0 < iz <= rz). */
00633 
00634         mix = ix = (i-(pGrid->is)) + pGrid->Disp[0];
00635         mjx = jx = (j-(pGrid->js)) + pGrid->Disp[1];
00636         mkx = kx = (k-(pGrid->ks)) + pGrid->Disp[2];
00637 
00638         n=0;
00639         while(mjx < 0){
00640           mjx += ry;
00641           n++;
00642         }
00643         if(n){
00644           mix -= n*rx;
00645         }
00646 
00647         n=0;
00648         while(mjx >= ry){
00649           mjx -= ry;
00650           n++;
00651         }
00652         if(n){
00653           mix += n*rx;
00654         }
00655 
00656         n=0;
00657         while(mkx < 0){
00658           mkx += rz;
00659           n++;
00660         }
00661         if(n){
00662           mix -= n*rx;
00663         }
00664 
00665         n=0;
00666         while(mkx >= rz){
00667           mkx -= rz;
00668           n++;
00669         }
00670         if(n){
00671           mix += n*rx;
00672         }
00673 
00674         if(mix < qa_min_ix){ /* This cell is in the left state */
00675           pGrid->U[k][j][i] = ql;
00676 #ifdef MHD
00677           pGrid->B1i[k][j][i] = ql.B1c;
00678           pGrid->B2i[k][j][i] = ql.B2c;
00679           pGrid->B3i[k][j][i] = ql.B3c;
00680 #endif
00681         }
00682         else if(mix >= qa_max_ix){ /* This cell is in the right state */
00683           pGrid->U[k][j][i] = qr;
00684 #ifdef MHD
00685           pGrid->B1i[k][j][i] = qr.B1c;
00686           pGrid->B2i[k][j][i] = qr.B2c;
00687           pGrid->B3i[k][j][i] = qr.B3c;
00688 #endif
00689         }
00690         else{ /* Otherwise it's in the qa[][][] array */
00691           pGrid->U[k][j][i] = qa[mkx][mjx][mix - qa_min_ix];
00692 #ifdef MHD
00693           pGrid->B1i[k][j][i] = aB1i[mkx][mjx][mix - qa_min_ix];
00694           pGrid->B2i[k][j][i] = aB2i[mkx][mjx][mix - qa_min_ix];
00695           pGrid->B3i[k][j][i] = aB3i[mkx][mjx][mix - qa_min_ix];
00696 #endif
00697         }
00698       }
00699     }
00700   }
00701 
00702 /*--- Step 7 -------------------------------------------------------------------
00703  * set function pointers for BCs, and conclude
00704  */
00705 
00706   bvals_mhd_fun(pDomain, left_x1,  lx_bc);
00707   bvals_mhd_fun(pDomain, right_x1, rx_bc);
00708   bvals_mhd_fun(pDomain, left_x2,  ly_bc);
00709   bvals_mhd_fun(pDomain, right_x2, ry_bc);
00710   bvals_mhd_fun(pDomain, left_x3,  lz_bc);
00711   bvals_mhd_fun(pDomain, right_x3, rz_bc);
00712 
00713   free_3d_array((void***)sqa );  sqa  = NULL;
00714 #ifdef MHD
00715   free_3d_array((void***)sB1i);  sB1i = NULL;
00716   free_3d_array((void***)sB2i);  sB2i = NULL;
00717   free_3d_array((void***)sB3i);  sB3i = NULL;
00718 #endif /* MHD */
00719 
00720   return;
00721 }
00722 
00723 
00724 /*==============================================================================
00725  * PROBLEM USER FUNCTIONS:
00726  * problem_write_restart() - writes problem-specific user data to restart files
00727  * problem_read_restart()  - reads problem-specific user data from restart files
00728  * get_usr_expr()          - sets pointer to expression for special output data
00729  * get_usr_out_fun()       - returns a user defined output function pointer
00730  * get_usr_par_prop()      - returns a user defined particle selection function
00731  * Userwork_in_loop        - problem specific work IN     main loop
00732  * Userwork_after_loop     - problem specific work AFTER  main loop
00733  *----------------------------------------------------------------------------*/
00734 
00735 void problem_write_restart(MeshS *pM,FILE *fp)
00736 {
00737   return;
00738 }
00739 
00740 void problem_read_restart(MeshS *pM, FILE *fp)
00741 {
00742   return;
00743 }
00744 
00745 ConsFun_t get_usr_expr(const char *expr)
00746 {
00747   return NULL;
00748 }
00749 
00750 VOutFun_t get_usr_out_fun(const char *name){
00751   return NULL;
00752 }
00753 
00754 void Userwork_in_loop(MeshS *pM)
00755 {
00756   return;
00757 }
00758 
00759 void Userwork_after_loop(MeshS *pM)
00760 {
00761   return;
00762 }
00763 
00764 /*=========================== PRIVATE FUNCTIONS ==============================*/
00765 
00766 /*----------------------------------------------------------------------------*/
00767 /*! \fn static void lx_bc(GridS *pG)
00768  *  \brief Apply boundary condition in left-x direction
00769  */
00770 
00771 static void lx_bc(GridS *pG)
00772 {
00773   int i, j, k, mi, mj, mk, is, js, je, ks, ke;
00774   is = pG->is;
00775   js = pG->js; je = pG->je;
00776   ks = pG->ks; ke = pG->ke;
00777 
00778   for(i=is-1; i>=0; i--){ /* Do NOT Change this loop ordering! */
00779     for(k = 0; k <= ke+nghost; k++){
00780       for(j = 0; j <= je+nghost; j++){
00781 /* Test the nearest "2D" translation vectors */
00782         if(k - rz >= ks){
00783           mi = i + rx;
00784           mj = j;
00785           mk = k - rz;
00786         }
00787         else if(j - ry >= js){
00788           mi = i + rx;
00789           mj = j - ry;
00790           mk = k;
00791         }
00792         else continue;
00793 
00794         pG->U[k][j][i] = pG->U[mk][mj][mi];
00795 #ifdef MHD
00796         pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00797         pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00798         pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00799 #endif
00800       }
00801     }
00802   }
00803 
00804   return;
00805 }
00806 
00807 /*----------------------------------------------------------------------------*/
00808 /*! \fn static void rx_bc(GridS *pG)
00809  *  \brief Apply boundary condition in right-x direction
00810  */
00811 
00812 static void rx_bc(GridS *pG)
00813 {
00814   int i, j, k, mi, mj, mk, ie, js, je, ks, ke;
00815   ie = pG->ie;
00816   js = pG->js; je = pG->je;
00817   ks = pG->ks; ke = pG->ke;
00818 
00819   for(i=ie+1; i<=ie+nghost; i++){ /* Do NOT Change this loop ordering! */
00820     for(k = 0; k <= ke+nghost; k++){
00821       for(j = 0; j <= je+nghost; j++){
00822 /* Test the nearest "2D" translation vectors */
00823         if(k + rz <= ke){
00824           mi = i - rx;
00825           mj = j;
00826           mk = k + rz;
00827         }
00828         else if(j + ry <= je){
00829           mi = i;
00830           mj = j;
00831           mk = k;
00832         }
00833         else continue;
00834 
00835         pG->U[k][j][i] = pG->U[mk][mj][mi];
00836 #ifdef MHD
00837         if(i > ie+1) pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00838         pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00839         pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00840 #endif
00841       }
00842     }
00843   }
00844 
00845   return;
00846 }
00847 
00848 /*----------------------------------------------------------------------------*/
00849 /*! \fn static void ly_bc(GridS *pG) 
00850  *  \brief Apply boundary condition in left-y direction
00851  */
00852 
00853 static void ly_bc(GridS *pG)
00854 {
00855   int i, j, k, mi, mj, mk, is, ie, js, ks, ke;
00856   is = pG->is; ie = pG->ie;
00857   js = pG->js;
00858   ks = pG->ks; ke = pG->ke;
00859 
00860   for(j=js-1; j>=0; j--){
00861     for(k = 0; k <= ke+nghost; k++){
00862       for(i = 0; i <= ie+nghost; i++){
00863 /* Test the nearest "2D" translation vectors */
00864         if(i - rx >= is){
00865           mi = i - rx;
00866           mj = j + ry;
00867           mk = k;
00868         }
00869         else if(k - rz >= ks){
00870           mi = i;
00871           mj = j + ry;
00872           mk = k - rz;
00873         }
00874         else continue;
00875 
00876         pG->U[k][j][i] = pG->U[mk][mj][mi];
00877 #ifdef MHD
00878         pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00879         pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00880         pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00881 #endif
00882       }
00883     }
00884   }
00885 
00886   return;
00887 }
00888 
00889 /*----------------------------------------------------------------------------*/
00890 /*! \fn static void ry_bc(GridS *pG)
00891  *  \brief  Apply boundary condition in right-y direction
00892  */
00893 
00894 static void ry_bc(GridS *pG)
00895 {
00896   int i, j, k, mi, mj, mk, is, ie, je, ks, ke;
00897   is = pG->is; ie = pG->ie;
00898   je = pG->je;
00899   ks = pG->ks; ke = pG->ke;
00900 
00901   for(j=je+1; j<=je+nghost; j++){
00902     for(k = 0; k <= ke+nghost; k++){
00903       for(i = 0; i <= ie+nghost; i++){
00904 /* Test the nearest "2D" translation vectors */
00905         if(i + rx <= ie){
00906           mi = i + rx;
00907           mj = j - ry;
00908           mk = k;
00909         }
00910         else if(k + rz <= ke){
00911           mi = i;
00912           mj = j - ry;
00913           mk = k + rz;
00914         }
00915         else continue;
00916 
00917         pG->U[k][j][i] = pG->U[mk][mj][mi];
00918 #ifdef MHD
00919         pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00920         if(j > je+1) pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00921         pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00922 #endif
00923       }
00924     }
00925   }
00926 
00927   return;
00928 }
00929 
00930 /*----------------------------------------------------------------------------*/
00931 /*! \fn static void lz_bc(GridS *pG)
00932  *  \brief Apply boundary condition in left-z direction
00933  */
00934 
00935 static void lz_bc(GridS *pG)
00936 {
00937   int i, j, k, mi, mj, mk, is, ie, js, je, ks;
00938   is = pG->is; ie = pG->ie;
00939   js = pG->js; je = pG->je;
00940   ks = pG->ks;
00941 
00942   for(k=ks-1; k>=0; k--){
00943     for(j = 0; j <= je+nghost; j++){
00944       for(i = 0; i <= ie+nghost; i++){
00945 /* Test the nearest "2D" translation vectors */
00946         if(i - rx >= is){
00947           mi = i - rx;
00948           mj = j;
00949           mk = k + rz;
00950         }
00951         else if(j - ry >= js){
00952           mi = i;
00953           mj = j - ry;
00954           mk = k + rz;
00955         }
00956         else continue;
00957 
00958         pG->U[k][j][i] = pG->U[mk][mj][mi];
00959 #ifdef MHD
00960         pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00961         pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00962         pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00963 #endif
00964       }
00965     }
00966   }
00967 
00968   return;
00969 }
00970 
00971 /*---------------------------------------------------------------------------*/
00972 /*! \fn static void rz_bc(GridS *pG)
00973  *  \brief Apply boundary condition in right-z direction
00974  */
00975 
00976 static void rz_bc(GridS *pG)
00977 {
00978   int i, j, k, mi, mj, mk, is, ie, js, je, ke;
00979   is = pG->is; ie = pG->ie;
00980   js = pG->js; je = pG->je;
00981   ke = pG->ke;
00982 
00983   for(k=ke+1; k<=ke+nghost; k++){
00984     for(j = 0; j <= je+nghost; j++){
00985       for(i = 0; i <= ie+nghost; i++){
00986 /* Test the nearest "2D" translation vectors */
00987         if(i + rx <= ie){
00988           mi = i + rx;
00989           mj = j;
00990           mk = k - rz;
00991         }
00992         else if(j + ry <= je){
00993           mi = i;
00994           mj = j + ry;
00995           mk = k - rz;
00996         }
00997         else continue;
00998 
00999         pG->U[k][j][i] = pG->U[mk][mj][mi];
01000 #ifdef MHD
01001         pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
01002         pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
01003         if(k > ke+1) pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
01004 #endif
01005       }
01006     }
01007   }
01008 
01009   return;
01010 }
01011 
01012 /*-----------------------------------------------------------------------------
01013  * Ax,Ay,Az: x,y,z-components of vector potential
01014  */
01015 
01016 #ifdef MHD
01017 /*! \fn static Real Ax(const Real x, const Real y, const Real z)
01018  *  \brief s-component of vector potential */
01019 static Real Ax(const Real x, const Real y, const Real z)
01020 {
01021   Real x1 = x*cos_a2*cos_a3 + y*cos_a2*sin_a3 + z*sin_a2;
01022   Real A2 =  x1*(x1 < 0.0 ? Bzl : Bzr);
01023   Real A3 = -x1*(x1 < 0.0 ? Byl : Byr);
01024 
01025   return -A2*sin_a3 - A3*sin_a2*cos_a3;
01026 }
01027 
01028 /*! \fn static Real Ay(const Real x, const Real y, const Real z)
01029  *  \brief y-component of vector potential */
01030 static Real Ay(const Real x, const Real y, const Real z)
01031 {
01032   Real x1 = x*cos_a2*cos_a3 + y*cos_a2*sin_a3 + z*sin_a2;
01033   Real A2 =  x1*(x1 < 0.0 ? Bzl : Bzr);
01034   Real A3 = -x1*(x1 < 0.0 ? Byl : Byr);
01035 
01036   return A2*cos_a3 - A3*sin_a2*sin_a3;
01037 }
01038 
01039 /*! \fn static Real Az(const Real x, const Real y, const Real z) 
01040  *  \brief z-component of vector potential */
01041 static Real Az(const Real x, const Real y, const Real z)
01042 {
01043   Real x1 = x*cos_a2*cos_a3 + y*cos_a2*sin_a3 + z*sin_a2;
01044 /* Real A2 =  x1*(x1 < 0.0 ? Bzl : Bzr); */
01045   Real A3 = -x1*(x1 < 0.0 ? Byl : Byr);
01046 
01047   return A3*cos_a2;
01048 }
01049 #endif /* MHD */