Athena: prob/shk_cloud.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file shk_cloud.c
00004  *  \brief Problem generator for shock-cloud problem; planar shock impacting
00005  *   a single spherical cloud.
00006  *
00007  * PURPOSE: Problem generator for shock-cloud problem; planar shock impacting
00008  *   a single spherical cloud.  Input parameters are:
00009  *    - problem/Mach   = Mach number of incident shock
00010  *    - problem/drat   = density ratio of cloud to ambient
00011  *    - problem/beta   = ratio of Pgas/Pmag
00012  *    - problem/iprob  = integer flag to determine problem
00013  *
00014  *   The cloud radius is fixed at 1.0.  The center of the coordinate system
00015  *   defines the center of the cloud, and should be in the middle of the cloud.
00016  *   The shock is initially at x1=-2.0.  A typical grid domain should span
00017  *   x1 in [-3.0,7.0] , y and z in [-2.5,2.5] (see input file in /tst)
00018  *   Various test cases are possible:
00019  *   - (iprob=1): B parallel to shock normal
00020  *   - (iprob=2): B perpendicular to shock normal -- NOT YET IMPLEMENTED
00021  *
00022  *   If the code is configured with nscalars>0, the cloud material is labeled
00023  *   with U[k][j][i].s[0]=1.                                                  
00024  *
00025  * PRIVATE FUNCTION PROTOTYPES:
00026  * shk_cloud_iib() - fixes BCs on L-x1 (left edge) of grid to postshock flow. */
00027 /*============================================================================*/
00028 
00029 #include <math.h>
00030 #include <stdio.h>
00031 #include "defs.h"
00032 #include "athena.h"
00033 #include "globals.h"
00034 #include "prototypes.h"
00035 
00036 /* postshock flow variables are shared with IIB function */
00037 
00038 static Real dl,pl,ul;
00039 #ifdef MHD
00040 static Real bxl,byl,bzl;
00041 #endif /* MHD */
00042 
00043 /*==============================================================================
00044  * PRIVATE FUNCTION PROTOTYPES:
00045  * shk_cloud_iib() - fixes BCs on L-x1 (left edge) of grid to postshock flow.
00046  *============================================================================*/
00047 
00048 void shk_cloud_iib(GridS *pGrid);
00049 
00050 /*=========================== PUBLIC FUNCTIONS ===============================*/
00051 /*----------------------------------------------------------------------------*/
00052 /* problem:  */
00053 
00054 void problem(DomainS *pDomain)
00055 {
00056   GridS *pGrid = pDomain->Grid;
00057   int i=0,j=0,k=0;
00058   int is,ie,js,je,ks,ke,iprob;
00059   Real jump1, jump2, jump3;
00060   Real x1,x2,x3,diag,rad,xshock;
00061   Real Mach,drat;
00062 #ifdef MHD
00063   Real beta,bxr,byr,bzr;
00064 #endif /* MHD */
00065   Real dr,pr,ur;
00066 
00067 /* Read input parameters */
00068 
00069   xshock = -2.0;
00070   rad    = 1.0;
00071   Mach = par_getd("problem","Mach");
00072   drat = par_getd("problem","drat");
00073   iprob = par_geti("problem","iprob");
00074   iprob = 1;
00075 #ifdef MHD
00076   beta = par_getd("problem","beta");
00077 #endif
00078   
00079 /* Set paramters in ambient medium ("R-state" for shock) */
00080 
00081   dr = 1.0;
00082   pr = 1.0/Gamma;
00083   ur = 0.0;
00084 
00085 /* Uses Rankine Hugoniot relations for adiabatic gas to initialize problem */
00086 
00087   jump1 = (Gamma + 1.0)/(Gamma_1 + 2.0/(Mach*Mach));
00088   jump2 = (2.0*Gamma*Mach*Mach - Gamma_1)/(Gamma + 1.0);
00089   jump3 = 2.0*(1.0 - 1.0/(Mach*Mach))/(Gamma + 1.0);
00090 
00091   dl = dr*jump1;
00092   pl = pr*jump2;
00093 #ifdef ISOTHERMAL
00094   ul = ur + jump3*Mach*Iso_csound;
00095 #else
00096   ul = ur + jump3*Mach*sqrt(Gamma*pr/dr);
00097 #endif
00098 
00099 /* Initialize the grid */
00100 
00101   is = pGrid->is;  ie = pGrid->ie;
00102   js = pGrid->js;  je = pGrid->je;
00103   ks = pGrid->ks;  ke = pGrid->ke;
00104 
00105   for (k=ks; k<=ke; k++) {
00106     for (j=js; j<=je; j++) {
00107       for (i=is; i<=ie; i++) {
00108         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00109         diag = sqrt(x1*x1 + x2*x2 + x3*x3);
00110 
00111 /* postshock flow */
00112         if(x1 < xshock) {
00113           pGrid->U[k][j][i].d  = dl;
00114           pGrid->U[k][j][i].M1 = ul*dl;
00115           pGrid->U[k][j][i].M2 = 0.0;
00116           pGrid->U[k][j][i].M3 = 0.0;
00117 #ifdef MHD
00118           pGrid->B1i[k][j][i] = bxl;
00119           pGrid->B2i[k][j][i] = byl;
00120           pGrid->B3i[k][j][i] = bzl;
00121           pGrid->U[k][j][i].B1c = bxl;
00122           pGrid->U[k][j][i].B2c = byl;
00123           pGrid->U[k][j][i].B3c = bzl;
00124 #endif
00125 #ifdef ADIABATIC
00126           pGrid->U[k][j][i].E = pl/Gamma_1 
00127 #ifdef MHD
00128             + 0.5*(bxl*bxl + byl*byl + bzl*bzl)
00129 #endif
00130             + 0.5*dl*(ul*ul);
00131 #endif
00132 #if (NSCALARS > 0)
00133           pGrid->U[k][j][i].s[0] = 0.0;
00134 #endif
00135 
00136 /* preshock ambient gas */
00137         } else {
00138           pGrid->U[k][j][i].d  = dr;
00139           pGrid->U[k][j][i].M1 = ur*dr;
00140           pGrid->U[k][j][i].M2 = 0.0;
00141           pGrid->U[k][j][i].M3 = 0.0;
00142 #ifdef MHD
00143           pGrid->B1i[k][j][i] = bxr;
00144           pGrid->B2i[k][j][i] = byr;
00145           pGrid->B3i[k][j][i] = bzr;
00146           pGrid->U[k][j][i].B1c = bxr;
00147           pGrid->U[k][j][i].B2c = byr;
00148           pGrid->U[k][j][i].B3c = bzr;
00149 #endif
00150 #ifdef ADIABATIC
00151           pGrid->U[k][j][i].E = pr/Gamma_1
00152 #ifdef MHD
00153             + 0.5*(bxr*bxr + byr*byr + bzr*bzr)
00154 #endif
00155             + 0.5*dr*(ur*ur);
00156 #endif
00157 #if (NSCALARS > 0)
00158           pGrid->U[k][j][i].s[0] = 0.0;
00159 #endif
00160         }
00161 
00162 /* cloud interior */
00163         if (diag < rad) {
00164           pGrid->U[k][j][i].d  = dr*drat;
00165           pGrid->U[k][j][i].M1 = ur*dr*drat;
00166           pGrid->U[k][j][i].M2 = 0.0;
00167           pGrid->U[k][j][i].M3 = 0.0;
00168 #ifdef MHD
00169           pGrid->B1i[k][j][i] = bxr;
00170           pGrid->B2i[k][j][i] = byr;
00171           pGrid->B3i[k][j][i] = bzr;
00172           pGrid->U[k][j][i].B1c = bxr;
00173           pGrid->U[k][j][i].B2c = byr;
00174           pGrid->U[k][j][i].B3c = bzr;
00175 #endif
00176 #ifdef ADIABATIC
00177           pGrid->U[k][j][i].E = pr/Gamma_1
00178 #ifdef MHD
00179             + 0.5*(bxr*bxr + byr*byr + bzr*bzr)
00180 #endif
00181             + 0.5*dr*drat*(ur*ur);
00182 #endif
00183 #if (NSCALARS > 0)
00184           pGrid->U[k][j][i].s[0] = drat;
00185 #endif
00186         }
00187       }
00188     }
00189   }
00190 
00191 /* boundary conditions on interface B */
00192 
00193 #ifdef MHD
00194   i = ie+1;
00195   for (k=ks; k<=ke; k++) {
00196     for (j=js; j<=je; j++) {
00197       pGrid->B1i[k][j][i] = bxr;
00198     }
00199   }
00200   j = je+1;
00201   for (k=ks; k<=ke; k++) {
00202     for (i=is; i<=ie; i++) {
00203       cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00204       if (x1 < xshock) {pGrid->B2i[k][j][i] = byl;}
00205       else {pGrid->B2i[k][j][i] = byr;}
00206     }
00207   }
00208   if (ke > ks) {
00209     k = ke+1;
00210     for (j=js; j<=je; j++) {
00211       for (i=is; i<=ie; i++) {
00212         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00213         if (x1 < xshock) {pGrid->B3i[k][j][i] = bzl;}
00214         else {pGrid->B3i[k][j][i] = bzr;}
00215       }
00216     }
00217   }
00218 #endif
00219 
00220 /* Set IIB value function pointer */
00221 
00222   if (pDomain->Disp[0] == 0) bvals_mhd_fun(pDomain,left_x1,shk_cloud_iib);
00223 
00224   return;
00225 }
00226 
00227 /*==============================================================================
00228  * PROBLEM USER FUNCTIONS:
00229  * problem_write_restart() - writes problem-specific user data to restart files
00230  * problem_read_restart()  - reads problem-specific user data from restart files
00231  * get_usr_expr()          - sets pointer to expression for special output data
00232  * get_usr_out_fun()       - returns a user defined output function pointer
00233  * get_usr_par_prop()      - returns a user defined particle selection function
00234  * Userwork_in_loop        - problem specific work IN     main loop
00235  * Userwork_after_loop     - problem specific work AFTER  main loop
00236  * color()   - returns first passively advected scalar s[0]
00237  *----------------------------------------------------------------------------*/
00238 
00239 void problem_write_restart(MeshS *pM, FILE *fp)
00240 {
00241   return;
00242 }
00243 
00244 void problem_read_restart(MeshS *pM, FILE *fp)
00245 {
00246   return;
00247 }
00248 
00249 #if (NSCALARS > 0)
00250 static Real color(const GridS *pG, const int i, const int j, const int k)
00251 {
00252   return pG->U[k][j][i].s[0]/pG->U[k][j][i].d;
00253 }
00254 #endif
00255 
00256 ConsFun_t get_usr_expr(const char *expr)
00257 {
00258 #if (NSCALARS > 0)
00259   if(strcmp(expr,"color")==0) return color;
00260 #endif
00261 
00262   return NULL;
00263 }
00264 
00265 VOutFun_t get_usr_out_fun(const char *name){
00266   return NULL;
00267 }
00268 
00269 void Userwork_in_loop(MeshS *pM)
00270 {
00271   return;
00272 }
00273 
00274 void Userwork_after_loop(MeshS *pM)
00275 {
00276   return;
00277 }
00278 
00279 /*=========================== PRIVATE FUNCTIONS ==============================*/
00280 
00281 /*----------------------------------------------------------------------------*/
00282 /*! \fn void shk_cloud_iib(GridS *pGrid)
00283  *  \brief Sets boundary condition on left X boundary (iib) 
00284  *
00285  * Note quantities at this boundary are held fixed at the downstream state
00286  */
00287 
00288 void shk_cloud_iib(GridS *pGrid)
00289 {
00290   int i=0,j=0,k=0;
00291   int js,je,ks,ke;
00292 
00293   js = pGrid->js; je = pGrid->je;
00294   ks = pGrid->ks; ke = pGrid->ke;
00295 
00296   for (k=ks; k<=ke; k++) {
00297     for (j=js; j<=je; j++) {
00298       for (i=1; i<=nghost; i++) {
00299         pGrid->U[k][j][i].d  = dl;
00300         pGrid->U[k][j][i].M1 = ul*dl;
00301         pGrid->U[k][j][i].M2 = 0.0;
00302         pGrid->U[k][j][i].M3 = 0.0;
00303 #ifdef MHD
00304         pGrid->B1i[k][j][i] = bxl;
00305         pGrid->B2i[k][j][i] = byl;
00306         pGrid->B3i[k][j][i] = bzl;
00307         pGrid->U[k][j][i].B1c = bxl;
00308         pGrid->U[k][j][i].B2c = byl;
00309         pGrid->U[k][j][i].B3c = bzl;
00310 #endif
00311 #ifdef ADIABATIC
00312         pGrid->U[k][j][i].E = pl/Gamma_1
00313 #ifdef MHD
00314           + 0.5*(bxl*bxl + byl*byl + bzl*bzl)
00315 #endif
00316           + 0.5*dl*(ul*ul);
00317 #endif
00318 #if (NSCALARS > 0)
00319         pGrid->U[k][j][i].s[0] = 0.0;
00320 #endif
00321       }
00322     }
00323   }
00324 }