00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <math.h>
00023 #include <stdio.h>
00024 #include <stdlib.h>
00025 #include "defs.h"
00026 #include "athena.h"
00027 #include "globals.h"
00028 #include "prototypes.h"
00029
00030 #ifndef MHD
00031 #error : The cpaw1d test only works for mhd.
00032 #endif
00033
00034
00035 static ConsS *RootSoln=NULL;
00036
00037
00038
00039
00040 void problem(DomainS *pDomain)
00041 {
00042 GridS *pGrid=(pDomain->Grid);
00043 int i, is = pGrid->is, ie = pGrid->ie;
00044 int j, js = pGrid->js;
00045 int k, ks = pGrid->ks;
00046 ConsS *Soln;
00047 int dir;
00048 Real x1,x2,x3,cs,sn,b_par,b_perp,fac,lambda,k_par,v_par,v_perp,den,pres;
00049 #ifdef RESISTIVITY
00050 Real v_A, kva, omega_h, omega_l, omega_r;
00051 #endif
00052
00053 if ((Soln = (ConsS*)malloc(((ie-is+1)+2*nghost)*sizeof(ConsS))) == NULL)
00054 ath_error("[problem] Error initializing Soln array");
00055
00056 if (pDomain->Level == 0) {
00057 if ((RootSoln = (ConsS*)malloc(((ie-is+1)+2*nghost)*sizeof(ConsS)))
00058 == NULL) ath_error("[problem] Error initializing RootSoln array");
00059 }
00060
00061 if (pGrid->Nx[1] > 1 || pGrid->Nx[2] > 1) {
00062 ath_error("[cpaw1d] grid must be 1D");
00063 }
00064
00065
00066
00067 lambda = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00068 k_par = 2.0*PI/lambda;
00069
00070 b_par = par_getd("problem","b_par");
00071 den = 1.0;
00072 b_perp = par_getd("problem","b_perp");
00073 v_perp = b_perp/sqrt((double)den);
00074
00075 dir = par_geti_def("problem","dir",1);
00076 if (dir == 1)
00077 fac = 1.0;
00078 else
00079 fac = -1.0;
00080
00081 #ifdef RESISTIVITY
00082 Q_Hall = par_getd("problem","Q_H");
00083 d_ind = 0.0;
00084 v_A = b_par/sqrt((double)den);
00085 if (Q_Hall > 0.0)
00086 {
00087 kva = k_par*v_A;
00088 omega_h = 1.0/Q_Hall;
00089
00090 omega_r = 0.5*SQR(kva)/omega_h*(sqrt(1.0+SQR(2.0*omega_h/kva)) + 1.0);
00091 omega_l = 0.5*SQR(kva)/omega_h*(sqrt(1.0+SQR(2.0*omega_h/kva)) - 1.0);
00092
00093 if (dir == 1)
00094 v_perp = v_perp * kva / omega_r;
00095 else
00096 v_perp = v_perp * kva / omega_l;
00097 }
00098 #endif
00099
00100
00101 pres = par_getd("problem","pres");
00102 v_par = par_getd("problem","v_par");
00103
00104
00105
00106 for (i=is; i<=ie+1; i++) {
00107 cc_pos(pGrid,i,js,ks,&x1,&x2,&x3);
00108
00109 sn = sin(k_par*x1);
00110 cs = cos(k_par*x1);
00111
00112 Soln[i].d = den;
00113
00114 Soln[i].M1 = den*v_par;
00115 Soln[i].M2 = -fac*den*v_perp*sn;
00116 Soln[i].M3 = -den*v_perp*cs;
00117
00118 Soln[i].B1c = b_par;
00119 Soln[i].B2c = fac*b_perp*sn;
00120 Soln[i].B3c = b_perp*cs;
00121 #ifndef ISOTHERMAL
00122 Soln[i].E = pres/Gamma_1
00123 + 0.5*den*(v_par*v_par + v_perp*v_perp)
00124 + 0.5*(b_par*b_par + b_perp*b_perp);
00125 #endif
00126 }
00127
00128
00129
00130 for (i=is; i<=ie+1; i++) {
00131 pGrid->U[ks][js][i].d = Soln[i].d;
00132 pGrid->U[ks][js][i].M1 = Soln[i].M1;
00133 pGrid->U[ks][js][i].M2 = Soln[i].M2;
00134 pGrid->U[ks][js][i].M3 = Soln[i].M3;
00135
00136 pGrid->U[ks][js][i].B1c = pGrid->B1i[ks][js][i] = Soln[i].B1c;
00137 pGrid->U[ks][js][i].B2c = pGrid->B2i[ks][js][i] = Soln[i].B2c;
00138 pGrid->U[ks][js][i].B3c = pGrid->B3i[ks][js][i] = Soln[i].B3c;
00139 #ifndef ISOTHERMAL
00140 pGrid->U[ks][js][i].E = Soln[i].E;
00141 #endif
00142 }
00143
00144
00145
00146 if (pDomain->Level == 0) {
00147 for (i=is; i<=ie+1; i++) {
00148 RootSoln[i].d = Soln[i].d;
00149 RootSoln[i].M1 = Soln[i].M1;
00150 RootSoln[i].M2 = Soln[i].M2;
00151 RootSoln[i].M3 = Soln[i].M3;
00152
00153 RootSoln[i].B1c = Soln[i].B1c;
00154 RootSoln[i].B2c = Soln[i].B2c;
00155 RootSoln[i].B3c = Soln[i].B3c;
00156 #ifndef ISOTHERMAL
00157 RootSoln[i].E = Soln[i].E;
00158 #endif
00159 }
00160 }
00161
00162 free(Soln);
00163 return;
00164 }
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 void problem_write_restart(MeshS *pM, FILE *fp)
00178 {
00179 return;
00180 }
00181
00182 void problem_read_restart(MeshS *pM, FILE *fp)
00183 {
00184 return;
00185 }
00186
00187 ConsFun_t get_usr_expr(const char *expr)
00188 {
00189 return NULL;
00190 }
00191
00192 VOutFun_t get_usr_out_fun(const char *name){
00193 return NULL;
00194 }
00195
00196 #ifdef RESISTIVITY
00197
00198
00199
00200 void get_eta_user(GridS *pG, int i, int j, int k,
00201 Real *eta_O, Real *eta_H, Real *eta_A)
00202 {
00203 *eta_O = 0.0;
00204 *eta_H = 0.0;
00205 *eta_A = 0.0;
00206
00207 return;
00208 }
00209 #endif
00210
00211
00212 void Userwork_in_loop(MeshS *pM)
00213 {
00214 }
00215
00216
00217
00218
00219
00220
00221
00222 void Userwork_after_loop(MeshS *pM)
00223 {
00224 GridS *pGrid;
00225 int i=0,is,ie,js,ks,Nx1;
00226 Real rms_error=0.0;
00227 ConsS error;
00228 FILE *fp;
00229 char *fname;
00230
00231 error.d = 0.0;
00232 error.M1 = 0.0;
00233 error.M2 = 0.0;
00234 error.M3 = 0.0;
00235 error.B1c = 0.0;
00236 error.B2c = 0.0;
00237 error.B3c = 0.0;
00238 #ifndef ISOTHERMAL
00239 error.E = 0.0;
00240 #endif
00241
00242
00243
00244 pGrid=pM->Domain[0][0].Grid;
00245 is = pGrid->is; ie = pGrid->ie;
00246 js = pGrid->js;
00247 ks = pGrid->ks;
00248 Nx1 = (ie-is+1);
00249
00250
00251
00252 for (i=is; i<=ie; i++) {
00253 error.d += fabs(pGrid->U[ks][js][i].d - RootSoln[i].d );
00254 error.M1 += fabs(pGrid->U[ks][js][i].M1 - RootSoln[i].M1 );
00255 error.M2 += fabs(pGrid->U[ks][js][i].M2 - RootSoln[i].M2 );
00256 error.M3 += fabs(pGrid->U[ks][js][i].M3 - RootSoln[i].M3 );
00257 error.B1c += fabs(pGrid->U[ks][js][i].B1c - RootSoln[i].B1c);
00258 error.B2c += fabs(pGrid->U[ks][js][i].B2c - RootSoln[i].B2c);
00259 error.B3c += fabs(pGrid->U[ks][js][i].B3c - RootSoln[i].B3c);
00260 #ifndef ISOTHERMAL
00261 error.E += fabs(pGrid->U[ks][js][i].E - RootSoln[i].E );
00262 #endif
00263 }
00264
00265
00266
00267 rms_error = SQR(error.d) + SQR(error.M1) + SQR(error.M2) + SQR(error.M3);
00268 rms_error += SQR(error.B1c) + SQR(error.B2c) + SQR(error.B3c);
00269 #ifndef ISOTHERMAL
00270 rms_error += SQR(error.E);
00271 #endif
00272 rms_error = sqrt(rms_error)/(double)Nx1;
00273
00274
00275
00276 fname = "cpaw1d-errors.dat";
00277
00278 if((fp=fopen(fname,"r")) != NULL){
00279 if((fp = freopen(fname,"a",fp)) == NULL){
00280 ath_error("[Userwork_after_loop]: Unable to reopen file.\n");
00281 return;
00282 }
00283 }
00284
00285 else{
00286 if((fp = fopen(fname,"w")) == NULL){
00287 ath_error("[Userwork_after_loop]: Unable to open file.\n");
00288 return;
00289 }
00290
00291 fprintf(fp,"# Nx1 RMS-Error d M1 M2 M3");
00292 #ifndef ISOTHERMAL
00293 fprintf(fp," E");
00294 #endif
00295 fprintf(fp," B1c B2c B3c");
00296 fprintf(fp,"\n#\n");
00297 }
00298
00299 fprintf(fp,"%d %e %e %e %e %e",Nx1,rms_error,
00300 (error.d/(double)Nx1),
00301 (error.M1/(double)Nx1),
00302 (error.M2/(double)Nx1),
00303 (error.M3/(double)Nx1));
00304 #ifndef ISOTHERMAL
00305 fprintf(fp," %e",(error.E/(double)Nx1));
00306 #endif
00307 fprintf(fp," %e %e %e",
00308 (error.B1c/(double)Nx1),
00309 (error.B2c/(double)Nx1),
00310 (error.B3c/(double)Nx1));
00311 fprintf(fp,"\n");
00312
00313 fclose(fp);
00314
00315 return;
00316 }