00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <math.h>
00030 #include <stdio.h>
00031 #include "defs.h"
00032 #include "athena.h"
00033 #include "globals.h"
00034 #include "prototypes.h"
00035
00036
00037
00038 static Real dl,pl,ul;
00039 #ifdef MHD
00040 static Real bxl,byl,bzl;
00041 #endif
00042
00043
00044
00045
00046
00047
00048 void shk_cloud_iib(GridS *pGrid);
00049
00050
00051
00052
00053
00054 void problem(DomainS *pDomain)
00055 {
00056 GridS *pGrid = pDomain->Grid;
00057 int i=0,j=0,k=0;
00058 int is,ie,js,je,ks,ke,iprob;
00059 Real jump1, jump2, jump3;
00060 Real x1,x2,x3,diag,rad,xshock;
00061 Real Mach,drat;
00062 #ifdef MHD
00063 Real beta,bxr,byr,bzr;
00064 #endif
00065 Real dr,pr,ur;
00066
00067
00068
00069 xshock = -2.0;
00070 rad = 1.0;
00071 Mach = par_getd("problem","Mach");
00072 drat = par_getd("problem","drat");
00073 iprob = par_geti("problem","iprob");
00074 iprob = 1;
00075 #ifdef MHD
00076 beta = par_getd("problem","beta");
00077 #endif
00078
00079
00080
00081 dr = 1.0;
00082 pr = 1.0/Gamma;
00083 ur = 0.0;
00084
00085
00086
00087 jump1 = (Gamma + 1.0)/(Gamma_1 + 2.0/(Mach*Mach));
00088 jump2 = (2.0*Gamma*Mach*Mach - Gamma_1)/(Gamma + 1.0);
00089 jump3 = 2.0*(1.0 - 1.0/(Mach*Mach))/(Gamma + 1.0);
00090
00091 dl = dr*jump1;
00092 pl = pr*jump2;
00093 #ifdef ISOTHERMAL
00094 ul = ur + jump3*Mach*Iso_csound;
00095 #else
00096 ul = ur + jump3*Mach*sqrt(Gamma*pr/dr);
00097 #endif
00098
00099
00100
00101 is = pGrid->is; ie = pGrid->ie;
00102 js = pGrid->js; je = pGrid->je;
00103 ks = pGrid->ks; ke = pGrid->ke;
00104
00105 for (k=ks; k<=ke; k++) {
00106 for (j=js; j<=je; j++) {
00107 for (i=is; i<=ie; i++) {
00108 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00109 diag = sqrt(x1*x1 + x2*x2 + x3*x3);
00110
00111
00112 if(x1 < xshock) {
00113 pGrid->U[k][j][i].d = dl;
00114 pGrid->U[k][j][i].M1 = ul*dl;
00115 pGrid->U[k][j][i].M2 = 0.0;
00116 pGrid->U[k][j][i].M3 = 0.0;
00117 #ifdef MHD
00118 pGrid->B1i[k][j][i] = bxl;
00119 pGrid->B2i[k][j][i] = byl;
00120 pGrid->B3i[k][j][i] = bzl;
00121 pGrid->U[k][j][i].B1c = bxl;
00122 pGrid->U[k][j][i].B2c = byl;
00123 pGrid->U[k][j][i].B3c = bzl;
00124 #endif
00125 #ifdef ADIABATIC
00126 pGrid->U[k][j][i].E = pl/Gamma_1
00127 #ifdef MHD
00128 + 0.5*(bxl*bxl + byl*byl + bzl*bzl)
00129 #endif
00130 + 0.5*dl*(ul*ul);
00131 #endif
00132 #if (NSCALARS > 0)
00133 pGrid->U[k][j][i].s[0] = 0.0;
00134 #endif
00135
00136
00137 } else {
00138 pGrid->U[k][j][i].d = dr;
00139 pGrid->U[k][j][i].M1 = ur*dr;
00140 pGrid->U[k][j][i].M2 = 0.0;
00141 pGrid->U[k][j][i].M3 = 0.0;
00142 #ifdef MHD
00143 pGrid->B1i[k][j][i] = bxr;
00144 pGrid->B2i[k][j][i] = byr;
00145 pGrid->B3i[k][j][i] = bzr;
00146 pGrid->U[k][j][i].B1c = bxr;
00147 pGrid->U[k][j][i].B2c = byr;
00148 pGrid->U[k][j][i].B3c = bzr;
00149 #endif
00150 #ifdef ADIABATIC
00151 pGrid->U[k][j][i].E = pr/Gamma_1
00152 #ifdef MHD
00153 + 0.5*(bxr*bxr + byr*byr + bzr*bzr)
00154 #endif
00155 + 0.5*dr*(ur*ur);
00156 #endif
00157 #if (NSCALARS > 0)
00158 pGrid->U[k][j][i].s[0] = 0.0;
00159 #endif
00160 }
00161
00162
00163 if (diag < rad) {
00164 pGrid->U[k][j][i].d = dr*drat;
00165 pGrid->U[k][j][i].M1 = ur*dr*drat;
00166 pGrid->U[k][j][i].M2 = 0.0;
00167 pGrid->U[k][j][i].M3 = 0.0;
00168 #ifdef MHD
00169 pGrid->B1i[k][j][i] = bxr;
00170 pGrid->B2i[k][j][i] = byr;
00171 pGrid->B3i[k][j][i] = bzr;
00172 pGrid->U[k][j][i].B1c = bxr;
00173 pGrid->U[k][j][i].B2c = byr;
00174 pGrid->U[k][j][i].B3c = bzr;
00175 #endif
00176 #ifdef ADIABATIC
00177 pGrid->U[k][j][i].E = pr/Gamma_1
00178 #ifdef MHD
00179 + 0.5*(bxr*bxr + byr*byr + bzr*bzr)
00180 #endif
00181 + 0.5*dr*drat*(ur*ur);
00182 #endif
00183 #if (NSCALARS > 0)
00184 pGrid->U[k][j][i].s[0] = drat;
00185 #endif
00186 }
00187 }
00188 }
00189 }
00190
00191
00192
00193 #ifdef MHD
00194 i = ie+1;
00195 for (k=ks; k<=ke; k++) {
00196 for (j=js; j<=je; j++) {
00197 pGrid->B1i[k][j][i] = bxr;
00198 }
00199 }
00200 j = je+1;
00201 for (k=ks; k<=ke; k++) {
00202 for (i=is; i<=ie; i++) {
00203 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00204 if (x1 < xshock) {pGrid->B2i[k][j][i] = byl;}
00205 else {pGrid->B2i[k][j][i] = byr;}
00206 }
00207 }
00208 if (ke > ks) {
00209 k = ke+1;
00210 for (j=js; j<=je; j++) {
00211 for (i=is; i<=ie; i++) {
00212 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00213 if (x1 < xshock) {pGrid->B3i[k][j][i] = bzl;}
00214 else {pGrid->B3i[k][j][i] = bzr;}
00215 }
00216 }
00217 }
00218 #endif
00219
00220
00221
00222 if (pDomain->Disp[0] == 0) bvals_mhd_fun(pDomain,left_x1,shk_cloud_iib);
00223
00224 return;
00225 }
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239 void problem_write_restart(MeshS *pM, FILE *fp)
00240 {
00241 return;
00242 }
00243
00244 void problem_read_restart(MeshS *pM, FILE *fp)
00245 {
00246 return;
00247 }
00248
00249 #if (NSCALARS > 0)
00250 static Real color(const GridS *pG, const int i, const int j, const int k)
00251 {
00252 return pG->U[k][j][i].s[0]/pG->U[k][j][i].d;
00253 }
00254 #endif
00255
00256 ConsFun_t get_usr_expr(const char *expr)
00257 {
00258 #if (NSCALARS > 0)
00259 if(strcmp(expr,"color")==0) return color;
00260 #endif
00261
00262 return NULL;
00263 }
00264
00265 VOutFun_t get_usr_out_fun(const char *name){
00266 return NULL;
00267 }
00268
00269 void Userwork_in_loop(MeshS *pM)
00270 {
00271 return;
00272 }
00273
00274 void Userwork_after_loop(MeshS *pM)
00275 {
00276 return;
00277 }
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288 void shk_cloud_iib(GridS *pGrid)
00289 {
00290 int i=0,j=0,k=0;
00291 int js,je,ks,ke;
00292
00293 js = pGrid->js; je = pGrid->je;
00294 ks = pGrid->ks; ke = pGrid->ke;
00295
00296 for (k=ks; k<=ke; k++) {
00297 for (j=js; j<=je; j++) {
00298 for (i=1; i<=nghost; i++) {
00299 pGrid->U[k][j][i].d = dl;
00300 pGrid->U[k][j][i].M1 = ul*dl;
00301 pGrid->U[k][j][i].M2 = 0.0;
00302 pGrid->U[k][j][i].M3 = 0.0;
00303 #ifdef MHD
00304 pGrid->B1i[k][j][i] = bxl;
00305 pGrid->B2i[k][j][i] = byl;
00306 pGrid->B3i[k][j][i] = bzl;
00307 pGrid->U[k][j][i].B1c = bxl;
00308 pGrid->U[k][j][i].B2c = byl;
00309 pGrid->U[k][j][i].B3c = bzl;
00310 #endif
00311 #ifdef ADIABATIC
00312 pGrid->U[k][j][i].E = pl/Gamma_1
00313 #ifdef MHD
00314 + 0.5*(bxl*bxl + byl*byl + bzl*bzl)
00315 #endif
00316 + 0.5*dl*(ul*ul);
00317 #endif
00318 #if (NSCALARS > 0)
00319 pGrid->U[k][j][i].s[0] = 0.0;
00320 #endif
00321 }
00322 }
00323 }
00324 }