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