Athena: prob/cshock1d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file cshock1d.c
00004  *  \brief Problem generator for 1-D standing C-type shock test.
00005  *
00006  * PURPOSE: Problem generator for 1-D standing C-type shock test. Only works for
00007  *   grid-aligned shock, the shock is along the x1 direction.
00008  *
00009  * REFERENCE: Mac Low, M-M et al., "Incorporation of Ambipolar Diffusion into
00010  *   the Zeus Magnetohydrodynamics Code", 1995, ApJ, 442, 726                 */
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 #ifndef MHD
00022 #error: The cshock1d test only works for mhd.
00023 #endif
00024 
00025 #ifndef ISOTHERMAL
00026 #error: The cshock1d test only works for isothermal equation of state.
00027 #endif
00028 
00029 #ifndef RESISTIVITY
00030 #error: The cshock1d test only works with RESISTIVITY.
00031 #endif
00032 
00033 /*--------------- Private Functions --------------------*/
00034 /* RK4 integrator for the semi-analytical solution */
00035 
00036 Real RK4(Real D, Real A, Real M, Real theta, Real h);
00037 
00038 Real Dprime(Real D, Real A, Real M, Real theta);
00039 
00040 /* Solution in the root domain */
00041 static ConsS *RootSoln=NULL;
00042 
00043 /*----------------------------------------------------------------------------*/
00044 /* problem:   */
00045 
00046 void problem(DomainS *pDomain)
00047 {
00048   GridS *pGrid=(pDomain->Grid);
00049   ConsS *Soln;
00050   int i, is = pGrid->is, ie = pGrid->ie;
00051   int j, js = pGrid->js, je = pGrid->je;
00052   int k, ks = pGrid->ks, ke = pGrid->ke;
00053   int i0, nx1, N_L;
00054   Real x1min,x1max,Lx,x1,x2,x3,xs,xe,h,x01,x02;
00055   Real Mach,Alfv,theta,d0,v0,vA,B0,Bx0,By0;
00056   Real Ls,Ns,D01,D02,myD;
00057 
00058   nx1 = (ie-is)+1 + 2*nghost;
00059 
00060 /* allocate memory for solution on this level.  If this is root level
00061  * also allocate memory for RootSoln */
00062   if ((Soln = (ConsS*)calloc_1d_array(nx1,sizeof(ConsS)))==NULL)
00063     ath_error("[problem]: Error allocating memory for Soln\n");
00064   if (pDomain->Level == 0){
00065     if ((RootSoln = (ConsS*)calloc_1d_array(nx1,sizeof(ConsS)))
00066       == NULL) ath_error("[problem]: Error alloc memory for RootSoln\n");
00067   }
00068 
00069 /* Model parameters */
00070   Mach = par_getd("problem","Mach"); /* Mach number */
00071   Alfv = par_getd("problem","Alfv"); /* Alfven Mach number */
00072 
00073   theta = par_getd("problem","theta"); /* magnetic field obliquety (deg) */
00074   theta = theta * PI / 180.0;
00075 
00076 /* upstream quantities (in the shock frame) */
00077   d0 = 1.0;
00078   v0 = Mach * Iso_csound;
00079 
00080   vA = (Mach/Alfv) * Iso_csound;
00081   B0 = sqrt(SQR(vA)*d0);
00082 
00083   Bx0 = B0 * cos(theta);
00084   By0 = B0 * sin(theta);
00085 
00086 /* root domain info */
00087   x1min = pDomain->RootMinX[0];
00088   x1max = pDomain->RootMaxX[0];
00089   Lx = x1max - x1min;
00090 
00091 /* Ambipolar diffusion coefficient 
00092  * ambipolar diffusion length scale is FIXED to be 1 in code unit! */
00093 
00094   Q_AD = 1.0/vA;            /* N.B. AD diffusion coefficient is set this way! */
00095   d_ind= 0.0;               /* constant ion density */
00096 
00097 /* info for semi-analytic calculation */
00098   Ls  = par_getd_def("problem","Ls",20.0);  /* Estimated thickness of the
00099                                              * C-shock (in unit of L_A) */
00100   Ns  = par_getd_def("problem","Ns",5e3);   /* number of cells in Ls in the
00101                                              * semi-analytic calculation */
00102   if (Ls > Lx)
00103     ath_error("Domain size in x1 is shorter than the C-Shosck thickness!\n");
00104 
00105   xs = x1min + 0.5*(Lx - Ls);
00106   xe = xs + Ls;
00107   h  = (xe - xs)/Ns;
00108 
00109 /* Find the 1D solution on this grid.
00110  * Semi-analytic solution is only applied in the middle of the root domain.
00111  * Have to initialize different segments separately. */
00112 
00113   /* Upstream region */
00114   i = is;
00115   cc_pos(pGrid,i,js,ks,&x1,&x2,&x3); 
00116   while (x1 < xs){
00117 
00118     Soln[i].d  = d0;
00119     Soln[i].M1 = d0 * v0;
00120     Soln[i].M2 = 0.0;
00121     Soln[i].M3 = 0.0;
00122     Soln[i].B1c= Bx0;
00123     Soln[i].B2c= By0;
00124     Soln[i].B3c= 0.0;
00125 
00126     i += 1;
00127     cc_pos(pGrid,i,js,ks,&x1,&x2,&x3);
00128   }
00129 
00130   /* C-shock region */
00131   x01 = xs;             x02 = x01+h;
00132   D01 = d0+1.0e-6;
00133   while (x02 <= xe)
00134   {
00135     D02 = RK4(D01, Alfv, Mach, theta, h);
00136     if ((x1 >= x01) && (x1 < x02))
00137     {
00138       myD = (D01*(x02-x1)+D02*(x1-x01))/h;
00139       Soln[i].d  = myD;
00140       Soln[i].B1c= Bx0;
00141       Soln[i].B2c= sqrt(SQR(By0)+2.0*SQR(Alfv*B0)*
00142                        (myD-1)*(1.0/myD-SQR(1.0/Mach)));
00143       Soln[i].B3c= 0.0;
00144       Soln[i].M1 = d0 * v0;
00145       Soln[i].M2 = myD* SQR(vA)/v0*cos(theta)*(Soln[i].B2c/B0-sin(theta));
00146       Soln[i].M3 = 0.0; 
00147 
00148       i += 1;
00149       cc_pos(pGrid,i,js,ks,&x1,&x2,&x3);
00150     }
00151     x01 = x02;
00152     x02 += h;
00153     D01 = D02;
00154   }
00155   i0 = i-1;
00156 
00157   /* Downstream region */
00158   while (i <= ie+1)
00159   {
00160     Soln[i] = Soln[i0];
00161     i += 1;
00162   }
00163 
00164   /* backup the root domain solution */
00165   if (pDomain->Level == 0){
00166     for (i=is; i<=ie+1; i++){
00167       RootSoln[i].d  = Soln[i].d;
00168       RootSoln[i].M1 = Soln[i].M1;
00169       RootSoln[i].M2 = Soln[i].M2;
00170       RootSoln[i].M3 = Soln[i].M3;
00171       RootSoln[i].B1c= Soln[i].B1c;
00172       RootSoln[i].B2c= Soln[i].B2c;
00173       RootSoln[i].B3c= Soln[i].B3c;
00174     }
00175   }
00176 
00177   /* set initial conditions */
00178   ie = ie+1;
00179   if (pGrid->Nx[1] > 1) je = je+1;
00180   if (pGrid->Nx[2] > 1) ke = ke+1;
00181 
00182   for (k=ks; k<=ke; k++) {
00183   for (j=js; j<=je; j++) {
00184   for (i=is; i<=ie; i++) {
00185 
00186     pGrid->U[k][j][i].d  = Soln[i].d;
00187     pGrid->U[k][j][i].M1 = d0 * v0;
00188     pGrid->U[k][j][i].M2 = Soln[i].M2;
00189     pGrid->U[k][j][i].M3 = 0.0;
00190 
00191     pGrid->U[k][j][i].B1c = pGrid->B1i[k][j][i] = Bx0;
00192     pGrid->U[k][j][i].B2c = pGrid->B2i[k][j][i] = Soln[i].B2c;
00193     pGrid->U[k][j][i].B3c = pGrid->B3i[k][j][i] = 0.0;
00194   }}}
00195 
00196   return;
00197 }
00198 
00199 /*==============================================================================
00200  * PROBLEM USER FUNCTIONS:
00201  * problem_write_restart() - writes problem-specific user data to restart files
00202  * problem_read_restart()  - reads problem-specific user data from restart files
00203  * get_usr_expr()          - sets pointer to expression for special output data
00204  * get_usr_out_fun()       - returns a user defined output function pointer
00205  * get_usr_par_prop()      - returns a user defined particle selection function
00206  * Userwork_in_loop        - problem specific work IN     main loop
00207  * Userwork_after_loop     - problem specific work AFTER  main loop
00208  *----------------------------------------------------------------------------*/
00209     
00210 void problem_write_restart(MeshS *pM, FILE *fp)
00211 {
00212   return;
00213 }
00214   
00215 void problem_read_restart(MeshS *pM, FILE *fp)
00216 {
00217   return;
00218 } 
00219     
00220 ConsFun_t get_usr_expr(const char *expr)
00221 {     
00222   return NULL;
00223 }     
00224 
00225 VOutFun_t get_usr_out_fun(const char *name){
00226   return NULL;
00227 }     
00228 
00229 #ifdef RESISTIVITY
00230 void get_eta_user(GridS *pG, int i, int j, int k,
00231                              Real *eta_O, Real *eta_H, Real *eta_A)
00232 { 
00233   *eta_O = 0.0;
00234   *eta_H = 0.0;
00235   *eta_A = 0.0;
00236 
00237   return;
00238 }
00239 #endif
00240 
00241 void Userwork_in_loop(MeshS *pM)
00242 {
00243 }
00244 
00245 /*---------------------------------------------------------------------------
00246  * Userwork_after_loop: computes L1-error in CPAW,
00247  * ASSUMING WAVE HAS PROPAGATED AN INTEGER NUMBER OF PERIODS
00248  * Must set parameters in input file appropriately so that this is true
00249  */
00250 
00251 void Userwork_after_loop(MeshS *pM)
00252 {
00253 }
00254 
00255 /*----------------------------------------------------------------------------*/
00256 /* Semi-analytical C-shock solution */
00257 
00258 /*! \fn Real RK4(Real D, Real A, Real M, Real theta, Real h)
00259  *  \brief 4th order Runge-Kutta integrator */
00260 Real RK4(Real D, Real A, Real M, Real theta, Real h)
00261 {
00262   Real k1, k2, k3, k4;
00263 
00264   k1 = Dprime(D,A,M,theta);
00265   k2 = Dprime(D+0.5*h*k1,A,M,theta);
00266   k3 = Dprime(D+0.5*h*k2,A,M,theta);
00267   k4 = Dprime(D+h*k3,A,M,theta);
00268 
00269   return D + h/6.0*(k1+2.0*k2+2.0*k3+k4);
00270 
00271 }
00272 
00273 /*! \fn Real Dprime(Real D, Real A, Real M, Real theta)
00274  *  \brief dD/dt */
00275 Real Dprime(Real D, Real A, Real M, Real theta)
00276 {
00277   Real sintheta, costheta, sintheta2, costheta2;
00278   Real M21, b, b2;
00279 
00280   sintheta = sin(theta);
00281   costheta = cos(theta);
00282   sintheta2 = SQR(sintheta);
00283   costheta2 = SQR(costheta);
00284   M21 = 1.0/SQR(M);
00285 
00286   b2 = sintheta2 + 2*SQR(A)*(D-1.0)*(1.0/D-M21);
00287   b  = sqrt(b2);
00288 
00289   return b/A*(b-D*((b-sintheta)/SQR(A)*costheta2+sintheta))/(b2+costheta2)/(1/SQR(D)-M21);
00290 }
00291