00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009 #include <math.h>
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012 #include <string.h>
00013 #include "defs.h"
00014 #include "athena.h"
00015 #include "globals.h"
00016 #include "prototypes.h"
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 static Real br0,omega0,vz0,rho0,pgas0,a,x2min,x2max,phi0,x1save,x3save;
00032 Real d(const Real x1, const Real x2, const Real x3);
00033 Real M2(const Real x1, const Real x2, const Real x3);
00034 Real B1(const Real x1, const Real x2, const Real x3);
00035 Real B1i(const Real x2);
00036 Real Pgas(const Real x1, const Real x2, const Real x3);
00037 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00038 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00039 void cylbr_ix1(GridS *pG);
00040 void cylbr_ox1(GridS *pG);
00041 static ConsS ***RootSoln=NULL;
00042
00043
00044
00045
00046 void problem(DomainS *pDomain)
00047 {
00048 GridS *pG = pDomain->Grid;
00049 int i,j,k;
00050 int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00051 int nx1,nx2,nx3;
00052 Real x1,x2,x3,x1a,x2a,x2b;
00053 Real Eint,Ekin,Emag;
00054
00055 is = pG->is; ie = pG->ie; nx1 = ie-is+1;
00056 js = pG->js; je = pG->je; nx2 = je-js+1;
00057 ks = pG->ks; ke = pG->ke; nx3 = ke-ks+1;
00058
00059 il = is-nghost*(nx1>1); iu = ie+nghost*(nx1>1); nx1 = iu-il+1;
00060 jl = js-nghost*(nx2>1); ju = je+nghost*(nx2>1); nx2 = ju-jl+1;
00061 kl = ks-nghost*(nx3>1); ku = ke+nghost*(nx3>1); nx3 = ku-kl+1;
00062
00063 #ifndef CYLINDRICAL
00064 ath_error("[cylbr]: This problem only works in cylindrical!\n");
00065 #endif
00066
00067 if (nx1==1) {
00068 ath_error("[cylbr]: This problem can only be run in 2D or 3D!\n");
00069 }
00070 else if (nx2==1 && nx3>1) {
00071 ath_error("[cylbr]: Only (R,phi) can be used in 2D!\n");
00072 }
00073
00074 x2min = par_getd("grid","x2min");
00075 x2max = par_getd("grid","x2max");
00076 a = 2.0*PI/(x2max-x2min);
00077 phi0 = a*(omega0*pG->time + x2min);
00078
00079 omega0 = par_getd("problem", "omega0");
00080 vz0 = par_getd("problem", "vz0");
00081 br0 = par_getd("problem", "br0");
00082 rho0 = par_getd("problem", "rho0");
00083 pgas0 = par_getd("problem", "pgas0");
00084
00085
00086 if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS))) == NULL)
00087 ath_error("[cylbr]: Error allocating memory for solution\n");
00088
00089 for (k=kl; k<=ku; k++) {
00090 for (j=jl; j<=ju; j++) {
00091 for (i=il; i<=iu; i++) {
00092 x1a = x1 - 0.5*pG->dx1;
00093 x2a = x2 - 0.5*pG->dx2;
00094 x2b = x2 + 0.5*pG->dx2;
00095 x1save = x1a;
00096 x3save = x3;
00097
00098 pG->U[k][j][i].d = avg2d(d,pG,i,j,k);
00099 pG->U[k][j][i].M2 = avg2d(M2,pG,i,j,k);
00100 pG->U[k][j][i].M3 = pG->U[k][j][i].d*vz0;
00101 pG->U[k][j][i].B1c = avg2d(B1,pG,i,j,k);
00102 pG->B1i[k][j][i] = qsimp(B1i,x2a,x2b)/pG->dx2;
00103
00104 Eint = avg2d(Pgas,pG,i,j,k)/Gamma_1;
00105 Emag = 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c));
00106 Ekin = 0.5*(SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2) + SQR(pG->U[k][j][i].M3))/pG->U[k][j][i].d;
00107 pG->U[k][j][i].E = Eint + Emag + Ekin;
00108
00109
00110 RootSoln[k][j][i] = pG->U[k][j][i];
00111 }
00112 }
00113 }
00114
00115 StaticGravPot = grav_pot;
00116 x1GravAcc = grav_acc;
00117 bvals_mhd_fun(pDomain,left_x1,cylbr_ix1);
00118 bvals_mhd_fun(pDomain,right_x1,cylbr_ox1);
00119
00120 return;
00121 }
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 void problem_write_restart(MeshS *pM, FILE *fp)
00134 {
00135 return;
00136 }
00137
00138 void problem_read_restart(MeshS *pM, FILE *fp)
00139 {
00140 return;
00141 }
00142
00143 ConsFun_t get_usr_expr(const char *expr)
00144 {
00145 return NULL;
00146 }
00147
00148 VOutFun_t get_usr_out_fun(const char *name){
00149 return NULL;
00150 }
00151
00152 void Userwork_in_loop(MeshS *pM)
00153 {
00154 printf("Max divB = %1.10e\n", compute_div_b(pM->Domain[0][0].Grid));
00155 }
00156
00157 void Userwork_after_loop(MeshS *pM)
00158 {
00159 compute_l1_error("CylBR", pM, RootSoln, 1);
00160 }
00161
00162
00163
00164
00165
00166 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00167 return 0.5*SQR(x1*omega0) - 0.5*SQR(br0)/SQR(x1);
00168 }
00169
00170
00171
00172 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00173 return x1*SQR(omega0) + SQR(br0)/pow(x1,3);
00174 }
00175
00176
00177
00178 Real d(const Real x1, const Real x2, const Real x3) {
00179 return rho0 + SQR(sin(a*x2-phi0));
00180 }
00181
00182
00183
00184 Real M2(const Real x1, const Real x2, const Real x3) {
00185 return (rho0 + SQR(sin(a*x2-phi0)))*omega0*x1;
00186 }
00187
00188
00189
00190 Real B1(const Real x1, const Real x2, const Real x3) {
00191 return br0*(cos(a*x2-phi0))/x1;
00192 }
00193
00194
00195
00196 Real B1i(const Real x2) {
00197 return B1(x1save,x2,x3save);
00198 }
00199
00200
00201
00202 Real Pgas(const Real x1, const Real x2, const Real x3) {
00203 return pgas0 + 0.5*SQR(br0)*(rho0 + SQR(sin(a*x2-phi0)))/SQR(x1);
00204 }
00205
00206
00207
00208
00209
00210
00211
00212 void cylbr_ix1(GridS *pG)
00213 {
00214 int is = pG->is;
00215 int js = pG->js, je = pG->je;
00216 int ks = pG->ks, ke = pG->ke;
00217 int i,j,k;
00218 Real x1,x2,x3,x1a,x2a,x2b;
00219 Real Eint,Emag,Ekin;
00220
00221 phi0 = a*(omega0*pG->time + x2min);
00222 for (k=ks; k<=ke; k++) {
00223 for (j=js; j<=je; j++) {
00224 for (i=1; i<=nghost; i++) {
00225 cc_pos(pG,is-i,j,k,&x1,&x2,&x3);
00226 x1a = x1 - 0.5*pG->dx1;
00227 x2a = x2 - 0.5*pG->dx2;
00228 x2b = x2 + 0.5*pG->dx2;
00229 x1save = x1a;
00230 x3save = x3;
00231
00232 pG->U[k][j][is-i].d = avg2d(d,pG,is-i,j,k);
00233 pG->U[k][j][is-i].M2 = avg2d(M2,pG,is-i,j,k);
00234 pG->U[k][j][is-i].M3 = pG->U[k][j][is-i].d*vz0;
00235 pG->U[k][j][is-i].B1c = avg2d(B1,pG,is-i,j,k);
00236 pG->B1i[k][j][is-i] = qsimp(B1i,x2a,x2b)/pG->dx2;
00237
00238 Eint = avg2d(Pgas,pG,is-i,j,k)/Gamma_1;
00239 Emag = 0.5*(SQR(pG->U[k][j][is-i].B1c) + SQR(pG->U[k][j][is-i].B2c) + SQR(pG->U[k][j][is-i].B3c));
00240 Ekin = 0.5*(SQR(pG->U[k][j][is-i].M1) + SQR(pG->U[k][j][is-i].M2) + SQR(pG->U[k][j][is-i].M3))/pG->U[k][j][is-i].d;
00241 pG->U[k][j][is-i].E = Eint + Emag + Ekin;
00242 }
00243 }
00244 }
00245
00246 return;
00247 }
00248
00249
00250
00251
00252
00253 void cylbr_ox1(GridS *pG)
00254 {
00255 int ie = pG->ie;
00256 int js = pG->js, je = pG->je;
00257 int ks = pG->ks, ke = pG->ke;
00258 int i,j,k;
00259 Real x1,x2,x3,x1a,x2a,x2b;
00260 Real Eint,Emag,Ekin;
00261
00262 phi0 = a*(omega0*pG->time + x2min);
00263 for (k=ks; k<=ke; k++) {
00264 for (j=js; j<=je; j++) {
00265 for (i=1; i<=nghost; i++) {
00266 cc_pos(pG,ie+i,j,k,&x1,&x2,&x3);
00267 x1a = x1 - 0.5*pG->dx1;
00268 x2a = x2 - 0.5*pG->dx2;
00269 x2b = x2 + 0.5*pG->dx2;
00270 x1save = x1a;
00271 x3save = x3;
00272
00273 pG->U[k][j][ie+i].d = avg2d(d,pG,ie+i,j,k);
00274 pG->U[k][j][ie+i].M2 = avg2d(M2,pG,ie+i,j,k);
00275 pG->U[k][j][ie+i].M3 = pG->U[k][j][ie+i].d*vz0;
00276 pG->U[k][j][ie+i].B1c = avg2d(B1,pG,ie+i,j,k);
00277 pG->B1i[k][j][ie+i] = qsimp(B1i,x2a,x2b)/pG->dx2;
00278
00279 Eint = avg2d(Pgas,pG,ie+i,j,k)/Gamma_1;
00280 Emag = 0.5*(SQR(pG->U[k][j][ie+i].B1c) + SQR(pG->U[k][j][ie+i].B2c) + SQR(pG->U[k][j][ie+i].B3c));
00281 Ekin = 0.5*(SQR(pG->U[k][j][ie+i].M1) + SQR(pG->U[k][j][ie+i].M2) + SQR(pG->U[k][j][ie+i].M3))/pG->U[k][j][ie+i].d;
00282 pG->U[k][j][ie+i].E = Eint + Emag + Ekin;
00283 }
00284 }
00285 }
00286
00287 return;
00288 }