53 int coords[3], dst1, dst2;
54 long int count, buf_size;
55 static double *qL, *qR;
58 print1 (
"! SB_SetBoundaryVar: wrong boundary\n");
77 if (grid[
JDIR].nproc > 1) SB_ShiftBoundaryVar(U, box, side, t, grid);
88 for (k = box->
kb; k <= box->ke; k++){
89 for (i = box->
ib; i <= box->ie; i++){
96 for (k = box->
kb; k <= box->ke; k++){
97 for (i = box->
ib; i <= box->ie; i++){
113 void SB_ShiftBoundaryVar(
double ***
q,
RBox *box,
int side,
double t,
Grid *grid)
131 int i,
j,
k, ngh_x, ngh_y;
132 int nx_buf1, ny_buf1, nz_buf1;
133 int nx_buf2, ny_buf2, nz_buf2;
134 long int count, buf1_size, buf2_size;
137 int dst1, dst2, src1, src2;
138 int Delta_j, shift_coordy_mod, shift_coordy_div;
141 static double *snd_buf, *rcv_buf;
146 static MPI_Comm cartcomm;
164 if (snd_buf == NULL){
166 #if defined STAGGERED_MHD && DIMENSIONS == 3
177 buf1.
ib = buf2.
ib = box->
ib;
178 buf1.
ie = buf2.
ie = box->
ie;
179 buf1.
kb = buf2.
kb = box->
kb;
180 buf1.
ke = buf2.
ke = box->
ke;
201 shift_coordy_mod = Delta_j%
NX2;
202 shift_coordy_div = Delta_j/
NX2;
206 ny_buf1 =
NX2 - (shift_coordy_mod) + ngh_y;
207 ny_buf2 = ngh_y + (shift_coordy_mod);
209 buf1.
jb = 0; buf1.
je = ny_buf1 - 1;
210 buf2.
jb = 0; buf2.
je = ny_buf2 - 1;
212 for (i = 0; i <
DIMENSIONS; i++) coords[i] = grid[i].rank_coord;
213 coords[
JDIR] += shift_coordy_div;
214 MPI_Cart_rank(cartcomm, coords, &dst1);
217 MPI_Cart_rank(cartcomm, coords, &dst2);
219 for (i = 0; i <
DIMENSIONS; i++) coords[i] = grid[i].rank_coord;
220 coords[
JDIR] -= shift_coordy_div;
221 MPI_Cart_rank(cartcomm, coords, &src1);
224 MPI_Cart_rank(cartcomm, coords, &src2);
226 }
else if (side ==
X1_END){
228 ny_buf1 = shift_coordy_mod + ngh_y;
229 ny_buf2 = ngh_y + (
NX2 - shift_coordy_mod);
231 buf1.
jb = 0; buf1.
je = ny_buf1 - 1;
232 buf2.
jb = 0; buf2.
je = ny_buf2 - 1;
234 for (i = 0; i <
DIMENSIONS; i++) coords[i] = grid[i].rank_coord;
235 coords[
JDIR] -= shift_coordy_div;
236 MPI_Cart_rank(cartcomm, coords, &dst2);
239 MPI_Cart_rank(cartcomm, coords, &dst1);
241 for (i = 0; i <
DIMENSIONS; i++) coords[i] = grid[i].rank_coord;
242 coords[
JDIR] += shift_coordy_div;
243 MPI_Cart_rank(cartcomm, coords, &src2);
246 MPI_Cart_rank(cartcomm, coords, &src1);
249 print1 (
"! SB_ShiftBoundaryVar: wrong boundary\n");
255 nx_buf1 = buf1.
ie - buf1.
ib + 1;
256 nx_buf2 = buf2.
ie - buf2.
ib + 1;
258 ny_buf1 = buf1.
je - buf1.
jb + 1;
259 ny_buf2 = buf2.
je - buf2.
jb + 1;
261 nz_buf1 = buf1.
ke - buf1.
kb + 1;
262 nz_buf2 = buf2.
ke - buf2.
kb + 1;
266 buf1_size = nz_buf1*ny_buf1*nx_buf1;
267 buf2_size = nz_buf2*ny_buf2*nx_buf2;
275 BOX_LOOP(pbuf1, k, j, i) snd_buf[count++] = q[k][
JBEG+j][i];
276 BOX_LOOP(pbuf2, k, j, i) snd_buf[count++] = q[k][
JEND-ny_buf2+1+j][i];
280 MPI_Sendrecv(snd_buf, buf1_size, MPI_DOUBLE, dst1, 1,
281 rcv_buf, buf1_size, MPI_DOUBLE, src1, 1,
284 MPI_Sendrecv(snd_buf + buf1_size, buf2_size, MPI_DOUBLE, dst2, 2,
285 rcv_buf + buf1_size, buf2_size, MPI_DOUBLE, src2, 2,
291 BOX_LOOP(pbuf1, k, j, i) q[k][j+ny_buf2][i] = rcv_buf[count++];
292 BOX_LOOP(pbuf2, k, j, i) q[k][j][i] = rcv_buf[count++];
298 int nghL,
int nghR,
Grid *grid)
319 long int count, buf_size;
320 static double *snd_buf1, *snd_buf2, *rcv_buf1, *rcv_buf2;
321 static int dst1, dst2;
322 static MPI_Comm cartcomm;
337 if (grid[
JDIR].nproc == 1){
346 box->
jb = jb0; box->
je = je0;
360 if (snd_buf1 == NULL){
376 for (i = 0; i <
DIMENSIONS; i++) coords[i] = grid[i].rank_coord;
378 MPI_Cart_rank(cartcomm, coords, &dst1);
381 MPI_Cart_rank(cartcomm, coords, &dst2);
400 box->
je =
JBEG + nghR - 1;
401 buf_size = (box->
ke - box->
kb + 1)
402 *(box->
je - box->
jb + 1)
403 *(box->
ie - box->
ib + 1);
404 BOX_LOOP(box, k, j, i) snd_buf1[count++] = U[
k][
j][
i];
406 MPI_Sendrecv(snd_buf1, buf_size, MPI_DOUBLE, dst1, 1,
407 rcv_buf2, buf_size, MPI_DOUBLE, dst2, 1,
413 box->
jb =
JEND - nghL + 1;
415 buf_size = (box->
ke - box->
kb + 1)
416 *(box->
je - box->
jb + 1)
417 *(box->
ie - box->
ib + 1);
418 BOX_LOOP(box, k, j, i) snd_buf2[count++] = U[
k][
j][
i];
420 MPI_Sendrecv(snd_buf2, buf_size, MPI_DOUBLE, dst2, 2,
421 rcv_buf1, buf_size, MPI_DOUBLE, dst1, 2,
429 BOX_LOOP(box, k, j, i) U[
k][
j][
i] = rcv_buf2[count++];
434 BOX_LOOP(box, k, j, i) U[
k][
j][
i] = rcv_buf1[count++];
438 box->
jb = jb0; box->
je = je0;
463 int j, Delta_j, nghost;
466 double Delta_L, Delta_y, dyfrac;
468 double *q_from, *q_to;
469 static double *dq, *dql, *flux;
484 }
else if (side ==
X1_END){
491 scrh = Delta_L/Delta_y;
496 dyfrac = fabs(scrh) - floor(fabs(scrh));
501 if (grid[
JDIR].nproc > 1) Delta_j = 0;
506 for (j = 1; j <=
JEND + nghost; j++) {
507 dq[
j] = q_from[
j] - q_from[j - 1];
512 for (j = 0; j <=
JEND + nghost; j++) dql[j] = 0.0;
515 for (j =
JBEG - 1; j <=
JEND; j++) flux[j] = q_from[j];
517 for (j =
JBEG - 1; j <=
JEND; j++) flux[j] = q_from[j + 1];
527 for (j =
JBEG - 1; j <=
JEND; j++) {
528 flux[
j] = q_from[
j] + 0.5*dql[
j]*(1.0 - dyfrac);
530 }
else if (side ==
X1_END){
531 for (j =
JBEG - 1; j <=
JEND; j++) {
532 flux[
j] = q_from[j + 1] - 0.5*dql[j + 1]*(1.0 - dyfrac);
540 static double *qp, *qm;
541 double ap, am, P0, P1, P2;
546 if ( (
JBEG - 3) < 0) {
547 printf (
"! SB_ShearingInterp: the value of SB_ORDER (%d)\n",
549 printf (
"! is not consistent with underlying algorithm\n");
555 for (j =
JBEG - 2; j <=
JEND + 2; j++) dql[j] =
MC(dq[j+1], dq[j]);
556 for (j =
JBEG - 1; j <=
JEND + 1; j++){
557 qp[
j] = 0.5*(q_from[
j] + q_from[j+1]) - (dql[j+1] - dql[j])/6.0;
558 qm[
j] = 0.5*(q_from[
j] + q_from[j-1]) - (dql[j] - dql[j-1])/6.0;
560 ap = qp[
j] - q_from[
j];
561 am = qm[
j] - q_from[
j];
566 if (fabs(ap) >= 2.0*fabs(am)) ap = -2.0*am;
567 if (fabs(am) >= 2.0*fabs(ap)) am = -2.0*ap;
569 qp[
j] = q_from[
j] + ap;
570 qm[
j] = q_from[
j] + am;
575 P0 = 1.5*q_from[
j] - 0.25*(qp[
j] + qm[
j]);
577 P2 = 3.0*(qp[
j] - 2.0*q_from[
j] + qm[
j]);
578 flux[
j] = P0 + 0.5*P1*(1.0 - dyfrac) +
579 P2*(0.25 - 0.5*dyfrac + dyfrac*dyfrac/3.0);
581 }
else if (side ==
X1_END) {
583 P0 = 1.5*q_from[j + 1] - 0.25*(qp[j + 1] + qm[j + 1]);
584 P1 = qp[j + 1] - qm[j + 1];
585 P2 = 3.0*(qp[j + 1] - 2.0*q_from[j + 1] + qm[j + 1]);
586 flux[
j] = P0 - 0.5*P1*(1.0 - dyfrac) +
587 P2*(0.25 - 0.5*dyfrac + dyfrac*dyfrac/3.0);
593 print (
"! SB_ORDER should lie between 1 and 3\n");
621 if (side ==
X1_END) dyfrac *= -1.0;
623 if (grid[
JDIR].nproc == 1){
627 q_to[
j] = q_from[jp] - dyfrac*(flux[jp] - flux[jm]);
633 q_to[
j] = q_from[jp] - dyfrac*(flux[jp] - flux[jm]);
650 if (jj < JBEG || jj >
JEND){
651 print (
" ! j index out of range in JSHIFT %d, %d\n", jc, jj);
#define X1_BEG
Boundary region at X1 beg.
int jb
Lower corner index in the x2 direction.
void print1(const char *fmt,...)
#define BOX_LOOP(B, k, j, i)
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
int kb
Lower corner index in the x3 direction.
#define X1_END
Boundary region at X1 end.
int ib
Lower corner index in the x1 direction.
int AL_Get_cart_comm(int, MPI_Comm *)
long int NX2
Number of interior zones in the X2 directions (boundaries excluded) for the local processor...
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
int nghost
Number of ghost zones.
double sb_vy
Velocity offset (>0), in SB_Boundary().
void print(const char *fmt,...)
double g_domBeg[3]
Lower limits of the computational domain.
#define ARRAY_1D(nx, type)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
int ie
Upper corner index in the x1 direction.
int ke
Upper corner index in the x3 direction.
int je
Upper corner index in the x2 direction.
double g_domEnd[3]
Upper limits of the computational domain.
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
#define QUIT_PLUTO(e_code)
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...