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 ang_mom,c_infty,lambda_s,vz0;
00026 static int iprob;
00027 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00028 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00029 Real myfunc(Real x, Real v);
00030 static ConsS ***RootSoln=NULL;
00031
00032
00033
00034
00035 void problem(DomainS *pDomain)
00036 {
00037 GridS *pG = pDomain->Grid;
00038 int i,j,k;
00039 int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00040 int nx1,nx2,nx3;
00041 Real x1,x2,x3;
00042 Real xs,vs,v,pgas0,pgas,alpha,beta,a,b,converged;
00043
00044 is = pG->is; ie = pG->ie; nx1 = ie-is+1;
00045 js = pG->js; je = pG->je; nx2 = je-js+1;
00046 ks = pG->ks; ke = pG->ke; nx3 = ke-ks+1;
00047
00048 il = is-nghost*(nx1>1); iu = ie+nghost*(nx1>1); nx1 = iu-il+1;
00049 jl = js-nghost*(nx2>1); ju = je+nghost*(nx2>1); nx2 = ju-jl+1;
00050 kl = ks-nghost*(nx3>1); ku = ke+nghost*(nx3>1); nx3 = ku-kl+1;
00051
00052 #ifdef MHD
00053 ath_error("[cylwindrot]: This problem only works in hydro!\n");
00054 #endif
00055
00056 #ifndef CYLINDRICAL
00057 ath_error("[cylwindrot]: This problem only works in cylindrical!\n");
00058 #endif
00059
00060 if (nx1==1) {
00061 ath_error("[cylwindrot]: This problem can only be run in 2D or 3D!\n");
00062 }
00063 else if (nx2==1 && nx3>1) {
00064 ath_error("[cylwindrot]: Only (R,phi) can be used in 2D!\n");
00065 }
00066
00067
00068 if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS))) == NULL)
00069 ath_error("[cylwindrot]: Error allocating memory for solution\n");
00070
00071 ang_mom = par_getd("problem","ang_mom");
00072 c_infty = par_getd("problem","c_infty");
00073 vz0 = par_getd("problem","vz0");
00074 iprob = par_geti("problem","iprob");
00075 printf("gamma = %f,\t ang_mom = %f,\t c_infty = %f\n", Gamma, ang_mom, c_infty);
00076
00077 beta = 2.0*Gamma_1/(Gamma+1.0);
00078 xs = (3.0-Gamma+sqrt(SQR(Gamma-3.0)-16.0*SQR(ang_mom)))/4.0;
00079 lambda_s = 1.0/Gamma_1*pow(xs,beta)+pow(xs,beta-1.0)-0.5*SQR(ang_mom)*pow(xs,beta-2.0);
00080 lambda_s = pow(lambda_s/(0.5+1.0/Gamma_1),1.0/beta);
00081 vs = c_infty*pow(lambda_s/xs,0.5*beta);
00082 printf("xs = %13.10f,\t lambda_s = %13.10f,\t vs = %13.10f\n", xs, lambda_s, vs);
00083
00084
00085 for (i=il; i<=iu; i++) {
00086 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00087 memset(&(pG->U[ks][js][i]),0.0,sizeof(ConsS));
00088 vs = pow(lambda_s/x1,0.5*beta);
00089
00090 switch(iprob) {
00091 case 1:
00092 if (x1 < xs) {
00093 a = TINY_NUMBER; b = vs;
00094 }
00095 else {
00096 a = vs; b = HUGE_NUMBER;
00097 }
00098 break;
00099 case 2:
00100 if (x1 < xs) {
00101 a = vs; b = HUGE_NUMBER;
00102 }
00103 else {
00104 a = TINY_NUMBER; b = vs;
00105 }
00106 break;
00107 default: ath_error("[cylwindrot]: Not an accepted problem number!\n");
00108 }
00109
00110 converged = bisection(myfunc,a,b,x1,&v);
00111 if (!converged) ath_error("[cylwindrot]: Bisection did not converge!\n");
00112
00113 pG->U[ks][js][i].d = lambda_s/(x1*v);
00114 pG->U[ks][js][i].M1 = lambda_s/x1;
00115 if (iprob==2)
00116 pG->U[ks][js][i].M1 *= -1.0;
00117 pG->U[ks][js][i].M2 = pG->U[ks][js][i].d*ang_mom/x1;
00118 pG->U[ks][js][i].M3 = pG->U[ks][js][i].d*vz0;
00119
00120
00121 #ifndef ISOTHERMAL
00122 pgas0 = 1.0/Gamma;
00123 pgas = pgas0*pow(pG->U[ks][js][i].d,Gamma);
00124 pG->U[ks][js][i].E = pgas/Gamma_1
00125 + 0.5*(SQR(pG->U[ks][js][i].M1) + SQR(pG->U[ks][js][i].M2) + SQR(pG->U[ks][js][i].M3))/pG->U[ks][js][i].d;
00126 #endif
00127 }
00128
00129
00130 for (k=kl; k<=ku; k++) {
00131 for (j=jl; j<=ju; j++) {
00132 for (i=il; i<=iu; i++) {
00133 pG->U[k][j][i] = pG->U[ks][js][i];
00134 RootSoln[k][j][i] = pG->U[ks][js][i];
00135 }
00136 }
00137 }
00138
00139 StaticGravPot = grav_pot;
00140 x1GravAcc = grav_acc;
00141 bvals_mhd_fun(pDomain,left_x1,do_nothing_bc);
00142 bvals_mhd_fun(pDomain,right_x1,do_nothing_bc);
00143
00144 return;
00145 }
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 void problem_write_restart(MeshS *pM, FILE *fp)
00158 {
00159 return;
00160 }
00161
00162 void problem_read_restart(MeshS *pM, FILE *fp)
00163 {
00164 return;
00165 }
00166
00167 ConsFun_t get_usr_expr(const char *expr)
00168 {
00169 return NULL;
00170 }
00171
00172 VOutFun_t get_usr_out_fun(const char *name){
00173 return NULL;
00174 }
00175
00176 void Userwork_in_loop(MeshS *pM)
00177 {
00178 }
00179
00180 void Userwork_after_loop(MeshS *pM)
00181 {
00182 compute_l1_error("CylWindRot", pM, RootSoln, 1);
00183 }
00184
00185
00186
00187
00188
00189
00190 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00191 return -SQR(c_infty)/x1;
00192 }
00193
00194
00195
00196 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00197 return SQR(c_infty/x1);
00198 }
00199
00200
00201
00202
00203
00204
00205 Real myfunc(Real x, Real v)
00206 {
00207 return Gamma_1*(1/x + 1/Gamma_1 - 0.5*(SQR(v/c_infty)+SQR(ang_mom/x)))*pow(v*x/c_infty,Gamma_1) - pow(lambda_s,Gamma_1);
00208 }