Go to the documentation of this file.00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <math.h>
00013 #include <stdio.h>
00014 #include <stdlib.h>
00015 #include "defs.h"
00016 #include "athena.h"
00017 #include "globals.h"
00018 #include "prototypes.h"
00019
00020
00021
00022
00023 void problem(DomainS *pDomain)
00024 {
00025 GridS *pGrid=(pDomain->Grid);
00026 Prim1DS W;
00027 Cons1DS U1d;
00028 int i, is = pGrid->is, ie = pGrid->ie;
00029 int j, js = pGrid->js, je = pGrid->je;
00030 int k, ks = pGrid->ks, ke = pGrid->ke;
00031 Real pressure,drat,prat,rad,pa,da,x1,x2,x3;
00032 Real b0=0.0,Bx=0.0,rin;
00033 double theta;
00034
00035 rin = par_getd("problem","radius");
00036 pa = par_getd("problem","pamb");
00037 da = par_getd_def("problem","damb",1.0);
00038 drat = par_getd_def("problem","drat",1.0);
00039 prat = par_getd("problem","prat");
00040 #ifdef MHD
00041 b0 = par_getd("problem","b0");
00042 theta = (PI/180.0)*par_getd("problem","angle");
00043 #endif
00044
00045
00046
00047 W.d = da;
00048 W.Vx = 0.0;
00049 W.Vy = 0.0;
00050 W.Vz = 0.0;
00051 #ifdef MHD
00052 Bx = b0*cos(theta);
00053 W.By = b0*sin(theta);
00054 W.Bz = 0.0;
00055 #endif
00056
00057 for (k=ks; k<=ke; k++) {
00058 for (j=js; j<=je; j++) {
00059 for (i=is; i<=ie; i++) {
00060 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00061 rad = sqrt(x1*x1 + x2*x2 + x3*x3);
00062 #ifndef ISOTHERMAL
00063 W.P = pa;
00064 if (rad < rin) W.P = prat*pa;
00065 #endif
00066 W.d = da;
00067 if (rad < rin) W.d = drat*da;
00068
00069 U1d = Prim1D_to_Cons1D(&(W),&Bx);
00070
00071 pGrid->U[k][j][i].d = U1d.d;
00072 pGrid->U[k][j][i].M1 = U1d.Mx;
00073 pGrid->U[k][j][i].M2 = U1d.My;
00074 pGrid->U[k][j][i].M3 = U1d.Mz;
00075 #ifndef ISOTHERMAL
00076 pGrid->U[k][j][i].E = U1d.E;
00077 #endif
00078 #ifdef MHD
00079 pGrid->B1i[k][j][i] = Bx;
00080 pGrid->B2i[k][j][i] = U1d.By;
00081 pGrid->B3i[k][j][i] = U1d.Bz;
00082 pGrid->U[k][j][i].B1c = Bx;
00083 pGrid->U[k][j][i].B2c = U1d.By;
00084 pGrid->U[k][j][i].B3c = U1d.Bz;
00085 if (i == ie && ie > is) pGrid->B1i[k][j][i+1] = Bx;
00086 if (j == je && je > js) pGrid->B2i[k][j+1][i] = U1d.By;
00087 if (k == ke && ke > ks) pGrid->B3i[k+1][j][i] = U1d.Bz;
00088 #endif
00089 }
00090 }
00091 }
00092 }
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 void problem_write_restart(MeshS *pM, FILE *fp)
00106 {
00107 return;
00108 }
00109
00110 void problem_read_restart(MeshS *pM, FILE *fp)
00111 {
00112 return;
00113 }
00114
00115 ConsFun_t get_usr_expr(const char *expr)
00116 {
00117 return NULL;
00118 }
00119
00120 VOutFun_t get_usr_out_fun(const char *name){
00121 return NULL;
00122 }
00123
00124 void Userwork_in_loop(MeshS *pM)
00125 {
00126 }
00127
00128 void Userwork_after_loop(MeshS *pM)
00129 {
00130 }