00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009 #include <math.h>
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012 #include <string.h>
00013 #include "defs.h"
00014 #include "athena.h"
00015 #include "globals.h"
00016 #include "prototypes.h"
00017
00018
00019
00020
00021
00022
00023
00024 static Real omega0,rho0,bz0,vz0,Pgas0,amp,R0,phi0,z0,rad,alpha,x2min,x2max;
00025 static int iprob;
00026 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00027 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00028 Real d(const Real x1, const Real x2, const Real x3);
00029 Real M2(const Real x1, const Real x2, const Real x3);
00030 void cyladvect_ix1(GridS *pG);
00031 void cyladvect_ox1(GridS *pG);
00032 static ConsS ***RootSoln=NULL;
00033
00034
00035
00036
00037 void problem(DomainS *pDomain)
00038 {
00039 GridS *pG = pDomain->Grid;
00040 int i,j,k;
00041 int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00042 int nx1,nx2,nx3;
00043 Real x1,x2,x3;
00044 Real Eint,Emag,Ekin;
00045 int mask=0;
00046
00047 is = pG->is; ie = pG->ie; nx1 = ie-is+1;
00048 js = pG->js; je = pG->je; nx2 = je-js+1;
00049 ks = pG->ks; ke = pG->ke; nx3 = ke-ks+1;
00050
00051 il = is-nghost*(nx1>1); iu = ie+nghost*(nx1>1); nx1 = iu-il+1;
00052 jl = js-nghost*(nx2>1); ju = je+nghost*(nx2>1); nx2 = ju-jl+1;
00053 kl = ks-nghost*(nx3>1); ku = ke+nghost*(nx3>1); nx3 = ku-kl+1;
00054
00055 #ifndef CYLINDRICAL
00056 ath_error("[cyladvect]: This problem only works in cylindrical!\n");
00057 #endif
00058
00059 if (nx1==1) {
00060 ath_error("[cyladvect]: This problem can only be run in 2D or 3D!\n");
00061 }
00062 else if (nx2==1 && nx3>1) {
00063 ath_error("[cyladvect]: Only (R,phi) can be used in 2D!\n");
00064 }
00065
00066
00067 if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS))) == NULL)
00068 ath_error("[cyladvect]: Error allocating memory for solution\n");
00069
00070
00071 x2min = par_getd("domain1", "x2min");
00072 x2max = par_getd("domain1", "x2max");
00073 omega0 = par_getd("problem", "omega0");
00074 rho0 = par_getd("problem", "rho0");
00075 bz0 = par_getd("problem", "bz0");
00076 vz0 = par_getd("problem", "vz0");
00077 Pgas0 = par_getd("problem", "Pgas0");
00078 amp = par_getd("problem", "amp");
00079 R0 = par_getd("problem", "R0");
00080 phi0 = par_getd("problem", "phi0");
00081 z0 = par_getd("problem", "z0");
00082 rad = par_getd("problem", "rad");
00083 iprob = par_geti("problem", "iprob");
00084
00085
00086 alpha = 2.0*PI/(x2max-x2min);
00087
00088
00089 for (k=kl; k<=ku; k++) {
00090 for (j=jl; j<=ju; j++) {
00091 for (i=il; i<=iu; i++) {
00092 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00093 memset(&(pG->U[k][j][i]),0.0,sizeof(ConsS));
00094
00095 switch (iprob) {
00096 case 1:
00097 mask = (int)((x1>=R0-rad) && (x1<=R0+rad) && (x2>=phi0-rad/R0) && (x2<=phi0+rad/R0));
00098 pG->U[k][j][i].d = rho0*(1.0 + amp*mask);
00099 pG->U[k][j][i].M2 = pG->U[k][j][i].d*x1vc(pG,i)*omega0;
00100 #ifdef MHD
00101 pG->U[k][j][i].B3c = bz0*mask;
00102 pG->B3i[k][j][i] = bz0*mask;
00103 #endif
00104 break;
00105 case 2:
00106 case 3:
00107 pG->U[k][j][i].d = avg2d(d, pG,i,j,k);
00108 pG->U[k][j][i].M2 = avg2d(M2,pG,i,j,k);
00109 #ifdef MHD
00110 pG->U[k][j][i].B3c = bz0;
00111 pG->B3i[k][j][i] = bz0;
00112 #endif
00113 break;
00114 default:
00115 ath_error("[cyladvect]: Not an accepted problem number\n");
00116 }
00117
00118 pG->U[k][j][i].M3 = pG->U[k][j][i].d*vz0;
00119
00120 #ifndef ISOTHERMAL
00121
00122 Emag = 0.0;
00123 #ifdef MHD
00124 Emag = 0.5*SQR(pG->U[k][j][i].B3c);
00125 #endif
00126 Eint = (Pgas0-Emag)/Gamma_1;
00127 Ekin = 0.5*(SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2) + SQR(pG->U[k][j][i].M3))/pG->U[k][j][i].d;
00128
00129 pG->U[k][j][i].E = Eint + Emag + Ekin;
00130 #endif
00131
00132
00133 RootSoln[k][j][i] = pG->U[k][j][i];
00134 }
00135 }
00136 }
00137
00138 StaticGravPot = grav_pot;
00139 x1GravAcc = grav_acc;
00140 if (iprob==2) {
00141 bvals_mhd_fun(pDomain, left_x1, cyladvect_ix1);
00142 bvals_mhd_fun(pDomain, right_x1, cyladvect_ox1);
00143 } else {
00144 bvals_mhd_fun(pDomain, left_x1, do_nothing_bc);
00145 bvals_mhd_fun(pDomain, right_x1, do_nothing_bc);
00146 }
00147
00148 return;
00149 }
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 void problem_write_restart(MeshS *pM, FILE *fp)
00163 {
00164 return;
00165 }
00166
00167 void problem_read_restart(MeshS *pM, FILE *fp)
00168 {
00169 return;
00170 }
00171
00172 ConsFun_t get_usr_expr(const char *expr)
00173 {
00174 return NULL;
00175 }
00176
00177 VOutFun_t get_usr_out_fun(const char *name){
00178 return NULL;
00179 }
00180
00181 void Userwork_in_loop(MeshS *pM)
00182 {
00183 #ifdef MHD
00184 printf("Max divB = %1.10e\n", compute_div_b(pM->Domain[0][0].Grid));
00185 #endif
00186 }
00187
00188 void Userwork_after_loop(MeshS *pM)
00189 {
00190 compute_l1_error("CylAdvect", pM, RootSoln, 0);
00191 }
00192
00193
00194
00195
00196
00197 Real d(const Real x1, const Real x2, const Real x3) {
00198 Real x,y,z,x0,y0,r;
00199
00200 if (iprob==2) {
00201 return rho0*(1.0 + amp*sin(alpha*x2-phi0));
00202 }
00203 else if (iprob==3) {
00204 x = x1*cos(x2); y = x1*sin(x2); z = x3;
00205 x0 = R0*cos(phi0); y0 = R0*sin(phi0);
00206 r = sqrt(SQR(x-x0) + SQR(y-y0) + SQR(z-z0));
00207 return rho0*(1.0 + exp(-0.5*SQR(3.0*r/rad)));
00208 }
00209
00210 return rho0;
00211 }
00212
00213
00214
00215 Real M2(const Real x1, const Real x2, const Real x3) {
00216 return d(x1,x2,x3)*omega0*x1;
00217 }
00218
00219
00220
00221 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00222 return 0.5*SQR(x1*omega0);
00223 }
00224
00225
00226
00227 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00228 return x1*SQR(omega0);
00229 }
00230
00231
00232
00233
00234
00235
00236
00237 void cyladvect_ix1(GridS *pG)
00238 {
00239 int is = pG->is;
00240 int js = pG->js, je = pG->je;
00241 int ks = pG->ks, ke = pG->ke;
00242 int i,j,k;
00243 Real Eint,Emag,Ekin;
00244
00245 phi0 = alpha*(omega0*pG->time + x2min);
00246 for (k=ks; k<=ke; k++) {
00247 for (j=js; j<=je; j++) {
00248 for (i=1; i<=nghost; i++) {
00249 pG->U[k][j][is-i].d = avg2d(d,pG,is-i,j,k);
00250 pG->U[k][j][is-i].M2 = avg2d(M2,pG,is-i,j,k);
00251 pG->U[k][j][is-i].M3 = pG->U[k][j][is-i].d*vz0;
00252 #ifdef MHD
00253 pG->U[k][j][is-i].B3c = bz0;
00254 pG->B3c[k][j][is-i] = bz0;
00255 #endif
00256
00257 #ifndef ISOTHERMAL
00258 Emag = 0.0;
00259 #ifdef MHD
00260 Emag = 0.5*SQR(pG->U[k][j][is-i].B3c);
00261 #endif
00262 Eint = (Pgas0-Emag)/Gamma_1;
00263 Ekin = 0.5*(SQR(pG->U[k][j][is-i].M1) + SQR(pG->U[k][j][is-i].M2) + SQR(pG->U[k][j][is-i].M3))/pG->U[k][j][is-i].d;
00264
00265 pG->U[k][j][is-i].E = Eint + Emag + Ekin;
00266 #endif
00267 }
00268 }
00269 }
00270
00271 return;
00272 }
00273
00274
00275
00276
00277
00278
00279 void cyladvect_ox1(GridS *pG)
00280 {
00281 int ie = pG->ie;
00282 int js = pG->js, je = pG->je;
00283 int ks = pG->ks, ke = pG->ke;
00284 int i,j,k;
00285 Real Eint,Emag,Ekin;
00286
00287 phi0 = alpha*(omega0*pG->time + x2min);
00288 for (k=ks; k<=ke; k++) {
00289 for (j=js; j<=je; j++) {
00290 for (i=1; i<=nghost; i++) {
00291 pG->U[k][j][ie+i].d = avg2d(d,pG,ie+i,j,k);
00292 pG->U[k][j][ie+i].M2 = avg2d(M2,pG,ie+i,j,k);
00293 pG->U[k][j][ie+i].M3 = pG->U[k][j][ie+i].d*vz0;
00294 #ifdef MHD
00295 pG->U[k][j][ie+i].B3c = bz0;
00296 pG->B3c[k][j][ie+i] = bz0;
00297 #endif
00298
00299 #ifndef ISOTHERMAL
00300 Emag = 0.0;
00301 #ifdef MHD
00302 Emag = 0.5*SQR(pG->U[k][j][ie+i].B3c);
00303 #endif
00304 Eint = (Pgas0-Emag)/Gamma_1;
00305 Ekin = 0.5*(SQR(pG->U[k][j][ie+i].M1) + SQR(pG->U[k][j][ie+i].M2) + SQR(pG->U[k][j][ie+i].M3))/pG->U[k][j][ie+i].d;
00306
00307 pG->U[k][j][ie+i].E = Eint + Emag + Ekin;
00308 #endif
00309 }
00310 }
00311 }
00312
00313 return;
00314 }