00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <math.h>
00021 #include <stdlib.h>
00022 #include <stdio.h>
00023 #include "defs.h"
00024 #include "athena.h"
00025 #include "globals.h"
00026 #include "prototypes.h"
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 void noh3d_oib(GridS *pGrid);
00037 void noh3d_ojb(GridS *pGrid);
00038 void noh3d_okb(GridS *pGrid);
00039
00040 #ifdef MHD
00041 #error : This is not a MHD problem.
00042 #endif
00043
00044
00045
00046
00047
00048 void problem(DomainS *pD)
00049 {
00050 GridS *pGrid=(pD->Grid);
00051 int i,j,k;
00052 Real x1,x2,x3,r;
00053
00054
00055
00056 for (k=pGrid->ks; k<=pGrid->ke; k++) {
00057 for (j=pGrid->js; j<=pGrid->je; j++) {
00058 for (i=pGrid->is; i<=pGrid->ie; i++) {
00059 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00060 if (pGrid->Nx[2] > 1) {
00061 r = sqrt(x1*x1 + x2*x2 + x3*x3);
00062 } else {
00063 r = sqrt(x1*x1 + x2*x2);
00064 }
00065 pGrid->U[k][j][i].d = 1.0;
00066 pGrid->U[k][j][i].M1 = -x1/r;
00067 pGrid->U[k][j][i].M2 = -x2/r;
00068 if (pGrid->Nx[2] > 1) {
00069 pGrid->U[k][j][i].M3 = -x3/r;
00070 } else {
00071 pGrid->U[k][j][i].M3 = 0.0;
00072 }
00073 pGrid->U[k][j][i].E = 1.0e-6/Gamma_1 + 0.5;
00074 }
00075 }
00076 }
00077
00078
00079
00080 bvals_mhd_fun(pD, right_x1,noh3d_oib);
00081 bvals_mhd_fun(pD, right_x2,noh3d_ojb);
00082 if (pGrid->Nx[2] > 1) bvals_mhd_fun(pD, right_x3,noh3d_okb);
00083
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 return NULL;
00110 }
00111
00112 VOutFun_t get_usr_out_fun(const char *name){
00113 return NULL;
00114 }
00115
00116 void Userwork_in_loop(MeshS *pM)
00117 {
00118 }
00119
00120 void Userwork_after_loop(MeshS *pM)
00121 {
00122 }
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 void noh3d_oib(GridS *pGrid)
00134 {
00135 int i, ie = pGrid->ie;
00136 int j, jl = pGrid->js - nghost, ju = pGrid->je + nghost;
00137 int k, kl, ku;
00138 Real x1,x2,x3,r,d0,f_t;
00139
00140 if (pGrid->Nx[2] > 1) {
00141 kl = pGrid->ks - nghost;
00142 ku = pGrid->ke + nghost;
00143 } else {
00144 kl = pGrid->ks;
00145 ku = pGrid->ke;
00146 }
00147
00148 for (k=kl; k<=ku; k++) {
00149 for (j=jl; j<=ju; j++) {
00150 for (i=ie+1; i<=ie+nghost; i++) {
00151 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00152 if (pGrid->Nx[2] > 1) {
00153 r = sqrt(x1*x1 + x2*x2 + x3*x3);
00154 f_t = (1.0 + pGrid->time/r)*(1.0 + pGrid->time/r);
00155 } else {
00156 r = sqrt(x1*x1 + x2*x2);
00157 f_t = (1.0 + pGrid->time/r);
00158 }
00159 d0 = 1.0*f_t;
00160
00161 pGrid->U[k][j][i].d = d0;
00162 pGrid->U[k][j][i].M1 = -x1*d0/r;
00163 pGrid->U[k][j][i].M2 = -x2*d0/r;
00164 if (pGrid->Nx[2] > 1) {
00165 pGrid->U[k][j][i].M3 = -x3*d0/r;
00166 pGrid->U[k][j][i].E = 1.0e-6*pow(f_t,1.0+Gamma)/Gamma_1 + 0.5*d0;
00167 } else {
00168 pGrid->U[k][j][i].M3 = 0.0;
00169 pGrid->U[k][j][i].E = 1.0e-6/Gamma_1 + 0.5*d0;
00170 }
00171 }
00172 }
00173 }
00174 }
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184 void noh3d_ojb(GridS *pGrid)
00185 {
00186 int j, je = pGrid->je;
00187 int i, il = pGrid->is - nghost, iu = pGrid->ie + nghost;
00188 int k, kl, ku;
00189 Real x1,x2,x3,r,d0,f_t;
00190
00191 if (pGrid->Nx[2] > 1) {
00192 kl = pGrid->ks - nghost;
00193 ku = pGrid->ke + nghost;
00194 } else {
00195 kl = pGrid->ks;
00196 ku = pGrid->ke;
00197 }
00198
00199 for (k=kl; k<=ku; k++) {
00200 for (j=je+1; j<=je+nghost; j++) {
00201 for (i=il; i<=iu; i++) {
00202 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00203 if (pGrid->Nx[2] > 1) {
00204 r = sqrt(x1*x1 + x2*x2 + x3*x3);
00205 f_t = (1.0 + pGrid->time/r)*(1.0 + pGrid->time/r);
00206 } else {
00207 r = sqrt(x1*x1 + x2*x2);
00208 f_t = (1.0 + pGrid->time/r);
00209 }
00210 d0 = 1.0*f_t;
00211
00212 pGrid->U[k][j][i].d = d0;
00213 pGrid->U[k][j][i].M1 = -x1*d0/r;
00214 pGrid->U[k][j][i].M2 = -x2*d0/r;
00215 if (pGrid->Nx[2] > 1) {
00216 pGrid->U[k][j][i].M3 = -x3*d0/r;
00217 pGrid->U[k][j][i].E = 1.0e-6*pow(f_t,1.0+Gamma)/Gamma_1 + 0.5*d0;
00218 } else {
00219 pGrid->U[k][j][i].M3 = 0.0;
00220 pGrid->U[k][j][i].E = 1.0e-6/Gamma_1 + 0.5*d0;
00221 }
00222 }
00223 }
00224 }
00225 }
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235 void noh3d_okb(GridS *pGrid)
00236 {
00237 int i, il = pGrid->is - nghost, iu = pGrid->ie + nghost;
00238 int j, jl = pGrid->js - nghost, ju = pGrid->je + nghost;
00239 int k, ke = pGrid->ke;
00240 Real x1,x2,x3,r,d0,f_t;
00241
00242 for (k=ke+1; k<=ke+nghost; k++) {
00243 for (j=jl; j<=ju; j++) {
00244 for (i=il; i<=iu; i++) {
00245 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00246 r = sqrt(x1*x1 + x2*x2 + x3*x3);
00247 f_t = (1.0 + pGrid->time/r)*(1.0 + pGrid->time/r);
00248 d0 = 1.0*f_t;
00249
00250 pGrid->U[k][j][i].d = d0;
00251 pGrid->U[k][j][i].M1 = -x1*d0/r;
00252 pGrid->U[k][j][i].M2 = -x2*d0/r;
00253 pGrid->U[k][j][i].M3 = -x3*d0/r;
00254 pGrid->U[k][j][i].E = 1.0e-6*pow(f_t,1.0+Gamma)/Gamma_1 + 0.5*d0;
00255 }
00256 }
00257 }
00258 }