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 static Real theta,omega,eta,E,vz;
00026 const Real rho_A = 1.0;
00027 const Real R_A = 1.0;
00028 const Real GM = 1.0;
00029 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00030 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00031 Real myfunc(const Real x, const Real y);
00032 static ConsS ***RootSoln=NULL;
00033
00034
00035
00036
00037 void problem(DomainS *pDomain)
00038 {
00039 GridS *pG = pDomain->Grid;
00040 int i,j,k,n,converged;
00041 int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00042 int nx1, nx2, nx3;
00043 Real x1, x2, x3;
00044 Real a,b,c,d,xmin,xmax,ymin,ymax;
00045 Real x,y,xslow,yslow,xfast,yfast;
00046 Real R0,R1,R2,rho,Mdot,K,Omega,Pgas,beta,vR,BR,vphi,Bphi;
00047 ConsS *Wind=NULL;
00048 Real *pU=NULL,*pUl=NULL,*pUr=NULL;
00049 Real lsf,rsf;
00050
00051 is = pG->is; ie = pG->ie; nx1 = ie-is+1;
00052 js = pG->js; je = pG->je; nx2 = je-js+1;
00053 ks = pG->ks; ke = pG->ke; nx3 = ke-ks+1;
00054
00055 il = is-nghost*(nx1>1); iu = ie+nghost*(nx1>1); nx1 = iu-il+1;
00056 jl = js-nghost*(nx2>1); ju = je+nghost*(nx2>1); nx2 = ju-jl+1;
00057 kl = ks-nghost*(nx3>1); ku = ke+nghost*(nx3>1); nx3 = ku-kl+1;
00058
00059 #ifndef CYLINDRICAL
00060 ath_error("[cylwindrotb]: This problem only works in cylindrical!\n");
00061 #endif
00062 #ifndef MHD
00063 ath_error("[cylwindrotb]: This problem only works in MHD!\n");
00064 #endif
00065
00066 if (nx1==1) {
00067 ath_error("[cylwindrotb]: Only R can be used in 1D!\n");
00068 }
00069 else if (nx2==1 && nx3>1) {
00070 ath_error("[cylwindrotb]: Only (R,phi) can be used in 2D!\n");
00071 }
00072
00073
00074 if ((Wind = (ConsS*)calloc_1d_array(nx1+1,sizeof(ConsS))) == NULL)
00075 ath_error("[cylwindrotb]: Error allocating memory\n");
00076
00077
00078 if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS))) == NULL)
00079 ath_error("[cylwindrotb]: Error allocating memory\n");
00080
00081 theta = par_getd("problem","theta");
00082 omega = par_getd("problem","omega");
00083 vz = par_getd("problem","vz");
00084
00085
00086
00087 xslow = 0.5243264128;
00088 yslow = 2.4985859152;
00089 xfast = 1.6383327831;
00090 yfast = 0.5373957134;
00091 E = 7.8744739104;
00092 eta = 2.3608500383;
00093
00094 xmin = par_getd("domain1","x1min")/R_A;
00095 xmax = par_getd("domain1","x1max")/R_A;
00096 ymin = 0.45/rho_A;
00097 ymax = 2.6/rho_A;
00098
00099 printf("theta = %f,\t omega = %f,\t eta = %f,\t E = %f\n", theta,omega,eta,E);
00100 printf("xslow = %f,\t yslow = %f,\t xfast = %f,\t yfast = %f\n", xslow,yslow,xfast,yfast);
00101 printf("xmin = %f,\t ymin = %f,\t xmax = %f,\t ymax = %f\n", xmin,ymin,xmax,ymax);
00102
00103
00104
00105 for (i=il; i<=iu+1; i++) {
00106 memset(&(Wind[i]),0.0,sizeof(ConsS));
00107 cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00108
00109
00110 R0 = x1 - 0.5*pG->dx1;
00111 x = R0/R_A;
00112
00113
00114 if (x < xslow) {
00115 sign_change(myfunc,yslow,10.0*ymax,x,&a,&b);
00116 sign_change(myfunc,b,10.0*ymax,x,&a,&b);
00117 } else if (x < 1.0) {
00118 sign_change(myfunc,1.0+TINY_NUMBER,yslow,x,&a,&b);
00119 } else if (x < xfast) {
00120 sign_change(myfunc,yfast,1.0-TINY_NUMBER,x,&a,&b);
00121 if (!sign_change(myfunc,b,1.0-TINY_NUMBER,x,&a,&b)) {
00122 a = yfast;
00123 b = 1.0-TINY_NUMBER;
00124 }
00125 } else {
00126 sign_change(myfunc,0.5*ymin,yfast,x,&a,&b);
00127 }
00128
00129
00130 converged = bisection(myfunc,a,b,x,&y);
00131 if(!converged) {
00132 ath_error("[cylwindrotb]: Bisection did not converge!\n");
00133 }
00134
00135
00136 rho = rho_A*y;
00137 Mdot = sqrt(R_A*SQR(rho_A)*GM*eta);
00138 Omega = sqrt((GM*omega)/pow(R_A,3));
00139 K = (GM*theta)/(Gamma*pow(rho_A,Gamma_1)*R_A);
00140 Pgas = K*pow(rho,Gamma);
00141 vR = Mdot/(R0*rho);
00142 beta = sqrt(1.0/rho_A);
00143 BR = beta*rho*vR;
00144 vphi = R0*Omega*(1.0/SQR(x)-y)/(1.0-y);
00145 Bphi = beta*rho*(vphi-R0*Omega);
00146
00147 Wind[i].d = rho;
00148 Wind[i].M1 = rho*vR;
00149 Wind[i].M2 = rho*vphi;
00150 Wind[i].M3 = rho*vz;
00151 Wind[i].B1c = BR;
00152 Wind[i].B2c = Bphi;
00153 Wind[i].B3c = 0.0;
00154 Wind[i].E = Pgas/Gamma_1
00155 + 0.5*(SQR(Wind[i].B1c) + SQR(Wind[i].B2c) + SQR(Wind[i].B3c))
00156 + 0.5*(SQR(Wind[i].M1 ) + SQR(Wind[i].M2 ) + SQR(Wind[i].M3 ))/Wind[i].d;
00157 }
00158
00159
00160 for (i=il; i<=iu; i++) {
00161 memset(&(pG->U[ks][js][i]),0.0,sizeof(ConsS));
00162 cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00163 lsf = (x1 - 0.5*pG->dx1)/x1;
00164 rsf = (x1 + 0.5*pG->dx1)/x1;
00165
00166 pU = (Real*)&(pG->U[ks][js][i]);
00167 pUl = (Real*)&(Wind[i]);
00168 pUr = (Real*)&(Wind[i+1]);
00169 for (n=0; n<NWAVE; n++) {
00170 pU[n] = 0.5*(lsf*pUl[n] + rsf*pUr[n]);
00171 }
00172
00173 pG->B1i[ks][js][i] = Wind[i].B1c;
00174 pG->B2i[ks][js][i] = Wind[i].B2c;
00175 pG->B3i[ks][js][i] = Wind[i].B3c;
00176 }
00177
00178
00179 for (k=kl; k<=ku; k++) {
00180 for (j=jl; j<=ju; j++) {
00181 for (i=il; i<=iu; i++) {
00182 pG->U[k][j][i] = pG->U[ks][js][i];
00183 RootSoln[k][j][i] = pG->U[ks][js][i];
00184 }
00185 }
00186 }
00187
00188 StaticGravPot = grav_pot;
00189 x1GravAcc = grav_acc;
00190 bvals_mhd_fun(pDomain,left_x1,do_nothing_bc);
00191 bvals_mhd_fun(pDomain,right_x1,do_nothing_bc);
00192
00193 free_1d_array((void *)Wind);
00194
00195 return;
00196 }
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208 void problem_write_restart(MeshS *pM, FILE *fp)
00209 {
00210 return;
00211 }
00212
00213 void problem_read_restart(MeshS *pM, FILE *fp)
00214 {
00215 return;
00216 }
00217
00218 ConsFun_t get_usr_expr(const char *expr)
00219 {
00220 return NULL;
00221 }
00222
00223 VOutFun_t get_usr_out_fun(const char *name){
00224 return NULL;
00225 }
00226
00227 void Userwork_in_loop(MeshS *pM)
00228 {
00229 }
00230
00231 void Userwork_after_loop(MeshS *pM)
00232 {
00233 compute_l1_error("CylWindRotB", pM, RootSoln, 1);
00234 }
00235
00236
00237
00238
00239
00240 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00241 return -GM/x1;
00242 }
00243
00244
00245
00246 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00247 return GM/SQR(x1);
00248 }
00249
00250
00251
00252
00253
00254
00255
00256
00257 Real myfunc(const Real x, const Real y)
00258 {
00259 return eta/(2.0*pow(x,2)*pow(y,2)) + (theta/Gamma_1)*pow(y,Gamma_1)
00260 - 1.0/x + 0.5*omega*(pow(x-1.0/x,2)/pow(y-1,2) - pow(x,2)) - E;
00261 }