Athena: prob/cpaw2d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file cpaw2d.c
00004  *  \brief Problem generator for 2-D circularly polarized Alfven wave (CPAW) 
00005  *  test.
00006  *
00007  * PURPOSE: Problem generator for 2-D circularly polarized Alfven wave (CPAW)
00008  *   test.  Works for any arbitrary wavevector in the x1-x2 plane.  The wave is
00009  *   defined with reference to a coordinate system (x,y,z) with transformation
00010  *   rules to the code coordinate system (x1,x2,x3)
00011  *   -  x =  x1*cos(alpha) + x2*sin(alpha)
00012  *   -  y = -x1*sin(alpha) + x2*cos(alpha)
00013  *   -  z = x3
00014  *   The magnetic field is given by:
00015  *   - B_x = b_par
00016  *   - B_y = b_perp*sin(k*x)
00017  *   - B_z = b_perp*cos(k*x)   where k = 2.0*PI/lambda
00018  *
00019  *   Can be used for either standing (problem/v_par=1.0) or travelling
00020  *   (problem/v_par=0.0) waves.
00021  *
00022  * USERWORK_AFTER_LOOP function computes L1 error norm in solution by comparing
00023  *   to initial conditions.  Problem must be evolved for an integer number of
00024  *   wave periods for this to work.
00025  *
00026  * REFERENCE: G. Toth,  "The div(B)=0 constraint in shock capturing MHD codes",
00027  *   JCP, 161, 605 (2000)                                                     */
00028 /*============================================================================*/
00029 
00030 #include <math.h>
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include "defs.h"
00034 #include "athena.h"
00035 #include "globals.h"
00036 #include "prototypes.h"
00037 
00038 #ifndef MHD
00039 #error : The cpaw2d test only works for MHD.
00040 #endif
00041 
00042 
00043 /* Initial solution, shared with Userwork_after_loop to compute L1 error.
00044  * Vector potential A3 defined in prvate function below.  B_par, etc. must
00045  * be defined as globals to be used by A3() */
00046  
00047 static ConsS **RootSoln=NULL;
00048 static Real A3(const Real x1, const Real x2);
00049 Real fac, sin_a, cos_a, b_par, b_perp;
00050 Real k_par;
00051 
00052 /*----------------------------------------------------------------------------*/
00053 /* problem:   */
00054 
00055 void problem(DomainS *pDomain)
00056 {
00057   GridS *pGrid = pDomain->Grid;
00058   ConsS **Soln;
00059   int i, is = pGrid->is, ie = pGrid->ie;
00060   int j, js = pGrid->js, je = pGrid->je;
00061   int k, ks = pGrid->ks, ke = pGrid->ke;
00062   int nx1, nx2;
00063   int dir;
00064   Real angle;    /* Angle the wave direction makes with the x1-direction */
00065   Real x1size,x2size,x1,x2,x3,cs,sn;
00066   Real v_par, v_perp, den, pres;
00067   Real lambda; /* Wavelength */
00068 #ifdef RESISTIVITY
00069   Real v_A, kva, omega_h, omega_l, omega_r;
00070 #endif
00071  
00072   nx1 = (ie - is + 1) + 2*nghost;
00073   nx2 = (je - js + 1) + 2*nghost;
00074 
00075   if (pGrid->Nx[1] == 1) {
00076     ath_error("[problem] Grid must be 2D");
00077   }
00078 
00079   if ((Soln = (ConsS**)calloc_2d_array(nx2,nx1,sizeof(ConsS))) == NULL)
00080     ath_error("[problem]: Error allocating memory for Soln\n");
00081 
00082   if (pDomain->Level == 0){
00083     if ((RootSoln =(ConsS**)calloc_2d_array(nx2,nx1,sizeof(ConsS)))==NULL)
00084       ath_error("[problem]: Error allocating memory for RootSoln\n");
00085   }
00086 
00087 /* An angle =  0.0 is a wave aligned with the x1-direction. */
00088 /* An angle = 90.0 is a wave aligned with the x2-direction. */
00089 
00090   angle = par_getd("problem","angle");
00091   x1size = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00092   x2size = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00093 
00094 /* Compute the sin and cos of the angle and the wavelength. */
00095 /* Put one wavelength in the grid */
00096 
00097   if (angle == 0.0) {
00098     sin_a = 0.0;
00099     cos_a = 1.0;
00100     lambda = x1size;
00101 ath_pout(0,"Yes!\n");
00102   }
00103   else if (angle == 90.0) {
00104     sin_a = 1.0;
00105     cos_a = 0.0;
00106     lambda = x2size;
00107   }
00108   else {
00109 
00110 /* We put 1 wavelength in each direction.  Hence the wavelength
00111  *      lambda = (pDomain->RootMaxX[0] - pDomain->RootMinX[0])*cos_a
00112  *  AND lambda = (pDomain->RootMaxX[1] - pDomain->RootMinX[1])*sin_a;
00113  *  are both satisfied. */
00114 
00115     if(x1size == x2size){
00116       cos_a = sin_a = sqrt(0.5);
00117     }
00118     else{
00119       angle = atan((double)(x1size/x2size));
00120       sin_a = sin(angle);
00121       cos_a = cos(angle);
00122     }
00123 /* Use the larger angle to determine the wavelength */
00124     if (cos_a >= sin_a) {
00125       lambda = x1size*cos_a;
00126     } else {
00127       lambda = x2size*sin_a;
00128     }
00129   }
00130 
00131 /* Initialize k_parallel */
00132 
00133   k_par = 2.0*PI/lambda;
00134   b_par = par_getd("problem","b_par");
00135   den = 1.0;
00136 
00137   ath_pout(0,"va_parallel = %g\nlambda = %g\n",b_par/sqrt(den),lambda);
00138 
00139   b_perp = par_getd("problem","b_perp");
00140   v_perp = b_perp/sqrt((double)den);
00141 
00142   dir    = par_geti_def("problem","dir",1); /* right(1)/left(2) polarization */
00143   if (dir == 1) /* right polarization */
00144     fac = 1.0;
00145   else          /* left polarization */
00146     fac = -1.0;
00147 
00148 #ifdef RESISTIVITY
00149   Q_Hall = par_getd("problem","Q_H");
00150   d_ind  = 0.0;
00151   v_A    = b_par/sqrt((double)den);
00152   if (Q_Hall > 0.0)
00153   {
00154     kva     = k_par*v_A;
00155     omega_h = 1.0/Q_Hall;
00156 
00157     omega_r = 0.5*SQR(kva)/omega_h*(sqrt(1.0+SQR(2.0*omega_h/kva)) + 1.0);
00158     omega_l = 0.5*SQR(kva)/omega_h*(sqrt(1.0+SQR(2.0*omega_h/kva)) - 1.0);
00159 
00160     if (dir == 1) /* right polarization (whistler wave) */
00161       v_perp = v_perp * kva / omega_r;
00162     else          /* left polarization */
00163       v_perp = v_perp * kva / omega_l;
00164   }
00165 #endif
00166 
00167 /* The gas pressure and parallel velocity are free parameters. */
00168 
00169   pres = par_getd("problem","pres");
00170   v_par = par_getd("problem","v_par");
00171 
00172 /* Use the vector potential to initialize the interface magnetic fields
00173  * The iterface fields are located at the left grid cell face normal */
00174 
00175   for (k=ks; k<=ke; k++) {
00176     for (j=js; j<=je+1; j++) {
00177       for (i=is; i<=ie+1; i++) {
00178         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00179         cs = cos(k_par*(x1*cos_a + x2*sin_a));
00180 
00181         x1 -= 0.5*pGrid->dx1;
00182         x2 -= 0.5*pGrid->dx2;
00183 
00184         pGrid->B1i[k][j][i] = -(A3(x1,(x2+pGrid->dx2)) - A3(x1,x2))/pGrid->dx2;
00185         pGrid->B2i[k][j][i] =  (A3((x1+pGrid->dx1),x2) - A3(x1,x2))/pGrid->dx1;
00186         pGrid->B3i[k][j][i] = b_perp*cs;
00187       }
00188     }
00189   }
00190   if (pGrid->Nx[2] > 1) {
00191     for (j=js; j<=je+1; j++) {
00192       for (i=is; i<=ie+1; i++) {
00193         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00194         cs = cos(k_par*(x1*cos_a + x2*sin_a));
00195         pGrid->B3i[ke+1][j][i] = b_perp*cs;
00196       }
00197     }
00198   }
00199 
00200 /* Now initialize the cell centered quantities */
00201 
00202   for (k=ks; k<=ke; k++) {
00203     for (j=js; j<=je; j++) {
00204       for (i=is; i<=ie; i++) {
00205         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00206 
00207         sn = sin(k_par*(x1*cos_a + x2*sin_a)) * fac;
00208         cs = cos(k_par*(x1*cos_a + x2*sin_a));
00209 
00210         Soln[j][i].d  = den;
00211         Soln[j][i].M1 = den*(v_par*cos_a + v_perp*sn*sin_a);
00212         Soln[j][i].M2 = den*(v_par*sin_a - v_perp*sn*cos_a);
00213         Soln[j][i].M3 = -den*v_perp*cs;
00214         pGrid->U[k][j][i].d  = Soln[j][i].d;
00215         pGrid->U[k][j][i].M1 = Soln[j][i].M1;
00216         pGrid->U[k][j][i].M2 = Soln[j][i].M2;
00217         pGrid->U[k][j][i].M3 = Soln[j][i].M3;
00218 
00219         Soln[j][i].B1c = 0.5*(pGrid->B1i[k][j][i] + pGrid->B1i[k][j][i+1]);
00220         Soln[j][i].B2c = 0.5*(pGrid->B2i[k][j][i] + pGrid->B2i[k][j+1][i]);
00221         Soln[j][i].B3c = b_perp*cs;
00222         pGrid->U[k][j][i].B1c = Soln[j][i].B1c;
00223         pGrid->U[k][j][i].B2c = Soln[j][i].B2c;
00224         pGrid->U[k][j][i].B3c = Soln[j][i].B3c;
00225 
00226 #ifndef ISOTHERMAL
00227         Soln[j][i].E = pres/Gamma_1 + 0.5*(SQR(pGrid->U[k][j][i].B1c) +
00228                  SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c) )
00229           + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2) +
00230                  SQR(pGrid->U[k][j][i].M3) )/den;
00231         pGrid->U[k][j][i].E = Soln[j][i].E;
00232 #endif
00233       }
00234     }
00235   }
00236 
00237 /* save solution on root grid */
00238 
00239   if (pDomain->Level == 0) {
00240     for (j=js; j<=je; j++) {
00241     for (i=is; i<=ie; i++) {
00242       RootSoln[j][i].d  = Soln[j][i].d ;
00243       RootSoln[j][i].M1 = Soln[j][i].M1;
00244       RootSoln[j][i].M2 = Soln[j][i].M2;
00245       RootSoln[j][i].M3 = Soln[j][i].M3;
00246 #ifndef ISOTHERMAL
00247       RootSoln[j][i].E  = Soln[j][i].E ;
00248 #endif /* ISOTHERMAL */
00249 #ifdef MHD
00250       RootSoln[j][i].B1c = Soln[j][i].B1c;
00251       RootSoln[j][i].B2c = Soln[j][i].B2c;
00252       RootSoln[j][i].B3c = Soln[j][i].B3c;
00253 #endif
00254 #if (NSCALARS > 0)
00255       for (n=0; n<NSCALARS; n++)
00256         RootSoln[j][i].s[n] = Soln[j][i].s[n];
00257 #endif
00258     }}
00259   }
00260 
00261   return;
00262 }
00263 
00264 /*==============================================================================
00265  * PROBLEM USER FUNCTIONS:
00266  * problem_write_restart() - writes problem-specific user data to restart files
00267  * problem_read_restart()  - reads problem-specific user data from restart files
00268  * get_usr_expr()          - sets pointer to expression for special output data
00269  * get_usr_out_fun()       - returns a user defined output function pointer
00270  * get_usr_par_prop()      - returns a user defined particle selection function
00271  * Userwork_in_loop        - problem specific work IN     main loop
00272  * Userwork_after_loop     - problem specific work AFTER  main loop
00273  * A3() - computes vector potential to initialize fields
00274  *----------------------------------------------------------------------------*/
00275 
00276 void problem_write_restart(MeshS *pM, FILE *fp)
00277 {
00278   return;
00279 }
00280 
00281 void problem_read_restart(MeshS *pM, FILE *fp)
00282 {
00283   return;
00284 }
00285 
00286 ConsFun_t get_usr_expr(const char *expr)
00287 {
00288   return NULL;
00289 }
00290 
00291 VOutFun_t get_usr_out_fun(const char *name){
00292   return NULL;
00293 }
00294 
00295 #ifdef RESISTIVITY
00296 void get_eta_user(GridS *pG, int i, int j, int k,
00297                              Real *eta_O, Real *eta_H, Real *eta_A)
00298 {   
00299   *eta_O = 0.0;
00300   *eta_H = 0.0;
00301   *eta_A = 0.0;
00302   
00303   return;
00304 }
00305 #endif
00306 
00307 void Userwork_in_loop(MeshS *pM)
00308 {
00309 }
00310 
00311 /*---------------------------------------------------------------------------
00312  * Userwork_after_loop: computes L1-error in CPAW,
00313  * ASSUMING WAVE HAS PROPAGATED AN INTEGER NUMBER OF PERIODS
00314  * Must set parameters in input file appropriately so that this is true
00315  */
00316 
00317 void Userwork_after_loop(MeshS *pM)
00318 {
00319   GridS *pGrid;
00320   int i,j,is,ie,js,je,ks,Nx1,Nx2;
00321   Real rms_error=0.0;
00322   ConsS error;
00323   FILE *fp;
00324   char *fname;
00325   error.d = 0.0;
00326   error.M1 = 0.0;
00327   error.M2 = 0.0;
00328   error.M3 = 0.0;
00329   error.B1c = 0.0;
00330   error.B2c = 0.0;
00331   error.B3c = 0.0;
00332 #ifndef ISOTHERMAL
00333   error.E = 0.0;
00334 #endif /* ISOTHERMAL */
00335 
00336 /* Compute error only on root Grid, which is in Domain[0][0] */
00337 
00338   pGrid=pM->Domain[0][0].Grid;
00339   if (pGrid == NULL) return;
00340 
00341   is = pGrid->is; ie = pGrid->ie;
00342   js = pGrid->js; je = pGrid->je;
00343   ks = pGrid->ks;
00344   Nx1 = (ie-is+1);
00345   Nx2 = (je-js+1);
00346 
00347 /* compute L1 error in each variable, and rms total error */
00348 
00349   for (j=js; j<=je; j++) {
00350     for (i=is; i<=ie; i++) {
00351       error.d   += fabs(pGrid->U[ks][j][i].d   - RootSoln[j][i].d  );
00352       error.M1  += fabs(pGrid->U[ks][j][i].M1  - RootSoln[j][i].M1 );
00353       error.M2  += fabs(pGrid->U[ks][j][i].M2  - RootSoln[j][i].M2 );
00354       error.M3  += fabs(pGrid->U[ks][j][i].M3  - RootSoln[j][i].M3 );
00355       error.B1c += fabs(pGrid->U[ks][j][i].B1c - RootSoln[j][i].B1c);
00356       error.B2c += fabs(pGrid->U[ks][j][i].B2c - RootSoln[j][i].B2c);
00357       error.B3c += fabs(pGrid->U[ks][j][i].B3c - RootSoln[j][i].B3c);
00358 #ifndef ISOTHERMAL
00359       error.E   += fabs(pGrid->U[ks][j][i].E   - RootSoln[j][i].E  );
00360 #endif /* ISOTHERMAL */
00361     }
00362   }
00363 
00364 /* Compute RMS error over all variables */
00365 
00366   rms_error = SQR(error.d) + SQR(error.M1) + SQR(error.M2) + SQR(error.M3);
00367   rms_error += SQR(error.B1c) + SQR(error.B2c) + SQR(error.B3c);
00368 #ifndef ISOTHERMAL
00369   rms_error += SQR(error.E);
00370 #endif /* ISOTHERMAL */
00371   rms_error = sqrt(rms_error)/(double)(Nx1*Nx2);
00372 
00373 /* Print error to file "cpaw2d-errors.dat" */
00374 
00375   fname = "cpaw2d-errors.dat";
00376 /* The file exists -- reopen the file in append mode */
00377   if((fp=fopen(fname,"r")) != NULL){
00378     if((fp = freopen(fname,"a",fp)) == NULL){
00379       ath_error("[Userwork_after_loop]: Unable to reopen file.\n");
00380       return;
00381     }
00382   }
00383 /* The file does not exist -- open the file in write mode */
00384   else{
00385     if((fp = fopen(fname,"w")) == NULL){
00386       ath_error("[Userwork_after_loop]: Unable to open file.\n");
00387       return;
00388     }
00389 /* write out some header information */
00390     fprintf(fp,"# Nx1   Nx2  RMS-Error  d  M1  M2  M3");
00391 #ifndef ISOTHERMAL
00392     fprintf(fp,"  E");
00393 #endif /* ISOTHERMAL */
00394     fprintf(fp,"  B1c  B2c  B3c");
00395     fprintf(fp,"\n#\n");
00396   }
00397 
00398   fprintf(fp,"%d  %d %e  %e  %e  %e  %e",Nx1,Nx2,rms_error,
00399           (error.d/(double)(Nx1*Nx2)),
00400           (error.M1/(double)(Nx1*Nx2)),
00401           (error.M2/(double)(Nx1*Nx2)),
00402           (error.M3/(double)(Nx1*Nx2)));
00403 #ifndef ISOTHERMAL
00404   fprintf(fp,"  %e",(error.E/(double)(Nx1*Nx2)));
00405 #endif /* ISOTHERMAL */
00406   fprintf(fp,"  %e  %e  %e",
00407           (error.B1c/(double)(Nx1*Nx2)),
00408           (error.B2c/(double)(Nx1*Nx2)),
00409           (error.B3c/(double)(Nx1*Nx2)));
00410   fprintf(fp,"\n");
00411 
00412   fclose(fp);
00413 
00414   return;
00415 }
00416 
00417 /*---------------------------------------------------------------------------*/
00418 /*! \fn static Real A3(const Real x1, const Real x2)
00419  *  \brief Define a scalar potential A3 such that:
00420  *   - B_x = - $\partial A3 / \partial y$
00421  *   - B_y =   $\partial A3 / \partial x$
00422  *   Then A3 is given in the function below.  */
00423 
00424 static Real A3(const Real x1, const Real x2)
00425 {
00426   return b_par*(x1*sin_a - x2*cos_a) 
00427      - (fac*b_perp/k_par)*cos(k_par*(x1*cos_a + x2*sin_a));
00428 }