00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007 #include <math.h>
00008 #include <stdio.h>
00009 #include <stdlib.h>
00010 #include "defs.h"
00011 #include "athena.h"
00012 #include "globals.h"
00013 #include "prototypes.h"
00014
00015
00016
00017
00018
00019
00020 void jet_iib(GridS *pGrid);
00021
00022
00023
00024 static Real rjet,Bxjet;
00025 static Prim1DS Wjet;
00026 static Cons1DS Ujet;
00027 static Real x1_mid,x2_mid,x3_mid;
00028
00029
00030
00031
00032 void problem(DomainS *pDomain){
00033 GridS *pGrid=(pDomain->Grid);
00034 int i, is = pGrid->is, ie = pGrid->ie;
00035 int j, js = pGrid->js, je = pGrid->je;
00036 int k, ks = pGrid->ks, ke = pGrid->ke;
00037 int il,iu,jl,ju,kl,ku;
00038 Prim1DS W;
00039 Cons1DS U;
00040 Real x1_min, x1_max;
00041 Real x2_min, x2_max;
00042 Real x3_min, x3_max;
00043 Real Bx=0.0;
00044 Bxjet = 0.0;
00045
00046
00047
00048 W.d = par_getd("problem", "d");
00049 W.P = par_getd("problem", "p");
00050 W.Vx = par_getd("problem", "vx");
00051 W.Vy = par_getd("problem", "vy");
00052 W.Vz = par_getd("problem", "vz");
00053 #ifdef MHD
00054 Bx = par_getd("problem", "bx");
00055 W.By = par_getd("problem", "by");
00056 W.Bz = par_getd("problem", "bz");
00057 #endif
00058
00059 Wjet.d = par_getd("problem", "djet");
00060 Wjet.P = par_getd("problem", "pjet");
00061 Wjet.Vx = par_getd("problem", "vxjet");
00062 Wjet.Vy = par_getd("problem", "vyjet");
00063 Wjet.Vz = par_getd("problem", "vzjet");
00064 #ifdef MHD
00065 Bxjet = par_getd("problem", "bxjet");
00066 Wjet.By = par_getd("problem", "byjet");
00067 Wjet.Bz = par_getd("problem", "bzjet");
00068 #endif
00069
00070 rjet = par_getd("problem", "rjet");
00071
00072 U = Prim1D_to_Cons1D(&W, &Bx);
00073 Ujet = Prim1D_to_Cons1D(&Wjet, &Bxjet);
00074
00075 x1_min = pDomain->RootMinX[0];
00076 x1_max = pDomain->RootMaxX[0];
00077 x2_min = pDomain->RootMinX[1];
00078 x2_max = pDomain->RootMaxX[1];
00079 x3_min = pDomain->RootMinX[2];
00080 x3_max = pDomain->RootMaxX[2];
00081
00082 x1_mid = 0.5 * (x1_max + x1_min);
00083 x2_mid = 0.5 * (x2_max + x2_min);
00084 x3_mid = 0.5 * (x3_max + x3_min);
00085
00086
00087 iu = ie + nghost;
00088 il = is - nghost;
00089 ju = je + nghost;
00090 jl = js - nghost;
00091 if(pGrid->Nx[2] > 1){
00092 ku = pGrid->ke + nghost;
00093 kl = pGrid->ks - nghost;
00094 }
00095 else{
00096 ku = pGrid->ke;
00097 kl = pGrid->ks;
00098 }
00099
00100
00101
00102 for(k=kl; k<=ku; k++){
00103 for(j=jl; j<=ju; j++){
00104 for(i=il; i<=iu; i++){
00105 pGrid->U[k][j][i].d = U.d;
00106 pGrid->U[k][j][i].M1 = U.Mx;
00107 pGrid->U[k][j][i].M2 = U.My;
00108 pGrid->U[k][j][i].M3 = U.Mz;
00109 pGrid->U[k][j][i].E = U.E;
00110 #ifdef MHD
00111 pGrid->U[k][j][i].B1c = Bx;
00112 pGrid->U[k][j][i].B2c = U.By;
00113 pGrid->U[k][j][i].B3c = U.Bz;
00114 pGrid->B1i[k][j][i] = Bx;
00115 pGrid->B2i[k][j][i] = U.By;
00116 pGrid->B3i[k][j][i] = U.Bz;
00117 #endif
00118 }
00119 }
00120 }
00121
00122
00123
00124 if (pDomain->Disp[0] == 0) bvals_mhd_fun(pDomain,left_x1,jet_iib);
00125
00126 }
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 void problem_write_restart(MeshS *pM, FILE *fp)
00141 {
00142 return;
00143 }
00144
00145 void problem_read_restart(MeshS *pM, FILE *fp)
00146 {
00147 int nl,nd;
00148
00149 Wjet.d = par_getd("problem", "djet");
00150 Wjet.P = par_getd("problem", "pjet");
00151 Wjet.Vx = par_getd("problem", "vxjet");
00152 Wjet.Vy = par_getd("problem", "vyjet");
00153 Wjet.Vz = par_getd("problem", "vzjet");
00154 #ifdef MHD
00155 Bxjet = par_getd("problem", "bxjet");
00156 Wjet.By = par_getd("problem", "byjet");
00157 Wjet.Bz = par_getd("problem", "bzjet");
00158 #endif
00159
00160 rjet = par_getd("problem", "rjet");
00161 Ujet = Prim1D_to_Cons1D(&Wjet, &Bxjet);
00162
00163 for (nl=0; nl<pM->NLevels; nl++){
00164 for (nd=0; nd<pM->DomainsPerLevel[nl]; nd++){
00165 if (pM->Domain[nl][nd].Grid != NULL) {
00166 if (pM->Domain[nl][nd].Disp[0] == 0)
00167 bvals_mhd_fun(&(pM->Domain[nl][nd]),left_x1,jet_iib);
00168 }
00169 }}
00170
00171 return;
00172 }
00173
00174 ConsFun_t get_usr_expr(const char *expr)
00175 {
00176 return NULL;
00177 }
00178
00179 VOutFun_t get_usr_out_fun(const char *name){
00180 return NULL;
00181 }
00182
00183 void Userwork_in_loop(MeshS *pM)
00184 {
00185 return;
00186 }
00187
00188 void Userwork_after_loop(MeshS *pM)
00189 {
00190 return;
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 void jet_iib(GridS *pGrid){
00203 int i, is = pGrid->is;
00204 int j, js = pGrid->js, je = pGrid->je;
00205 int k, ks = pGrid->ks, ke = pGrid->ke;
00206 Real rad,x1,x2,x3;
00207
00208 for(k=ks; k<=ke; k++){
00209 for(j=js; j<=je; j++){
00210 for(i=1; i<=nghost; i++){
00211 cc_pos(pGrid,(is-i),j,k,&x1,&x2,&x3);
00212 rad = sqrt(SQR(x2 - x2_mid) + SQR(x3 - x3_mid));
00213
00214 if(rad <= rjet){
00215 pGrid->U[k][j][is-i].d = Ujet.d;
00216 pGrid->U[k][j][is-i].M1 = Ujet.Mx;
00217 pGrid->U[k][j][is-i].M2 = Ujet.My;
00218 pGrid->U[k][j][is-i].M3 = Ujet.Mz;
00219 pGrid->U[k][j][is-i].E = Ujet.E;
00220 #ifdef MHD
00221 pGrid->U[k][j][is-i].B1c = Bxjet;
00222 pGrid->U[k][j][is-i].B2c = Ujet.By;
00223 pGrid->U[k][j][is-i].B3c = Ujet.Bz;
00224 pGrid->B1i[k][j][is-i] = Bxjet;
00225 pGrid->B2i[k][j][is-i] = Ujet.By;
00226 pGrid->B3i[k][j][is-i] = Ujet.Bz;
00227 #endif
00228 } else{
00229 pGrid->U[k][j][is-i] = pGrid->U[k][j][is+(i-1)];
00230 pGrid->U[k][j][is-i].M1 = -pGrid->U[k][j][is-i].M1;
00231 #ifdef MHD
00232 pGrid->U[k][j][is-i].B2c = -pGrid->U[k][j][is-i].B2c;
00233 pGrid->U[k][j][is-i].B3c = -pGrid->U[k][j][is-i].B3c;
00234 pGrid->B1i[k][j][is-i] = pGrid->B1i[k][j][is+i];
00235 pGrid->B2i[k][j][is-i] = -pGrid->B2i[k][j][is+(i-1)];
00236 pGrid->B3i[k][j][is-i] = -pGrid->B3i[k][j][is+(i-1)];
00237 #endif
00238 }
00239 }
00240 }
00241 }
00242 }