Athena: prob/shkset2d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file shkset2d.c
00004  *  \brief Sets up shock at angle to grid to test multidimensional algorithm.
00005  *
00006  * PURPOSE: Sets up shock at angle to grid to test multidimensional algorithm.
00007  *   The grid angle atan(Ly/Lx) is fixed to be atan(0.5), or atan(1), and 
00008  *   Nx1/Nx2 must be the same ratio as Lx/Ly.  Uses the angle of the shock to
00009  *   remap ghost cells to the equivalent active grid cells, which requires
00010  *   that Nx1>32, using special function pointers.  The shock is initialized
00011  *   with reference to a coordinate system (x,y,z) with transformation rules to
00012  *   the code coordinate system (x1,x2,x3)
00013  *   -  x =  x1*cos(alpha) + x2*sin(alpha)
00014  *   -  y = -x1*sin(alpha) + x2*cos(alpha)
00015  *   -  z = x3
00016 
00017  *   This inverts to:
00018  *   -  x1 = x*cos(alpha) - y*sin(alpha)
00019  *   -  x2 = x*sin(alpha) + y*cos(alpha)
00020  *   -  x3 = z                                                            
00021  *
00022  * PRIVATE FUNCTION PROTOTYPES:
00023  * - shkset2d_iib() - sets BCs on L-x1 (left edge) of grid.
00024  * - shkset2d_oib() - sets BCs on R-x1 (right edge) of grid.
00025  * - shkset2d_ijb() - sets BCs on L-x2 (bottom edge) of grid.
00026  * - shkset2d_ojb() - sets BCs on R-x2 (top edge) of grid.                    */
00027 /*============================================================================*/
00028 
00029 #include <math.h>
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include "defs.h"
00033 #include "athena.h"
00034 #include "globals.h"
00035 #include "prototypes.h"
00036 
00037 /*==============================================================================
00038  * PRIVATE FUNCTION PROTOTYPES:
00039  * shkset2d_iib() - sets BCs on L-x1 (left edge) of grid.
00040  * shkset2d_oib() - sets BCs on R-x1 (right edge) of grid.
00041  * shkset2d_ijb() - sets BCs on L-x2 (bottom edge) of grid.
00042  * shkset2d_ojb() - sets BCs on R-x2 (top edge) of grid.
00043  *============================================================================*/
00044 
00045 void shkset2d_iib(GridS *pGrid);
00046 void shkset2d_oib(GridS *pGrid);
00047 void shkset2d_ijb(GridS *pGrid);
00048 void shkset2d_ojb(GridS *pGrid);
00049 
00050 /* Make size of box and dimension of unit cell (r1 x r2) static globals so they
00051  * can be accessed by boundary value functions */
00052 static Real Lx,Ly;
00053 static int r1,r2;
00054 
00055 /*----------------------------------------------------------------------------*/
00056 /* problem:   */
00057 
00058 void problem(DomainS *pDomain)
00059 {
00060   GridS *pGrid = pDomain->Grid;
00061   int i, is = pGrid->is, ie = pGrid->ie;
00062   int j, js = pGrid->js, je = pGrid->je;
00063   int k, ks = pGrid->ks, ke = pGrid->ke;
00064   int kl,ku,irefine,ir,ix1,ix2,nx1,nx2,gcd;
00065   Real angle, sin_a, cos_a; /* Angle the shock makes with the x1-direction */
00066   Real rootdx1, rootdx2;
00067   Prim1DS Wl, Wr;
00068   Cons1DS Ul, Ur;
00069   ConsS ql, qr;
00070   Real Bxl=0.0,Bxr=0.0;
00071   div_t id;   /* structure containing remainder and quotient */
00072 
00073 /* Following are used to compute volume of cell crossed by initial interface
00074  * that is assigned to left/right states */
00075   int dll, dlr, drr, drl;
00076   Real afl_lx, afr_lx, afl_rx, afr_rx;
00077   Real afl_ly, afr_ly, afl_ry, afr_ry;
00078   Real vfl, vfr, B1r, B2r;
00079 
00080   if (pGrid->Nx[1] == 1)
00081     ath_error("[shkset2d]: This problem can only be run in 2D or 3D\n");
00082 
00083   if (pGrid->Nx[2] > 1){
00084     ku = pGrid->ke + nghost;
00085     kl = pGrid->ks - nghost;
00086   } else {
00087     ku = pGrid->ke;
00088     kl = pGrid->ks;
00089   }
00090 
00091 /* Find number of cells on root grid */
00092 
00093   irefine = 1;
00094   for (ir=1;ir<=pDomain->Level;ir++) irefine *= 2;
00095 
00096   Lx = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00097   Ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00098 
00099   rootdx1 = pGrid->dx1*((double)(irefine)); 
00100   rootdx2 = pGrid->dx2*((double)(irefine)); 
00101 
00102   nx1 = (int)(Lx/rootdx1);
00103   nx2 = (int)(Ly/rootdx2);
00104 
00105 /* Compute greatest common divisor of nx1,nx2.  The size of the "unit cell"
00106  * is nx1/gcd by nx2/gcd */
00107 
00108   if((gcd = ath_gcd(nx1,nx2)) < 10)
00109     ath_error("[shkset2d]: Greatest Common Divisor (nx1,nx2) = %d\n",gcd);
00110 
00111   id = div(nx1,gcd);
00112   r1 = id.quot;
00113   if(id.rem != 0)
00114     ath_error("[shkset2d]: GCD failed, Remainder of %d / %d is %d\n",
00115               nx1,gcd,id.rem);
00116 
00117   id = div(nx2,gcd);
00118   r2 = id.quot;
00119   if(id.rem != 0)
00120     ath_error("[shkset2d]: GCD failed, Remainder of %d / %d is %d\n",
00121               nx2,gcd,id.rem);
00122 
00123   ath_pout(1,"The unit cell is (%d,1,%d) grid cells in size\n",r1,r2);
00124 
00125 /* Compute angle initial interface makes to the grid */
00126 
00127   if(Lx == Ly){
00128     cos_a = sin_a = sqrt(0.5);
00129   }
00130   else{
00131     angle = atan((double)(Lx/Ly));
00132     sin_a = sin(angle);
00133     cos_a = cos(angle);
00134   }
00135 
00136 /* Parse left state read from input file: dl,pl,ul,vl,wl,bxl,byl,bzl */
00137 
00138   Wl.d = par_getd("problem","dl");
00139 #ifdef ADIABATIC
00140   Wl.P = par_getd("problem","pl");
00141 #endif
00142   Wl.Vx = par_getd("problem","v1l");
00143   Wl.Vy = par_getd("problem","v2l");
00144   Wl.Vz = par_getd("problem","v3l");
00145 #ifdef MHD
00146   Bxl = par_getd("problem","b1l");
00147   Wl.By = par_getd("problem","b2l");
00148   Wl.Bz = par_getd("problem","b3l");
00149 #endif
00150 
00151 /* Parse right state read from input file: dr,pr,ur,vr,wr,bxr,byr,bzr */
00152 
00153   Wr.d = par_getd("problem","dr");
00154 #ifdef ADIABATIC
00155   Wr.P = par_getd("problem","pr");
00156 #endif
00157   Wr.Vx = par_getd("problem","v1r");
00158   Wr.Vy = par_getd("problem","v2r");
00159   Wr.Vz = par_getd("problem","v3r");
00160 #ifdef MHD
00161   Bxr = par_getd("problem","b1r");
00162   Wr.By = par_getd("problem","b2r");
00163   Wr.Bz = par_getd("problem","b3r");
00164   if (Bxr != Bxl) ath_error(0,"[shkset2d] L/R values of Bx not the same\n");
00165 #endif
00166 
00167   Ul = Prim1D_to_Cons1D(&Wl, &Bxl);
00168   Ur = Prim1D_to_Cons1D(&Wr, &Bxr);
00169 
00170 /* Initialize ql rotated to the (x1,x2,x3) coordinate system */
00171   ql.d   = Ul.d;
00172   ql.M1  = Ul.Mx*cos_a - Ul.My*sin_a;
00173   ql.M2  = Ul.Mx*sin_a + Ul.My*cos_a;
00174   ql.M3  = Ul.Mz;
00175 #ifdef MHD
00176   ql.B1c = Bxl*cos_a - Ul.By*sin_a;
00177   ql.B2c = Bxl*sin_a + Ul.By*cos_a;
00178   ql.B3c = Ul.Bz;
00179 #endif
00180 #ifdef ADIABATIC
00181   ql.E   = Ul.E;
00182 #endif
00183 
00184 /* Initialize qr rotated to the (x1,x2,x3) coordinate system */
00185   qr.d   = Ur.d;
00186   qr.M1  = Ur.Mx*cos_a - Ur.My*sin_a;
00187   qr.M2  = Ur.Mx*sin_a + Ur.My*cos_a;
00188   qr.M3  = Ur.Mz;
00189 #ifdef MHD
00190   qr.B1c = Bxr*cos_a - Ur.By*sin_a;
00191   qr.B2c = Bxr*sin_a + Ur.By*cos_a;
00192   qr.B3c = Ur.Bz;
00193 #endif
00194 #ifdef ADIABATIC
00195   qr.E   = Ur.E;
00196 #endif
00197 
00198 /* Initialize the grid */
00199 
00200   for (k=kl; k<=ku; k++) {
00201     for (j=0; j<=je+nghost; j++) {
00202       ix2 = j + pGrid->Disp[1];
00203       for (i=0; i<=ie+nghost; i++) {
00204         ix1 = i + pGrid->Disp[0];
00205 
00206 /* cell is completely in the left state */
00207         if((drr = r2*(ix1) + r1*(ix2) - gcd*r1*r2) <= 0){
00208           pGrid->U[k][j][i] = ql;
00209 #ifdef MHD
00210           pGrid->B1i[k][j][i] = ql.B1c;
00211           pGrid->B2i[k][j][i] = ql.B2c;
00212           pGrid->B3i[k][j][i] = ql.B3c;
00213 #endif /* MHD */
00214         }
00215 /* cell is completely in the right state */
00216         else if((dll = r2*(ix1-1) + r1*(ix2-1) - gcd*r1*r2) >= 0){
00217           pGrid->U[k][j][i] = qr;
00218 #ifdef MHD
00219           pGrid->B1i[k][j][i] = qr.B1c;
00220           pGrid->B2i[k][j][i] = qr.B2c;
00221           pGrid->B3i[k][j][i] = qr.B3c;
00222 #endif /* MHD */
00223         }
00224 /* The more complicated case of a cell  split by the interface boundary */
00225         else{
00226           dlr = r2*(ix1-1) + r1*(ix2) - gcd*r1*r2;
00227 
00228           if(dlr < 0){ /* The boundary hits the right y-face */
00229             afl_lx = 1.0;
00230             afr_lx = 0.0;
00231             afl_ry = (Real)(-dlr)/(Real)(r2);
00232             afr_ry = 1.0 - afl_ry;
00233           }
00234           else if(dlr > 0){ /* The boundary hits the left x-face */
00235             afl_lx = (Real)(-dll)/(Real)(r1);
00236             afr_lx = 1.0 - afl_lx;
00237             afl_ry = 0.0;
00238             afr_ry = 1.0;
00239           }
00240           else{ /* dlr == 0.0, The boundary hits the grid cell corner */
00241             afl_lx = 1.0;
00242             afr_lx = 0.0;
00243             afl_ry = 0.0;
00244             afr_ry = 1.0;
00245           }
00246 
00247           drl = r2*(ix1) + r1*(ix2-1) - gcd*r1*r2;
00248 
00249           if(drl < 0){ /* The boundary hits the right x-face */
00250             afl_rx = (Real)(-drl)/(Real)(r1);
00251             afr_rx = 1.0 - afl_rx;
00252             afl_ly = 1.0;
00253             afr_ly = 0.0;
00254           }
00255           else if(drl > 0){ /* The boundary hits the left y-face */
00256             afl_rx = 0.0;
00257             afr_rx = 1.0;
00258             afl_ly = (Real)(-dll)/(Real)(r2);
00259             afr_ly = 1.0 - afl_ly;
00260           }
00261           else{ /* drl == 0.0, The boundary hits the grid cell corner */
00262             afl_rx = 0.0;
00263             afr_rx = 1.0;
00264             afl_ly = 1.0;
00265             afr_ly = 0.0;
00266           }
00267 
00268 /* The boundary hits both x-interfaces */
00269           if(dlr > 0 && drl < 0){ 
00270             vfl = 0.5*(afl_lx + afl_rx);
00271             vfr = 1.0 - vfl;
00272           }
00273 /* The boundary hits both y-interfaces */
00274           else if(dlr < 0 && drl > 0){ 
00275             vfl = 0.5*(afl_ly + afl_ry);
00276             vfr = 1.0 - vfl;
00277           }
00278 /* The boundary hits both grid cell corners */
00279           else if(dlr == 0 && drl == 0){ 
00280             vfl = vfr = 0.5;
00281           }
00282 /* The boundary hits the left x- and left y-interface */
00283           else if(dlr > 0 && drl > 0){
00284             vfl = 0.5*afl_lx*afl_ly;
00285             vfr = 1.0 - vfl;
00286           }
00287 /* dlr<0 && drl<0:  The boundary hits the right x- and right y-interface */
00288           else{ 
00289             vfr = 0.5*afr_rx*afr_ry;
00290             vfl = 1.0 - vfr;
00291           }
00292 
00293 /* Initialize the x- and y-interface magnetic fields */
00294 #ifdef MHD
00295           pGrid->B1i[k][j][i] = afl_lx*ql.B1c + afr_lx*qr.B1c;
00296           B1r              = afl_rx*ql.B1c + afr_rx*qr.B1c;
00297 
00298           pGrid->B2i[k][j][i] = afl_ly*ql.B2c + afr_ly*qr.B2c;
00299           B2r              = afl_ry*ql.B2c + afr_ry*qr.B2c;
00300 
00301           pGrid->B3i[k][j][i] = vfl*ql.B3c + vfr*qr.B3c;
00302 #endif /* MHD */
00303 
00304 /* Initialize the volume averaged quantities */
00305           pGrid->U[k][j][i].d  = vfl*ql.d + vfr*qr.d;
00306           pGrid->U[k][j][i].M1 = vfl*ql.M1 + vfr*qr.M1;
00307           pGrid->U[k][j][i].M2 = vfl*ql.M2 + vfr*qr.M2;
00308           pGrid->U[k][j][i].M3 = vfl*ql.M3 + vfr*qr.M3;
00309 #ifdef MHD
00310           pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i] + B1r);
00311           pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i] + B2r);
00312           pGrid->U[k][j][i].B3c = vfl*ql.B3c + vfr*qr.B3c;
00313 #endif /* MHD */
00314 #ifndef ISOTHERMAL
00315           pGrid->U[k][j][i].E  = vfl*ql.E + vfr*qr.E;
00316 #endif
00317         }
00318       }
00319     }
00320   }
00321 
00322 /* Set boundary value function pointers */
00323 
00324   bvals_mhd_fun(pDomain, left_x1,  shkset2d_iib);
00325   bvals_mhd_fun(pDomain, left_x2,  shkset2d_ijb);
00326   bvals_mhd_fun(pDomain, right_x1, shkset2d_oib);
00327   bvals_mhd_fun(pDomain, right_x2, shkset2d_ojb);
00328 
00329   return;
00330 }
00331 
00332 /*==============================================================================
00333  * PROBLEM USER FUNCTIONS:
00334  * problem_write_restart() - writes problem-specific user data to restart files
00335  * problem_read_restart()  - reads problem-specific user data from restart files
00336  * get_usr_expr()          - sets pointer to expression for special output data
00337  * get_usr_out_fun()       - returns a user defined output function pointer
00338  * get_usr_par_prop()      - returns a user defined particle selection function
00339  * Userwork_in_loop        - problem specific work IN     main loop
00340  * Userwork_after_loop     - problem specific work AFTER  main loop
00341  *----------------------------------------------------------------------------*/
00342 
00343 void problem_write_restart(MeshS *pM, FILE *fp)
00344 {
00345   return;
00346 }
00347 
00348 void problem_read_restart(MeshS *pM, FILE *fp)
00349 {
00350   return;
00351 }
00352 
00353 ConsFun_t get_usr_expr(const char *expr)
00354 {
00355   return NULL;
00356 }
00357 
00358 VOutFun_t get_usr_out_fun(const char *name){
00359   return NULL;
00360 }
00361 
00362 void Userwork_in_loop(MeshS *pM)
00363 {
00364 }
00365 
00366 void Userwork_after_loop(MeshS *pM)
00367 {
00368 }
00369 
00370 /*=========================== PRIVATE FUNCTIONS ==============================*/
00371 
00372 /*----------------------------------------------------------------------------*/
00373 /*! \fn void shkset2d_iib(GridS *pGrid)
00374  *  \brief Sets ghost zones using the nearest computational grid
00375  * cells implied by the size of the unit cell (r1xr2).
00376  */
00377 
00378 void shkset2d_iib(GridS *pGrid)
00379 {
00380   const int is = pGrid->is;
00381   int i, j, k, ju, jl, kl, ku; /* j-upper, j-lower */
00382 
00383   if (pGrid->Nx[1] > 1){
00384     ju = pGrid->je + nghost;
00385     jl = pGrid->js - nghost + r2;
00386   } else {
00387     ju = pGrid->je;
00388     jl = pGrid->js;
00389   }
00390 
00391   if (pGrid->Nx[2] > 1){
00392     ku = pGrid->ke + nghost;
00393     kl = pGrid->ks - nghost;
00394   } else {
00395     ku = pGrid->ke;
00396     kl = pGrid->ks;
00397   }
00398 
00399   for (k=kl; k<=ku; k++) {
00400     for (i=1; i<=nghost; i++) { /* Do NOT Change this loop ordering! */
00401       for (j=jl; j<=ju; j++) {
00402         pGrid->U  [k][j][is-i] = pGrid->U  [k][j-r2][is-i+r1];
00403 #ifdef MHD
00404         pGrid->B1i[k][j][is-i] = pGrid->B1i[k][j-r2][is-i+r1];
00405         pGrid->B2i[k][j][is-i] = pGrid->B2i[k][j-r2][is-i+r1];
00406         pGrid->B3i[k][j][is-i] = pGrid->B3i[k][j-r2][is-i+r1];
00407 #endif
00408       }
00409     }
00410   }
00411   return;
00412 }
00413 
00414 /*----------------------------------------------------------------------------*/
00415 /*! \fn void shkset2d_oib(GridS *pGrid)
00416  *  \brief Sets ghost zones using the nearest computational grid
00417  * cells implied by the size of the unit cell (r1xr2).
00418  */
00419 
00420 void shkset2d_oib(GridS *pGrid)
00421 {
00422   const int ie = pGrid->ie;
00423   int i, j, k, ju, jl, kl, ku; /* j-upper, j-lower */
00424 
00425   if (pGrid->Nx[1] > 1){
00426     ju = pGrid->je + nghost - r2;
00427     jl = pGrid->js - nghost;
00428   } else {
00429     ju = pGrid->je;
00430     jl = pGrid->js;
00431   }
00432 
00433   if (pGrid->Nx[2] > 1){
00434     ku = pGrid->ke + nghost;
00435     kl = pGrid->ks - nghost;
00436   } else {
00437     ku = pGrid->ke;
00438     kl = pGrid->ks;
00439   }
00440 
00441 /* Note that i=ie+1 is not a boundary condition for the interface field B1i */
00442 
00443   for (k=kl; k<=ku; k++) {
00444     for (i=1; i<=nghost; i++) { /* Do NOT Change this loop ordering! */
00445       for (j=jl; j<=ju; j++) {
00446         pGrid->U[k][j][ie+i] = pGrid->U[k][j+r2][ie+i-r1];
00447 #ifdef MHD
00448         if(i>1) pGrid->B1i[k][j][ie+i] = pGrid->B1i[k][j+r2][ie+i-r1];
00449         pGrid->B2i[k][j][ie+i] = pGrid->B2i[k][j+r2][ie+i-r1];
00450         pGrid->B3i[k][j][ie+i] = pGrid->B3i[k][j+r2][ie+i-r1];
00451 #endif
00452       }
00453     }
00454   }
00455   return;
00456 }
00457 
00458 /*----------------------------------------------------------------------------*/
00459 /*! \fn void shkset2d_ijb(GridS *pGrid)
00460  *  \brief Sets ghost zones using the nearest computational grid
00461  * cells implied by the size of the unit cell (r1xr2).
00462  */
00463 
00464 void shkset2d_ijb(GridS *pGrid)
00465 {
00466   const int js = pGrid->js;
00467   int i, j, k, iu, il, kl, ku; /* i-upper, i-lower */
00468 
00469   iu = pGrid->ie + nghost;
00470   il = pGrid->is - nghost + r1;
00471 
00472   if (pGrid->Nx[2] > 1){
00473     ku = pGrid->ke + nghost;
00474     kl = pGrid->ks - nghost;
00475   } else {
00476     ku = pGrid->ke;
00477     kl = pGrid->ks;
00478   }
00479 
00480   for (k=kl; k<=ku; k++) {
00481     for (j=1; j<=nghost; j++) {
00482       for (i=il; i<=iu; i++) {
00483         pGrid->U  [k][js-j][i] = pGrid->U  [k][js-j+r2][i-r1];
00484 #ifdef MHD
00485         pGrid->B1i[k][js-j][i] = pGrid->B1i[k][js-j+r2][i-r1];
00486         pGrid->B2i[k][js-j][i] = pGrid->B2i[k][js-j+r2][i-r1];
00487         pGrid->B3i[k][js-j][i] = pGrid->B3i[k][js-j+r2][i-r1];
00488 #endif
00489       }
00490     }
00491   }
00492   return;
00493 }
00494 
00495 /*----------------------------------------------------------------------------*/
00496 /*! \fn void shkset2d_ojb(GridS *pGrid)
00497  *  \brief Sets ghost zones using the nearest computational grid
00498  * cells implied by the size of the unit cell (r1xr2).
00499  */
00500 
00501 void shkset2d_ojb(GridS *pGrid)
00502 {
00503   const int je = pGrid->je;
00504   int i, j, k, iu, il, kl, ku; /* i-upper, i-lower */
00505 
00506   iu = pGrid->ie + nghost - r1;
00507   il = pGrid->is - nghost;
00508 
00509   if (pGrid->Nx[2] > 1){
00510     ku = pGrid->ke + nghost;
00511     kl = pGrid->ks - nghost;
00512   } else {
00513     ku = pGrid->ke;
00514     kl = pGrid->ks;
00515   }
00516 
00517 /* Note that j=je+1 is not a boundary condition for the interface field B2i */
00518 
00519   for (k=kl; k<=ku; k++) {
00520     for (j=1; j<=nghost; j++) {
00521       for (i=il; i<=iu; i++) {
00522         pGrid->U[k][je+j][i] = pGrid->U[k][je+j-r2][i+r1];
00523 #ifdef MHD
00524         pGrid->B1i[k][je+j][i] = pGrid->B1i[k][je+j-r2][i+r1];
00525         if(j>1) pGrid->B2i[k][je+j][i] = pGrid->B2i[k][je+j-r2][i+r1];
00526         pGrid->B3i[k][je+j][i] = pGrid->B3i[k][je+j-r2][i+r1];
00527 #endif
00528       }
00529     }
00530   }
00531   return;
00532 }