Go to the documentation of this file.00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <math.h>
00012 #include <stdlib.h>
00013 #include <stdio.h>
00014 #include <string.h>
00015 #include "defs.h"
00016 #include "athena.h"
00017 #include "globals.h"
00018 #include "prototypes.h"
00019
00020 #ifdef MHD
00021 #error : The lw_implode problem only works for hydro.
00022 #endif
00023
00024
00025 static Real expr_diff_d(const GridS *pG, const int i, const int j, const int k);
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 ir,irefine,nx2;
00037 Real d_in,p_in,d_out,p_out,Ly,rootdx2;
00038
00039
00040 if (pGrid->Nx[0] <= 1 || pGrid->Nx[1] <= 1) {
00041 ath_error("[problem]: This problem requires Nx1 > 1, Nx2 > 1\n");
00042 }
00043
00044 d_in = par_getd("problem","d_in");
00045 p_in = par_getd("problem","p_in");
00046
00047 d_out = par_getd("problem","d_out");
00048 p_out = par_getd("problem","p_out");
00049
00050
00051
00052 irefine = 1;
00053 for (ir=1;ir<=pDomain->Level;ir++) irefine *= 2;
00054
00055 Ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00056 rootdx2 = pGrid->dx2*((double)(irefine));
00057 nx2 = (int)(Ly/rootdx2);
00058 nx2 /= 2;
00059
00060
00061 for (k=ks; k<=ke; k++) {
00062 for (j=js; j<=je; j++) {
00063 for (i=is; i<=ie; i++) {
00064 pGrid->U[k][j][i].M1 = 0.0;
00065 pGrid->U[k][j][i].M2 = 0.0;
00066 pGrid->U[k][j][i].M3 = 0.0;
00067
00068 if(((j-js + pDomain->Disp[1])+(i-is + pDomain->Disp[0])) > (nx2*irefine)) {
00069 pGrid->U[k][j][i].d = d_out;
00070 #ifndef ISOTHERMAL
00071 pGrid->U[k][j][i].E = p_out/Gamma_1;
00072 #endif
00073 } else {
00074 pGrid->U[k][j][i].d = d_in;
00075 #ifndef ISOTHERMAL
00076 pGrid->U[k][j][i].E = p_in/Gamma_1;
00077 #endif
00078 }
00079 }
00080 }
00081 }
00082
00083 return;
00084 }
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 void problem_write_restart(MeshS *pM, FILE *fp)
00098 {
00099 return;
00100 }
00101
00102 void problem_read_restart(MeshS *pM, FILE *fp)
00103 {
00104 return;
00105 }
00106
00107 ConsFun_t get_usr_expr(const char *expr)
00108 {
00109 if(strcmp(expr,"diff_d")==0) return expr_diff_d;
00110 return NULL;
00111 }
00112
00113 VOutFun_t get_usr_out_fun(const char *name){
00114 return NULL;
00115 }
00116
00117 void Userwork_in_loop(MeshS *pM)
00118 {
00119 }
00120
00121 void Userwork_after_loop(MeshS *pM)
00122 {
00123 }
00124
00125
00126
00127 static Real expr_diff_d(const GridS *pG, const int i, const int j, const int k)
00128 {
00129 return (pG->U[k][j][i].d - pG->U[k][i][j].d);
00130 }