Athena: prob/kh.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file kh.c
00004  *  \brief Problem generator for KH instability. 
00005  *
00006  * PURPOSE: Problem generator for KH instability.  Sets up two versions:
00007  * - iprob=1: slip surface with random perturbations
00008  * - iprob=2: tanh profile at interface, with single-mode perturbation        */
00009 /*============================================================================*/
00010 
00011 #include <float.h>
00012 #include <math.h>
00013 #include <stdio.h>
00014 #include <stdlib.h>
00015 #include "defs.h"
00016 #include "athena.h"
00017 #include "globals.h"
00018 #include "prototypes.h"
00019 
00020 /*==============================================================================
00021  * PRIVATE FUNCTION PROTOTYPES:
00022  * ran2()
00023  *============================================================================*/
00024 
00025 static double ran2(long int *idum);
00026 
00027 /*=========================== PUBLIC FUNCTIONS ===============================*/
00028 /*----------------------------------------------------------------------------*/
00029 /* problem:  */
00030 
00031 void problem(DomainS *pDomain)
00032 {
00033   GridS *pGrid = pDomain->Grid;
00034   int i=0,j=0,k=0;
00035   int is,ie,js,je,ks,ke,iprob;
00036   Real amp,drat,vflow,b0,a,sigma,x1,x2,x3;
00037   long int iseed = -1;
00038 
00039   is = pGrid->is; ie = pGrid->ie;
00040   js = pGrid->js; je = pGrid->je;
00041   ks = pGrid->ks; ke = pGrid->ke;
00042 
00043 /* Read problem parameters */
00044 
00045   iprob = par_geti("problem","iprob");
00046   vflow = par_getd("problem","vflow");
00047   drat = par_getd("problem","drat");
00048   amp = par_getd("problem","amp");
00049 #ifdef MHD
00050   b0  = par_getd("problem","b0");
00051 #endif
00052 
00053 /* iprob=1.  Two uniform streams moving at +/- vflow, random perturbations */
00054 
00055   if (iprob == 1) {
00056     for (k=ks; k<=ke; k++) {
00057       for (j=js; j<=je; j++) {
00058         for (i=is; i<=ie; i++) {
00059           cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00060           pGrid->U[k][j][i].d = 1.0;
00061           pGrid->U[k][j][i].M1 = vflow + amp*(ran2(&iseed) - 0.5);
00062           pGrid->U[k][j][i].M2 = amp*(ran2(&iseed) - 0.5);
00063           pGrid->U[k][j][i].M3 = 0.0;
00064           if (fabs(x2) < 0.25) {
00065             pGrid->U[k][j][i].d = drat;
00066             pGrid->U[k][j][i].M1 = -drat*(vflow + amp*(ran2(&iseed) - 0.5));
00067             pGrid->U[k][j][i].M2 = drat*amp*(ran2(&iseed) - 0.5);
00068           }
00069 /* Pressure scaled to give a sound speed of 1 with gamma=1.4 */
00070 #ifndef BAROTROPIC
00071           pGrid->U[k][j][i].E = 2.5/Gamma_1
00072              + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00073              + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00074 #endif /* BAROTROPIC */
00075 #ifdef MHD
00076           pGrid->B1i[k][j][i] = b0;
00077           pGrid->U[k][j][i].B1c = b0;
00078 #ifndef BAROTROPIC
00079           pGrid->U[k][j][i].E += 0.5*b0*b0;
00080 #endif /* BAROTROPIC */
00081 #endif /* MHD */
00082         }
00083 #ifdef MHD
00084       pGrid->B1i[k][j][ie+1] = b0;
00085 #endif
00086       }
00087     }
00088   }
00089 
00090 /* iprob=2.  Test suggested by E. Zweibel, based on Ryu & Jones.
00091  * Two uniform density flows with single mode perturbation
00092  */
00093 
00094   if (iprob == 2) {
00095     a = 0.05;
00096     sigma = 0.2;
00097     for (k=ks; k<=ke; k++) {
00098       for (j=js; j<=je; j++) {
00099         for (i=is; i<=ie; i++) {
00100           cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00101           pGrid->U[k][j][i].d = 1.0;
00102           pGrid->U[k][j][i].M1 = vflow*tanh(x2/a);
00103           pGrid->U[k][j][i].M2 = amp*sin(2.0*PI*x1)*exp(-(x2*x2)/(sigma*sigma));
00104           pGrid->U[k][j][i].M3 = 0.0;
00105 #ifndef BAROTROPIC
00106           pGrid->U[k][j][i].E = 1.0/Gamma_1
00107              + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00108              + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00109 #endif /* BAROTROPIC */
00110 #ifdef MHD
00111           pGrid->B1i[k][j][i] = b0;
00112           pGrid->U[k][j][i].B1c = b0;
00113 #ifndef BAROTROPIC
00114           pGrid->U[k][j][i].E += 0.5*b0*b0;
00115 #endif /* BAROTROPIC */
00116 #endif /* MHD */
00117 /* Use passive scalar to keep track of the fluids, since densities are same */
00118 #if (NSCALARS > 0)
00119           pGrid->U[k][j][i].s[0] = 0.0;
00120           if (x2 > 0) pGrid->U[k][j][i].s[0] = 1.0;
00121 #endif
00122         }
00123 #ifdef MHD
00124       pGrid->B1i[k][j][ie+1] = b0;
00125 #endif
00126       }
00127     }
00128   }
00129 
00130 /* With viscosity and/or resistivity, read eta_R and nu_V */
00131 
00132 #ifdef OHMIC
00133   eta_Ohm = par_getd("problem","eta");
00134 #endif
00135 #ifdef NAVIER_STOKES
00136   nu_V = par_getd("problem","nu");
00137 #endif
00138 #ifdef BRAGINSKII
00139   nu_V = par_getd("problem","nu");
00140 #endif
00141 
00142 }
00143 
00144 /*==============================================================================
00145  * PROBLEM USER FUNCTIONS:
00146  * problem_write_restart() - writes problem-specific user data to restart files
00147  * problem_read_restart()  - reads problem-specific user data from restart files
00148  * get_usr_expr()          - sets pointer to expression for special output data
00149  * get_usr_out_fun()       - returns a user defined output function pointer
00150  * get_usr_par_prop()      - returns a user defined particle selection function
00151  * Userwork_in_loop        - problem specific work IN     main loop
00152  * Userwork_after_loop     - problem specific work AFTER  main loop
00153  *----------------------------------------------------------------------------*/
00154 
00155 void problem_write_restart(MeshS *pM, FILE *fp)
00156 {
00157   return;
00158 }
00159 
00160 void problem_read_restart(MeshS *pM, FILE *fp)
00161 {
00162 #ifdef OHMIC
00163   eta_Ohm = par_getd("problem","eta");
00164 #endif
00165 #ifdef NAVIER_STOKES
00166   nu_V = par_getd("problem","nu");
00167 #endif
00168 #ifdef BRAGINSKII
00169   nu_V = par_getd("problem","nu");
00170 #endif
00171   return;
00172 }
00173 
00174 #if (NSCALARS > 0)
00175 /*! \fn static Real color(const GridS *pG, const int i, const int j,const int k)
00176  *  \brief Returns first passively advected scalar s[0] */
00177 static Real color(const GridS *pG, const int i, const int j, const int k)
00178 {
00179   return pG->U[k][j][i].s[0]/pG->U[k][j][i].d;
00180 }
00181 #endif
00182 
00183 ConsFun_t get_usr_expr(const char *expr)
00184 {
00185 #if (NSCALARS > 0)
00186   if(strcmp(expr,"color")==0) return color;
00187 #endif
00188   return NULL;
00189 }
00190 
00191 VOutFun_t get_usr_out_fun(const char *name){
00192   return NULL;
00193 }
00194 
00195 void Userwork_in_loop(MeshS *pM)
00196 {
00197   return;
00198 }
00199 
00200 void Userwork_after_loop(MeshS *pM)
00201 {
00202   return;
00203 }
00204 
00205 /*=========================== PRIVATE FUNCTIONS ==============================*/
00206 /*------------------------------------------------------------------------------
00207  */
00208 
00209 #define IM1 2147483563
00210 #define IM2 2147483399
00211 #define AM (1.0/IM1)
00212 #define IMM1 (IM1-1)
00213 #define IA1 40014
00214 #define IA2 40692
00215 #define IQ1 53668
00216 #define IQ2 52774
00217 #define IR1 12211
00218 #define IR2 3791
00219 #define NTAB 32
00220 #define NDIV (1+IMM1/NTAB)
00221 #define RNMX (1.0-DBL_EPSILON)
00222 
00223 /*! \fn double ran2(long int *idum)
00224  *  \brief  Extracted from the Numerical Recipes in C (version 2) code. Modified
00225  *   to use doubles instead of floats. -- T. A. Gardiner -- Aug. 12, 2003
00226  *
00227  * Long period (> 2 x 10^{18}) random number generator of L'Ecuyer
00228  * with Bays-Durham shuffle and added safeguards.  Returns a uniform
00229  * random deviate between 0.0 and 1.0 (exclusive of the endpoint
00230  * values).  Call with idum = a negative integer to initialize;
00231  * thereafter, do not alter idum between successive deviates in a
00232  * sequence.  RNMX should appriximate the largest floating point value
00233  * that is less than 1. 
00234 
00235  */
00236 double ran2(long int *idum){
00237   int j;
00238   long int k;
00239   static long int idum2=123456789;
00240   static long int iy=0;
00241   static long int iv[NTAB];
00242   double temp;
00243 
00244   if (*idum <= 0) { /* Initialize */
00245     if (-(*idum) < 1) *idum=1; /* Be sure to prevent idum = 0 */
00246     else *idum = -(*idum);
00247     idum2=(*idum);
00248     for (j=NTAB+7;j>=0;j--) { /* Load the shuffle table (after 8 warm-ups) */
00249       k=(*idum)/IQ1;
00250       *idum=IA1*(*idum-k*IQ1)-k*IR1;
00251       if (*idum < 0) *idum += IM1;
00252       if (j < NTAB) iv[j] = *idum;
00253     }
00254     iy=iv[0];
00255   }
00256   k=(*idum)/IQ1;                 /* Start here when not initializing */
00257   *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without */
00258   if (*idum < 0) *idum += IM1;   /* overflows by Schrage's method */
00259   k=idum2/IQ2;
00260   idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise */
00261   if (idum2 < 0) idum2 += IM2;
00262   j=(int)(iy/NDIV);              /* Will be in the range 0...NTAB-1 */
00263   iy=iv[j]-idum2;                /* Here idum is shuffled, idum and idum2 */
00264   iv[j] = *idum;                 /* are combined to generate output */
00265   if (iy < 1) iy += IMM1;
00266   if ((temp=AM*iy) > RNMX) return RNMX; /* No endpoint values */
00267   else return temp;
00268 }
00269 
00270 #undef IM1
00271 #undef IM2
00272 #undef AM
00273 #undef IMM1
00274 #undef IA1
00275 #undef IA2
00276 #undef IQ1
00277 #undef IQ2
00278 #undef IR1
00279 #undef IR2
00280 #undef NTAB
00281 #undef NDIV
00282 #undef RNMX