00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <math.h>
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include "defs.h"
00033 #include "athena.h"
00034 #include "globals.h"
00035 #include "prototypes.h"
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 void shkset2d_iib(GridS *pGrid);
00046 void shkset2d_oib(GridS *pGrid);
00047 void shkset2d_ijb(GridS *pGrid);
00048 void shkset2d_ojb(GridS *pGrid);
00049
00050
00051
00052 static Real Lx,Ly;
00053 static int r1,r2;
00054
00055
00056
00057
00058 void problem(DomainS *pDomain)
00059 {
00060 GridS *pGrid = pDomain->Grid;
00061 int i, is = pGrid->is, ie = pGrid->ie;
00062 int j, js = pGrid->js, je = pGrid->je;
00063 int k, ks = pGrid->ks, ke = pGrid->ke;
00064 int kl,ku,irefine,ir,ix1,ix2,nx1,nx2,gcd;
00065 Real angle, sin_a, cos_a;
00066 Real rootdx1, rootdx2;
00067 Prim1DS Wl, Wr;
00068 Cons1DS Ul, Ur;
00069 ConsS ql, qr;
00070 Real Bxl=0.0,Bxr=0.0;
00071 div_t id;
00072
00073
00074
00075 int dll, dlr, drr, drl;
00076 Real afl_lx, afr_lx, afl_rx, afr_rx;
00077 Real afl_ly, afr_ly, afl_ry, afr_ry;
00078 Real vfl, vfr, B1r, B2r;
00079
00080 if (pGrid->Nx[1] == 1)
00081 ath_error("[shkset2d]: This problem can only be run in 2D or 3D\n");
00082
00083 if (pGrid->Nx[2] > 1){
00084 ku = pGrid->ke + nghost;
00085 kl = pGrid->ks - nghost;
00086 } else {
00087 ku = pGrid->ke;
00088 kl = pGrid->ks;
00089 }
00090
00091
00092
00093 irefine = 1;
00094 for (ir=1;ir<=pDomain->Level;ir++) irefine *= 2;
00095
00096 Lx = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00097 Ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00098
00099 rootdx1 = pGrid->dx1*((double)(irefine));
00100 rootdx2 = pGrid->dx2*((double)(irefine));
00101
00102 nx1 = (int)(Lx/rootdx1);
00103 nx2 = (int)(Ly/rootdx2);
00104
00105
00106
00107
00108 if((gcd = ath_gcd(nx1,nx2)) < 10)
00109 ath_error("[shkset2d]: Greatest Common Divisor (nx1,nx2) = %d\n",gcd);
00110
00111 id = div(nx1,gcd);
00112 r1 = id.quot;
00113 if(id.rem != 0)
00114 ath_error("[shkset2d]: GCD failed, Remainder of %d / %d is %d\n",
00115 nx1,gcd,id.rem);
00116
00117 id = div(nx2,gcd);
00118 r2 = id.quot;
00119 if(id.rem != 0)
00120 ath_error("[shkset2d]: GCD failed, Remainder of %d / %d is %d\n",
00121 nx2,gcd,id.rem);
00122
00123 ath_pout(1,"The unit cell is (%d,1,%d) grid cells in size\n",r1,r2);
00124
00125
00126
00127 if(Lx == Ly){
00128 cos_a = sin_a = sqrt(0.5);
00129 }
00130 else{
00131 angle = atan((double)(Lx/Ly));
00132 sin_a = sin(angle);
00133 cos_a = cos(angle);
00134 }
00135
00136
00137
00138 Wl.d = par_getd("problem","dl");
00139 #ifdef ADIABATIC
00140 Wl.P = par_getd("problem","pl");
00141 #endif
00142 Wl.Vx = par_getd("problem","v1l");
00143 Wl.Vy = par_getd("problem","v2l");
00144 Wl.Vz = par_getd("problem","v3l");
00145 #ifdef MHD
00146 Bxl = par_getd("problem","b1l");
00147 Wl.By = par_getd("problem","b2l");
00148 Wl.Bz = par_getd("problem","b3l");
00149 #endif
00150
00151
00152
00153 Wr.d = par_getd("problem","dr");
00154 #ifdef ADIABATIC
00155 Wr.P = par_getd("problem","pr");
00156 #endif
00157 Wr.Vx = par_getd("problem","v1r");
00158 Wr.Vy = par_getd("problem","v2r");
00159 Wr.Vz = par_getd("problem","v3r");
00160 #ifdef MHD
00161 Bxr = par_getd("problem","b1r");
00162 Wr.By = par_getd("problem","b2r");
00163 Wr.Bz = par_getd("problem","b3r");
00164 if (Bxr != Bxl) ath_error(0,"[shkset2d] L/R values of Bx not the same\n");
00165 #endif
00166
00167 Ul = Prim1D_to_Cons1D(&Wl, &Bxl);
00168 Ur = Prim1D_to_Cons1D(&Wr, &Bxr);
00169
00170
00171 ql.d = Ul.d;
00172 ql.M1 = Ul.Mx*cos_a - Ul.My*sin_a;
00173 ql.M2 = Ul.Mx*sin_a + Ul.My*cos_a;
00174 ql.M3 = Ul.Mz;
00175 #ifdef MHD
00176 ql.B1c = Bxl*cos_a - Ul.By*sin_a;
00177 ql.B2c = Bxl*sin_a + Ul.By*cos_a;
00178 ql.B3c = Ul.Bz;
00179 #endif
00180 #ifdef ADIABATIC
00181 ql.E = Ul.E;
00182 #endif
00183
00184
00185 qr.d = Ur.d;
00186 qr.M1 = Ur.Mx*cos_a - Ur.My*sin_a;
00187 qr.M2 = Ur.Mx*sin_a + Ur.My*cos_a;
00188 qr.M3 = Ur.Mz;
00189 #ifdef MHD
00190 qr.B1c = Bxr*cos_a - Ur.By*sin_a;
00191 qr.B2c = Bxr*sin_a + Ur.By*cos_a;
00192 qr.B3c = Ur.Bz;
00193 #endif
00194 #ifdef ADIABATIC
00195 qr.E = Ur.E;
00196 #endif
00197
00198
00199
00200 for (k=kl; k<=ku; k++) {
00201 for (j=0; j<=je+nghost; j++) {
00202 ix2 = j + pGrid->Disp[1];
00203 for (i=0; i<=ie+nghost; i++) {
00204 ix1 = i + pGrid->Disp[0];
00205
00206
00207 if((drr = r2*(ix1) + r1*(ix2) - gcd*r1*r2) <= 0){
00208 pGrid->U[k][j][i] = ql;
00209 #ifdef MHD
00210 pGrid->B1i[k][j][i] = ql.B1c;
00211 pGrid->B2i[k][j][i] = ql.B2c;
00212 pGrid->B3i[k][j][i] = ql.B3c;
00213 #endif
00214 }
00215
00216 else if((dll = r2*(ix1-1) + r1*(ix2-1) - gcd*r1*r2) >= 0){
00217 pGrid->U[k][j][i] = qr;
00218 #ifdef MHD
00219 pGrid->B1i[k][j][i] = qr.B1c;
00220 pGrid->B2i[k][j][i] = qr.B2c;
00221 pGrid->B3i[k][j][i] = qr.B3c;
00222 #endif
00223 }
00224
00225 else{
00226 dlr = r2*(ix1-1) + r1*(ix2) - gcd*r1*r2;
00227
00228 if(dlr < 0){
00229 afl_lx = 1.0;
00230 afr_lx = 0.0;
00231 afl_ry = (Real)(-dlr)/(Real)(r2);
00232 afr_ry = 1.0 - afl_ry;
00233 }
00234 else if(dlr > 0){
00235 afl_lx = (Real)(-dll)/(Real)(r1);
00236 afr_lx = 1.0 - afl_lx;
00237 afl_ry = 0.0;
00238 afr_ry = 1.0;
00239 }
00240 else{
00241 afl_lx = 1.0;
00242 afr_lx = 0.0;
00243 afl_ry = 0.0;
00244 afr_ry = 1.0;
00245 }
00246
00247 drl = r2*(ix1) + r1*(ix2-1) - gcd*r1*r2;
00248
00249 if(drl < 0){
00250 afl_rx = (Real)(-drl)/(Real)(r1);
00251 afr_rx = 1.0 - afl_rx;
00252 afl_ly = 1.0;
00253 afr_ly = 0.0;
00254 }
00255 else if(drl > 0){
00256 afl_rx = 0.0;
00257 afr_rx = 1.0;
00258 afl_ly = (Real)(-dll)/(Real)(r2);
00259 afr_ly = 1.0 - afl_ly;
00260 }
00261 else{
00262 afl_rx = 0.0;
00263 afr_rx = 1.0;
00264 afl_ly = 1.0;
00265 afr_ly = 0.0;
00266 }
00267
00268
00269 if(dlr > 0 && drl < 0){
00270 vfl = 0.5*(afl_lx + afl_rx);
00271 vfr = 1.0 - vfl;
00272 }
00273
00274 else if(dlr < 0 && drl > 0){
00275 vfl = 0.5*(afl_ly + afl_ry);
00276 vfr = 1.0 - vfl;
00277 }
00278
00279 else if(dlr == 0 && drl == 0){
00280 vfl = vfr = 0.5;
00281 }
00282
00283 else if(dlr > 0 && drl > 0){
00284 vfl = 0.5*afl_lx*afl_ly;
00285 vfr = 1.0 - vfl;
00286 }
00287
00288 else{
00289 vfr = 0.5*afr_rx*afr_ry;
00290 vfl = 1.0 - vfr;
00291 }
00292
00293
00294 #ifdef MHD
00295 pGrid->B1i[k][j][i] = afl_lx*ql.B1c + afr_lx*qr.B1c;
00296 B1r = afl_rx*ql.B1c + afr_rx*qr.B1c;
00297
00298 pGrid->B2i[k][j][i] = afl_ly*ql.B2c + afr_ly*qr.B2c;
00299 B2r = afl_ry*ql.B2c + afr_ry*qr.B2c;
00300
00301 pGrid->B3i[k][j][i] = vfl*ql.B3c + vfr*qr.B3c;
00302 #endif
00303
00304
00305 pGrid->U[k][j][i].d = vfl*ql.d + vfr*qr.d;
00306 pGrid->U[k][j][i].M1 = vfl*ql.M1 + vfr*qr.M1;
00307 pGrid->U[k][j][i].M2 = vfl*ql.M2 + vfr*qr.M2;
00308 pGrid->U[k][j][i].M3 = vfl*ql.M3 + vfr*qr.M3;
00309 #ifdef MHD
00310 pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i] + B1r);
00311 pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i] + B2r);
00312 pGrid->U[k][j][i].B3c = vfl*ql.B3c + vfr*qr.B3c;
00313 #endif
00314 #ifndef ISOTHERMAL
00315 pGrid->U[k][j][i].E = vfl*ql.E + vfr*qr.E;
00316 #endif
00317 }
00318 }
00319 }
00320 }
00321
00322
00323
00324 bvals_mhd_fun(pDomain, left_x1, shkset2d_iib);
00325 bvals_mhd_fun(pDomain, left_x2, shkset2d_ijb);
00326 bvals_mhd_fun(pDomain, right_x1, shkset2d_oib);
00327 bvals_mhd_fun(pDomain, right_x2, shkset2d_ojb);
00328
00329 return;
00330 }
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343 void problem_write_restart(MeshS *pM, FILE *fp)
00344 {
00345 return;
00346 }
00347
00348 void problem_read_restart(MeshS *pM, FILE *fp)
00349 {
00350 return;
00351 }
00352
00353 ConsFun_t get_usr_expr(const char *expr)
00354 {
00355 return NULL;
00356 }
00357
00358 VOutFun_t get_usr_out_fun(const char *name){
00359 return NULL;
00360 }
00361
00362 void Userwork_in_loop(MeshS *pM)
00363 {
00364 }
00365
00366 void Userwork_after_loop(MeshS *pM)
00367 {
00368 }
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 void shkset2d_iib(GridS *pGrid)
00379 {
00380 const int is = pGrid->is;
00381 int i, j, k, ju, jl, kl, ku;
00382
00383 if (pGrid->Nx[1] > 1){
00384 ju = pGrid->je + nghost;
00385 jl = pGrid->js - nghost + r2;
00386 } else {
00387 ju = pGrid->je;
00388 jl = pGrid->js;
00389 }
00390
00391 if (pGrid->Nx[2] > 1){
00392 ku = pGrid->ke + nghost;
00393 kl = pGrid->ks - nghost;
00394 } else {
00395 ku = pGrid->ke;
00396 kl = pGrid->ks;
00397 }
00398
00399 for (k=kl; k<=ku; k++) {
00400 for (i=1; i<=nghost; i++) {
00401 for (j=jl; j<=ju; j++) {
00402 pGrid->U [k][j][is-i] = pGrid->U [k][j-r2][is-i+r1];
00403 #ifdef MHD
00404 pGrid->B1i[k][j][is-i] = pGrid->B1i[k][j-r2][is-i+r1];
00405 pGrid->B2i[k][j][is-i] = pGrid->B2i[k][j-r2][is-i+r1];
00406 pGrid->B3i[k][j][is-i] = pGrid->B3i[k][j-r2][is-i+r1];
00407 #endif
00408 }
00409 }
00410 }
00411 return;
00412 }
00413
00414
00415
00416
00417
00418
00419
00420 void shkset2d_oib(GridS *pGrid)
00421 {
00422 const int ie = pGrid->ie;
00423 int i, j, k, ju, jl, kl, ku;
00424
00425 if (pGrid->Nx[1] > 1){
00426 ju = pGrid->je + nghost - r2;
00427 jl = pGrid->js - nghost;
00428 } else {
00429 ju = pGrid->je;
00430 jl = pGrid->js;
00431 }
00432
00433 if (pGrid->Nx[2] > 1){
00434 ku = pGrid->ke + nghost;
00435 kl = pGrid->ks - nghost;
00436 } else {
00437 ku = pGrid->ke;
00438 kl = pGrid->ks;
00439 }
00440
00441
00442
00443 for (k=kl; k<=ku; k++) {
00444 for (i=1; i<=nghost; i++) {
00445 for (j=jl; j<=ju; j++) {
00446 pGrid->U[k][j][ie+i] = pGrid->U[k][j+r2][ie+i-r1];
00447 #ifdef MHD
00448 if(i>1) pGrid->B1i[k][j][ie+i] = pGrid->B1i[k][j+r2][ie+i-r1];
00449 pGrid->B2i[k][j][ie+i] = pGrid->B2i[k][j+r2][ie+i-r1];
00450 pGrid->B3i[k][j][ie+i] = pGrid->B3i[k][j+r2][ie+i-r1];
00451 #endif
00452 }
00453 }
00454 }
00455 return;
00456 }
00457
00458
00459
00460
00461
00462
00463
00464 void shkset2d_ijb(GridS *pGrid)
00465 {
00466 const int js = pGrid->js;
00467 int i, j, k, iu, il, kl, ku;
00468
00469 iu = pGrid->ie + nghost;
00470 il = pGrid->is - nghost + r1;
00471
00472 if (pGrid->Nx[2] > 1){
00473 ku = pGrid->ke + nghost;
00474 kl = pGrid->ks - nghost;
00475 } else {
00476 ku = pGrid->ke;
00477 kl = pGrid->ks;
00478 }
00479
00480 for (k=kl; k<=ku; k++) {
00481 for (j=1; j<=nghost; j++) {
00482 for (i=il; i<=iu; i++) {
00483 pGrid->U [k][js-j][i] = pGrid->U [k][js-j+r2][i-r1];
00484 #ifdef MHD
00485 pGrid->B1i[k][js-j][i] = pGrid->B1i[k][js-j+r2][i-r1];
00486 pGrid->B2i[k][js-j][i] = pGrid->B2i[k][js-j+r2][i-r1];
00487 pGrid->B3i[k][js-j][i] = pGrid->B3i[k][js-j+r2][i-r1];
00488 #endif
00489 }
00490 }
00491 }
00492 return;
00493 }
00494
00495
00496
00497
00498
00499
00500
00501 void shkset2d_ojb(GridS *pGrid)
00502 {
00503 const int je = pGrid->je;
00504 int i, j, k, iu, il, kl, ku;
00505
00506 iu = pGrid->ie + nghost - r1;
00507 il = pGrid->is - nghost;
00508
00509 if (pGrid->Nx[2] > 1){
00510 ku = pGrid->ke + nghost;
00511 kl = pGrid->ks - nghost;
00512 } else {
00513 ku = pGrid->ke;
00514 kl = pGrid->ks;
00515 }
00516
00517
00518
00519 for (k=kl; k<=ku; k++) {
00520 for (j=1; j<=nghost; j++) {
00521 for (i=il; i<=iu; i++) {
00522 pGrid->U[k][je+j][i] = pGrid->U[k][je+j-r2][i+r1];
00523 #ifdef MHD
00524 pGrid->B1i[k][je+j][i] = pGrid->B1i[k][je+j-r2][i+r1];
00525 if(j>1) pGrid->B2i[k][je+j][i] = pGrid->B2i[k][je+j-r2][i+r1];
00526 pGrid->B3i[k][je+j][i] = pGrid->B3i[k][je+j-r2][i+r1];
00527 #endif
00528 }
00529 }
00530 }
00531 return;
00532 }