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