Athena: prob/linear_wave1d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file linear_wave1d.c
00004  *  \brief Problem generator for plane-parallel, grid-aligned linear wave tests.
00005  *
00006  * PURPOSE: Problem generator for plane-parallel, grid-aligned linear wave
00007  *   tests.  If wavevector is in x2 (x3) direction, then the grid must be
00008  *   2D (3D). Only works for wavevector parallel to grid.  Tests in 2D and 3D
00009  *   with wavevector at an arbitrary angle are initialized with different
00010  *   functions (linear_wave2d/3d).
00011  *
00012  *   Can be used for either standing (problem/vflow=1.0) or travelling
00013  *   (problem/vflow=0.0) waves.
00014  *
00015  *   Configure --with-gravity=fft to check Jeans stability of plane waves
00016  *   propagating parallel to grid.
00017  *
00018  *   Configure --enable-resistivity and/or --with-viscosity=ns to check damping
00019  *   of linear waves by resistivity/ambipolar diffusion and/or viscosity
00020  *
00021  * USERWORK_AFTER_LOOP function computes L1 error norm in solution by comparing
00022  *   to initial conditions.  Problem must be evolved for an integer number of
00023  *   wave periods for this to work (only works for ideal MHD).                */
00024 /*============================================================================*/
00025 
00026 #include <math.h>
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include "defs.h"
00030 #include "athena.h"
00031 #include "globals.h"
00032 #include "prototypes.h"
00033 
00034 /* Initial solution on root level, shared with Userwork_after_loop  */
00035 static ConsS ***RootSoln=NULL;
00036 static int wave_flag;
00037 
00038 
00039 /*----------------------------------------------------------------------------*/
00040 /* problem:   */
00041 
00042 void problem(DomainS *pDomain)
00043 {
00044   GridS *pGrid=(pDomain->Grid);
00045   ConsS ***Soln;
00046   int i=0,j=0,k=0;
00047   int is,ie,js,je,ks,ke,n,m,nx1,nx2,nx3,wave_dir;
00048   Real amp,vflow;
00049   Real d0,p0,u0,v0,w0,h0;
00050   Real x1,x2,x3,r,ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00051 #ifdef MHD
00052   Real bx0,by0,bz0,xfact,yfact;
00053 #endif /* MHD */
00054   is = pGrid->is; ie = pGrid->ie;
00055   js = pGrid->js; je = pGrid->je;
00056   ks = pGrid->ks; ke = pGrid->ke;
00057   nx1 = (ie-is)+1 + 2*nghost;
00058   nx2 = (je-js)+1 + 2*nghost;
00059   nx3 = (ke-ks)+1 + 2*nghost;
00060 
00061 /* Read initial conditions  */
00062 
00063   wave_flag = par_geti("problem","wave_flag");
00064   amp = par_getd("problem","amp");
00065   vflow = par_getd("problem","vflow");
00066   wave_dir = par_getd("problem","wave_dir");
00067 
00068 /* allocate memory for solution on this level.  If this is root level
00069  * also allocate memory for RootSoln */
00070 
00071   if ((Soln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))==NULL)
00072     ath_error("[problem]: Error allocating memory for Soln\n");
00073   if (pDomain->Level == 0){
00074     if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))
00075       == NULL) ath_error("[problem]: Error alloc memory for RootSoln\n");
00076   }
00077 
00078 /* Get eigenmatrix, where the quantities u0 and bx0 are parallel to the
00079  * wavevector, and v0,w0,by0,bz0 are perpendicular.  Background state chosen
00080  * carefully so wave speeds are well-separated: Va=1, Vf=2, Vs=0.5 */
00081 
00082   d0 = 1.0;
00083 #ifndef ISOTHERMAL
00084   p0 = 1.0/Gamma;
00085   u0 = vflow*sqrt(Gamma*p0/d0);
00086 #else
00087   u0 = vflow*Iso_csound;
00088 #endif
00089   v0 = 0.0;
00090   w0 = 0.0;
00091 #ifdef MHD
00092   bx0 = 1.0;
00093   by0 = sqrt(2.0);
00094   bz0 = 0.5;
00095   xfact = 0.0;
00096   yfact = 1.0;
00097 #endif
00098 
00099   for (n=0; n<NWAVE; n++) {
00100     for (m=0; m<NWAVE; m++) {
00101       rem[n][m] = 0.0;
00102       lem[n][m] = 0.0;
00103     }
00104   }
00105 
00106 #ifdef HYDRO
00107 #ifdef ISOTHERMAL
00108   esys_roe_iso_hyd(u0,v0,w0,   ev,rem,lem);
00109 #else
00110   h0 = ((p0/Gamma_1 + 0.5*d0*(u0*u0+v0*v0+w0*w0)) + p0)/d0;
00111   esys_roe_adb_hyd(u0,v0,w0,h0,ev,rem,lem);
00112   ath_pout(0,"Ux - Cs = %e, %e\n",ev[0],rem[0][wave_flag]);
00113   ath_pout(0,"Ux      = %e, %e\n",ev[1],rem[1][wave_flag]);
00114   ath_pout(0,"Ux + Cs = %e, %e\n",ev[4],rem[4][wave_flag]);
00115 #endif /* ISOTHERMAL */
00116 #endif /* HYDRO */
00117 
00118 #ifdef MHD
00119 #if defined(ISOTHERMAL)
00120   esys_roe_iso_mhd(d0,u0,v0,w0,bx0,by0,bz0,xfact,yfact,ev,rem,lem);
00121   ath_pout(0,"Ux - Cf = %e, %e\n",ev[0],rem[0][wave_flag]);
00122   ath_pout(0,"Ux - Ca = %e, %e\n",ev[1],rem[1][wave_flag]);
00123   ath_pout(0,"Ux - Cs = %e, %e\n",ev[2],rem[2][wave_flag]);
00124   ath_pout(0,"Ux + Cs = %e, %e\n",ev[3],rem[3][wave_flag]);
00125   ath_pout(0,"Ux + Ca = %e, %e\n",ev[4],rem[4][wave_flag]);
00126   ath_pout(0,"Ux + Cf = %e, %e\n",ev[5],rem[5][wave_flag]);
00127 #else
00128   h0 = ((p0/Gamma_1+0.5*(bx0*bx0+by0*by0+bz0*bz0)+0.5*d0*(u0*u0+v0*v0+w0*w0))
00129                + (p0+0.5*(bx0*bx0+by0*by0+bz0*bz0)))/d0;
00130   esys_roe_adb_mhd(d0,u0,v0,w0,h0,bx0,by0,bz0,xfact,yfact,ev,rem,lem);
00131   ath_pout(0,"Ux - Cf = %e, %e\n",ev[0],rem[0][wave_flag]);
00132   ath_pout(0,"Ux - Ca = %e, %e\n",ev[1],rem[1][wave_flag]);
00133   ath_pout(0,"Ux - Cs = %e, %e\n",ev[2],rem[2][wave_flag]);
00134   ath_pout(0,"Ux      = %e, %e\n",ev[3],rem[3][wave_flag]);
00135   ath_pout(0,"Ux + Cs = %e, %e\n",ev[4],rem[4][wave_flag]);
00136   ath_pout(0,"Ux + Ca = %e, %e\n",ev[5],rem[5][wave_flag]);
00137   ath_pout(0,"Ux + Cf = %e, %e\n",ev[6],rem[6][wave_flag]);
00138 #endif /* ISOTHERMAL */
00139 #endif /* MHD */
00140 
00141 /* Now initialize 1D solution vector  */
00142 
00143   for (k=ks; k<=ke; k++) {
00144   for (j=js; j<=je; j++) {
00145   for (i=is; i<=ie; i++) {
00146 
00147 /* Set background state */
00148 
00149     cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00150     Soln[k][j][i].d = d0;
00151 #ifndef ISOTHERMAL
00152     Soln[k][j][i].E = p0/Gamma_1 + 0.5*d0*u0*u0;
00153 #ifdef MHD
00154     Soln[k][j][i].E += 0.5*(bx0*bx0 + by0*by0 + bz0*bz0);
00155 #endif /* MHD */
00156 #endif /* ISOTHERMAL */
00157 
00158 /* Select appropriate solution based on direction of wavevector */
00159 
00160     switch(wave_dir){
00161 
00162 /* wave in x1-direction */
00163 
00164     case 1:
00165       Soln[k][j][i].d += amp*sin(2.0*PI*x1)*rem[0][wave_flag];
00166 #ifndef ISOTHERMAL
00167       Soln[k][j][i].E += amp*sin(2.0*PI*x1)*rem[4][wave_flag];
00168 #endif /* ISOTHERMAL */
00169       Soln[k][j][i].M1 = d0*vflow;
00170       Soln[k][j][i].M2 = 0.0;
00171       Soln[k][j][i].M3 = 0.0;
00172       Soln[k][j][i].M1 += amp*sin(2.0*PI*x1)*rem[1][wave_flag];
00173       Soln[k][j][i].M2 += amp*sin(2.0*PI*x1)*rem[2][wave_flag];
00174       Soln[k][j][i].M3 += amp*sin(2.0*PI*x1)*rem[3][wave_flag];
00175 #ifdef MHD
00176       Soln[k][j][i].B1c = bx0;
00177       Soln[k][j][i].B2c = by0 + amp*sin(2.0*PI*x1)*rem[NWAVE-2][wave_flag];
00178       Soln[k][j][i].B3c = bz0 + amp*sin(2.0*PI*x1)*rem[NWAVE-1][wave_flag];
00179 #endif /* MHD */
00180 #if (NSCALARS > 0)
00181       for (n=0; n<NSCALARS; n++)
00182         Soln[k][j][i].s[n] = amp*(1.0 + sin(2.0*PI*x1));
00183 #endif
00184       break;
00185 
00186 /* Wave in x2-direction */
00187 
00188     case 2:
00189       Soln[k][j][i].d += amp*sin(2.0*PI*x2)*rem[0][wave_flag];
00190 #ifndef ISOTHERMAL
00191       Soln[k][j][i].E += amp*sin(2.0*PI*x2)*rem[4][wave_flag];
00192 #endif /* ISOTHERMAL */
00193       Soln[k][j][i].M1 = 0.0;
00194       Soln[k][j][i].M2 = d0*vflow;
00195       Soln[k][j][i].M3 = 0.0;
00196       Soln[k][j][i].M1 += amp*sin(2.0*PI*x2)*rem[3][wave_flag];
00197       Soln[k][j][i].M2 += amp*sin(2.0*PI*x2)*rem[1][wave_flag];
00198       Soln[k][j][i].M3 += amp*sin(2.0*PI*x2)*rem[2][wave_flag];
00199 #ifdef MHD
00200       Soln[k][j][i].B1c = bz0 + amp*sin(2.0*PI*x2)*rem[NWAVE-1][wave_flag];
00201       Soln[k][j][i].B2c = bx0;
00202       Soln[k][j][i].B3c = by0 + amp*sin(2.0*PI*x2)*rem[NWAVE-2][wave_flag];
00203 #endif /* MHD */
00204 #if (NSCALARS > 0)
00205       for (n=0; n<NSCALARS; n++)
00206         Soln[k][j][i].s[n] = amp*(1.0 + sin(2.0*PI*x2));
00207 #endif
00208       break;
00209 
00210 /* Wave in x3-direction */
00211 
00212     case 3:
00213       Soln[k][j][i].d += amp*sin(2.0*PI*x3)*rem[0][wave_flag];
00214 #ifndef ISOTHERMAL
00215       Soln[k][j][i].E += amp*sin(2.0*PI*x3)*rem[4][wave_flag];
00216 #endif /* ISOTHERMAL */
00217       Soln[k][j][i].M1 = 0.0;
00218       Soln[k][j][i].M2 = 0.0;
00219       Soln[k][j][i].M3 = d0*vflow;
00220       Soln[k][j][i].M1 += amp*sin(2.0*PI*x3)*rem[2][wave_flag];
00221       Soln[k][j][i].M2 += amp*sin(2.0*PI*x3)*rem[3][wave_flag];
00222       Soln[k][j][i].M3 += amp*sin(2.0*PI*x3)*rem[1][wave_flag];
00223 #ifdef MHD
00224       Soln[k][j][i].B1c = by0 + amp*sin(2.0*PI*x3)*rem[NWAVE-2][wave_flag];
00225       Soln[k][j][i].B2c = bz0 + amp*sin(2.0*PI*x3)*rem[NWAVE-1][wave_flag];
00226       Soln[k][j][i].B3c = bx0;
00227 #endif /* MHD */
00228 #if (NSCALARS > 0)
00229       for (n=0; n<NSCALARS; n++)
00230         Soln[k][j][i].s[n] = amp*(1.0 + sin(2.0*PI*x3));
00231 #endif
00232       break;
00233     default:
00234       ath_error("[linear_wave1d]: wave direction %d not allowed\n",wave_dir);
00235     }
00236   }}}
00237 
00238 /* Now set initial conditions to wave solution */ 
00239 
00240   for (k=ks; k<=ke; k++) {
00241   for (j=js; j<=je; j++) {
00242   for (i=is; i<=ie; i++) {
00243     pGrid->U[k][j][i].d = Soln[k][j][i].d;
00244 #ifndef ISOTHERMAL
00245     pGrid->U[k][j][i].E = Soln[k][j][i].E;
00246 #endif /* ISOTHERMAL */
00247     pGrid->U[k][j][i].M1 = Soln[k][j][i].M1;
00248     pGrid->U[k][j][i].M2 = Soln[k][j][i].M2;
00249     pGrid->U[k][j][i].M3 = Soln[k][j][i].M3;
00250 #ifdef MHD
00251     pGrid->U[k][j][i].B1c = Soln[k][j][i].B1c;
00252     pGrid->U[k][j][i].B2c = Soln[k][j][i].B2c;
00253     pGrid->U[k][j][i].B3c = Soln[k][j][i].B3c;
00254     pGrid->B1i[k][j][i] = Soln[k][j][i].B1c;
00255     pGrid->B2i[k][j][i] = Soln[k][j][i].B2c;
00256     pGrid->B3i[k][j][i] = Soln[k][j][i].B3c;
00257 #endif /* MHD */
00258 #if (NSCALARS > 0)
00259       for (n=0; n<NSCALARS; n++)
00260         pGrid->U[k][j][i].s[n] = Soln[k][j][i].s[n]; 
00261 #endif
00262   }}}
00263 #ifdef MHD
00264   for (k=ks; k<=ke; k++) {
00265     for (j=js; j<=je; j++) {
00266       pGrid->B1i[k][j][ie+1] = pGrid->B1i[k][j][ie];
00267     }
00268   }
00269   if (pGrid->Nx[1] > 1) {
00270     for (k=ks; k<=ke; k++) {
00271       for (i=is; i<=ie; i++) {
00272         pGrid->B2i[k][je+1][i] = pGrid->B2i[k][je][i];
00273       }
00274     }
00275   }
00276   if (pGrid->Nx[2] > 1) {
00277     for (j=js; j<=je; j++) {
00278       for (i=is; i<=ie; i++) {
00279         pGrid->B3i[ke+1][j][i] = pGrid->B3i[ke][j][i];
00280       }
00281     }
00282   }
00283 #endif /* MHD */
00284 
00285 /* For self-gravitating problems, read 4\piG and compute mean density */
00286 
00287 #ifdef SELF_GRAVITY
00288   four_pi_G = par_getd("problem","four_pi_G");
00289   grav_mean_rho = d0;
00290 #endif /* SELF_GRAVITY */
00291 
00292 /* With viscosity and/or resistivity, read eta_R and nu_V */
00293 
00294 #ifdef RESISTIVITY
00295   eta_Ohm = par_getd("problem","eta_O");
00296   Q_AD    = par_getd_def("problem","Q_AD",0.0);
00297   d_ind   = par_getd_def("problem","d_ind",0.0);
00298 #endif
00299 #ifdef NAVIER_STOKES
00300   nu_V = par_getd("problem","nu");
00301 #endif
00302 
00303 /* save solution on root grid */
00304 
00305   if (pDomain->Level == 0) {
00306     for (k=ks; k<=ke; k++) {
00307     for (j=js; j<=je; j++) {
00308     for (i=is; i<=ie; i++) {
00309       RootSoln[k][j][i].d  = Soln[k][j][i].d ;
00310       RootSoln[k][j][i].M1 = Soln[k][j][i].M1;
00311       RootSoln[k][j][i].M2 = Soln[k][j][i].M2;
00312       RootSoln[k][j][i].M3 = Soln[k][j][i].M3;
00313 #ifndef ISOTHERMAL
00314       RootSoln[k][j][i].E  = Soln[k][j][i].E ;
00315 #endif /* ISOTHERMAL */
00316 #ifdef MHD
00317       RootSoln[k][j][i].B1c = Soln[k][j][i].B1c;
00318       RootSoln[k][j][i].B2c = Soln[k][j][i].B2c;
00319       RootSoln[k][j][i].B3c = Soln[k][j][i].B3c;
00320 #endif
00321 #if (NSCALARS > 0)
00322       for (n=0; n<NSCALARS; n++) 
00323         RootSoln[k][j][i].s[n] = Soln[k][j][i].s[n];
00324 #endif
00325     }}}
00326   }
00327 
00328   free_3d_array(Soln);
00329 
00330   return;
00331 }
00332 
00333 /*==============================================================================
00334  * PROBLEM USER FUNCTIONS:
00335  * problem_write_restart() - writes problem-specific user data to restart files
00336  * problem_read_restart()  - reads problem-specific user data from restart files
00337  * get_usr_expr()          - sets pointer to expression for special output data
00338  * get_usr_out_fun()       - returns a user defined output function pointer
00339  * get_usr_par_prop()      - returns a user defined particle selection function
00340  * Userwork_in_loop        - problem specific work IN     main loop
00341  * Userwork_after_loop     - problem specific work AFTER  main loop
00342  *----------------------------------------------------------------------------*/
00343 
00344 void problem_write_restart(MeshS *pM, FILE *fp)
00345 {
00346   return;
00347 }
00348 
00349 void problem_read_restart(MeshS *pM, FILE *fp)
00350 {
00351 #ifdef SELF_GRAVITY
00352   Real d0 = 1.0;
00353   four_pi_G = par_getd("problem","four_pi_G");
00354   grav_mean_rho = d0;
00355 #endif /* SELF_GRAVITY */
00356 #ifdef RESISTIVITY
00357   eta_Ohm = par_getd_def("problem","eta_O",0.0);
00358   Q_Hall  = par_getd_def("problem","Q_H",0.0);
00359   Q_AD    = par_getd_def("problem","Q_AD",0.0);
00360 #endif
00361 #ifdef NAVIER_STOKES
00362   nu_V = par_getd("problem","nu");
00363 #endif
00364   return;
00365 }
00366 
00367 ConsFun_t get_usr_expr(const char *expr)
00368 {
00369   return NULL;
00370 }
00371 
00372 VOutFun_t get_usr_out_fun(const char *name){
00373   return NULL;
00374 }
00375 
00376 #ifdef RESISTIVITY
00377 /*! \fn void get_eta_user(GridS *pG, int i, int j, int k,
00378  *                           Real *eta_O, Real *eta_H, Real *eta_A)
00379  *  \brief Get user defined resistivity. */
00380 void get_eta_user(GridS *pG, int i, int j, int k,
00381                              Real *eta_O, Real *eta_H, Real *eta_A)
00382 {
00383   *eta_O = 0.0;
00384   *eta_H = 0.0;
00385   *eta_A = 0.0;
00386 
00387   return;
00388 }
00389 #endif
00390 
00391 void Userwork_in_loop(MeshS *pM)
00392 {
00393 }
00394 
00395 /*------------------------------------------------------------------------------
00396  * Userwork_after_loop: computes L1-error in linear waves,
00397  * ASSUMING WAVE HAS PROPAGATED AN INTEGER NUMBER OF PERIODS
00398  * Must set parameters (tlim, wavespeeds, box size) in input file appropriately
00399  * so that this is true
00400  * With nested grids, computes L1 error only on root grid (level=0).
00401  */
00402 
00403 void Userwork_after_loop(MeshS *pM)
00404 {
00405   GridS *pGrid;
00406   int i=0,j=0,k=0;
00407   int is,ie,js,je,ks,ke;
00408 #if (NSCALARS > 0)
00409    int n;
00410 #endif
00411 
00412   Real rms_error=0.0;
00413   ConsS error,total_error;
00414   FILE *fp;
00415   char *fname;
00416   int Nx1, Nx2, Nx3, count;
00417 
00418   total_error.d = 0.0;
00419   total_error.M1 = 0.0;
00420   total_error.M2 = 0.0;
00421   total_error.M3 = 0.0;
00422 #ifdef MHD
00423   total_error.B1c = 0.0;
00424   total_error.B2c = 0.0;
00425   total_error.B3c = 0.0;
00426 #endif /* MHD */
00427 #ifndef ISOTHERMAL
00428   total_error.E = 0.0;
00429 #endif /* ISOTHERMAL */
00430 #if (NSCALARS > 0)
00431   for (n=0; n<NSCALARS; n++) total_error.s[n] = 0.0;
00432 #endif
00433 
00434 /* Compute error only on root Grid, which is in Domain[0][0] */
00435 
00436   pGrid=pM->Domain[0][0].Grid;
00437   if (pGrid == NULL) return;
00438 
00439 /* compute L1 error in each variable, and rms total error */
00440 
00441   is = pGrid->is; ie = pGrid->ie;
00442   js = pGrid->js; je = pGrid->je;
00443   ks = pGrid->ks; ke = pGrid->ke;
00444   for (k=ks; k<=ke; k++) {
00445   for (j=js; j<=je; j++) {
00446     error.d = 0.0;
00447     error.M1 = 0.0;
00448     error.M2 = 0.0;
00449     error.M3 = 0.0;
00450 #ifdef MHD
00451     error.B1c = 0.0;
00452     error.B2c = 0.0;
00453     error.B3c = 0.0;
00454 #endif /* MHD */
00455 #ifndef ISOTHERMAL
00456     error.E = 0.0;
00457 #endif /* ISOTHERMAL */
00458 #if (NSCALARS > 0)
00459     for (n=0; n<NSCALARS; n++) error.s[n] = 0.0;
00460 #endif
00461 
00462     for (i=is; i<=ie; i++) {
00463       error.d   += fabs(pGrid->U[k][j][i].d   - RootSoln[k][j][i].d );
00464       error.M1  += fabs(pGrid->U[k][j][i].M1  - RootSoln[k][j][i].M1);
00465       error.M2  += fabs(pGrid->U[k][j][i].M2  - RootSoln[k][j][i].M2);
00466       error.M3  += fabs(pGrid->U[k][j][i].M3  - RootSoln[k][j][i].M3); 
00467 #ifdef MHD
00468       error.B1c += fabs(pGrid->U[k][j][i].B1c - RootSoln[k][j][i].B1c);
00469       error.B2c += fabs(pGrid->U[k][j][i].B2c - RootSoln[k][j][i].B2c);
00470       error.B3c += fabs(pGrid->U[k][j][i].B3c - RootSoln[k][j][i].B3c);
00471 #endif /* MHD */
00472 #ifndef ISOTHERMAL
00473       error.E   += fabs(pGrid->U[k][j][i].E   - RootSoln[k][j][i].E );
00474 #endif /* ISOTHERMAL */
00475 #if (NSCALARS > 0)
00476       for (n=0; n<NSCALARS; n++) 
00477         error.s[n] += fabs(pGrid->U[k][j][i].s[n] - RootSoln[k][j][i].s[n]);;
00478 #endif
00479     }
00480 
00481     total_error.d += error.d;
00482     total_error.M1 += error.M1;
00483     total_error.M2 += error.M2;
00484     total_error.M3 += error.M3;
00485 #ifdef MHD
00486     total_error.B1c += error.B1c;
00487     total_error.B2c += error.B2c;
00488     total_error.B3c += error.B3c;
00489 #endif /* MHD */
00490 #ifndef ISOTHERMAL
00491     total_error.E += error.E;
00492 #endif /* ISOTHERMAL */
00493 #if (NSCALARS > 0)
00494     for (n=0; n<NSCALARS; n++) total_error.s[n] += error.s[n];
00495 #endif
00496   }}
00497 
00498   Nx1 = ie - is + 1;
00499   Nx2 = je - js + 1;
00500   Nx3 = ke - ks + 1;
00501 
00502   count = Nx1*Nx2*Nx3;
00503 
00504 /* Compute RMS error over all variables, and print out */
00505 
00506   rms_error = SQR(total_error.d) + SQR(total_error.M1) + SQR(total_error.M2)
00507                 + SQR(total_error.M3);
00508 #ifdef MHD
00509   rms_error += SQR(total_error.B1c) + SQR(total_error.B2c) 
00510                + SQR(total_error.B3c);
00511 #endif /* MHD */
00512 #ifndef ISOTHERMAL
00513   rms_error += SQR(total_error.E);
00514 #endif /* ISOTHERMAL */
00515   rms_error = sqrt(rms_error)/(double)count;
00516 
00517 /* Print error to file "LinWave-errors.#.dat", where #=wave_flag  */
00518 
00519   fname = ath_fname(NULL,"LinWave-errors",NULL,NULL,1,wave_flag,NULL,"dat");
00520 /* The file exists -- reopen the file in append mode */
00521   if((fp=fopen(fname,"r")) != NULL){
00522     if((fp = freopen(fname,"a",fp)) == NULL){
00523       ath_error("[Userwork_after_loop]: Unable to reopen file.\n");
00524       return;
00525     }
00526   }
00527 /* The file does not exist -- open the file in write mode */
00528   else{
00529     if((fp = fopen(fname,"w")) == NULL){
00530       ath_error("[Userwork_after_loop]: Unable to open file.\n");
00531       return;
00532     }
00533 /* Now write out some header information */
00534     fprintf(fp,"# Nx1  Nx2  Nx3  RMS-Error  d  M1  M2  M3");
00535 #ifndef ISOTHERMAL
00536     fprintf(fp,"  E");
00537 #endif /* ISOTHERMAL */
00538 #ifdef MHD
00539     fprintf(fp,"  B1c  B2c  B3c");
00540 #endif /* MHD */
00541 #if (NSCALARS > 0)
00542     for (n=0; n<NSCALARS; n++) {
00543       fprintf(fp,"  S[ %d ]",n);
00544     }
00545 #endif
00546 
00547     fprintf(fp,"\n#\n");
00548   }
00549  
00550   fprintf(fp,"%d  %d  %d  %e",Nx1,Nx2,Nx3,rms_error);
00551 
00552   fprintf(fp,"  %e  %e  %e  %e",
00553           (total_error.d/(double) count),
00554           (total_error.M1/(double)count),
00555           (total_error.M2/(double)count),
00556           (total_error.M3/(double)count));
00557 #ifndef ISOTHERMAL
00558   fprintf(fp,"  %e",(total_error.E/(double)count));
00559 #endif /* ISOTHERMAL */
00560 #ifdef MHD
00561   fprintf(fp,"  %e  %e  %e",
00562           (total_error.B1c/(double)count),
00563           (total_error.B2c/(double)count),
00564           (total_error.B3c/(double)count));
00565 #endif /* MHD */
00566 #if (NSCALARS > 0)
00567     for (n=0; n<NSCALARS; n++) {
00568       fprintf(fp,"  %e",total_error.s[n]/(double)count);
00569     }
00570 #endif
00571   fprintf(fp,"\n");
00572 
00573   fclose(fp);
00574 
00575   return;
00576 }