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
00023
00024
00025
00026
00027
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
00044
00045
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
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;
00065 Real x1size,x2size,x1,x2,x3,cs,sn;
00066 Real v_par, v_perp, den, pres;
00067 Real lambda;
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
00088
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
00095
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
00111
00112
00113
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
00124 if (cos_a >= sin_a) {
00125 lambda = x1size*cos_a;
00126 } else {
00127 lambda = x2size*sin_a;
00128 }
00129 }
00130
00131
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);
00143 if (dir == 1)
00144 fac = 1.0;
00145 else
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)
00161 v_perp = v_perp * kva / omega_r;
00162 else
00163 v_perp = v_perp * kva / omega_l;
00164 }
00165 #endif
00166
00167
00168
00169 pres = par_getd("problem","pres");
00170 v_par = par_getd("problem","v_par");
00171
00172
00173
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
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
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
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
00266
00267
00268
00269
00270
00271
00272
00273
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
00313
00314
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
00335
00336
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
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
00361 }
00362 }
00363
00364
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
00371 rms_error = sqrt(rms_error)/(double)(Nx1*Nx2);
00372
00373
00374
00375 fname = "cpaw2d-errors.dat";
00376
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
00384 else{
00385 if((fp = fopen(fname,"w")) == NULL){
00386 ath_error("[Userwork_after_loop]: Unable to open file.\n");
00387 return;
00388 }
00389
00390 fprintf(fp,"# Nx1 Nx2 RMS-Error d M1 M2 M3");
00391 #ifndef ISOTHERMAL
00392 fprintf(fp," E");
00393 #endif
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
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
00419
00420
00421
00422
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 }