Athena: prob/dmr.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file dmr.c
00004  *  \brief Problem generator for double Mach reflection test.
00005  *
00006  * PURPOSE: Problem generator for double Mach reflection test.  Only works for
00007  *   genuinely 2D problems in X1-X2 plane.
00008  *
00009  * REFERENCE: P. Woodward & P. Colella, "The numerical simulation of 
00010  *   two-dimensional fluid flow with strong shocks", JCP, 54, 115, sect. IVc. */
00011 /*============================================================================*/
00012 
00013 #include <math.h>
00014 #include <stdio.h>
00015 #include <stdlib.h>
00016 #include "defs.h"
00017 #include "athena.h"
00018 #include "globals.h"
00019 #include "prototypes.h"
00020 
00021 #ifdef ISOTHERMAL
00022 #error : The dmr test only works for adiabatic EOS.
00023 #endif /* ISOTHERMAL */
00024 #ifndef HYDRO
00025 #error : The dmr test only works for hydro.
00026 #endif /* HYDRO */
00027 
00028 /*==============================================================================
00029  * PRIVATE FUNCTION PROTOTYPES:
00030  * dmrbv_iib() - sets BCs on L-x1 (left edge) of grid.  
00031  * dmrbv_ijb() - sets BCs on L-x2 (bottom edge) of grid.  
00032  * dmrbv_ojb() - sets BCs on R-x2 (top edge) of grid.  
00033  *============================================================================*/
00034 
00035 void dmrbv_iib(GridS *pGrid);
00036 void dmrbv_ijb(GridS *pGrid);
00037 void dmrbv_ojb(GridS *pGrid);
00038 
00039 /*=========================== PUBLIC FUNCTIONS ===============================*/
00040 /*----------------------------------------------------------------------------*/
00041 /* problem:  */
00042 
00043 void problem(DomainS *pDomain)
00044 {
00045   GridS *pGrid = pDomain->Grid;
00046   int i=0,j=0;
00047   int is,ie,js,je,ks;
00048   Real d0,e0,u0,v0,x1_shock,x1,x2,x3;
00049 
00050   is = pGrid->is; ie = pGrid->ie;
00051   js = pGrid->js; je = pGrid->je;
00052   ks = pGrid->ks;
00053   if (pGrid->Nx[0] == 1 || pGrid->Nx[1] == 1) {
00054     ath_error("[dmr]: this test only works with Nx1 & Nx2 > 1\n");
00055   }
00056   if (pGrid->Nx[2] > 1) {
00057     ath_error("[dmr]: this test only works for 2D problems, with Nx3=1\n");
00058   }
00059 
00060 /* Initialize shock using parameters defined in Woodward & Colella */
00061 
00062   d0 = 8.0;
00063   e0 = 291.25;
00064   u0 =  8.25*sqrt(3.0)/2.0;
00065   v0 = -8.25*0.5;
00066   for (j=js; j<=je; j++) {
00067     for (i=is; i<=ie; i++) {
00068       cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00069       x1_shock = 0.1666666666 + x2/sqrt((double)3.0);
00070 /* upstream conditions */
00071       pGrid->U[ks][j][i].d = 1.4;
00072       pGrid->U[ks][j][i].E = 2.5;
00073       pGrid->U[ks][j][i].M1 = 0.0;
00074       pGrid->U[ks][j][i].M2 = 0.0;
00075 /* downstream conditions */
00076       if (x1 < x1_shock) {
00077         pGrid->U[ks][j][i].d = d0;
00078         pGrid->U[ks][j][i].E = e0 + 0.5*d0*(u0*u0+v0*v0);
00079         pGrid->U[ks][j][i].M1 = d0*u0;
00080         pGrid->U[ks][j][i].M2 = d0*v0;
00081       }
00082     }
00083   }
00084 
00085 /* Set boundary value function pointers */
00086 
00087   if (pDomain->Disp[0] == 0) bvals_mhd_fun(pDomain, left_x1,  dmrbv_iib);
00088   if (pDomain->Disp[1] == 0) bvals_mhd_fun(pDomain, left_x2,  dmrbv_ijb);
00089   if (pDomain->MaxX[1] == pDomain->RootMaxX[1])
00090     bvals_mhd_fun(pDomain, right_x2, dmrbv_ojb);
00091 }
00092 
00093 /*==============================================================================
00094  * PROBLEM USER FUNCTIONS:
00095  * problem_write_restart() - writes problem-specific user data to restart files
00096  * problem_read_restart()  - reads problem-specific user data from restart files
00097  * get_usr_expr()          - sets pointer to expression for special output data
00098  * get_usr_out_fun()       - returns a user defined output function pointer
00099  * get_usr_par_prop()      - returns a user defined particle selection function
00100  * Userwork_in_loop        - problem specific work IN     main loop
00101  * Userwork_after_loop     - problem specific work AFTER  main loop
00102  *----------------------------------------------------------------------------*/
00103 
00104 void problem_write_restart(MeshS *pM, FILE *fp)
00105 {
00106   return;
00107 }
00108 
00109 void problem_read_restart(MeshS *pM, FILE *fp)
00110 {
00111   return;
00112 }
00113 
00114 ConsFun_t get_usr_expr(const char *expr)
00115 {
00116   return NULL;
00117 }
00118 
00119 VOutFun_t get_usr_out_fun(const char *name){
00120   return NULL;
00121 }
00122 
00123 void Userwork_in_loop(MeshS *pM)
00124 {
00125 }
00126 
00127 void Userwork_after_loop(MeshS *pM)
00128 {
00129 }
00130 
00131 /*=========================== PRIVATE FUNCTIONS ==============================*/
00132 
00133 /*----------------------------------------------------------------------------*/
00134 /*! \fn void dmrbv_iib(GridS *pGrid)
00135  *  \brief Sets boundary condition on left X boundary (iib) for dmr test
00136  *
00137  * Note quantities at this boundary are held fixed at the downstream state
00138  */
00139 
00140 void dmrbv_iib(GridS *pGrid)
00141 {
00142 int i=0,j=0;
00143 int is,ie,js,je,ks,jl,ju;
00144 Real d0,e0,u0,v0;
00145 
00146   d0 = 8.0;
00147   e0 = 291.25;
00148   u0 =  8.25*sqrt(3.0)/2.0;
00149   v0 = -8.25*0.5;
00150 
00151   is = pGrid->is; ie = pGrid->ie;
00152   js = pGrid->js; je = pGrid->je;
00153   ks = pGrid->ks;
00154   ju = pGrid->je + nghost;
00155   jl = pGrid->js - nghost;
00156 
00157   for (j=jl; j<=ju;  j++) {
00158     for (i=1;  i<=nghost;  i++) {
00159       pGrid->U[ks][j][is-i].d  = d0;
00160       pGrid->U[ks][j][is-i].M1 = d0*u0;
00161       pGrid->U[ks][j][is-i].M2 = d0*v0;
00162       pGrid->U[ks][j][is-i].E  = e0 + 0.5*d0*(u0*u0+v0*v0);
00163     }
00164   }
00165 }
00166 
00167 /*----------------------------------------------------------------------------*/
00168 /*! \fn void dmrbv_ijb(GridS *pGrid)
00169  *  \brief  Sets boundary condition on lower Y boundary (ijb) for dmr test.
00170  *
00171  * Note quantaties at this boundary are held fixed at the downstream state for
00172  * x1 < 0.16666666, and are reflected for x1 > 0.16666666
00173  */
00174 
00175 void dmrbv_ijb(GridS *pGrid)
00176 {
00177 int i=0,j=0;
00178 int is,ie,js,je,ks,il,iu;
00179 Real d0,e0,u0,v0,x1,x2,x3;
00180 
00181   d0 = 8.0;
00182   e0 = 291.25;
00183   u0 =  8.25*sqrt(3.0)/2.0;
00184   v0 = -8.25*0.5;
00185 
00186   is = pGrid->is; ie = pGrid->ie;
00187   js = pGrid->js; je = pGrid->je;
00188   ks = pGrid->ks;
00189   iu = pGrid->ie + nghost;
00190   il = pGrid->is - nghost;
00191 
00192   for (j=1;  j<=nghost;  j++) {
00193     for (i=il; i<=iu; i++) {
00194       cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00195       if (x1 < 0.1666666666) {
00196 /* fixed at downstream state */
00197         pGrid->U[ks][js-j][i].d  = d0;
00198         pGrid->U[ks][js-j][i].M1 = d0*u0;
00199         pGrid->U[ks][js-j][i].M2 = d0*v0;
00200         pGrid->U[ks][js-j][i].E  = e0 + 0.5*d0*(u0*u0+v0*v0);
00201       } else {
00202 /* reflected */
00203         pGrid->U[ks][js-j][i].d  = pGrid->U[ks][js+(j-1)][i].d;
00204         pGrid->U[ks][js-j][i].M1 = pGrid->U[ks][js+(j-1)][i].M1;
00205         pGrid->U[ks][js-j][i].M2 = -pGrid->U[ks][js+(j-1)][i].M2;
00206         pGrid->U[ks][js-j][i].E  = pGrid->U[ks][js+(j-1)][i].E;
00207       }
00208     }
00209   }
00210 }
00211 
00212 /*----------------------------------------------------------------------------*/
00213 /*! \fn void dmrbv_ojb(GridS *pGrid)
00214  *  \brief Sets TIME-DEPENDENT boundary condition on upper Y boundary (ojb)
00215  * for dmr test.  
00216  *
00217  * Quantaties at this boundary are held fixed at the downstream
00218  * state for x1 < 0.16666666+v1_shock*time, and at the upstream state for
00219  * x1 > 0.16666666+v1_shock*time
00220  */
00221 
00222 void dmrbv_ojb(GridS *pGrid)
00223 {
00224 int i=0,j=0;
00225 int is,ie,js,je,ks,il,iu;
00226 Real d0,e0,u0,v0,x1_shock,x1,x2,x3;
00227 
00228   d0 = 8.0;
00229   e0 = 291.25;
00230   u0 =  8.25*sqrt(3.0)/2.0;
00231   v0 = -8.25*0.5;
00232   x1_shock = 0.1666666666 + (1.0 + 20.0*pGrid->time)/sqrt(3.0);
00233 
00234   is = pGrid->is; ie = pGrid->ie;
00235   js = pGrid->js; je = pGrid->je;
00236   ks = pGrid->ks;
00237   iu = pGrid->ie + nghost;
00238   il = pGrid->is - nghost;
00239 
00240   for (j=1;  j<=nghost;  j++) {
00241     for (i=il; i<=iu; i++) {
00242       cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00243       if (x1 < x1_shock) {
00244 /* fixed at downstream state */
00245         pGrid->U[ks][je+j][i].d  = d0;
00246         pGrid->U[ks][je+j][i].M1 = d0*u0;
00247         pGrid->U[ks][je+j][i].M2 = d0*v0;
00248         pGrid->U[ks][je+j][i].E  = e0 + 0.5*d0*(u0*u0+v0*v0);
00249       } else {
00250 /* fixed at upstream state */
00251         pGrid->U[ks][je+j][i].d  = 1.4;
00252         pGrid->U[ks][je+j][i].M1 = 0.0;
00253         pGrid->U[ks][je+j][i].M2 = 0.0;
00254         pGrid->U[ks][je+j][i].E  = 2.5;
00255       }
00256     }
00257   }
00258 }