00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008 #include <math.h>
00009 #include <stdio.h>
00010 #include <stdlib.h>
00011 #include <string.h>
00012 #include "defs.h"
00013 #include "athena.h"
00014 #include "globals.h"
00015 #include "prototypes.h"
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 static Real omega0,rho0;
00026 static int iprob;
00027 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00028 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00029 Real M2(const Real x1, const Real x2, const Real x3);
00030 static ConsS ***RootSoln=NULL;
00031
00032
00033
00034
00035 void problem(DomainS *pDomain)
00036 {
00037 GridS *pG = pDomain->Grid;
00038 int i,j,k;
00039 int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00040 int nx1,nx2,nx3;
00041 Real x1,x2,x3;
00042 Real Eint,Emag,Ekin,vphi0,vz0,pgas0;
00043
00044 is = pG->is; ie = pG->ie; nx1 = ie-is+1;
00045 js = pG->js; je = pG->je; nx2 = je-js+1;
00046 ks = pG->ks; ke = pG->ke; nx3 = ke-ks+1;
00047
00048 il = is-nghost*(nx1>1); iu = ie+nghost*(nx1>1); nx1 = iu-il+1;
00049 jl = js-nghost*(nx2>1); ju = je+nghost*(nx2>1); nx2 = ju-jl+1;
00050 kl = ks-nghost*(nx3>1); ku = ke+nghost*(nx3>1); nx3 = ku-kl+1;
00051
00052 #ifndef MHD
00053 ath_error("[cylbphi]: This problem only works in MHD!\n");
00054 #endif
00055 #ifndef CYLINDRICAL
00056 ath_error("[cylbphi]: This problem only works in cylindrical!\n");
00057 #endif
00058
00059 if (nx1==1) {
00060 ath_error("[cylbphi]: This problem can only be run in 2D or 3D!\n");
00061 }
00062 else if (nx2==1 && nx3>1) {
00063 ath_error("[cylbphi]: Only (R,phi) can be used in 2D!\n");
00064 }
00065
00066 omega0 = par_getd("problem", "omega0");
00067 vz0 = par_getd("problem", "vz0");
00068 bphi0 = par_getd("problem", "bphi0");
00069 rho0 = par_getd("problem", "rho0");
00070 pgas0 = par_getd("problem", "pgas0");
00071 iprob = par_geti("problem", "iprob");
00072
00073
00074
00075 if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS))) == NULL)
00076 ath_error("[cylbphi]: Error allocating memory for solution!\n");
00077
00078
00079
00080
00081
00082
00083
00084 for (k=kl; k<=ku; k++) {
00085 for (j=jl; j<=ju; j++) {
00086 for (i=il; i<=iu; i++) {
00087 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00088 memset(&(pG->U[k][j][i]),0.0,sizeof(Gas));
00089
00090 pG->U[k][j][i].d = rho0;
00091 pG->U[k][j][i].M1 = 0.0;
00092 pG->U[k][j][i].M2 = avg1d(M2,pG,i,j,k);
00093 pG->U[k][j][i].M3 = pG->U[k][j][i].d*vz0;
00094 switch (iprob) {
00095 case 1:
00096 pG->B2i[k][j][i] = bphi0;
00097 pG->U[k][j][i].B2c = bphi0;
00098 break;
00099 case 2:
00100 pG->B2i[k][j][i] = bphi0/x1;
00101 pG->U[k][j][i].B2c = bphi0/x1;
00102 break;
00103 default: printf("[cylbphi]: Not an accepted problem number\n");
00104 }
00105
00106
00107 Eint = pgas0/Gamma_1;
00108 Emag = 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c)
00109 + SQR(pG->U[k][j][i].B3c));
00110 Ekin = 0.5*(SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2)
00111 + SQR(pG->U[k][j][i].M3))/pG->U[k][j][i].d;
00112 pG->U[k][j][i].E = Eint + Emag + Ekin;
00113
00114
00115 RootSoln[k][j][i] = pG->U[ks][j][i];
00116 }
00117 }
00118 }
00119
00120 StaticGravPot = grav_pot;
00121 x1GravAcc = grav_acc;
00122 bvals_mhd_fun(pDomain,left_x1,do_nothing_bc);
00123 bvals_mhd_fun(pDomain,right_x1,do_nothing_bc);
00124
00125 return;
00126 }
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 void problem_write_restart(MeshS *pM, FILE *fp)
00139 {
00140 return;
00141 }
00142
00143 void problem_read_restart(MeshS *pM, FILE *fp)
00144 {
00145 return;
00146 }
00147
00148 ConsFun_t get_usr_expr(const char *expr)
00149 {
00150 return NULL;
00151 }
00152
00153 VOutFun_t get_usr_out_fun(const char *name){
00154 return NULL;
00155 }
00156
00157 void Userwork_in_loop(MeshS *pM)
00158 {
00159
00160 }
00161
00162 void Userwork_after_loop(MeshS *pM)
00163 {
00164 compute_l1_error("CylBPhi", pM, RootSoln, 1);
00165 }
00166
00167
00168
00169
00170
00171 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00172 switch (iprob) {
00173 case 1: return 0.5*SQR(x1*omega0) - (SQR(bphi0)/rho0)*log(x1);
00174 break;
00175 case 2: return 0.5*SQR(x1*omega0);
00176 break;
00177 default: return 0.0;
00178 }
00179 }
00180
00181
00182
00183 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00184 switch (iprob) {
00185 case 1: return x1*SQR(omega0) - SQR(bphi0)/(rho0*x1);
00186 break;
00187 case 2: return x1*SQR(omega0);
00188 break;
00189 default: return 0.0;
00190 }
00191 }
00192
00193
00194
00195 Real M2(const Real x1, const Real x2, const Real x3) {
00196 return rho0*omega0*x1;
00197 }