00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <math.h>
00021 #include <float.h>
00022 #include "../defs.h"
00023 #include "../athena.h"
00024 #include "../globals.h"
00025 #include "prototypes.h"
00026 #include "../prototypes.h"
00027
00028 #ifdef SELF_GRAVITY_USING_MULTIGRID
00029
00030 #ifdef STATIC_MESH_REFINEMENT
00031 #error self gravity with multigrid not yet implemented to work with SMR
00032 #endif
00033
00034 #ifdef MPI_PARALLEL
00035 #error self gravity with multigrid not yet implemented to work with MPI
00036 #endif
00037
00038 #ifdef MPI_PARALLEL
00039 static double *send_buf=NULL, *recv_buf=NULL;
00040 #endif
00041
00042
00043
00044
00045 typedef struct MGrid_s{
00046 Real ***rhs,***Phi;
00047 Real dx1,dx2,dx3;
00048 int Nx1,Nx2,Nx3;
00049 int is,ie;
00050 int js,je;
00051 int ks,ke;
00052 int rx1_id, lx1_id;
00053 int rx2_id, lx2_id;
00054 int rx3_id, lx3_id;
00055 }MGrid;
00056
00057
00058 Real ***error;
00059
00060
00061
00062
00063
00064
00065 void multig_3d(MGrid *pMG);
00066 void Jacobi(MGrid *pMG);
00067 void Restriction_3d(MGrid *pMG_fine, MGrid *pMG_coarse);
00068 void Prolongation_3d(MGrid *pMG_coarse, MGrid *pMG_fine);
00069
00070 #ifdef MPI_PARALLEL
00071 void set_mg_bvals(MGrid *pMG);
00072 void swap_mg_ix1(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq);
00073 void swap_mg_ox1(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq);
00074 void swap_mg_ix2(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq);
00075 void swap_mg_ox2(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq);
00076 void swap_mg_ix3(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq);
00077 void swap_mg_ox3(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq);
00078 #endif
00079
00080
00081
00082
00083
00084
00085
00086
00087 void selfg_multig_1d(DomainS *pD)
00088 {
00089 ath_error("[selfg_multig_1d] 1D multigrid not implemented yet\n");
00090 return;
00091 }
00092
00093
00094
00095
00096
00097
00098 void selfg_multig_2d(DomainS *pD)
00099 {
00100 ath_error("[selfg_multig_2d] 2D multigrid not implemented yet\n");
00101 return;
00102 }
00103
00104
00105
00106
00107
00108
00109
00110 void selfg_multig_3d(DomainS *pD)
00111 {
00112 GridS *pG = (pD->Grid);
00113 int i, is = pG->is, ie = pG->ie;
00114 int j, js = pG->js, je = pG->je;
00115 int k, ks = pG->ks, ke = pG->ke;
00116 int Nx1z, Nx2z, Nx3z;
00117 MGrid Root_grid;
00118 Real mass = 0.0, tmass, dVol, rad, x1, x2, x3;
00119 Real Grav_const = four_pi_G/(4.0*PI);
00120 #ifdef MPI_PARALLEL
00121 Real mpi_err;
00122 #endif
00123
00124
00125
00126 for (k=ks-nghost; k<=ke+nghost; k++){
00127 for (j=js-nghost; j<=je+nghost; j++){
00128 for (i=is-nghost; i<=ie+nghost; i++){
00129 pG->Phi_old[k][j][i] = pG->Phi[k][j][i];
00130 }
00131 }
00132 }
00133
00134
00135
00136 dVol = pG->dx1*pG->dx2*pG->dx3;
00137 for (k=ks; k<=ke; k++){
00138 for (j=js; j<=je; j++){
00139 for (i=is; i<=ie; i++){
00140 mass += pG->U[k][j][i].d*dVol;
00141 }
00142 }
00143 }
00144
00145 #ifdef MPI_PARALLEL
00146 mpi_err = MPI_Allreduce(&mass, &tmass,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
00147 if (mpi_err) ath_error("[selfg_multigrid]: MPI_Reduce returned err = %d\n",
00148 mpi_err);
00149 #else
00150 tmass = mass;
00151 #endif
00152
00153
00154
00155 if (pG->lx1_id < 0) {
00156 for (k=ks; k<=ke; k++) {
00157 for (j=js; j<=je; j++) {
00158 for (i=1; i<=nghost; i++){
00159 cc_pos(pG,is-i,j,k,&x1,&x2,&x3);
00160 rad = sqrt(x1*x1 + x2*x2 + x3*x3);
00161 pG->Phi[k][j][is-i] = -Grav_const*tmass/rad;
00162 }
00163 }
00164 }
00165 }
00166
00167 if (pG->rx1_id < 0) {
00168 for (k=ks; k<=ke; k++) {
00169 for (j=js; j<=je; j++) {
00170 for (i=1; i<=nghost; i++){
00171 cc_pos(pG,ie+i,j,k,&x1,&x2,&x3);
00172 rad = sqrt(x1*x1 + x2*x2 + x3*x3);
00173 pG->Phi[k][j][ie+i] = -Grav_const*tmass/rad;
00174 }
00175 }
00176 }
00177 }
00178
00179
00180
00181 if (pG->lx2_id < 0) {
00182 for (k=ks; k<=ke; k++){
00183 for (j=1; j<=nghost; j++){
00184 for (i=is-nghost; i<=ie+nghost; i++){
00185 cc_pos(pG,i,js-j,k,&x1,&x2,&x3);
00186 rad = sqrt(x1*x1 + x2*x2 + x3*x3);
00187 pG->Phi[k][js-j][i] = -Grav_const*tmass/rad;
00188 }
00189 }
00190 }
00191 }
00192
00193 if (pG->rx2_id < 0) {
00194 for (k=ks; k<=ke; k++){
00195 for (j=1; j<=nghost; j++){
00196 for (i=is-nghost; i<=ie+nghost; i++){
00197 cc_pos(pG,i,je+j,k,&x1,&x2,&x3);
00198 rad = sqrt(x1*x1 + x2*x2 + x3*x3);
00199 pG->Phi[k][je+j][i] = -Grav_const*tmass/rad;
00200 }
00201 }
00202 }
00203 }
00204
00205
00206
00207 if (pG->lx3_id < 0) {
00208 for (k=1; k<=nghost; k++){
00209 for (j=js-nghost; j<=je+nghost; j++){
00210 for (i=is-nghost; i<=ie+nghost; i++){
00211 cc_pos(pG,i,j,ks-k,&x1,&x2,&x3);
00212 rad = sqrt(x1*x1 + x2*x2 + x3*x3);
00213 pG->Phi[ks-k][j][i] = -Grav_const*tmass/rad;
00214 }
00215 }
00216 }
00217 }
00218
00219 if (pG->rx3_id < 0) {
00220 for (k=1; k<=nghost; k++){
00221 for (j=js-nghost; j<=je+nghost; j++){
00222 for (i=is-nghost; i<=ie+nghost; i++){
00223 cc_pos(pG,i,j,ke+k,&x1,&x2,&x3);
00224 rad = sqrt(x1*x1 + x2*x2 + x3*x3);
00225 pG->Phi[ke+k][j][i] = -Grav_const*tmass/rad;
00226 }
00227 }
00228 }
00229 }
00230
00231
00232
00233 Root_grid.Nx1 = pG->Nx[0];
00234 Root_grid.Nx2 = pG->Nx[1];
00235 Root_grid.Nx3 = pG->Nx[2];
00236 Root_grid.is = 1; Root_grid.ie = pG->Nx[0];
00237 Root_grid.js = 1; Root_grid.je = pG->Nx[1];
00238 Root_grid.ks = 1; Root_grid.ke = pG->Nx[2];
00239 Root_grid.dx1 = pG->dx1;
00240 Root_grid.dx2 = pG->dx2;
00241 Root_grid.dx3 = pG->dx3;
00242 Root_grid.rx1_id = pG->rx1_id; Root_grid.lx1_id = pG->lx1_id;
00243 Root_grid.rx2_id = pG->rx2_id; Root_grid.lx2_id = pG->lx2_id;
00244 Root_grid.rx3_id = pG->rx3_id; Root_grid.lx3_id = pG->lx3_id;
00245
00246
00247 Nx1z = pG->Nx[0] + 2;
00248 Nx2z = pG->Nx[1] + 2;
00249 Nx3z = pG->Nx[2] + 2;
00250 Root_grid.rhs = (Real ***) calloc_3d_array(Nx3z,Nx2z,Nx1z,sizeof(Real));
00251 Root_grid.Phi = (Real ***) calloc_3d_array(Nx3z,Nx2z,Nx1z,sizeof(Real));
00252 if (Root_grid.rhs == NULL) {
00253 ath_error("[selfg_by_multig_3d]: Error allocating memory\n");
00254 }
00255 if (Root_grid.Phi == NULL) {
00256 ath_error("[selfg_by_multig_3d]: Error allocating memory\n");
00257 }
00258
00259
00260 for (k=ks-1; k<=ke+1; k++){
00261 for (j=js-1; j<=je+1; j++){
00262 for (i=is-1; i<=ie+1; i++){
00263 Root_grid.rhs[k-ks+1][j-js+1][i-is+1] = four_pi_G*pG->U[k][j][i].d;
00264 Root_grid.Phi[k-ks+1][j-js+1][i-is+1] = pG->Phi[k][j][i];
00265 }
00266 }
00267 }
00268 #ifdef MPI_PARALLEL
00269 set_mg_bvals(&Root_grid);
00270 #endif
00271
00272
00273
00274 multig_3d(&Root_grid);
00275
00276
00277
00278
00279 for (k=ks; k<=ke; k++){
00280 for (j=js; j<=je; j++){
00281 for (i=is; i<=ie; i++){
00282 pG->Phi[k][j][i] = Root_grid.Phi[k-ks+1][j-js+1][i-is+1];
00283 }
00284 }
00285 }
00286
00287 free_3d_array(Root_grid.rhs);
00288 free_3d_array(Root_grid.Phi);
00289
00290 return;
00291 }
00292
00293
00294
00295
00296
00297
00298 void multig_3d(MGrid *pMG)
00299 {
00300 int i, is = pMG->is, ie = pMG->ie;
00301 int j, js = pMG->js, je = pMG->je;
00302 int k, ks = pMG->ks, ke = pMG->ke;
00303 int Nx1z, Nx2z, Nx3z;
00304 MGrid Coarse_grid;
00305
00306
00307
00308 if (pMG->Nx1==4 || pMG->Nx2==4 || pMG->Nx3==4)
00309 Jacobi(pMG);
00310
00311
00312
00313
00314 else {
00315 Jacobi(pMG);
00316
00317
00318
00319 Coarse_grid.Nx1 = pMG->Nx1/2;
00320 Coarse_grid.Nx2 = pMG->Nx2/2;
00321 Coarse_grid.Nx3 = pMG->Nx3/2;
00322 Coarse_grid.is = 1; Coarse_grid.ie = Coarse_grid.Nx1;
00323 Coarse_grid.js = 1; Coarse_grid.je = Coarse_grid.Nx2;
00324 Coarse_grid.ks = 1; Coarse_grid.ke = Coarse_grid.Nx3;
00325 Coarse_grid.dx1 = 2.0*pMG->dx1;
00326 Coarse_grid.dx2 = 2.0*pMG->dx2;
00327 Coarse_grid.dx3 = 2.0*pMG->dx3;
00328 Coarse_grid.rx1_id = pMG->rx1_id; Coarse_grid.lx1_id = pMG->lx1_id;
00329 Coarse_grid.rx2_id = pMG->rx2_id; Coarse_grid.lx2_id = pMG->lx2_id;
00330 Coarse_grid.rx3_id = pMG->rx3_id; Coarse_grid.lx3_id = pMG->lx3_id;
00331
00332
00333 Nx1z = (pMG->Nx1)/2 + 2;
00334 Nx2z = (pMG->Nx2)/2 + 2;
00335 Nx3z = (pMG->Nx3)/2 + 2;
00336 Coarse_grid.rhs = (Real ***) calloc_3d_array(Nx3z,Nx2z,Nx1z,sizeof(Real));
00337 Coarse_grid.Phi = (Real ***) calloc_3d_array(Nx3z,Nx2z,Nx1z,sizeof(Real));
00338 if (Coarse_grid.rhs == NULL) {
00339 ath_error("[multig_3d]: Error allocating memory for some level\n");
00340 }
00341 if (Coarse_grid.Phi == NULL) {
00342 ath_error("[multig_3d]: Error allocating memory for some level\n");
00343 }
00344
00345 Restriction_3d(pMG, &Coarse_grid);
00346 #ifdef MPI_PARALLEL
00347 set_mg_bvals(&Coarse_grid);
00348 #endif
00349
00350 multig_3d(&Coarse_grid);
00351
00352
00353
00354
00355
00356
00357
00358 Prolongation_3d(&Coarse_grid, pMG);
00359 free_3d_array(Coarse_grid.rhs);
00360 free_3d_array(Coarse_grid.Phi);
00361 #ifdef MPI_PARALLEL
00362 set_mg_bvals(pMG);
00363 #endif
00364
00365
00366
00367 Jacobi(pMG);
00368 }
00369
00370 return;
00371 }
00372
00373
00374
00375
00376
00377
00378
00379
00380 void Jacobi(MGrid *pMG)
00381 {
00382 int i, is = pMG->is, ie = pMG->ie;
00383 int j, js = pMG->js, je = pMG->je;
00384 int k, ks = pMG->ks, ke = pMG->ke;
00385 int n;
00386 Real dx1sq = (pMG->dx1*pMG->dx1);
00387 Real dx2sq = (pMG->dx2*pMG->dx2);
00388 Real dx3sq = (pMG->dx3*pMG->dx3);
00389
00390
00391
00392 if (pMG->Nx3 > 1) {
00393 for (n=0; n<=10; n++){
00394 for (k=ks; k<=ke; k++){
00395 for (j=js; j<=je; j++){
00396 for (i=is; i<=ie; i++){
00397 error[k][j][i] = -pMG->rhs[k][j][i];
00398 error[k][j][i] += (pMG->Phi[k][j][i+1] + pMG->Phi[k][j][i-1])/dx1sq;
00399 error[k][j][i] += (pMG->Phi[k][j+1][i] + pMG->Phi[k][j-1][i])/dx2sq;
00400 error[k][j][i] += (pMG->Phi[k+1][j][i] + pMG->Phi[k-1][j][i])/dx3sq;
00401 }
00402 }
00403 }
00404
00405 for (k=ks; k<=ke; k++){
00406 for (j=js; j<=je; j++){
00407 for (i=is; i<=ie; i++){
00408 pMG->Phi[k][j][i] = error[k][j][i]/(2.0/dx1sq + 2.0/dx2sq + 2.0/dx3sq);
00409 }
00410 }
00411 }
00412 #ifdef MPI_PARALLEL
00413 set_mg_bvals(pMG);
00414 #endif
00415 }
00416 }
00417
00418
00419
00420 if (pMG->Nx3 == 1) {
00421 for (n=0; n<=10; n++){
00422 for (j=js; j<=je; j++){
00423 for (i=is; i<=ie; i++){
00424 error[ks][j][i] = -pMG->rhs[ks][j][i];
00425 error[ks][j][i] += (pMG->Phi[ks][j][i+1] + pMG->Phi[ks][j][i-1])/dx1sq;
00426 error[ks][j][i] += (pMG->Phi[ks][j+1][i] + pMG->Phi[ks][j-1][i])/dx2sq;
00427 }
00428 }
00429
00430 for (j=js; j<=je; j++){
00431 for (i=is; i<=ie; i++){
00432 pMG->Phi[ks][j][i] = error[ks][j][i]/(2.0/dx1sq + 2.0/dx2sq);
00433 }
00434 }
00435 #ifdef MPI_PARALLEL
00436 set_mg_bvals(pMG);
00437 #endif
00438 }
00439 }
00440
00441 return;
00442 }
00443
00444
00445
00446
00447
00448
00449 void Restriction_3d(MGrid *pMG_fine, MGrid *pMG_coarse)
00450 {
00451 int i, is = pMG_fine->is, ie = pMG_fine->ie;
00452 int j, js = pMG_fine->js, je = pMG_fine->je;
00453 int k, ks = pMG_fine->ks, ke = pMG_fine->ke;
00454 Real dx1sq = (pMG_fine->dx1*pMG_fine->dx1);
00455 Real dx2sq = (pMG_fine->dx2*pMG_fine->dx2);
00456 Real dx3sq = (pMG_fine->dx3*pMG_fine->dx3);
00457
00458 for (k=ks; k<=ke; k++){
00459 for (j=js; j<=je; j++){
00460 for (i=is; i<=ie; i++){
00461 error[k][j][i] = pMG_fine->rhs[k][j][i];
00462 error[k][j][i] -= (pMG_fine->Phi[k][j][i+1] + pMG_fine->Phi[k][j][i-1]
00463 - 2.0*pMG_fine->Phi[k][j][i]) / dx1sq;
00464 error[k][j][i] -= (pMG_fine->Phi[k][j+1][i] + pMG_fine->Phi[k][j-1][i]
00465 - 2.0*pMG_fine->Phi[k][j][i]) / dx2sq;
00466 error[k][j][i] -= (pMG_fine->Phi[k+1][j][i] + pMG_fine->Phi[k-1][j][i]
00467 - 2.0*pMG_fine->Phi[k][j][i]) / dx3sq;
00468 }
00469 }
00470 }
00471
00472 for(k=ks; k<=pMG_coarse->ke; k++){
00473 for (j=js; j<=pMG_coarse->je; j++){
00474 for (i=is; i<=pMG_coarse->ie; i++){
00475 pMG_coarse->Phi[k][j][i] =
00476 (pMG_fine->Phi[2*k ][2*j ][2*i] + pMG_fine->Phi[2*k ][2*j ][2*i-1]
00477 + pMG_fine->Phi[2*k ][2*j-1][2*i] + pMG_fine->Phi[2*k ][2*j-1][2*i-1]
00478 + pMG_fine->Phi[2*k-1][2*j ][2*i] + pMG_fine->Phi[2*k-1][2*j ][2*i-1]
00479 + pMG_fine->Phi[2*k-1][2*j-1][2*i] + pMG_fine->Phi[2*k-1][2*j-1][2*i-1])/8.0 ;
00480 pMG_coarse->rhs[k][j][i] =
00481 (error[2*k ][2*j ][2*i] + error[2*k ][2*j ][2*i-1]
00482 + error[2*k ][2*j-1][2*i] + error[2*k ][2*j-1][2*i-1]
00483 + error[2*k-1][2*j ][2*i] + error[2*k-1][2*j ][2*i-1]
00484 + error[2*k-1][2*j-1][2*i] + error[2*k-1][2*j-1][2*i-1])/8.0 ;
00485 }
00486 }
00487 }
00488
00489 return;
00490 }
00491
00492
00493
00494
00495
00496 void Prolongation_3d(MGrid *pMG_coarse, MGrid *pMG_fine)
00497 {
00498 int i, is = pMG_coarse->is, ie = pMG_coarse->ie;
00499 int j, js = pMG_coarse->js, je = pMG_coarse->je;
00500 int k, ks = pMG_coarse->ks, ke = pMG_coarse->ke;
00501
00502 for (k=ks; k<=ke; k++){
00503 for (j=js; j<=je; j++){
00504 for (i=is; i<=ie; i++){
00505 pMG_fine->Phi[2*k ][2*j ][2*i ] += 0.75*pMG_coarse->Phi[k ][j ][i ];
00506 pMG_fine->Phi[2*k ][2*j ][2*i ] += 0.25*pMG_coarse->Phi[k+1][j+1][i+1];
00507
00508 pMG_fine->Phi[2*k ][2*j ][2*i-1] += 0.75*pMG_coarse->Phi[k ][j ][i ];
00509 pMG_fine->Phi[2*k ][2*j ][2*i-1] += 0.25*pMG_coarse->Phi[k+1][j+1][i-1];
00510
00511 pMG_fine->Phi[2*k ][2*j-1][2*i ] += 0.75*pMG_coarse->Phi[k ][j ][i ];
00512 pMG_fine->Phi[2*k ][2*j-1][2*i ] += 0.25*pMG_coarse->Phi[k+1][j-1][i+1];
00513
00514 pMG_fine->Phi[2*k-1][2*j ][2*i ] += 0.75*pMG_coarse->Phi[k ][j ][i ];
00515 pMG_fine->Phi[2*k-1][2*j ][2*i ] += 0.25*pMG_coarse->Phi[k-1][j+1][i+1];
00516
00517 pMG_fine->Phi[2*k ][2*j-1][2*i-1] += 0.75*pMG_coarse->Phi[k ][j ][i ];
00518 pMG_fine->Phi[2*k ][2*j-1][2*i-1] += 0.25*pMG_coarse->Phi[k+1][j-1][i-1];
00519
00520 pMG_fine->Phi[2*k-1][2*j-1][2*i ] += 0.75*pMG_coarse->Phi[k ][j ][i ];
00521 pMG_fine->Phi[2*k-1][2*j-1][2*i ] += 0.25*pMG_coarse->Phi[k-1][j-1][i+1];
00522
00523 pMG_fine->Phi[2*k-1][2*j] [2*i-1] += 0.75*pMG_coarse->Phi[k ][j ][i ];
00524 pMG_fine->Phi[2*k-1][2*j] [2*i-1] += 0.25*pMG_coarse->Phi[k-1][j+1][i-1];
00525
00526 pMG_fine->Phi[2*k-1][2*j-1][2*i-1] += 0.75*pMG_coarse->Phi[k ][j ][i ];
00527 pMG_fine->Phi[2*k-1][2*j-1][2*i-1] += 0.25*pMG_coarse->Phi[k-1][j-1][i-1];
00528 }
00529 }}
00530
00531 return;
00532 }
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547 #ifdef MPI_PARALLEL
00548 void set_mg_bvals(MGrid *pMG)
00549 {
00550 int cnt3, cnt, err;
00551 MPI_Status stat;
00552 MPI_Request rq;
00553
00554
00555
00556
00557 cnt3 = 1;
00558 if (pMG->Nx3 > 1) cnt3 = pMG->Nx3;
00559 cnt = pMG->Nx2*cnt3;
00560
00561
00562 if (pMG->rx1_id >= 0 && pMG->lx1_id >= 0) {
00563
00564 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->lx1_id,
00565 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00566 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00567
00568 swap_mg_ox1(pMG,cnt,0,&rq);
00569 swap_mg_ix1(pMG,cnt,1,&rq);
00570
00571
00572 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->rx1_id,
00573 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00574 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00575
00576 swap_mg_ix1(pMG,cnt,0,&rq);
00577 swap_mg_ox1(pMG,cnt,1,&rq);
00578 }
00579
00580
00581 if (pMG->rx1_id >= 0 && pMG->lx1_id < 0) {
00582
00583 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->rx1_id,
00584 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00585 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00586
00587 swap_mg_ox1(pMG,cnt,0,&rq);
00588 swap_mg_ox1(pMG,cnt,1,&rq);
00589 }
00590
00591
00592 if (pMG->rx1_id < 0 && pMG->lx1_id >= 0) {
00593
00594 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->lx1_id,
00595 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00596 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00597
00598 swap_mg_ix1(pMG,cnt,0,&rq);
00599 swap_mg_ix1(pMG,cnt,1,&rq);
00600 }
00601
00602
00603
00604
00605 cnt3 = 1;
00606 if (pMG->Nx3 > 1) cnt3 = pMG->Nx3;
00607 cnt = (pMG->Nx1 + 2)*cnt3;
00608
00609
00610 if (pMG->rx2_id >= 0 && pMG->lx2_id >= 0) {
00611
00612 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->lx2_id,
00613 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00614 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00615
00616 swap_mg_ox2(pMG,cnt,0,&rq);
00617 swap_mg_ix2(pMG,cnt,1,&rq);
00618
00619
00620 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->rx2_id,
00621 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00622 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00623
00624 swap_mg_ix2(pMG,cnt,0,&rq);
00625 swap_mg_ox2(pMG,cnt,1,&rq);
00626 }
00627
00628
00629 if (pMG->rx2_id >= 0 && pMG->lx2_id < 0) {
00630
00631 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->rx2_id,
00632 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00633 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00634
00635 swap_mg_ox2(pMG,cnt,0,&rq);
00636 swap_mg_ox2(pMG,cnt,1,&rq);
00637 }
00638
00639
00640 if (pMG->rx2_id < 0 && pMG->lx2_id >= 0) {
00641
00642 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->lx2_id,
00643 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00644 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00645
00646 swap_mg_ix2(pMG,cnt,0,&rq);
00647 swap_mg_ix2(pMG,cnt,1,&rq);
00648 }
00649
00650
00651
00652
00653 if (pMG->Nx3 > 1){
00654
00655 cnt = (pMG->Nx1 + 2)*(pMG->Nx2 + 2);
00656
00657
00658 if (pMG->rx3_id >= 0 && pMG->lx3_id >= 0) {
00659
00660 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->lx3_id,
00661 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00662 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00663
00664 swap_mg_ox3(pMG,cnt,0,&rq);
00665 swap_mg_ix3(pMG,cnt,1,&rq);
00666
00667
00668 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->rx3_id,
00669 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00670 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00671
00672 swap_mg_ix3(pMG,cnt,0,&rq);
00673 swap_mg_ox3(pMG,cnt,1,&rq);
00674 }
00675
00676
00677 if (pMG->rx3_id >= 0 && pMG->lx3_id < 0) {
00678
00679 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->rx3_id,
00680 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00681 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00682
00683 swap_mg_ox3(pMG,cnt,0,&rq);
00684 swap_mg_ox3(pMG,cnt,1,&rq);
00685 }
00686
00687
00688 if (pMG->rx3_id < 0 && pMG->lx3_id >= 0) {
00689
00690 err = MPI_Irecv(recv_buf, cnt, MPI_DOUBLE, pMG->lx3_id,
00691 boundary_cells_tag, MPI_COMM_WORLD, &rq);
00692 if(err) ath_error("[set_mg_bvals]: MPI_Irecv error = %d\n",err);
00693
00694 swap_mg_ix3(pMG,cnt,0,&rq);
00695 swap_mg_ix3(pMG,cnt,1,&rq);
00696 }
00697 }
00698
00699 return;
00700 }
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710 void swap_mg_ix1(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00711 {
00712 int j,jl,ju,k,kl,ku,err;
00713 MPI_Status stat;
00714 double *psb = send_buf;
00715 double *prb = recv_buf;
00716
00717 jl = pMG->js;
00718 ju = pMG->je;
00719
00720 if(pMG->Nx3 > 1){
00721 kl = pMG->ks;
00722 ku = pMG->ke;
00723 } else {
00724 kl = ku = pMG->ks;
00725 }
00726
00727
00728
00729 if (swap_flag == 0) {
00730 for (k=kl; k<=ku; k++){
00731 for (j=jl; j<=ju; j++){
00732 *(psb++) = pMG->Phi[k][j][pMG->is];
00733 }
00734 }
00735
00736
00737 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pMG->lx1_id,
00738 boundary_cells_tag, MPI_COMM_WORLD);
00739 if(err) ath_error("[swap_mg_ix1]: MPI_Send error = %d\n",err);
00740 }
00741
00742
00743
00744 if (swap_flag == 1) {
00745
00746
00747 err = MPI_Wait(prq, &stat);
00748 if(err) ath_error("[swap_mg_ix1]: MPI_Wait error = %d\n",err);
00749
00750 for (k=kl; k<=ku; k++){
00751 for (j=jl; j<=ju; j++){
00752 pMG->Phi[k][j][pMG->is-1] = *(prb++);
00753 }
00754 }
00755 }
00756
00757 return;
00758 }
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768 void swap_mg_ox1(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00769 {
00770 int j,jl,ju,k,kl,ku,err;
00771 MPI_Status stat;
00772 double *psb = send_buf;
00773 double *prb = recv_buf;
00774
00775 jl = pMG->js;
00776 ju = pMG->je;
00777
00778 if(pMG->Nx3 > 1){
00779 kl = pMG->ks;
00780 ku = pMG->ke;
00781 } else {
00782 kl = ku = pMG->ks;
00783 }
00784
00785
00786
00787 if (swap_flag == 0) {
00788 for (k=kl; k<=ku; k++){
00789 for (j=jl; j<=ju; j++){
00790 *(psb++) = pMG->Phi[k][j][pMG->ie];
00791 }
00792 }
00793
00794
00795 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pMG->rx1_id,
00796 boundary_cells_tag, MPI_COMM_WORLD);
00797 if(err) ath_error("[swap_mg_ox1]: MPI_Send error = %d\n",err);
00798 }
00799
00800
00801
00802 if (swap_flag == 1) {
00803
00804
00805 err = MPI_Wait(prq, &stat);
00806 if(err) ath_error("[swap_mg_ox1]: MPI_Wait error = %d\n",err);
00807
00808 for (k=kl; k<=ku; k++){
00809 for (j=jl; j<=ju; j++){
00810 pMG->Phi[k][j][pMG->ie+1] = *(prb++);
00811 }
00812 }
00813 }
00814
00815 return;
00816 }
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826 void swap_mg_ix2(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00827 {
00828 int i,il,iu,k,kl,ku,err;
00829 MPI_Status stat;
00830 double *pd = send_buf;
00831
00832 il = pMG->is - 1;
00833 iu = pMG->ie + 1;
00834
00835 if(pMG->Nx3 > 1){
00836 kl = pMG->ks;
00837 ku = pMG->ke;
00838 } else {
00839 kl = ku = pMG->ks;
00840 }
00841
00842
00843
00844 if (swap_flag == 0) {
00845 for (k=kl; k<=ku; k++){
00846 for (i=il; i<=iu; i++){
00847 *(pd++) = pMG->Phi[k][pMG->js][i];
00848 }
00849 }
00850
00851
00852 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pMG->lx2_id,
00853 boundary_cells_tag, MPI_COMM_WORLD);
00854 if(err) ath_error("[swap_mg_ix2]: MPI_Send error = %d\n",err);
00855 }
00856
00857
00858
00859 if (swap_flag == 1) {
00860
00861
00862 err = MPI_Wait(prq, &stat);
00863 if(err) ath_error("[swap_mg_ix2]: MPI_Wait error = %d\n",err);
00864
00865 for (k=kl; k<=ku; k++){
00866 for (i=il; i<=iu; i++){
00867 pMG->Phi[k][pMG->js-1][i] = *(pd++);
00868 }
00869 }
00870 }
00871
00872 return;
00873 }
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883 void swap_mg_ox2(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00884 {
00885 int i,il,iu,k,kl,ku,err;
00886 MPI_Status stat;
00887 double *pd = send_buf;
00888
00889 il = pMG->is - 1;
00890 iu = pMG->ie + 1;
00891
00892 if(pMG->Nx3 > 1){
00893 kl = pMG->ks;
00894 ku = pMG->ke;
00895 } else {
00896 kl = ku = pMG->ks;
00897 }
00898
00899
00900
00901 if (swap_flag == 0) {
00902 for (k=kl; k<=ku; k++){
00903 for (i=il; i<=iu; i++){
00904 *(pd++) = pMG->Phi[k][pMG->je][i];
00905 }
00906 }
00907
00908
00909 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pMG->rx2_id,
00910 boundary_cells_tag, MPI_COMM_WORLD);
00911 if(err) ath_error("[swap_mg_ox2]: MPI_Send error = %d\n",err);
00912 }
00913
00914
00915
00916 if (swap_flag == 1) {
00917
00918
00919 err = MPI_Wait(prq, &stat);
00920 if(err) ath_error("[swap_mg_ox2]: MPI_Wait error = %d\n",err);
00921
00922 for (k=kl; k<=ku; k++){
00923 for (i=il; i<=iu; i++){
00924 pMG->Phi[k][pMG->je+1][i] = *(pd++);
00925 }
00926 }
00927 }
00928
00929 return;
00930 }
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940 void swap_mg_ix3(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00941 {
00942 int i,il,iu,j,jl,ju,err;
00943 MPI_Status stat;
00944 double *pd = send_buf;
00945
00946 il = pMG->is - 1;
00947 iu = pMG->ie + 1;
00948 jl = pMG->js - 1;
00949 ju = pMG->je + 1;
00950
00951
00952
00953 if (swap_flag == 0) {
00954 for (j=jl; j<=ju; j++){
00955 for (i=il; i<=iu; i++){
00956 *(pd++) = pMG->Phi[pMG->ks][j][i];
00957 }
00958 }
00959
00960
00961 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pMG->lx3_id,
00962 boundary_cells_tag, MPI_COMM_WORLD);
00963 if(err) ath_error("[swap_mg_ix3]: MPI_Send error = %d\n",err);
00964 }
00965
00966
00967
00968 if (swap_flag == 1) {
00969
00970
00971 err = MPI_Wait(prq, &stat);
00972 if(err) ath_error("[swap_mg_ix3]: MPI_Wait error = %d\n",err);
00973
00974 for (j=jl; j<=ju; j++){
00975 for (i=il; i<=iu; i++){
00976 pMG->Phi[pMG->ks-1][j][i] = *(pd++);
00977 }
00978 }
00979 }
00980
00981 return;
00982 }
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992 void swap_mg_ox3(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00993 {
00994 int i,il,iu,j,jl,ju,err;
00995 MPI_Status stat;
00996 double *pd = send_buf;
00997
00998 il = pMG->is - 1;
00999 iu = pMG->ie + 1;
01000 jl = pMG->js - 1;
01001 ju = pMG->je + 1;
01002
01003
01004
01005 if (swap_flag == 0) {
01006 for (j=jl; j<=ju; j++){
01007 for (i=il; i<=iu; i++){
01008 *(pd++) = pMG->Phi[pMG->ke][j][i];
01009 }
01010 }
01011
01012
01013 err = MPI_Send(send_buf, cnt, MPI_DOUBLE, pMG->rx3_id,
01014 boundary_cells_tag, MPI_COMM_WORLD);
01015 if(err) ath_error("[swap_mg_ox3]: MPI_Send error = %d\n",err);
01016 }
01017
01018
01019
01020 if (swap_flag == 1) {
01021
01022
01023 err = MPI_Wait(prq, &stat);
01024 if(err) ath_error("[swap_mg_ox3]: MPI_Wait error = %d\n",err);
01025
01026 for (j=jl; j<=ju; j++){
01027 for (i=il; i<=iu; i++){
01028 pMG->Phi[pMG->ke+1][j][i] = *(pd++);
01029 }
01030 }
01031 }
01032
01033 return;
01034 }
01035 #endif
01036
01037
01038
01039
01040
01041
01042
01043 void selfg_multig_3d_init(MeshS *pM)
01044 {
01045 #ifdef MPI_PARALLEL
01046 int i,j,k,x1cnt,x2cnt,x3cnt;
01047 int nx1t,nx2t,nx3t, size;
01048 int NGrid_x1, NGrid_x2, NGrid_x3;
01049 #endif
01050 int nmax,size1=0,size2=0,size3=0,nl,nd;
01051
01052
01053 for (nl=0; nl<(pM->NLevels); nl++){
01054 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
01055 if (pM->Domain[nl][nd].Grid != NULL) {
01056 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
01057 size1 = pM->Domain[nl][nd].Grid->Nx[0];
01058 }
01059 if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
01060 size2 = pM->Domain[nl][nd].Grid->Nx[1];
01061 }
01062 if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
01063 size3 = pM->Domain[nl][nd].Grid->Nx[2];
01064 }
01065 }
01066 }
01067 }
01068
01069 size1 += 2;
01070 size2 += 2;
01071 size3 += 2;
01072 error = (Real ***)calloc_3d_array(size3, size2, size1, sizeof(Real));
01073 if (error == NULL) {
01074 ath_error("[multig_3d]: Error allocating memory for error array\n");
01075 }
01076
01077
01078
01079
01080 #ifdef MPI_PARALLEL
01081 NGrid_x1 = par_geti("parallel","NGrid_x1");
01082 NGrid_x2 = par_geti("parallel","NGrid_x2");
01083 NGrid_x3 = par_geti("parallel","NGrid_x3");
01084
01085 x1cnt = x2cnt = x3cnt = 0;
01086
01087 for (k=0; k<NGrid_x3; k++){
01088 for (j=0; j<NGrid_x2; j++){
01089 for (i=0; i<NGrid_x1; i++){
01090 if(NGrid_x1 > 1){
01091 nx2t = pD->grid_block[k][j][i].jde - pD->grid_block[k][j][i].jds + 1;
01092 nx3t = pD->grid_block[k][j][i].kde - pD->grid_block[k][j][i].kds + 1;
01093
01094 x1cnt = nx2t*nx3t > x1cnt ? nx2t*nx3t : x1cnt;
01095 }
01096
01097 if(NGrid_x2 > 1){
01098 nx1t = pD->grid_block[k][j][i].ide - pD->grid_block[k][j][i].ids + 1;
01099 if(nx1t > 1) nx1t += 2;
01100 nx3t = pD->grid_block[k][j][i].kde - pD->grid_block[k][j][i].kds + 1;
01101
01102 x2cnt = nx1t*nx3t > x2cnt ? nx1t*nx3t : x2cnt;
01103 }
01104
01105
01106 if(NGrid_x3 > 1){
01107 nx1t = pD->grid_block[k][j][i].ide - pD->grid_block[k][j][i].ids + 1;
01108 if(nx1t > 1) nx1t += 2;
01109 nx2t = pD->grid_block[k][j][i].jde - pD->grid_block[k][j][i].jds + 1;
01110 if(nx2t > 1) nx2t += 2;
01111
01112 x3cnt = nx1t*nx2t > x3cnt ? nx1t*nx2t : x3cnt;
01113 }
01114 }
01115 }
01116 }
01117
01118 size = x1cnt > x2cnt ? x1cnt : x2cnt;
01119 size = x3cnt > size ? x3cnt : size;
01120
01121 if (size > 0) {
01122 if((send_buf = (double*)malloc(size*sizeof(double))) == NULL)
01123 ath_error("[selfg_by_multig_3d_init]: Failed to allocate send buffer\n");
01124
01125 if((recv_buf = (double*)malloc(size*sizeof(double))) == NULL)
01126 ath_error("[selfg_by_multig_3d_init]: Failed to allocate recv buffer\n");
01127 }
01128 #endif
01129 return;
01130 }
01131
01132 #endif