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 #ifdef ISOTHERMAL
00021 #error : The twoibw test only works for adiabatic EOS.
00022 #endif
00023 #ifndef HYDRO
00024 #error : The twoibw test only works for hydro.
00025 #endif
00026
00027
00028
00029
00030 void problem(DomainS *pDomain)
00031 {
00032 GridS *pGrid=(pDomain->Grid);
00033 int i, is = pGrid->is, ie = pGrid->ie;
00034 int j, js = pGrid->js, je = pGrid->je;
00035 int k, ks = pGrid->ks, ke = pGrid->ke;
00036 int shk_dir;
00037 Real x1,x2,x3;
00038
00039 shk_dir = par_geti("problem","shk_dir");
00040 if (shk_dir < 1 || shk_dir > 3) {
00041 ath_error("[problem]: shk_dir = %d must be either 1, 2 or 3\n",shk_dir);
00042 }
00043
00044 if(shk_dir == 1) {
00045
00046 for (k=ks; k<=ke; k++) {
00047 for (j=js; j<=je; j++) {
00048 for (i=is; i<=ie; i++) {
00049
00050 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00051
00052 pGrid->U[k][j][i].d = 1.0;
00053 pGrid->U[k][j][i].M1 = 0.0;
00054 pGrid->U[k][j][i].M2 = 0.0;
00055 pGrid->U[k][j][i].M3 = 0.0;
00056 if (x1 < 0.1) {
00057 pGrid->U[k][j][i].E = 1.0e3/Gamma_1 ;
00058 }
00059 else if (x1 > 0.9) {
00060 pGrid->U[k][j][i].E = 1.0e2/Gamma_1 ;
00061 }
00062 else {
00063 pGrid->U[k][j][i].E = 0.01/Gamma_1 ;
00064 }
00065 }
00066 }
00067 }
00068 }
00069 else if (shk_dir == 2) {
00070
00071 for (k=ks; k<=ke; k++) {
00072 for (j=js; j<=je; j++) {
00073 for (i=is; i<=ie; i++) {
00074
00075 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00076
00077 pGrid->U[k][j][i].d = 1.0;
00078 pGrid->U[k][j][i].M1 = 0.0;
00079 pGrid->U[k][j][i].M2 = 0.0;
00080 pGrid->U[k][j][i].M3 = 0.0;
00081 if (x2 < 0.1) {
00082 pGrid->U[k][j][i].E = 1.0e3/Gamma_1 ;
00083 }
00084 else if (x2 > 0.9) {
00085 pGrid->U[k][j][i].E = 1.0e2/Gamma_1 ;
00086 }
00087 else {
00088 pGrid->U[k][j][i].E = 0.01/Gamma_1 ;
00089 }
00090 }
00091 }
00092 }
00093 }
00094 else {
00095
00096 for (k=ks; k<=ke; k++) {
00097 for (j=js; j<=je; j++) {
00098 for (i=is; i<=ie; i++) {
00099
00100 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00101
00102 pGrid->U[k][j][i].d = 1.0;
00103 pGrid->U[k][j][i].M1 = 0.0;
00104 pGrid->U[k][j][i].M2 = 0.0;
00105 pGrid->U[k][j][i].M3 = 0.0;
00106 if (x3 < 0.1) {
00107 pGrid->U[k][j][i].E = 1.0e3/Gamma_1 ;
00108 }
00109 else if (x3 > 0.9) {
00110 pGrid->U[k][j][i].E = 1.0e2/Gamma_1 ;
00111 }
00112 else {
00113 pGrid->U[k][j][i].E = 0.01/Gamma_1 ;
00114 }
00115 }
00116 }
00117 }
00118 }
00119 return;
00120 }
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 void problem_write_restart(MeshS *pM, FILE *fp)
00134 {
00135 return;
00136 }
00137
00138 void problem_read_restart(MeshS *pM, FILE *fp)
00139 {
00140 return;
00141 }
00142
00143 ConsFun_t get_usr_expr(const char *expr)
00144 {
00145 return NULL;
00146 }
00147
00148 VOutFun_t get_usr_out_fun(const char *name){
00149 return NULL;
00150 }
00151
00152 void Userwork_in_loop(MeshS *pM)
00153 {
00154 }
00155
00156 void Userwork_after_loop(MeshS *pM)
00157 {
00158 }