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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080 #include <math.h>
00081 #include <stdio.h>
00082 #include <stdlib.h>
00083 #include "defs.h"
00084 #include "athena.h"
00085 #include "globals.h"
00086 #include "prototypes.h"
00087
00088 #ifdef STATIC_MESH_REFINEMENT
00089 #error The shkset3d test does not work with SMR.
00090 #endif
00091
00092
00093
00094
00095 static int rx, ry, rz;
00096 static Real ang_2, ang_3;
00097 static Real sin_a2, cos_a2, sin_a3, cos_a3;
00098
00099 static Real dl, vxl, vyl, vzl;
00100 static Real dr, vxr, vyr, vzr;
00101 #ifdef MHD
00102 static Real Bxl, Byl, Bzl, Bxr, Byr, Bzr;
00103 #endif
00104 #ifdef ADIABATIC
00105 static Real Pl, Pr;
00106 #endif
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121 static void lx_bc(GridS *pG);
00122 static void rx_bc(GridS *pG);
00123 static void ly_bc(GridS *pG);
00124 static void ry_bc(GridS *pG);
00125 static void lz_bc(GridS *pG);
00126 static void rz_bc(GridS *pG);
00127 static Real Ax(const Real x1, const Real x2, const Real x3);
00128 static Real Ay(const Real x1, const Real x2, const Real x3);
00129 static Real Az(const Real x1, const Real x2, const Real x3);
00130
00131
00132
00133
00134
00135 void problem(DomainS *pDomain)
00136 {
00137 GridS *pGrid = pDomain->Grid;
00138 int n;
00139 int i, j, k, ix, jx, kx, mix, mjx, mkx;
00140
00141 int Nx1, Nx2, Nx3;
00142 int sx, sy, sz;
00143 div_t id;
00144 double d_ix;
00145 Real xl, x, xr, dt_usr;
00146 Real lx1,lx2,lx3,cx1,cx2,cx3,rx1,rx2,rx3;
00147 Real sp0_x1, sp0_x2, sp0_x3;
00148
00149
00150
00151 int N, isqa, jsqa, ksqa;
00152 int scx, scy, scz;
00153 Real sdx, sdy, sdz;
00154 ConsS *pq, ***sqa=NULL;
00155 #ifdef MHD
00156 Real ***sB1i=NULL, ***sB2i=NULL, ***sB3i=NULL;
00157 #endif
00158
00159
00160 ConsS ql, qr;
00161
00162 ConsS ***qa=NULL;
00163 #ifdef MHD
00164 Real ***aB1i=NULL, ***aB2i=NULL, ***aB3i=NULL;
00165 #endif
00166 int qa_min_ix, qa_min_jx, qa_min_kx;
00167 int qa_max_ix, qa_max_jx, qa_max_kx;
00168
00169
00170
00171
00172
00173
00174
00175 Nx1 = pDomain->Nx[0];
00176 Nx2 = pDomain->Nx[1];
00177 Nx3 = pDomain->Nx[2];
00178 if(pDomain->Nx[2] <= 0)
00179 ath_error("[get_set_input]: This problem assumes a 3D domain\n");
00180
00181
00182
00183
00184
00185 rx = par_geti("problem","rx");
00186 ry = par_geti("problem","ry");
00187 rz = par_geti("problem","rz");
00188
00189 if(rx <= 0 || ry <= 0 || rz <= 0)
00190 ath_error("[get_set_input]: Only rx > 0, ry > 0 and rz > 0 is allowed\n");
00191
00192 id = div(Nx1,rx);
00193 sx = id.quot;
00194 if(id.rem != 0)
00195 ath_error("[get_set_input]: Choose Nx1 % rx = 0 instead of %d\n",Nx1%rx);
00196
00197 id = div(Nx2,ry);
00198 sy = id.quot;
00199 if(id.rem != 0)
00200 ath_error("[get_set_input]: Choose Nx2 % ry = 0 instead of %d\n",Nx2%ry);
00201
00202 id = div(Nx3,rz);
00203 sz = id.quot;
00204 if(id.rem != 0)
00205 ath_error("[get_set_input]: Choose Nx3 % rz = 0 instead of %d\n",Nx3%rz);
00206
00207 if(sy > sx || sz > sx)
00208 ath_error("Nx1/rx is assumed to be larger than the y- and z-dir. ratios\n");
00209
00210
00211
00212 ang_3 = atan((rx*pGrid->dx1)/(ry*pGrid->dx2));
00213 sin_a3 = sin(ang_3);
00214 cos_a3 = cos(ang_3);
00215
00216 ang_2 = atan((rx*pGrid->dx1*cos_a3)/(rz*pGrid->dx3));
00217 sin_a2 = sin(ang_2);
00218 cos_a2 = cos(ang_2);
00219
00220
00221
00222 dl = par_getd("problem","dl");
00223 vxl = par_getd("problem","vxl");
00224 vyl = par_getd("problem","vyl");
00225 vzl = par_getd("problem","vzl");
00226 #ifdef MHD
00227 Bxl = par_getd("problem","Bxl");
00228 Byl = par_getd("problem","Byl");
00229 Bzl = par_getd("problem","Bzl");
00230 #endif
00231 #ifdef ADIABATIC
00232 Pl = par_getd("problem","pl");
00233 #endif
00234
00235 dr = par_getd("problem","dr");
00236 vxr = par_getd("problem","vxr");
00237 vyr = par_getd("problem","vyr");
00238 vzr = par_getd("problem","vzr");
00239 #ifdef MHD
00240 Bxr = par_getd("problem","Bxr");
00241 Byr = par_getd("problem","Byr");
00242 Bzr = par_getd("problem","Bzr");
00243 #endif
00244 #ifdef ADIABATIC
00245 Pr = par_getd("problem","pr");
00246 #endif
00247
00248
00249
00250
00251 ql.d = dl;
00252 ql.M1 = dl*(vxl*cos_a2*cos_a3 - vyl*sin_a3 - vzl*sin_a2*cos_a3);
00253 ql.M2 = dl*(vxl*cos_a2*sin_a3 + vyl*cos_a3 - vzl*sin_a2*sin_a3);
00254 ql.M3 = dl*(vxl*sin_a2 + vzl*cos_a2);
00255
00256 #ifdef MHD
00257 ql.B1c = (Bxl*cos_a2*cos_a3 - Byl*sin_a3 - Bzl*sin_a2*cos_a3);
00258 ql.B2c = (Bxl*cos_a2*sin_a3 + Byl*cos_a3 - Bzl*sin_a2*sin_a3);
00259 ql.B3c = (Bxl*sin_a2 + Bzl*cos_a2);
00260 #endif
00261
00262 #ifdef ADIABATIC
00263 ql.E = Pl/Gamma_1
00264 #ifdef MHD
00265 + 0.5*(Bxl*Bxl + Byl*Byl + Bzl*Bzl)
00266 #endif
00267 + 0.5*dl*(vxl*vxl + vyl*vyl + vzl*vzl);
00268 #endif
00269
00270
00271 qr.d = dr;
00272 qr.M1 = dr*(vxr*cos_a2*cos_a3 - vyr*sin_a3 - vzr*sin_a2*cos_a3);
00273 qr.M2 = dr*(vxr*cos_a2*sin_a3 + vyr*cos_a3 - vzr*sin_a2*sin_a3);
00274 qr.M3 = dr*(vxr*sin_a2 + vzr*cos_a2);
00275
00276 #ifdef MHD
00277 qr.B1c = (Bxr*cos_a2*cos_a3 - Byr*sin_a3 - Bzr*sin_a2*cos_a3);
00278 qr.B2c = (Bxr*cos_a2*sin_a3 + Byr*cos_a3 - Bzr*sin_a2*sin_a3);
00279 qr.B3c = (Bxr*sin_a2 + Bzr*cos_a2);
00280 #endif
00281
00282 #ifdef ADIABATIC
00283 qr.E = Pr/Gamma_1
00284 #ifdef MHD
00285 + 0.5*(Bxr*Bxr + Byr*Byr + Bzr*Bzr)
00286 #endif
00287 + 0.5*dr*(vxr*vxr + vyr*vyr + vzr*vzr);
00288 #endif
00289
00290
00291
00292 if((qa = (ConsS***)calloc_3d_array(rz, ry, rx+rx, sizeof(ConsS)))==NULL)
00293 ath_error("[shkset3d]: Error allocating qa[%d][%d][%d]\n",rz,ry,rx+rx);
00294
00295 #ifdef MHD
00296 if((aB1i = (Real***)calloc_3d_array(rz, ry, rx+rx+1, sizeof(Real))) == NULL)
00297 ath_error("[shkset3d]: Error allocating aB1i[%d][%d][%d]\n",rz,ry,rx+rx+1);
00298
00299 if((aB2i = (Real***)calloc_3d_array(rz, ry+1, rx+rx, sizeof(Real))) == NULL)
00300 ath_error("[shkset3d]: Error allocating aB2i[%d][%d][%d]\n",rz,ry+1,rx+rx);
00301
00302 if((aB3i = (Real***)calloc_3d_array(rz+1, ry, rx+rx, sizeof(Real))) == NULL)
00303 ath_error("[shkset3d]: Error allocating aB3i[%d][%d][%d]\n",rz+1,ry,rx+rx);
00304 #endif
00305
00306
00307 qa_min_jx = 0; qa_max_jx = ry;
00308 qa_min_kx = 0; qa_max_kx = rz;
00309
00310
00311
00312
00313
00314 d_ix = -pGrid->MinX[0]/pGrid->dx1 - rx*(pGrid->MinX[1]/(ry*pGrid->dx2)
00315 + pGrid->MinX[2]/(rz*pGrid->dx3));
00316 qa_max_ix = ceil(d_ix);
00317 qa_min_ix = qa_max_ix - rx - rx;
00318
00319
00320
00321
00322
00323
00324 if(qa_max_ix - d_ix > 1.0e-12){
00325 ath_perr(-1,"[shkset3d]: qa_max_ix - d_ix = %g\n",
00326 qa_max_ix - d_ix);
00327 ath_error("[shkset3d]: Try setting the x2min and x3min = 0.0\n");
00328 }
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 N = rx*ry*rz;
00339 scx = ry*rz;
00340 scy = rx*rz;
00341 scz = rx*ry;
00342
00343 sdx = pGrid->dx1/scx;
00344 sdy = pGrid->dx2/scy;
00345 sdz = pGrid->dx3/scz;
00346
00347 if((sqa = (ConsS***)calloc_3d_array(N, N, N+N, sizeof(ConsS))) == NULL)
00348 ath_error("[shkset3d]: Error allocating sqa[%d][%d][%d]\n",N,N,N+N);
00349
00350 #ifdef MHD
00351 if((sB1i = (Real***)calloc_3d_array(N, N, N+N+1, sizeof(Real))) == NULL)
00352 ath_error("[shkset3d]: Error allocating sB1i[%d][%d][%d]\n",N,N,N+N+1);
00353
00354 if((sB2i = (Real***)calloc_3d_array(N, N+1, N+N, sizeof(Real))) == NULL)
00355 ath_error("[shkset3d]: Error allocating sB2i[%d][%d][%d]\n",N,N+1,N+N);
00356
00357 if((sB3i = (Real***)calloc_3d_array(N+1, N, N+N, sizeof(Real))) == NULL)
00358 ath_error("[shkset3d]: Error allocating sB3i[%d][%d][%d]\n",N+1,N,N+N);
00359 #endif
00360
00361
00362 sp0_x1 = pGrid->MinX[0] + qa_min_ix*pGrid->dx1;
00363 sp0_x2 = pGrid->MinX[1] + qa_min_jx*pGrid->dx2;
00364 sp0_x3 = pGrid->MinX[2] + qa_min_kx*pGrid->dx3;
00365
00366
00367
00368
00369
00370 #ifdef MHD
00371 for(k=0; k<=N; k++){
00372 for(j=0; j<=N; j++){
00373 for(i=0; i<=(N+N); i++){
00374
00375 rx1 = sp0_x1 + (i+1)*sdx;
00376 rx2 = sp0_x2 + (j+1)*sdy;
00377 rx3 = sp0_x3 + (k+1)*sdz;
00378
00379
00380 lx1 = sp0_x1 + i*sdx;
00381 lx2 = sp0_x2 + j*sdy;
00382 lx3 = sp0_x3 + k*sdz;
00383
00384
00385 cx1 = lx1 + 0.5*sdx;
00386 cx2 = lx2 + 0.5*sdy;
00387 cx3 = lx3 + 0.5*sdz;
00388
00389
00390 if(j<N && k<N){
00391
00392 xl = lx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + lx3*sin_a2;
00393 xr = lx1*cos_a2*cos_a3 + rx2*cos_a2*sin_a3 + rx3*sin_a2;
00394
00395 if(xl >= 0.0){
00396 sB1i[k][j][i] = qr.B1c;
00397 }
00398 else if(xr <= 0.0){
00399 sB1i[k][j][i] = ql.B1c;
00400 }
00401 else{
00402 sB1i[k][j][i] = Bxl*cos_a2*cos_a3 +
00403 (Az(lx1, rx2, cx3) - Az(lx1, lx2, cx3) )/sdy -
00404 (Ay(lx1, cx2, rx3) - Ay(lx1, cx2, lx3) )/sdz;
00405 }
00406 }
00407
00408
00409 if(i<(N+N) && k<N){
00410
00411 xl = lx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + lx3*sin_a2;
00412 xr = rx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + rx3*sin_a2;
00413
00414 if(xl >= 0.0){
00415 sB2i[k][j][i] = qr.B2c;
00416 }
00417 else if(xr <= 0.0){
00418 sB2i[k][j][i] = ql.B2c;
00419 }
00420 else{
00421 sB2i[k][j][i] = Bxl*cos_a2*sin_a3 +
00422 (Ax(cx1, lx2, rx3) - Ax(cx1, lx2, lx3) )/sdz -
00423 (Az(rx1, lx2, cx3) - Az(lx1, lx2, cx3) )/sdx;
00424 }
00425 }
00426
00427
00428 if(i<(N+N) && j<N){
00429
00430 xl = lx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + lx3*sin_a2;
00431 xr = rx1*cos_a2*cos_a3 + rx2*cos_a2*sin_a3 + lx3*sin_a2;
00432
00433 if(xl >= 0.0){
00434 sB3i[k][j][i] = qr.B3c;
00435 }
00436 else if(xr <= 0.0){
00437 sB3i[k][j][i] = ql.B3c;
00438 }
00439 else{
00440 sB3i[k][j][i] = Bxl*sin_a2 +
00441 (Ay(rx1, cx2, lx3) - Ay(lx1, cx2, lx3) )/sdx -
00442 (Ax(cx1, rx2, lx3) - Ax(cx1, lx2, lx3) )/sdy;
00443 }
00444 }
00445 }
00446 }
00447 }
00448 #endif
00449
00450
00451
00452
00453
00454
00455 for(k=0; k<N; k++){
00456 for(j=0; j<N; j++){
00457 for(i=0; i<(N+N); i++){
00458
00459 rx1 = sp0_x1 + (i+1)*sdx;
00460 rx2 = sp0_x2 + (j+1)*sdy;
00461 rx3 = sp0_x3 + (k+1)*sdz;
00462
00463 xr = rx1*cos_a2*cos_a3 + rx2*cos_a2*sin_a3 + rx3*sin_a2;
00464
00465
00466 lx1 = sp0_x1 + i*sdx;
00467 lx2 = sp0_x2 + j*sdy;
00468 lx3 = sp0_x3 + k*sdz;
00469
00470 xl = lx1*cos_a2*cos_a3 + lx2*cos_a2*sin_a3 + lx3*sin_a2;
00471
00472
00473 cx1 = lx1 + 0.5*sdx;
00474 cx2 = lx2 + 0.5*sdy;
00475 cx3 = lx3 + 0.5*sdz;
00476 x = cx1*cos_a2*cos_a3 + cx2*cos_a2*sin_a3 + cx3*sin_a2;
00477
00478 if(xr <= 0.0){
00479 sqa[k][j][i] = ql;
00480 continue;
00481 }
00482
00483 if(xl >= 0.0){
00484 sqa[k][j][i] = qr;
00485 continue;
00486 }
00487
00488 pq = &(sqa[k][j][i]);
00489
00490
00491 if(x < 0.0){
00492 pq->d = dl;
00493 pq->M1 = dl*(vxl*cos_a2*cos_a3 - vyl*sin_a3 - vzl*sin_a2*cos_a3);
00494 pq->M2 = dl*(vxl*cos_a2*sin_a3 + vyl*cos_a3 - vzl*sin_a2*sin_a3);
00495 pq->M3 = dl*(vxl*sin_a2 + vzl*cos_a2);
00496 #ifdef ADIABATIC
00497 pq->E = Pl/Gamma_1 + 0.5*dl*(vxl*vxl + vyl*vyl + vzl*vzl);
00498 #endif
00499 }
00500 else{
00501 pq->d = dr;
00502 pq->M1 = dr*(vxr*cos_a2*cos_a3 - vyr*sin_a3 - vzr*sin_a2*cos_a3);
00503 pq->M2 = dr*(vxr*cos_a2*sin_a3 + vyr*cos_a3 - vzr*sin_a2*sin_a3);
00504 pq->M3 = dr*(vxr*sin_a2 + vzr*cos_a2);
00505 #ifdef ADIABATIC
00506 pq->E = Pr/Gamma_1 + 0.5*dr*(vxr*vxr + vyr*vyr + vzr*vzr);
00507 #endif
00508 }
00509
00510
00511 #ifdef MHD
00512 pq->B1c = 0.5*(sB1i[k][j][i] + sB1i[k ][j ][i+1]);
00513 pq->B2c = 0.5*(sB2i[k][j][i] + sB2i[k ][j+1][i ]);
00514 pq->B3c = 0.5*(sB3i[k][j][i] + sB3i[k+1][j ][i ]);
00515
00516 #ifdef ADIABATIC
00517 pq->E += 0.5*(pq->B1c*pq->B1c + pq->B2c*pq->B2c + pq->B3c*pq->B3c);
00518 #endif
00519 #endif
00520 }
00521 }
00522 }
00523
00524
00525
00526
00527
00528
00529
00530 #ifdef MHD
00531 for(k=0; k<=rz; k++){
00532 for(j=0; j<=ry; j++){
00533 for(i=0; i<=(rx+rx); i++){
00534
00535
00536 if(j<ry && k<rz){
00537 aB1i[k][j][i] = 0.0;
00538 isqa = i*scx;
00539 for(ksqa=k*scz; ksqa<(k+1)*scz; ksqa++){
00540 for(jsqa=j*scy; jsqa<(j+1)*scy; jsqa++){
00541 aB1i[k][j][i] += sB1i[ksqa][jsqa][isqa];
00542 }
00543 }
00544 aB1i[k][j][i] /= (Real)(scy*scz);
00545 }
00546
00547
00548 if(i<(rx+rx) && k<rz){
00549 aB2i[k][j][i] = 0.0;
00550 jsqa = j*scy;
00551 for(ksqa=k*scz; ksqa<(k+1)*scz; ksqa++){
00552 for(isqa=i*scx; isqa<(i+1)*scx; isqa++){
00553 aB2i[k][j][i] += sB2i[ksqa][jsqa][isqa];
00554 }
00555 }
00556 aB2i[k][j][i] /= (Real)(scx*scz);
00557 }
00558
00559
00560 if(i<(rx+rx) && j<ry){
00561 aB3i[k][j][i] = 0.0;
00562 ksqa = k*scz;
00563 for(jsqa=j*scy; jsqa<(j+1)*scy; jsqa++){
00564 for(isqa=i*scx; isqa<(i+1)*scx; isqa++){
00565 aB3i[k][j][i] += sB3i[ksqa][jsqa][isqa];
00566 }
00567 }
00568 aB3i[k][j][i] /= (Real)(scx*scy);
00569 }
00570 }
00571 }
00572 }
00573 #endif
00574
00575
00576 for(k=0; k<rz; k++){
00577 for(j=0; j<ry; j++){
00578 for(i=0; i<(rx+rx); i++){
00579 qa[k][j][i].d = 0.0;
00580 qa[k][j][i].M1 = qa[k][j][i].M2 = qa[k][j][i].M3 = 0.0;
00581 #ifdef ADIABATIC
00582 qa[k][j][i].E = 0.0;
00583 #endif
00584
00585
00586
00587
00588 for(ksqa=k*scz; ksqa<(k+1)*scz; ksqa++){
00589 for(jsqa=j*scy; jsqa<(j+1)*scy; jsqa++){
00590 for(isqa=i*scx; isqa<(i+1)*scx; isqa++){
00591 qa[k][j][i].d += sqa[ksqa][jsqa][isqa].d;
00592 qa[k][j][i].M1 += sqa[ksqa][jsqa][isqa].M1;
00593 qa[k][j][i].M2 += sqa[ksqa][jsqa][isqa].M2;
00594 qa[k][j][i].M3 += sqa[ksqa][jsqa][isqa].M3;
00595 #ifdef ADIABATIC
00596 qa[k][j][i].E += sqa[ksqa][jsqa][isqa].E;
00597 #endif
00598 }
00599 }
00600 }
00601
00602
00603
00604 qa[k][j][i].d /= (Real)(N*N);
00605 qa[k][j][i].M1 /= (Real)(N*N);
00606 qa[k][j][i].M2 /= (Real)(N*N);
00607 qa[k][j][i].M3 /= (Real)(N*N);
00608 #ifdef ADIABATIC
00609 qa[k][j][i].E /= (Real)(N*N);
00610 #endif
00611
00612
00613 #ifdef MHD
00614 qa[k][j][i].B1c = 0.5*(aB1i[k][j][i] + aB1i[k ][j ][i+1]);
00615 qa[k][j][i].B2c = 0.5*(aB2i[k][j][i] + aB2i[k ][j+1][i ]);
00616 qa[k][j][i].B3c = 0.5*(aB3i[k][j][i] + aB3i[k+1][j ][i ]);
00617 #endif
00618 }
00619 }
00620 }
00621
00622
00623
00624
00625
00626 for (k=0; k<=(pGrid->ke)+nghost; k++) {
00627 for (j=0; j<=(pGrid->je)+nghost; j++) {
00628 for (i=0; i<=(pGrid->ie)+nghost; i++) {
00629
00630
00631
00632
00633
00634 mix = ix = (i-(pGrid->is)) + pGrid->Disp[0];
00635 mjx = jx = (j-(pGrid->js)) + pGrid->Disp[1];
00636 mkx = kx = (k-(pGrid->ks)) + pGrid->Disp[2];
00637
00638 n=0;
00639 while(mjx < 0){
00640 mjx += ry;
00641 n++;
00642 }
00643 if(n){
00644 mix -= n*rx;
00645 }
00646
00647 n=0;
00648 while(mjx >= ry){
00649 mjx -= ry;
00650 n++;
00651 }
00652 if(n){
00653 mix += n*rx;
00654 }
00655
00656 n=0;
00657 while(mkx < 0){
00658 mkx += rz;
00659 n++;
00660 }
00661 if(n){
00662 mix -= n*rx;
00663 }
00664
00665 n=0;
00666 while(mkx >= rz){
00667 mkx -= rz;
00668 n++;
00669 }
00670 if(n){
00671 mix += n*rx;
00672 }
00673
00674 if(mix < qa_min_ix){
00675 pGrid->U[k][j][i] = ql;
00676 #ifdef MHD
00677 pGrid->B1i[k][j][i] = ql.B1c;
00678 pGrid->B2i[k][j][i] = ql.B2c;
00679 pGrid->B3i[k][j][i] = ql.B3c;
00680 #endif
00681 }
00682 else if(mix >= qa_max_ix){
00683 pGrid->U[k][j][i] = qr;
00684 #ifdef MHD
00685 pGrid->B1i[k][j][i] = qr.B1c;
00686 pGrid->B2i[k][j][i] = qr.B2c;
00687 pGrid->B3i[k][j][i] = qr.B3c;
00688 #endif
00689 }
00690 else{
00691 pGrid->U[k][j][i] = qa[mkx][mjx][mix - qa_min_ix];
00692 #ifdef MHD
00693 pGrid->B1i[k][j][i] = aB1i[mkx][mjx][mix - qa_min_ix];
00694 pGrid->B2i[k][j][i] = aB2i[mkx][mjx][mix - qa_min_ix];
00695 pGrid->B3i[k][j][i] = aB3i[mkx][mjx][mix - qa_min_ix];
00696 #endif
00697 }
00698 }
00699 }
00700 }
00701
00702
00703
00704
00705
00706 bvals_mhd_fun(pDomain, left_x1, lx_bc);
00707 bvals_mhd_fun(pDomain, right_x1, rx_bc);
00708 bvals_mhd_fun(pDomain, left_x2, ly_bc);
00709 bvals_mhd_fun(pDomain, right_x2, ry_bc);
00710 bvals_mhd_fun(pDomain, left_x3, lz_bc);
00711 bvals_mhd_fun(pDomain, right_x3, rz_bc);
00712
00713 free_3d_array((void***)sqa ); sqa = NULL;
00714 #ifdef MHD
00715 free_3d_array((void***)sB1i); sB1i = NULL;
00716 free_3d_array((void***)sB2i); sB2i = NULL;
00717 free_3d_array((void***)sB3i); sB3i = NULL;
00718 #endif
00719
00720 return;
00721 }
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735 void problem_write_restart(MeshS *pM,FILE *fp)
00736 {
00737 return;
00738 }
00739
00740 void problem_read_restart(MeshS *pM, FILE *fp)
00741 {
00742 return;
00743 }
00744
00745 ConsFun_t get_usr_expr(const char *expr)
00746 {
00747 return NULL;
00748 }
00749
00750 VOutFun_t get_usr_out_fun(const char *name){
00751 return NULL;
00752 }
00753
00754 void Userwork_in_loop(MeshS *pM)
00755 {
00756 return;
00757 }
00758
00759 void Userwork_after_loop(MeshS *pM)
00760 {
00761 return;
00762 }
00763
00764
00765
00766
00767
00768
00769
00770
00771 static void lx_bc(GridS *pG)
00772 {
00773 int i, j, k, mi, mj, mk, is, js, je, ks, ke;
00774 is = pG->is;
00775 js = pG->js; je = pG->je;
00776 ks = pG->ks; ke = pG->ke;
00777
00778 for(i=is-1; i>=0; i--){
00779 for(k = 0; k <= ke+nghost; k++){
00780 for(j = 0; j <= je+nghost; j++){
00781
00782 if(k - rz >= ks){
00783 mi = i + rx;
00784 mj = j;
00785 mk = k - rz;
00786 }
00787 else if(j - ry >= js){
00788 mi = i + rx;
00789 mj = j - ry;
00790 mk = k;
00791 }
00792 else continue;
00793
00794 pG->U[k][j][i] = pG->U[mk][mj][mi];
00795 #ifdef MHD
00796 pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00797 pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00798 pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00799 #endif
00800 }
00801 }
00802 }
00803
00804 return;
00805 }
00806
00807
00808
00809
00810
00811
00812 static void rx_bc(GridS *pG)
00813 {
00814 int i, j, k, mi, mj, mk, ie, js, je, ks, ke;
00815 ie = pG->ie;
00816 js = pG->js; je = pG->je;
00817 ks = pG->ks; ke = pG->ke;
00818
00819 for(i=ie+1; i<=ie+nghost; i++){
00820 for(k = 0; k <= ke+nghost; k++){
00821 for(j = 0; j <= je+nghost; j++){
00822
00823 if(k + rz <= ke){
00824 mi = i - rx;
00825 mj = j;
00826 mk = k + rz;
00827 }
00828 else if(j + ry <= je){
00829 mi = i;
00830 mj = j;
00831 mk = k;
00832 }
00833 else continue;
00834
00835 pG->U[k][j][i] = pG->U[mk][mj][mi];
00836 #ifdef MHD
00837 if(i > ie+1) pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00838 pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00839 pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00840 #endif
00841 }
00842 }
00843 }
00844
00845 return;
00846 }
00847
00848
00849
00850
00851
00852
00853 static void ly_bc(GridS *pG)
00854 {
00855 int i, j, k, mi, mj, mk, is, ie, js, ks, ke;
00856 is = pG->is; ie = pG->ie;
00857 js = pG->js;
00858 ks = pG->ks; ke = pG->ke;
00859
00860 for(j=js-1; j>=0; j--){
00861 for(k = 0; k <= ke+nghost; k++){
00862 for(i = 0; i <= ie+nghost; i++){
00863
00864 if(i - rx >= is){
00865 mi = i - rx;
00866 mj = j + ry;
00867 mk = k;
00868 }
00869 else if(k - rz >= ks){
00870 mi = i;
00871 mj = j + ry;
00872 mk = k - rz;
00873 }
00874 else continue;
00875
00876 pG->U[k][j][i] = pG->U[mk][mj][mi];
00877 #ifdef MHD
00878 pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00879 pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00880 pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00881 #endif
00882 }
00883 }
00884 }
00885
00886 return;
00887 }
00888
00889
00890
00891
00892
00893
00894 static void ry_bc(GridS *pG)
00895 {
00896 int i, j, k, mi, mj, mk, is, ie, je, ks, ke;
00897 is = pG->is; ie = pG->ie;
00898 je = pG->je;
00899 ks = pG->ks; ke = pG->ke;
00900
00901 for(j=je+1; j<=je+nghost; j++){
00902 for(k = 0; k <= ke+nghost; k++){
00903 for(i = 0; i <= ie+nghost; i++){
00904
00905 if(i + rx <= ie){
00906 mi = i + rx;
00907 mj = j - ry;
00908 mk = k;
00909 }
00910 else if(k + rz <= ke){
00911 mi = i;
00912 mj = j - ry;
00913 mk = k + rz;
00914 }
00915 else continue;
00916
00917 pG->U[k][j][i] = pG->U[mk][mj][mi];
00918 #ifdef MHD
00919 pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00920 if(j > je+1) pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00921 pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00922 #endif
00923 }
00924 }
00925 }
00926
00927 return;
00928 }
00929
00930
00931
00932
00933
00934
00935 static void lz_bc(GridS *pG)
00936 {
00937 int i, j, k, mi, mj, mk, is, ie, js, je, ks;
00938 is = pG->is; ie = pG->ie;
00939 js = pG->js; je = pG->je;
00940 ks = pG->ks;
00941
00942 for(k=ks-1; k>=0; k--){
00943 for(j = 0; j <= je+nghost; j++){
00944 for(i = 0; i <= ie+nghost; i++){
00945
00946 if(i - rx >= is){
00947 mi = i - rx;
00948 mj = j;
00949 mk = k + rz;
00950 }
00951 else if(j - ry >= js){
00952 mi = i;
00953 mj = j - ry;
00954 mk = k + rz;
00955 }
00956 else continue;
00957
00958 pG->U[k][j][i] = pG->U[mk][mj][mi];
00959 #ifdef MHD
00960 pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
00961 pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
00962 pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
00963 #endif
00964 }
00965 }
00966 }
00967
00968 return;
00969 }
00970
00971
00972
00973
00974
00975
00976 static void rz_bc(GridS *pG)
00977 {
00978 int i, j, k, mi, mj, mk, is, ie, js, je, ke;
00979 is = pG->is; ie = pG->ie;
00980 js = pG->js; je = pG->je;
00981 ke = pG->ke;
00982
00983 for(k=ke+1; k<=ke+nghost; k++){
00984 for(j = 0; j <= je+nghost; j++){
00985 for(i = 0; i <= ie+nghost; i++){
00986
00987 if(i + rx <= ie){
00988 mi = i + rx;
00989 mj = j;
00990 mk = k - rz;
00991 }
00992 else if(j + ry <= je){
00993 mi = i;
00994 mj = j + ry;
00995 mk = k - rz;
00996 }
00997 else continue;
00998
00999 pG->U[k][j][i] = pG->U[mk][mj][mi];
01000 #ifdef MHD
01001 pG->B1i[k][j][i] = pG->B1i[mk][mj][mi];
01002 pG->B2i[k][j][i] = pG->B2i[mk][mj][mi];
01003 if(k > ke+1) pG->B3i[k][j][i] = pG->B3i[mk][mj][mi];
01004 #endif
01005 }
01006 }
01007 }
01008
01009 return;
01010 }
01011
01012
01013
01014
01015
01016 #ifdef MHD
01017
01018
01019 static Real Ax(const Real x, const Real y, const Real z)
01020 {
01021 Real x1 = x*cos_a2*cos_a3 + y*cos_a2*sin_a3 + z*sin_a2;
01022 Real A2 = x1*(x1 < 0.0 ? Bzl : Bzr);
01023 Real A3 = -x1*(x1 < 0.0 ? Byl : Byr);
01024
01025 return -A2*sin_a3 - A3*sin_a2*cos_a3;
01026 }
01027
01028
01029
01030 static Real Ay(const Real x, const Real y, const Real z)
01031 {
01032 Real x1 = x*cos_a2*cos_a3 + y*cos_a2*sin_a3 + z*sin_a2;
01033 Real A2 = x1*(x1 < 0.0 ? Bzl : Bzr);
01034 Real A3 = -x1*(x1 < 0.0 ? Byl : Byr);
01035
01036 return A2*cos_a3 - A3*sin_a2*sin_a3;
01037 }
01038
01039
01040
01041 static Real Az(const Real x, const Real y, const Real z)
01042 {
01043 Real x1 = x*cos_a2*cos_a3 + y*cos_a2*sin_a3 + z*sin_a2;
01044
01045 Real A3 = -x1*(x1 < 0.0 ? Byl : Byr);
01046
01047 return A3*cos_a2;
01048 }
01049 #endif