00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <math.h>
00014 #include <stdio.h>
00015 #include <stdlib.h>
00016 #include "defs.h"
00017 #include "athena.h"
00018 #include "globals.h"
00019 #include "prototypes.h"
00020
00021 #ifdef ISOTHERMAL
00022 #error : The dmr test only works for adiabatic EOS.
00023 #endif
00024 #ifndef HYDRO
00025 #error : The dmr test only works for hydro.
00026 #endif
00027
00028
00029
00030
00031
00032
00033
00034
00035 void dmrbv_iib(GridS *pGrid);
00036 void dmrbv_ijb(GridS *pGrid);
00037 void dmrbv_ojb(GridS *pGrid);
00038
00039
00040
00041
00042
00043 void problem(DomainS *pDomain)
00044 {
00045 GridS *pGrid = pDomain->Grid;
00046 int i=0,j=0;
00047 int is,ie,js,je,ks;
00048 Real d0,e0,u0,v0,x1_shock,x1,x2,x3;
00049
00050 is = pGrid->is; ie = pGrid->ie;
00051 js = pGrid->js; je = pGrid->je;
00052 ks = pGrid->ks;
00053 if (pGrid->Nx[0] == 1 || pGrid->Nx[1] == 1) {
00054 ath_error("[dmr]: this test only works with Nx1 & Nx2 > 1\n");
00055 }
00056 if (pGrid->Nx[2] > 1) {
00057 ath_error("[dmr]: this test only works for 2D problems, with Nx3=1\n");
00058 }
00059
00060
00061
00062 d0 = 8.0;
00063 e0 = 291.25;
00064 u0 = 8.25*sqrt(3.0)/2.0;
00065 v0 = -8.25*0.5;
00066 for (j=js; j<=je; j++) {
00067 for (i=is; i<=ie; i++) {
00068 cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00069 x1_shock = 0.1666666666 + x2/sqrt((double)3.0);
00070
00071 pGrid->U[ks][j][i].d = 1.4;
00072 pGrid->U[ks][j][i].E = 2.5;
00073 pGrid->U[ks][j][i].M1 = 0.0;
00074 pGrid->U[ks][j][i].M2 = 0.0;
00075
00076 if (x1 < x1_shock) {
00077 pGrid->U[ks][j][i].d = d0;
00078 pGrid->U[ks][j][i].E = e0 + 0.5*d0*(u0*u0+v0*v0);
00079 pGrid->U[ks][j][i].M1 = d0*u0;
00080 pGrid->U[ks][j][i].M2 = d0*v0;
00081 }
00082 }
00083 }
00084
00085
00086
00087 if (pDomain->Disp[0] == 0) bvals_mhd_fun(pDomain, left_x1, dmrbv_iib);
00088 if (pDomain->Disp[1] == 0) bvals_mhd_fun(pDomain, left_x2, dmrbv_ijb);
00089 if (pDomain->MaxX[1] == pDomain->RootMaxX[1])
00090 bvals_mhd_fun(pDomain, right_x2, dmrbv_ojb);
00091 }
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 void problem_write_restart(MeshS *pM, FILE *fp)
00105 {
00106 return;
00107 }
00108
00109 void problem_read_restart(MeshS *pM, FILE *fp)
00110 {
00111 return;
00112 }
00113
00114 ConsFun_t get_usr_expr(const char *expr)
00115 {
00116 return NULL;
00117 }
00118
00119 VOutFun_t get_usr_out_fun(const char *name){
00120 return NULL;
00121 }
00122
00123 void Userwork_in_loop(MeshS *pM)
00124 {
00125 }
00126
00127 void Userwork_after_loop(MeshS *pM)
00128 {
00129 }
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 void dmrbv_iib(GridS *pGrid)
00141 {
00142 int i=0,j=0;
00143 int is,ie,js,je,ks,jl,ju;
00144 Real d0,e0,u0,v0;
00145
00146 d0 = 8.0;
00147 e0 = 291.25;
00148 u0 = 8.25*sqrt(3.0)/2.0;
00149 v0 = -8.25*0.5;
00150
00151 is = pGrid->is; ie = pGrid->ie;
00152 js = pGrid->js; je = pGrid->je;
00153 ks = pGrid->ks;
00154 ju = pGrid->je + nghost;
00155 jl = pGrid->js - nghost;
00156
00157 for (j=jl; j<=ju; j++) {
00158 for (i=1; i<=nghost; i++) {
00159 pGrid->U[ks][j][is-i].d = d0;
00160 pGrid->U[ks][j][is-i].M1 = d0*u0;
00161 pGrid->U[ks][j][is-i].M2 = d0*v0;
00162 pGrid->U[ks][j][is-i].E = e0 + 0.5*d0*(u0*u0+v0*v0);
00163 }
00164 }
00165 }
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 void dmrbv_ijb(GridS *pGrid)
00176 {
00177 int i=0,j=0;
00178 int is,ie,js,je,ks,il,iu;
00179 Real d0,e0,u0,v0,x1,x2,x3;
00180
00181 d0 = 8.0;
00182 e0 = 291.25;
00183 u0 = 8.25*sqrt(3.0)/2.0;
00184 v0 = -8.25*0.5;
00185
00186 is = pGrid->is; ie = pGrid->ie;
00187 js = pGrid->js; je = pGrid->je;
00188 ks = pGrid->ks;
00189 iu = pGrid->ie + nghost;
00190 il = pGrid->is - nghost;
00191
00192 for (j=1; j<=nghost; j++) {
00193 for (i=il; i<=iu; i++) {
00194 cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00195 if (x1 < 0.1666666666) {
00196
00197 pGrid->U[ks][js-j][i].d = d0;
00198 pGrid->U[ks][js-j][i].M1 = d0*u0;
00199 pGrid->U[ks][js-j][i].M2 = d0*v0;
00200 pGrid->U[ks][js-j][i].E = e0 + 0.5*d0*(u0*u0+v0*v0);
00201 } else {
00202
00203 pGrid->U[ks][js-j][i].d = pGrid->U[ks][js+(j-1)][i].d;
00204 pGrid->U[ks][js-j][i].M1 = pGrid->U[ks][js+(j-1)][i].M1;
00205 pGrid->U[ks][js-j][i].M2 = -pGrid->U[ks][js+(j-1)][i].M2;
00206 pGrid->U[ks][js-j][i].E = pGrid->U[ks][js+(j-1)][i].E;
00207 }
00208 }
00209 }
00210 }
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222 void dmrbv_ojb(GridS *pGrid)
00223 {
00224 int i=0,j=0;
00225 int is,ie,js,je,ks,il,iu;
00226 Real d0,e0,u0,v0,x1_shock,x1,x2,x3;
00227
00228 d0 = 8.0;
00229 e0 = 291.25;
00230 u0 = 8.25*sqrt(3.0)/2.0;
00231 v0 = -8.25*0.5;
00232 x1_shock = 0.1666666666 + (1.0 + 20.0*pGrid->time)/sqrt(3.0);
00233
00234 is = pGrid->is; ie = pGrid->ie;
00235 js = pGrid->js; je = pGrid->je;
00236 ks = pGrid->ks;
00237 iu = pGrid->ie + nghost;
00238 il = pGrid->is - nghost;
00239
00240 for (j=1; j<=nghost; j++) {
00241 for (i=il; i<=iu; i++) {
00242 cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00243 if (x1 < x1_shock) {
00244
00245 pGrid->U[ks][je+j][i].d = d0;
00246 pGrid->U[ks][je+j][i].M1 = d0*u0;
00247 pGrid->U[ks][je+j][i].M2 = d0*v0;
00248 pGrid->U[ks][je+j][i].E = e0 + 0.5*d0*(u0*u0+v0*v0);
00249 } else {
00250
00251 pGrid->U[ks][je+j][i].d = 1.4;
00252 pGrid->U[ks][je+j][i].M1 = 0.0;
00253 pGrid->U[ks][je+j][i].M2 = 0.0;
00254 pGrid->U[ks][je+j][i].E = 2.5;
00255 }
00256 }
00257 }
00258 }