Athena: prob/linear_wave2d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file linear_wave2d.c
00004  *  \brief Problem generator for linear wave convergence tests in 2D. 
00005  *
00006  * PURPOSE: Problem generator for linear wave convergence tests in 2D.  In 2D,
00007  *   the angle the wave propagates to the grid is automatically computed
00008  *   to be tan^{-1} (X/Y), so that periodic boundary conditions can be used,
00009  *   and the size of the box is automatically rescaled so that the fast wave
00010  *   crossing time (across a diagonal) is one-half, the Alfven wave crossing
00011  *   time is one, and the slow wave crossing time is two.
00012  *
00013  *   Note angle=0 or 90 [grid aligned waves] is not allowed in this routine.
00014  *   Use linear_wave1d for grid aligned wave on 1D/2D/3D grids.
00015  *
00016  *   Can be used for either standing (problem/vflow=1.0) or travelling
00017  *   (problem/vflow=0.0) waves.
00018  *
00019  *   Configure --with-gravity=fft to check Jeans stability of plane
00020  *   waves propagating at an angle to the grid.
00021  *
00022  *   Configure --enable-resistivity and/or --with-viscosity=ns to check
00023  *   damping of linear waves by resistivity/ambipolar diffusion and/or viscosity
00024  *
00025  * USERWORK_AFTER_LOOP function computes L1 error norm in solution by comparing
00026  *   to initial conditions.  Problem must be evolved for an integer number of
00027  *   wave periods for this to work (only works for ideal MHD).
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 /* Initial solution, shared with Userwork_after_loop to compute L1 error */
00039 static ConsS ***RootSoln=NULL;
00040 static int wave_flag;
00041 
00042 /*----------------------------------------------------------------------------*/
00043 /* problem:   */
00044 
00045 void problem(DomainS *pDomain)
00046 {
00047   GridS *pGrid=(pDomain->Grid);
00048   ConsS ***Soln;
00049   int i=0,j=0,k=0;
00050   int is,ie,js,je,ks,ke,n,m,nx1,nx2,nx3,Nx1,Nx2;
00051   Real amp,vflow,angle;
00052   Real d0,p0,u0,v0,w0,h0;
00053   Real x1,x2,x3,r,ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00054   Real x1size,x2size,sin_a,cos_a,lambda;
00055 #ifdef MHD
00056   Real bx0,by0,bz0,xfact,yfact;
00057   Real **az;
00058 #endif /* MHD */
00059   is = pGrid->is; ie = pGrid->ie;
00060   js = pGrid->js; je = pGrid->je;
00061   ks = pGrid->ks; ke = pGrid->ke;
00062   nx1 = (ie-is)+1 + 2*nghost;
00063   nx2 = (je-js)+1 + 2*nghost;
00064   nx3 = (ke-ks)+1 + 2*nghost;
00065 
00066 /* NOTE: For parallel calculations Nx1 != nx1 and Nx2 != nx2 */
00067 
00068   Nx1 = pDomain->Nx[0];
00069   Nx2 = pDomain->Nx[1];
00070   if (Nx1 == 1 || Nx2 == 1) {
00071     ath_error("[problem]: this test only works with Nx1 & Nx2 > 1\n");
00072   }
00073 
00074 /* allocate memory for solution and vector potential */
00075 
00076 #ifdef MHD
00077   if ((az = (Real**)calloc_2d_array(nx2, nx1, sizeof(Real))) == NULL)
00078     ath_error("[problem]: Error allocating memory for \"az\"\n");
00079 #endif /* MHD */
00080 
00081   if ((Soln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))==NULL)
00082     ath_error("[problem]: Error allocating memory\n");
00083   if (pDomain->Level == 0){
00084     if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))
00085       == NULL) ath_error("[problem]: Error alloc memory for RootSoln\n");
00086   }
00087 
00088 /* Read initial conditions */
00089 
00090   wave_flag = par_geti("problem","wave_flag");
00091   amp = par_getd("problem","amp");
00092   vflow = par_getd("problem","vflow");
00093 
00094 /* Set angle, wavelength */
00095 
00096   x1size = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00097   x2size = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00098   angle = atan(x1size/x2size);
00099   sin_a = sin(angle);
00100   cos_a = cos(angle);
00101 
00102 /* Use the larger angle to determine the wavelength */
00103   if (cos_a >= sin_a) {
00104     lambda = x1size*cos_a;
00105   } else {
00106     lambda = x2size*sin_a;
00107   }
00108 
00109 /* Get eigenmatrix, where the quantities u0 and bx0 are parallel to the
00110  * wavevector, and v0,w0,by0,bz0 are perpendicular. */
00111 
00112   d0 = 1.0;
00113 #ifndef ISOTHERMAL
00114   p0 = 1.0/Gamma;
00115   u0 = vflow*sqrt(Gamma*p0/d0);
00116 #else
00117   u0 = vflow*Iso_csound;
00118 #endif
00119   v0 = 0.0;
00120   w0 = 0.0;
00121 #ifdef MHD
00122   bx0 = 1.0;
00123   by0 = sqrt(2.0);
00124   bz0 = 0.5;
00125   xfact = 0.0;
00126   yfact = 1.0;
00127 #endif
00128 
00129   for (n=0; n<NWAVE; n++) {
00130     for (m=0; m<NWAVE; m++) {
00131       rem[n][m] = 0.0;
00132       lem[n][m] = 0.0;
00133     }
00134   }
00135 
00136 #ifdef HYDRO
00137 #ifdef ISOTHERMAL
00138   esys_roe_iso_hyd(u0,v0,w0,   ev,rem,lem);
00139 #else
00140   h0 = ((p0/Gamma_1 + 0.5*d0*(u0*u0+v0*v0+w0*w0)) + p0)/d0;
00141   esys_roe_adb_hyd(u0,v0,w0,h0,ev,rem,lem);
00142   ath_pout(0,"Ux - Cs = %e, %e\n",ev[0],rem[0][wave_flag]);
00143   ath_pout(0,"Ux      = %e, %e\n",ev[1],rem[1][wave_flag]);
00144   ath_pout(0,"Ux + Cs = %e, %e\n",ev[4],rem[4][wave_flag]);
00145 #endif /* ISOTHERMAL */
00146 #endif /* HYDRO */
00147 
00148 #ifdef MHD
00149 #if defined(ISOTHERMAL)
00150   esys_roe_iso_mhd(d0,u0,v0,w0,bx0,by0,bz0,xfact,yfact,ev,rem,lem);
00151   ath_pout(0,"Ux - Cf = %e, %e\n",ev[0],rem[0][wave_flag]);
00152   ath_pout(0,"Ux - Ca = %e, %e\n",ev[1],rem[1][wave_flag]);
00153   ath_pout(0,"Ux - Cs = %e, %e\n",ev[2],rem[2][wave_flag]);
00154   ath_pout(0,"Ux + Cs = %e, %e\n",ev[3],rem[3][wave_flag]);
00155   ath_pout(0,"Ux + Ca = %e, %e\n",ev[4],rem[4][wave_flag]);
00156   ath_pout(0,"Ux + Cf = %e, %e\n",ev[5],rem[5][wave_flag]);
00157 #else
00158   h0 = ((p0/Gamma_1+0.5*(bx0*bx0+by0*by0+bz0*bz0)+0.5*d0*(u0*u0+v0*v0+w0*w0))
00159                + (p0+0.5*(bx0*bx0+by0*by0+bz0*bz0)))/d0;
00160   esys_roe_adb_mhd(d0,u0,v0,w0,h0,bx0,by0,bz0,xfact,yfact,ev,rem,lem);
00161   ath_pout(0,"Ux - Cf = %e, %e\n",ev[0],rem[0][wave_flag]);
00162   ath_pout(0,"Ux - Ca = %e, %e\n",ev[1],rem[1][wave_flag]);
00163   ath_pout(0,"Ux - Cs = %e, %e\n",ev[2],rem[2][wave_flag]);
00164   ath_pout(0,"Ux      = %e, %e\n",ev[3],rem[3][wave_flag]);
00165   ath_pout(0,"Ux + Cs = %e, %e\n",ev[4],rem[4][wave_flag]);
00166   ath_pout(0,"Ux + Ca = %e, %e\n",ev[5],rem[5][wave_flag]);
00167   ath_pout(0,"Ux + Cf = %e, %e\n",ev[6],rem[6][wave_flag]);
00168 #endif /* ISOTHERMAL */
00169 #endif /* MHD */
00170 
00171 /* Now initialize 2D solution */
00172 /* Fields are initialized using vector potential in 2D (except B3) */
00173 
00174 #ifdef MHD
00175   for (j=js; j<=je+1; j++) {
00176     for (i=is; i<=ie+1; i++) {
00177       cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00178       x1 = x1 - 0.5*pGrid->dx1;
00179       x2 = x2 - 0.5*pGrid->dx2;
00180       r = (x1*cos_a + x2*sin_a)/lambda;
00181 
00182       az[j][i] = -bx0*(x1*sin_a - x2*cos_a)
00183                 - by0*(x1*cos_a + x2*sin_a)
00184                 + amp*lambda*cos(2.0*PI*r)*rem[NWAVE-2][wave_flag]/(2.0*PI);
00185     }
00186   }
00187 #endif /* MHD */
00188 
00189 /* The initial solution is stored into global 2D arrays so the error norm can
00190  * be computed at the end of the run.  After calculating these arrays, the
00191  * conserved variables are set equal to them */
00192 
00193   for (k=ks; k<=ke; k++) {
00194   for (j=js; j<=je; j++) {
00195   for (i=is; i<=ie; i++) {
00196     cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00197     r = (x1*cos_a + x2*sin_a)/lambda;
00198 
00199     Soln[k][j][i].d = d0 + amp*sin(2.0*PI*r)*rem[0][wave_flag];
00200 
00201 #ifndef ISOTHERMAL
00202 #ifdef HYDRO
00203     Soln[k][j][i].E = p0/Gamma_1 + 0.5*d0*u0*u0 
00204                     + amp*sin(2.0*PI*r)*rem[4][wave_flag];
00205 #else
00206     Soln[k][j][i].E = p0/Gamma_1 + 0.5*d0*u0*u0 + 0.5*(bx0*bx0+by0*by0+bz0*bz0)
00207                     + amp*sin(2.0*PI*r)*rem[4][wave_flag];
00208 #endif /* HYDRO */
00209 #endif /* ISOTHERMAL */
00210 
00211     Soln[k][j][i].M1 = d0*vflow*cos_a
00212                     + amp*sin(2.0*PI*r)*rem[1][wave_flag]*cos_a
00213                     - amp*sin(2.0*PI*r)*rem[2][wave_flag]*sin_a;
00214     Soln[k][j][i].M2 = d0*vflow*sin_a
00215                     + amp*sin(2.0*PI*r)*rem[1][wave_flag]*sin_a
00216                     + amp*sin(2.0*PI*r)*rem[2][wave_flag]*cos_a;
00217     Soln[k][j][i].M3 = amp*sin(2.0*PI*r)*rem[3][wave_flag];
00218 
00219 #ifdef MHD
00220     Soln[k][j][i].B3c = bz0 + amp*sin(2.0*PI*r)*rem[NWAVE-1][wave_flag];
00221 #endif /* MHD */
00222 
00223 #if (NSCALARS > 0)
00224     for (n=0; n<NSCALARS; n++)
00225       Soln[k][j][i].s[n] = amp*(1.0 + sin(2.0*PI*r));
00226 #endif
00227 
00228   }}}
00229 
00230 #ifdef MHD
00231 /* Set face-centered fields */
00232 
00233   for (k=ks; k<=ke; k++) {
00234   for (j=js; j<=je; j++) {
00235   for (i=is; i<=ie+1; i++) {
00236     pGrid->B1i[k][j][i] = (az[j+1][i] - az[j][i])/pGrid->dx2;
00237   }}}
00238   for (k=ks; k<=ke; k++) {
00239   for (j=js; j<=je+1; j++) {
00240   for (i=is; i<=ie; i++) {
00241     pGrid->B2i[k][j][i] =-(az[j][i+1] - az[j][i])/pGrid->dx1;
00242   }}}
00243   for (k=ks; k<=ke; k++) {
00244   for (j=js; j<=je; j++) {
00245   for (i=is; i<=ie; i++) {
00246     pGrid->B3i[k][j][i] = Soln[k][j][i].B3c;
00247   }}}
00248   if (pGrid->Nx[2] > 1) {
00249     for (j=js; j<=je; j++) {
00250     for (i=is; i<=ie; i++) {
00251         pGrid->B3i[ke+1][j][i] = Soln[ke][j][i].B3c;
00252     }}
00253   }
00254 
00255 /* compute cell-centered fields */
00256 
00257   for (k=ks; k<=ke; k++) {
00258   for (j=js; j<=je; j++) {
00259   for (i=is; i<=ie; i++) {
00260      Soln[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i] + pGrid->B1i[k][j][i+1]);
00261      Soln[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i] + pGrid->B2i[k][j+1][i]);
00262   }}}
00263   free_2d_array(az);
00264 #endif /* MHD */
00265 
00266 /* Now set initial conditions to 2d wave solution */ 
00267 
00268   for (k=ks; k<=ke; k++) {
00269   for (j=js; j<=je; j++) {
00270   for (i=is; i<=ie; i++) {
00271     pGrid->U[k][j][i].d = Soln[k][j][i].d;
00272 #ifndef ISOTHERMAL
00273     pGrid->U[k][j][i].E = Soln[k][j][i].E;
00274 #endif /* ISOTHERMAL */
00275     pGrid->U[k][j][i].M1 = Soln[k][j][i].M1;
00276     pGrid->U[k][j][i].M2 = Soln[k][j][i].M2;
00277     pGrid->U[k][j][i].M3 = Soln[k][j][i].M3;
00278 #ifdef MHD
00279     pGrid->U[k][j][i].B1c = Soln[k][j][i].B1c;
00280     pGrid->U[k][j][i].B2c = Soln[k][j][i].B2c;
00281     pGrid->U[k][j][i].B3c = Soln[k][j][i].B3c;
00282 #endif /* MHD */
00283 #if (NSCALARS > 0)
00284       for (n=0; n<NSCALARS; n++)
00285         pGrid->U[k][j][i].s[n] = Soln[k][j][i].s[n];
00286 #endif
00287   }}}
00288 
00289 /* For self-gravitating problems, read 4\piG and compute mean density */
00290 
00291 #ifdef SELF_GRAVITY
00292   four_pi_G = par_getd("problem","four_pi_G");
00293   grav_mean_rho = d0;
00294 #endif /* SELF_GRAVITY */
00295 
00296 /* With viscosity and/or resistivity, read eta_R and nu_V */
00297 
00298 #ifdef RESISTIVITY 
00299   eta_Ohm = par_getd("problem","eta_O");
00300   Q_AD    = par_getd_def("problem","Q_AD",0.0);
00301   d_ind   = par_getd_def("problem","d_ind",0.0);
00302 #endif
00303 #ifdef NAVIER_STOKES
00304   nu_V = par_getd("problem","nu");
00305 #endif
00306 
00307 /* save solution on root grid */
00308 
00309   if (pDomain->Level == 0) {
00310     for (k=ks; k<=ke; k++) {
00311     for (j=js; j<=je; j++) {
00312     for (i=is; i<=ie; i++) {
00313       RootSoln[k][j][i].d  = Soln[k][j][i].d ;
00314       RootSoln[k][j][i].M1 = Soln[k][j][i].M1;
00315       RootSoln[k][j][i].M2 = Soln[k][j][i].M2;
00316       RootSoln[k][j][i].M3 = Soln[k][j][i].M3;
00317 #ifndef ISOTHERMAL
00318       RootSoln[k][j][i].E  = Soln[k][j][i].E ;
00319 #endif /* ISOTHERMAL */
00320 #ifdef MHD
00321       RootSoln[k][j][i].B1c = Soln[k][j][i].B1c;
00322       RootSoln[k][j][i].B2c = Soln[k][j][i].B2c;
00323       RootSoln[k][j][i].B3c = Soln[k][j][i].B3c;
00324 #endif
00325 #if (NSCALARS > 0)
00326       for (n=0; n<NSCALARS; n++)
00327         RootSoln[k][j][i].s[n] = Soln[k][j][i].s[n];
00328 #endif
00329     }}}
00330   }
00331 
00332   free_3d_array(Soln);
00333 
00334   return;
00335 }
00336 
00337 /*==============================================================================
00338  * PROBLEM USER FUNCTIONS:
00339  * problem_write_restart() - writes problem-specific user data to restart files
00340  * problem_read_restart()  - reads problem-specific user data from restart files
00341  * get_usr_expr()          - sets pointer to expression for special output data
00342  * get_usr_out_fun()       - returns a user defined output function pointer
00343  * get_usr_par_prop()      - returns a user defined particle selection function
00344  * Userwork_in_loop        - problem specific work IN     main loop
00345  * Userwork_after_loop     - problem specific work AFTER  main loop
00346  * color()   - returns first passively advected scalar s[0]
00347  *----------------------------------------------------------------------------*/
00348 
00349 void problem_write_restart(MeshS *pM, FILE *fp)
00350 {
00351   return;
00352 }
00353 
00354 void problem_read_restart(MeshS *pM, FILE *fp)
00355 {
00356   return;
00357 }
00358 
00359 #if (NSCALARS > 0)
00360 /*! \fn static Real color(const GridS *pG, const int i, const int j,const int k)
00361  *  \brief returns first passively advected scalar s[0] */
00362 
00363 static Real color(const GridS *pG, const int i, const int j, const int k)
00364 {
00365   return pG->U[k][j][i].s[0]/pG->U[k][j][i].d;
00366 }
00367 #endif
00368 
00369 ConsFun_t get_usr_expr(const char *expr)
00370 {
00371 #if (NSCALARS > 0)
00372   if(strcmp(expr,"color")==0) return color;
00373 #endif
00374   return NULL;
00375 }
00376 
00377 VOutFun_t get_usr_out_fun(const char *name){
00378   return NULL;
00379 }
00380 
00381 #ifdef RESISTIVITY
00382 void get_eta_user(GridS *pG, int i, int j, int k,
00383                              Real *eta_O, Real *eta_H, Real *eta_A)
00384 {
00385   *eta_O = 0.0;
00386   *eta_H = 0.0;
00387   *eta_A = 0.0;
00388 
00389   return;
00390 }
00391 #endif
00392 
00393 void Userwork_in_loop(MeshS *pM)
00394 {
00395 }
00396 
00397 /*---------------------------------------------------------------------------
00398  * Userwork_after_loop: computes L1-error in linear waves,
00399  * ASSUMING WAVE HAS PROPAGATED AN INTEGER NUMBER OF PERIODS
00400  * Must set parameters in input file appropriately so that this is true
00401  */
00402 
00403 void Userwork_after_loop(MeshS *pM)
00404 {
00405   GridS *pGrid;
00406   int i=0,j=0,k=0;
00407 #if (NSCALARS > 0)
00408    int n;
00409 #endif
00410   int is,ie,js,je,ks,ke;
00411   Real rms_error=0.0;
00412   ConsS error,total_error;
00413   FILE *fp;
00414   char *fname;
00415   int Nx1, Nx2, Nx3, count;
00416 #if defined MPI_PARALLEL
00417   double err[8+NSCALARS], tot_err[8+NSCALARS];
00418   int ierr,myID;
00419 #endif
00420 
00421   total_error.d = 0.0;
00422   total_error.M1 = 0.0;
00423   total_error.M2 = 0.0;
00424   total_error.M3 = 0.0;
00425 #ifdef MHD
00426   total_error.B1c = 0.0;
00427   total_error.B2c = 0.0;
00428   total_error.B3c = 0.0;
00429 #endif /* MHD */
00430 #ifndef ISOTHERMAL
00431   total_error.E = 0.0;
00432 #endif /* ISOTHERMAL */
00433 #if (NSCALARS > 0)
00434   for (n=0; n<NSCALARS; n++) total_error.s[n] = 0.0;
00435 #endif
00436 
00437 /* Compute error only on root Grid, which is in Domain[0][0] */
00438 
00439   pGrid=pM->Domain[0][0].Grid;
00440   if (pGrid == NULL) return;
00441 
00442 /* compute L1 error in each variable, and rms total error */
00443 
00444   is = pGrid->is; ie = pGrid->ie;
00445   js = pGrid->js; je = pGrid->je;
00446   ks = pGrid->ks; ke = pGrid->ke;
00447   for (k=ks; k<=ke; k++) {
00448   for (j=js; j<=je; j++) {
00449     error.d = 0.0;
00450     error.M1 = 0.0;
00451     error.M2 = 0.0;
00452     error.M3 = 0.0;
00453 #ifdef MHD
00454     error.B1c = 0.0;
00455     error.B2c = 0.0;
00456     error.B3c = 0.0;
00457 #endif /* MHD */
00458 #ifndef ISOTHERMAL
00459     error.E = 0.0;
00460 #endif /* ISOTHERMAL */
00461 #if (NSCALARS > 0)
00462     for (n=0; n<NSCALARS; n++) error.s[n] = 0.0;
00463 #endif
00464 
00465     for (i=is; i<=ie; i++) {
00466       error.d   += fabs(pGrid->U[k][j][i].d   - RootSoln[k][j][i].d );
00467       error.M1  += fabs(pGrid->U[k][j][i].M1  - RootSoln[k][j][i].M1);
00468       error.M2  += fabs(pGrid->U[k][j][i].M2  - RootSoln[k][j][i].M2);
00469       error.M3  += fabs(pGrid->U[k][j][i].M3  - RootSoln[k][j][i].M3); 
00470 #ifdef MHD
00471       error.B1c += fabs(pGrid->U[k][j][i].B1c - RootSoln[k][j][i].B1c);
00472       error.B2c += fabs(pGrid->U[k][j][i].B2c - RootSoln[k][j][i].B2c);
00473       error.B3c += fabs(pGrid->U[k][j][i].B3c - RootSoln[k][j][i].B3c);
00474 #endif /* MHD */
00475 #ifndef ISOTHERMAL
00476       error.E   += fabs(pGrid->U[k][j][i].E   - RootSoln[k][j][i].E );
00477 #endif /* ISOTHERMAL */
00478 #if (NSCALARS > 0)
00479       for (n=0; n<NSCALARS; n++)
00480         error.s[n] += fabs(pGrid->U[k][j][i].s[n] - RootSoln[k][j][i].s[n]);;
00481 #endif
00482     }
00483 
00484     total_error.d += error.d;
00485     total_error.M1 += error.M1;
00486     total_error.M2 += error.M2;
00487     total_error.M3 += error.M3;
00488 #ifdef MHD
00489     total_error.B1c += error.B1c;
00490     total_error.B2c += error.B2c;
00491     total_error.B3c += error.B3c;
00492 #endif /* MHD */
00493 #ifndef ISOTHERMAL
00494     total_error.E += error.E;
00495 #endif /* ISOTHERMAL */
00496 #if (NSCALARS > 0)
00497     for (n=0; n<NSCALARS; n++) total_error.s[n] += error.s[n];
00498 #endif
00499   }}
00500 
00501 #ifdef MPI_PARALLEL
00502   Nx1 = pM->Domain[0][0].Nx[0];
00503   Nx2 = pM->Domain[0][0].Nx[1];
00504   Nx3 = pM->Domain[0][0].Nx[2];
00505 #else
00506   Nx1 = ie - is + 1;
00507   Nx2 = je - js + 1;
00508   Nx3 = ke - ks + 1;
00509 #endif
00510   count = Nx1*Nx2*Nx3;
00511 
00512 #ifdef MPI_PARALLEL 
00513 /* Now we have to use an All_Reduce to get the total error over all the MPI
00514  * grids.  Begin by copying the error into the err[] array */
00515 
00516   err[0] = total_error.d;
00517   err[1] = total_error.M1;
00518   err[2] = total_error.M2;
00519   err[3] = total_error.M3;
00520 #ifdef MHD
00521   err[4] = total_error.B1c;
00522   err[5] = total_error.B2c;
00523   err[6] = total_error.B3c;
00524 #endif /* MHD */
00525 #ifndef ISOTHERMAL
00526   err[7] = total_error.E;
00527 #endif /* ISOTHERMAL */
00528 #if (NSCALARS > 0)
00529   for (n=0; n<NSCALARS; n++) err[8+n] = total_error.s[n];
00530 #endif
00531 
00532 /* Sum up the Computed Error */
00533   ierr = MPI_Reduce(err,tot_err,(8+NSCALARS),MPI_DOUBLE,MPI_SUM,0,
00534     pM->Domain[0][0].Comm_Domain);
00535 
00536 /* If I'm the parent, copy the sum back to the total_error variable */
00537 
00538   ierr = MPI_Comm_rank(pM->Domain[0][0].Comm_Domain, &myID);
00539   if(myID == 0){ /* I'm the parent */
00540     total_error.d   = tot_err[0];
00541     total_error.M1  = tot_err[1];
00542     total_error.M2  = tot_err[2];
00543     total_error.M3  = tot_err[3];
00544 #ifdef MHD
00545     total_error.B1c = tot_err[4];
00546     total_error.B2c = tot_err[5];
00547     total_error.B3c = tot_err[6];
00548 #endif /* MHD */
00549 #ifndef ISOTHERMAL
00550     total_error.E   = tot_err[7];
00551 #endif /* ISOTHERMAL */
00552 #if (NSCALARS > 0)
00553   for (n=0; n<NSCALARS; n++) total_error.s[n] = err[8+n];
00554 #endif
00555 
00556   }
00557   else return; /* The child grids do not do any of the following code */
00558 
00559 #endif /* MPI_PARALLEL */
00560 
00561 /* Compute RMS error over all variables, and print out */
00562 
00563   rms_error = SQR(total_error.d) + SQR(total_error.M1) + SQR(total_error.M2)
00564                 + SQR(total_error.M3);
00565 #ifdef MHD
00566   rms_error += SQR(total_error.B1c) + SQR(total_error.B2c) 
00567                + SQR(total_error.B3c);
00568 #endif /* MHD */
00569 #ifndef ISOTHERMAL
00570   rms_error += SQR(total_error.E);
00571 #endif /* ISOTHERMAL */
00572   rms_error = sqrt(rms_error)/(double)count;
00573 
00574 
00575 /* Print error to file "LinWave-errors.#.dat", where #=wave_flag  */
00576 
00577 #ifdef MPI_PARALLEL
00578   fname = ath_fname("../","LinWave-errors",NULL,NULL,1,wave_flag,NULL,"dat");
00579 #else
00580   fname = ath_fname(NULL,"LinWave-errors",NULL,NULL,1,wave_flag,NULL,"dat");
00581 #endif
00582 
00583 /* The file exists -- reopen the file in append mode */
00584   if((fp=fopen(fname,"r")) != NULL){
00585     if((fp = freopen(fname,"a",fp)) == NULL){
00586       ath_error("[Userwork_after_loop]: Unable to reopen file.\n");
00587       free(fname);
00588       return;
00589     }
00590   }
00591 /* The file does not exist -- open the file in write mode */
00592   else{
00593     if((fp = fopen(fname,"w")) == NULL){
00594       ath_error("[Userwork_after_loop]: Unable to open file.\n");
00595       free(fname);
00596       return;
00597     }
00598 /* Now write out some header information */
00599     fprintf(fp,"# Nx1  Nx2  Nx3  RMS-Error  d  M1  M2  M3");
00600 #ifndef ISOTHERMAL
00601     fprintf(fp,"  E");
00602 #endif /* ISOTHERMAL */
00603 #ifdef MHD
00604     fprintf(fp,"  B1c  B2c  B3c");
00605 #endif /* MHD */
00606 #if (NSCALARS > 0)
00607     for (n=0; n<NSCALARS; n++) {
00608       fprintf(fp,"  S[ %d ]",n);
00609     }
00610 #endif
00611     fprintf(fp,"\n#\n");
00612   }
00613 
00614   fprintf(fp,"%d  %d  %d  %e",Nx1,Nx2,Nx3,rms_error);
00615 
00616   fprintf(fp,"  %e  %e  %e  %e",
00617           (total_error.d/(double) count),
00618           (total_error.M1/(double)count),
00619           (total_error.M2/(double)count),
00620           (total_error.M3/(double)count));
00621 
00622 #ifndef ISOTHERMAL
00623   fprintf(fp,"  %e",(total_error.E/(double)count));
00624 #endif /* ISOTHERMAL */
00625 
00626 #ifdef MHD
00627   fprintf(fp,"  %e  %e  %e",
00628           (total_error.B1c/(double)count),
00629           (total_error.B2c/(double)count),
00630           (total_error.B3c/(double)count));
00631 #endif /* MHD */
00632 #if (NSCALARS > 0)
00633     for (n=0; n<NSCALARS; n++) {
00634       fprintf(fp,"  %e",total_error.s[n]/(double)count);
00635     }
00636 #endif
00637 
00638   fprintf(fp,"\n");
00639 
00640   fclose(fp);
00641   free(fname);
00642 
00643   return;
00644 }