Go to the documentation of this file.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 #include <math.h>
00025 #include <stdio.h>
00026 #include <stdlib.h>
00027 #include "defs.h"
00028 #include "athena.h"
00029 #include "globals.h"
00030 #include "prototypes.h"
00031
00032
00033 static Real Mach,dl,pl,ul,dr,pr,ur;
00034
00035
00036
00037
00038
00039
00040 static void initialize_states(void);
00041
00042
00043
00044
00045
00046 void problem(DomainS *pDomain)
00047 {
00048 GridS *pGrid = pDomain->Grid;
00049 int i,j,k,is,ie,js,je,ks,ke,shk_flag;
00050 Real amp,x1,x2,x3,xshock;
00051 div_t index;
00052
00053 is = pGrid->is;
00054 ie = pGrid->ie;
00055
00056 js = pGrid->js;
00057 je = pGrid->je;
00058
00059 ks = pGrid->ks;
00060 ke = pGrid->ke;
00061
00062
00063 Mach = par_getd("problem","Mach");
00064 amp = par_getd("problem","amp");
00065 shk_flag = par_getd("problem","shk_flag");
00066
00067
00068 dr = 1.0;
00069 #ifdef ADIABATIC
00070 pr = 1.0/Gamma;
00071 ur = Mach*sqrt(Gamma*pr/dr);
00072 #else
00073 ur = Mach*Iso_csound;
00074 #endif
00075
00076
00077 if (shk_flag == 0) {
00078 initialize_states();
00079 xshock = 0.5*(pDomain->RootMaxX[0] + pDomain->RootMinX[0]);
00080
00081
00082 } else {
00083 dl = dr;
00084 pl = pr;
00085 ul = ur;
00086 xshock = 0.9*pDomain->RootMaxX[0] + 0.1*pDomain->RootMinX[0];
00087 }
00088
00089
00090
00091
00092 for (k=ks; k<=ke; k++) {
00093 for (j=js-nghost; j<=je+nghost; j++) {
00094 for (i=is-nghost; i<=ie+nghost; i++) {
00095
00096 cc_pos(pGrid, i, j, k, &x1, &x2, &x3);
00097 if (x1 < xshock) {
00098 pGrid->U[k][j][i].d = dr;
00099 pGrid->U[k][j][i].M1 = dr*ur;
00100 pGrid->U[k][j][i].M2 = 0.0;
00101 pGrid->U[k][j][i].M3 = 0.0;
00102 #ifdef ADIABATIC
00103 pGrid->U[k][j][i].E = pr/Gamma_1 + 0.5*dr*ur*ur;
00104 #endif
00105 } else {
00106
00107 pGrid->U[k][j][i].d = dl;
00108 pGrid->U[k][j][i].M1 = dl*ul;
00109 pGrid->U[k][j][i].M2 = 0.0;
00110 pGrid->U[k][j][i].M3 = 0.0;
00111 #ifdef ADIABATIC
00112 pGrid->U[k][j][i].E = pl/Gamma_1 + 0.5*dl*ul*ul;
00113 #endif
00114 }
00115 }
00116
00117
00118
00119
00120
00121 for (i=is-nghost; i<=ie+nghost; i++) {
00122 cc_pos(pGrid, i, j, k, &x1, &x2, &x3);
00123 if (x1 < xshock) {
00124 index = div((i+j+k),2);
00125 if (index.rem == 0) {
00126 pGrid->U[k][j][i].M2 = amp;
00127 } else {
00128 pGrid->U[k][j][i].M2 = -amp;
00129 }
00130 }
00131 }
00132 }}
00133
00134 return;
00135 }
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 void problem_write_restart(MeshS *pM, FILE *fp)
00149 {
00150 return;
00151 }
00152
00153 void problem_read_restart(MeshS *pM, FILE *fp)
00154 {
00155 return;
00156 }
00157
00158 ConsFun_t get_usr_expr(const char *expr)
00159 {
00160 return NULL;
00161 }
00162
00163 VOutFun_t get_usr_out_fun(const char *name){
00164 return NULL;
00165 }
00166
00167 void Userwork_in_loop(MeshS *pM)
00168 {
00169 }
00170
00171 void Userwork_after_loop(MeshS *pM)
00172 {
00173 }
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 static void initialize_states(void)
00184 {
00185 Real jump1, jump2, jump3;
00186
00187 #ifdef ADIABATIC
00188 jump1 = (Gamma + 1.0)/(Gamma_1 + 2.0/(Mach*Mach));
00189 jump2 = (2.0*Gamma*Mach*Mach - Gamma_1)/(Gamma + 1.0);
00190 jump3 = 2.0*(1.0 - 1.0/(Mach*Mach))/(Gamma + 1.0);
00191
00192 dl = dr*jump1;
00193 pl = pr*jump2;
00194 ul = ur - jump3*Mach*sqrt(Gamma*pr/dr);
00195
00196
00197 ur = Mach*sqrt(Gamma*pr/dr);
00198 ul = ur/jump1;
00199 #endif
00200
00201 return;
00202 }