00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include <math.h>
00020 #include <stdio.h>
00021 #include <stdlib.h>
00022 #include <string.h>
00023 #include "defs.h"
00024 #include "athena.h"
00025 #include "globals.h"
00026 #include "prototypes.h"
00027
00028
00029
00030
00031
00032
00033
00034
00035 static Real omega0,rho0;
00036 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00037 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00038 Real M2(const Real x1, const Real x2, const Real x3);
00039
00040
00041
00042
00043
00044 void problem(DomainS *pDomain)
00045 {
00046 GridS *pG = pDomain->Grid;
00047 int i,j,k;
00048 int is,ie,js,je,ks,ke,nx1,nx2,nx3;
00049 int il,iu,jl,ju,kl,ku;
00050 Real r0,phi0,x0,y0,z0,angle,radius,prat,pamb,b0;
00051 Real x1,x2,x3,x2i;
00052 Real x,y,z,Eint,Emag,Ekin;
00053
00054 is = pG->is; ie = pG->ie; nx1 = ie-is+1;
00055 js = pG->js; je = pG->je; nx2 = je-js+1;
00056 ks = pG->ks; ke = pG->ke; nx3 = ke-ks+1;
00057
00058 il = is-nghost*(nx1>1); iu = ie+nghost*(nx1>1); nx1 = iu-il+1;
00059 jl = js-nghost*(nx2>1); ju = je+nghost*(nx2>1); nx2 = ju-jl+1;
00060 kl = ks-nghost*(nx3>1); ku = ke+nghost*(nx3>1); nx3 = ku-kl+1;
00061
00062 #ifndef CYLINDRICAL
00063 ath_error("[cylblast]: This problem only works in cylindrical!\n");
00064 #endif
00065
00066 if (nx1==1) {
00067 ath_error("[cylblast]: This problem can only be run in 2D or 3D!\n");
00068 }
00069 else if (nx2==1 && nx3>1) {
00070 ath_error("[cylblast]: Only (R,phi) can be used in 2D!\n");
00071 }
00072
00073
00074 radius = par_getd("problem","radius");
00075 pamb = par_getd("problem","pamb");
00076 prat = par_getd("problem","prat");
00077 rho0 = par_getd("problem","rho0");
00078 omega0 = par_getd("problem","omega0");
00079 b0 = par_getd("problem","b0");
00080
00081
00082 r0 = par_getd("problem","r0");
00083 phi0 = par_getd("problem","phi0");
00084 z0 = par_getd("problem","z0");
00085
00086
00087 angle = (PI/180.0)*par_getd("problem","angle");
00088
00089 x0 = r0*cos(phi0);
00090 y0 = r0*sin(phi0);
00091
00092
00093 for (k=kl; k<=ku; k++) {
00094 for (j=jl; j<=ju; j++) {
00095 for (i=il; i<=iu; i++) {
00096 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00097 x2i = x2 - 0.5*pG->dx2;
00098
00099 pG->U[k][j][i].d = rho0;
00100 pG->U[k][j][i].M1 = 0.0;
00101 pG->U[k][j][i].M2 = pG->U[k][j][i].d*x1*omega0;
00102 pG->U[k][j][i].M3 = 0.0;
00103 #ifdef MHD
00104
00105 pG->B1i[k][j][i] = b0*(cos(angle)*cos(x2)+sin(angle)*sin(x2));
00106 pG->B2i[k][j][i] = b0*(-cos(angle)*sin(x2i)+sin(angle)*cos(x2i));
00107 pG->B3i[k][j][i] = 0.0;
00108 pG->U[k][j][i].B1c = b0*(cos(angle)*cos(x2)+sin(angle)*sin(x2));
00109 pG->U[k][j][i].B2c = b0*(-cos(angle)*sin(x2)+sin(angle)*cos(x2));
00110 pG->U[k][j][i].B3c = 0.0;
00111 #endif
00112
00113 x = x1*cos(x2);
00114 y = x1*sin(x2);
00115 z = x3;
00116
00117
00118 #ifndef ISOTHERMAL
00119 Eint = pamb/Gamma_1;
00120 if (SQR(x-x0) + SQR(y-y0) + SQR(z-z0) < SQR(radius)) {
00121 Eint *= prat;
00122 }
00123 Emag = 0.0;
00124 #ifdef MHD
00125 Emag = 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c));
00126 #endif
00127 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;
00128 pG->U[k][j][i].E = Eint + Emag + Ekin;
00129 #endif
00130 }
00131 }
00132 }
00133
00134
00135 StaticGravPot = grav_pot;
00136 x1GravAcc = grav_acc;
00137 bvals_mhd_fun(pDomain,left_x1, do_nothing_bc);
00138 bvals_mhd_fun(pDomain,right_x1,do_nothing_bc);
00139 bvals_mhd_fun(pDomain,left_x2, do_nothing_bc);
00140 bvals_mhd_fun(pDomain,right_x2,do_nothing_bc);
00141
00142 return;
00143 }
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 void problem_write_restart(MeshS *pM, FILE *fp)
00156 {
00157 return;
00158 }
00159
00160 void problem_read_restart(MeshS *pM, FILE *fp)
00161 {
00162 return;
00163 }
00164
00165 ConsFun_t get_usr_expr(const char *expr)
00166 {
00167 return NULL;
00168 }
00169
00170 VOutFun_t get_usr_out_fun(const char *name){
00171 return NULL;
00172 }
00173
00174 void Userwork_in_loop(MeshS *pM)
00175 {
00176 }
00177
00178 void Userwork_after_loop(MeshS *pM)
00179 {
00180 }
00181
00182
00183
00184
00185
00186 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00187 return 0.5*SQR(x1*omega0);
00188 }
00189
00190
00191
00192 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00193 return x1*SQR(omega0);
00194 }
00195
00196
00197
00198 Real M2(const Real x1, const Real x2, const Real x3) {
00199 return rho0*omega0*x1;
00200 }