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 int i,il,iu,j,jl,ju,k,kl,ku;
00027 int shk_dir;
00028 Real x1,x2,x3;
00029 Prim1DS Wl, Wr;
00030 Cons1DS U1d, Ul, Ur;
00031 Real Bxl=0.0, Bxr=0.0;
00032
00033
00034
00035 Wl.d = par_getd("problem","dl");
00036 #ifdef ADIABATIC
00037 Wl.P = par_getd("problem","pl");
00038 #endif
00039 Wl.Vx = par_getd("problem","v1l");
00040 Wl.Vy = par_getd("problem","v2l");
00041 Wl.Vz = par_getd("problem","v3l");
00042 #ifdef MHD
00043 Bxl = par_getd("problem","b1l");
00044 Wl.By = par_getd("problem","b2l");
00045 Wl.Bz = par_getd("problem","b3l");
00046 #endif
00047 #if (NSCALARS > 0)
00048 Wl.r[0] = par_getd("problem","r[0]l");
00049 #endif
00050
00051
00052
00053 Wr.d = par_getd("problem","dr");
00054 #ifdef ADIABATIC
00055 Wr.P = par_getd("problem","pr");
00056 #endif
00057 Wr.Vx = par_getd("problem","v1r");
00058 Wr.Vy = par_getd("problem","v2r");
00059 Wr.Vz = par_getd("problem","v3r");
00060 #ifdef MHD
00061 Bxr = par_getd("problem","b1r");
00062 Wr.By = par_getd("problem","b2r");
00063 Wr.Bz = par_getd("problem","b3r");
00064 if (Bxr != Bxl) ath_error(0,"[shkset1d] L/R values of Bx not the same\n");
00065 #endif
00066 #if (NSCALARS > 0)
00067 Wr.r[0] = par_getd("problem","r[0]r");
00068 #endif
00069
00070 Ul = Prim1D_to_Cons1D(&Wl, &Bxl);
00071 Ur = Prim1D_to_Cons1D(&Wr, &Bxr);
00072
00073
00074 shk_dir = par_geti("problem","shk_dir");
00075 if (shk_dir != 1 && shk_dir != 2 && shk_dir != 3) {
00076 ath_error("[problem]: shk_dir = %d must be either 1,2 or 3\n",shk_dir);
00077 }
00078
00079
00080 iu = pGrid->ie + nghost;
00081 il = pGrid->is - nghost;
00082
00083 if (pGrid->Nx[1] > 1) {
00084 ju = pGrid->je + nghost;
00085 jl = pGrid->js - nghost;
00086 }
00087 else {
00088 ju = pGrid->je;
00089 jl = pGrid->js;
00090 }
00091
00092 if (pGrid->Nx[2] > 1) {
00093 ku = pGrid->ke + nghost;
00094 kl = pGrid->ks - nghost;
00095 }
00096 else {
00097 ku = pGrid->ke;
00098 kl = pGrid->ks;
00099 }
00100
00101
00102
00103
00104 switch(shk_dir) {
00105
00106 case 1:
00107 for (k=kl; k<=ku; k++) {
00108 for (j=jl; j<=ju; j++) {
00109 for (i=il; i<=iu; i++) {
00110 cc_pos(pGrid, i, j, k, &x1, &x2, &x3);
00111
00112
00113 if (x1 <= 0.0) {
00114 U1d = Ul;
00115 } else {
00116 U1d = Ur;
00117 }
00118
00119
00120 pGrid->U[k][j][i].d = U1d.d;
00121 pGrid->U[k][j][i].M1 = U1d.Mx;
00122 pGrid->U[k][j][i].M2 = U1d.My;
00123 pGrid->U[k][j][i].M3 = U1d.Mz;
00124 #ifdef MHD
00125 pGrid->B1i[k][j][i] = Bxl;
00126 pGrid->B2i[k][j][i] = U1d.By;
00127 pGrid->B3i[k][j][i] = U1d.Bz;
00128 pGrid->U[k][j][i].B1c = Bxl;
00129 pGrid->U[k][j][i].B2c = U1d.By;
00130 pGrid->U[k][j][i].B3c = U1d.Bz;
00131 #endif
00132 #ifdef ADIABATIC
00133 pGrid->U[k][j][i].E = U1d.E;
00134 #endif
00135 #if (NSCALARS > 0)
00136 pGrid->U[k][j][i].s[0] = U1d.s[0];
00137 #endif
00138 }
00139 }
00140 }
00141 break;
00142
00143
00144 case 2:
00145 for (k=kl; k<=ku; k++) {
00146 for (j=jl; j<=ju; j++) {
00147 for (i=il; i<=iu; i++) {
00148 cc_pos(pGrid, i, j, k, &x1, &x2, &x3);
00149
00150
00151 if (x2 <= 0.0) {
00152 U1d = Ul;
00153 } else {
00154 U1d = Ur;
00155 }
00156
00157
00158 pGrid->U[k][j][i].d = U1d.d;
00159 pGrid->U[k][j][i].M1 = U1d.Mz;
00160 pGrid->U[k][j][i].M2 = U1d.Mx;
00161 pGrid->U[k][j][i].M3 = U1d.My;
00162 #ifdef MHD
00163 pGrid->B1i[k][j][i] = U1d.Bz;
00164 pGrid->B2i[k][j][i] = Bxl;
00165 pGrid->B3i[k][j][i] = U1d.By;
00166 pGrid->U[k][j][i].B1c = U1d.Bz;
00167 pGrid->U[k][j][i].B2c = Bxl;
00168 pGrid->U[k][j][i].B3c = U1d.By;
00169 #endif
00170 #ifdef ADIABATIC
00171 pGrid->U[k][j][i].E = U1d.E;
00172 #endif
00173 #if (NSCALARS > 0)
00174 pGrid->U[k][j][i].s[0] = U1d.s[0];
00175 #endif
00176 }
00177 }
00178 }
00179 break;
00180
00181
00182 case 3:
00183 for (k=kl; k<=ku; k++) {
00184 for (j=jl; j<=ju; j++) {
00185 for (i=il; i<=iu; i++) {
00186 cc_pos(pGrid, i, j, k, &x1, &x2, &x3);
00187
00188
00189 if (x3 <= 0.0) {
00190 U1d = Ul;
00191 } else {
00192 U1d = Ur;
00193 }
00194
00195
00196 pGrid->U[k][j][i].d = U1d.d;
00197 pGrid->U[k][j][i].M1 = U1d.My;
00198 pGrid->U[k][j][i].M2 = U1d.Mz;
00199 pGrid->U[k][j][i].M3 = U1d.Mx;
00200 #ifdef MHD
00201 pGrid->B1i[k][j][i] = U1d.By;
00202 pGrid->B2i[k][j][i] = U1d.Bz;
00203 pGrid->B3i[k][j][i] = Bxl;
00204 pGrid->U[k][j][i].B1c = U1d.By;
00205 pGrid->U[k][j][i].B2c = U1d.Bz;
00206 pGrid->U[k][j][i].B3c = Bxl;
00207 #endif
00208 #ifdef ADIABATIC
00209 pGrid->U[k][j][i].E = U1d.E;
00210 #endif
00211 #if (NSCALARS > 0)
00212 pGrid->U[k][j][i].s[0] = U1d.s[0];
00213 #endif
00214 }
00215 }
00216 }
00217 break;
00218 default:
00219 ath_error("[shkset1d]: invalid shk_dir = %i\n",shk_dir);
00220 }
00221
00222 return;
00223 }
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 void problem_write_restart(MeshS *pM, FILE *fp)
00237 {
00238 return;
00239 }
00240
00241 void problem_read_restart(MeshS *pM, FILE *fp)
00242 {
00243 return;
00244 }
00245
00246 ConsFun_t get_usr_expr(const char *expr)
00247 {
00248 return NULL;
00249 }
00250
00251 VOutFun_t get_usr_out_fun(const char *name){
00252 return NULL;
00253 }
00254
00255 void Userwork_in_loop(MeshS *pM)
00256 {
00257 return;
00258 }
00259
00260 void Userwork_after_loop(MeshS *pM)
00261 {
00262 return;
00263 }