Athena: gravity/selfg_multigrid.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file selfg_multigrid.c
00004  *  \brief Contains functions to solve Poisson's equation for self-gravity in
00005  *   3D using multigrid.
00006  *
00007  *   These functions work for non-periodic domains.  A low-order multipole
00008  *   expansion is used to compute the potential on the boundaries.
00009  *
00010  * HISTORY:
00011  * - june-2007 - 2D and 3D solvers written by Irene Balmes
00012  * - july-2007 - routines incorporated into Athena by JMS and IB
00013  *
00014  * CONTAINS PUBLIC FUNCTIONS:
00015  * - selfg_by_multig_2d() - 2D Poisson solver using multigrid
00016  * - selfg_by_multig_3d() - 3D Poisson solver using multigrid
00017  * - selfg_by_multig_3d_init() - Initializes send/receive buffers for MPI */
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 /*! \struct MGrid
00043  *  \brief Holds RHS, potential, and information about grid
00044  * size for a given level in the multi-grid hierarchy  */
00045 typedef struct MGrid_s{
00046   Real ***rhs,***Phi;  /* RHS of elliptic equation, and solution */
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 /* 3D temporary array needed for restriction of errors  */
00058 Real ***error;
00059 
00060 /*==============================================================================
00061  * PRIVATE FUNCTION PROTOTYPES:
00062  *   multig_3d() -
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 /*=========================== PUBLIC FUNCTIONS ===============================*/
00082 /*----------------------------------------------------------------------------*/
00083 /*! \fn void selfg_multig_1d(DomainS *pD)
00084  *  \brief 1-D multigrid gravity
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 /*! \fn void selfg_multig_2d(DomainS *pD)
00095  *  \brief 2-D multigrid gravity
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 /*! \fn void selfg_multig_3d(DomainS *pD)
00106  *  \brief Do not use with periodic BCs, uses multipole expansion
00107  *   to compute potential at boundary
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 /* Copy current potential into old */
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 /* Compute solution at boundaries using monopole expansion */
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 /* MPI_PARALLEL */
00152 
00153 /*  Inner and outer x1 boundaries */
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 /*  Inner and outer x2 boundaries */
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 /*  Inner and outer x3 boundaries */
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 /* Initialize MGrid structure for top level (the root, or finest, level) */
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 /* There is only one ghost zone needed at each level, not nghost */
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 /* Initialize solution on root grid, including single ghost zone */
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 /* Compute new potential.  Note multig_3d calls itself recursively. */
00273 
00274   multig_3d(&Root_grid);
00275 
00276 /* copy solution for potential from MGrid into Grid structure.  Boundary
00277  * conditions for nghost ghost cells are set by set_bvals() call in main() */
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 /*! \fn void multig_3d(MGrid *pMG)
00295  *  \brief Functions needed for the multigrid solver in 3D
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 /* If we are down to 4 cells in any dimension do 10 iterations and return */
00307 
00308   if (pMG->Nx1==4 || pMG->Nx2==4 || pMG->Nx3==4)
00309     Jacobi(pMG);
00310 
00311 /* Else, do 10 iterations at this level, restrict to a coarser grid, and call
00312  * multig_3d again with this coarse grid */
00313 
00314   else { 
00315     Jacobi(pMG);
00316 
00317 /* Allocate and initialize MGrid at next coarsest level */
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 /* Again, only one ghost zone on each level */
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 /* The following code is first reached after 10 iterations at the coarsest
00353  * level.  We then prolongate, do 10 iterations, and return.  This will return
00354  * execution to this same spot for the next coarsest level, so we will
00355  * prolongate, do 10 iterations, return, and so on.
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 /* End with 10 iterations at the finest level */
00366 
00367     Jacobi(pMG);
00368   }
00369 
00370   return;
00371 }
00372 
00373 /*----------------------------------------------------------------------------*/
00374 /*! \fn void Jacobi(MGrid *pMG)
00375  *  \brief Jacobi iterations. 
00376  *
00377  *   Do not use with periodic BCs, uses multipole expansion
00378  *   to compute potential at boundary
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 /* Jacobi iterations in 3D */
00391 
00392   if (pMG->Nx3 > 1) {
00393     for (n=0; n<=10; n++){  /* hardwired to do 10 iterations */
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 /* Jacobi iterations in 2D (x-y plane) */
00419 
00420   if (pMG->Nx3 == 1) {
00421     for (n=0; n<=10; n++){  /* hardwired to do 10 iterations */
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 /*! \fn void Restriction_3d(MGrid *pMG_fine, MGrid *pMG_coarse) 
00446  *  \brief Averages fine grid solution onto coarse
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 /*! \fn void Prolongation_3d(MGrid *pMG_coarse, MGrid *pMG_fine)
00494  *  \brief Linear interpolation of coarse grid onto fine  
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 /*! \fn void set_mg_bvals(MGrid *pMG)
00536  *  \brief Sets BC for Jacobi iterates for MPI parallel jobs.
00537  *
00538  *   With self-gravity using multigrid, the boundary conditions at the edge of
00539  *   the Domain are held fixed.  So only ghostzones associated with internal
00540  *   boundaries between MPI grids need to be passed.
00541  *
00542  * This routine is largely a copy of set_bvals().
00543  * Order for updating boundary conditions must always be x1-x2-x3 in order to
00544  * fill the corner cells properly
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 /*--- Step 1. ------------------------------------------------------------------
00555  * Boundary Conditions in x1-direction */
00556 
00557   cnt3 = 1;
00558   if (pMG->Nx3 > 1) cnt3 = pMG->Nx3;
00559   cnt = pMG->Nx2*cnt3;
00560 
00561 /* MPI blocks to both left and right */
00562   if (pMG->rx1_id >= 0 && pMG->lx1_id >= 0) {
00563     /* Post a non-blocking receive for the input data from the left grid */
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);  /* send R */
00569     swap_mg_ix1(pMG,cnt,1,&rq);  /* listen L */
00570 
00571     /* Post a non-blocking receive for the input data from the right grid */
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);  /* send L */
00577     swap_mg_ox1(pMG,cnt,1,&rq);  /* listen R */
00578   }
00579 
00580 /* Physical boundary on left, MPI block on right */
00581   if (pMG->rx1_id >= 0 && pMG->lx1_id < 0) {
00582     /* Post a non-blocking receive for the input data from the right grid */
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);  /* send R */
00588     swap_mg_ox1(pMG,cnt,1,&rq);  /* listen R */
00589   }
00590 
00591 /* MPI block on left, Physical boundary on right */
00592   if (pMG->rx1_id < 0 && pMG->lx1_id >= 0) {
00593     /* Post a non-blocking receive for the input data from the left grid */
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);  /* send L */
00599     swap_mg_ix1(pMG,cnt,1,&rq);  /* listen L */
00600   }
00601 
00602 /*--- Step 2. ------------------------------------------------------------------
00603  * Boundary Conditions in x2-direction */
00604 
00605   cnt3 = 1;
00606   if (pMG->Nx3 > 1) cnt3 = pMG->Nx3;
00607   cnt = (pMG->Nx1 + 2)*cnt3;
00608 
00609 /* MPI blocks to both left and right */
00610   if (pMG->rx2_id >= 0 && pMG->lx2_id >= 0) {
00611     /* Post a non-blocking receive for the input data from the left grid */
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);  /* send R */
00617     swap_mg_ix2(pMG,cnt,1,&rq);  /* listen L */
00618 
00619     /* Post a non-blocking receive for the input data from the right grid */
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);  /* send L */
00625     swap_mg_ox2(pMG,cnt,1,&rq);  /* listen R */
00626   }
00627 
00628 /* Physical boundary on left, MPI block on right */
00629   if (pMG->rx2_id >= 0 && pMG->lx2_id < 0) {
00630     /* Post a non-blocking receive for the input data from the right grid */
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);  /* send R */
00636     swap_mg_ox2(pMG,cnt,1,&rq);  /* listen R */
00637   }
00638 
00639 /* MPI block on left, Physical boundary on right */
00640   if (pMG->rx2_id < 0 && pMG->lx2_id >= 0) {
00641     /* Post a non-blocking receive for the input data from the left grid */
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);  /* send L */
00647     swap_mg_ix2(pMG,cnt,1,&rq);  /* listen L */
00648   }
00649 
00650 /*--- Step 3. ------------------------------------------------------------------
00651  * Boundary Conditions in x3-direction */
00652 
00653   if (pMG->Nx3 > 1){
00654 
00655     cnt = (pMG->Nx1 + 2)*(pMG->Nx2 + 2);
00656 
00657 /* MPI blocks to both left and right */
00658     if (pMG->rx3_id >= 0 && pMG->lx3_id >= 0) {
00659       /* Post a non-blocking receive for the input data from the left grid */
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);  /* send R */
00665       swap_mg_ix3(pMG,cnt,1,&rq);  /* listen L */
00666 
00667       /* Post a non-blocking receive for the input data from the right grid */
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);  /* send L */
00673       swap_mg_ox3(pMG,cnt,1,&rq);  /* listen R */
00674     }
00675 
00676 /* Physical boundary on left, MPI block on right */
00677     if (pMG->rx3_id >= 0 && pMG->lx3_id < 0) {
00678       /* Post a non-blocking receive for the input data from the right grid */
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);  /* send R */
00684       swap_mg_ox3(pMG,cnt,1,&rq);  /* listen R */
00685     }
00686 
00687 /* MPI block on left, Physical boundary on right */
00688     if (pMG->rx3_id < 0 && pMG->lx3_id >= 0) {
00689       /* Post a non-blocking receive for the input data from the left grid */
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);  /* send L */
00695       swap_mg_ix3(pMG,cnt,1,&rq);  /* listen L */
00696     }
00697   }
00698 
00699   return;
00700 }
00701 
00702 /*----------------------------------------------------------------------------*/
00703 /*! \fn void swap_mg_ix1(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00704  *  \brief MPI_SWAP of boundary conditions, Inner x1 boundary
00705  *
00706  *   This function does either a send (swap_flag=0), or receive (swap_flag=1)
00707  *   Largely copied from set_bvals/send_ix1 and set_bvals/receive_ix1
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 /* Pack single row i=is of iterates into send buffer */
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     /* send contents of buffer to the neighboring grid on L-x1 */
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 /* Receive message and unpack single row i=is-1 of iteratas into ghost zonee */
00743 
00744   if (swap_flag == 1) {
00745 
00746     /* Wait to receive the input data from the left grid */
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 /*! \fn void swap_mg_ox1(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00762  *  \brief MPI_SWAP of boundary conditions, Outer x1 boundary
00763  *
00764  *   This function does either a send (swap_flag=0), or receive (swap_flag=1)
00765  *   Largely copied from set_bvals/send_ox1 and set_bvals/receive_ox1
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 /* Pack single row i=ie of iterates into send buffer */
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     /* send contents of buffer to the neighboring grid on R-x1 */
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 /* Receive message and unpack single row i=ie+1 of iterates into ghost zone */
00801 
00802   if (swap_flag == 1) {
00803 
00804     /* Wait to receive the input data from the right grid */
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 /*! \fn void swap_mg_ix2(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00820  *  \brief MPI_SWAP of boundary conditions, Inner x2 boundary
00821  *
00822  *   This function does either a send (swap_flag=0), or receive (swap_flag=1)
00823  *   Largely copied from set_bvals/send_ix2 and set_bvals/receive_ix2
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 /* Pack single row j=js of iterates into send buffer */
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     /* send contents of buffer to the neighboring grid on L-x2 */
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 /* Receive message and unpack single row j=js-1 of iterates into ghost zone */
00858 
00859   if (swap_flag == 1) {
00860 
00861     /* Wait to receive the input data from the left grid */
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 /*! \fn void swap_mg_ox2(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00877  *  \brief MPI_SWAP of boundary conditions, Outer x2 boundary
00878  *
00879  *   This function does either a send (swap_flag=0), or receive (swap_flag=1)
00880  *   Largely copied from set_bvals/send_ix2 and set_bvals/receive_ix2
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 /* Pack single row j=je of iterates into send buffer */
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     /* send contents of buffer to the neighboring grid on R-x2 */
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 /* Receive message and unpack single row j=je+1 of iterates into ghost zone */
00915 
00916   if (swap_flag == 1) {
00917 
00918     /* Wait to receive the input data from the right grid */
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 /*! \fn void swap_mg_ix3(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00934  *  \brief MPI_SWAP of boundary conditions, Inner x3 boundary
00935  *
00936  *   This function does either a send (swap_flag=0), or receive (swap_flag=1)
00937  *   Largely copied from set_bvals/send_ix3 and set_bvals/receive_ix3
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 /* Pack single row k=ks of iterates into send buffer */
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     /* send contents of buffer to the neighboring grid on L-x3 */
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 /* Receive message and unpack single row k=ks-1 of iterates into ghost zone */
00967 
00968   if (swap_flag == 1) {
00969 
00970     /* Wait to receive the input data from the left grid */
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 /*! \fn void swap_mg_ox3(MGrid *pMG, int cnt, int swap_flag, MPI_Request *prq)
00986  *  \brief MPI_SWAP of boundary conditions, Outer x3 boundary
00987  *
00988  *   This function does either a send (swap_flag=0), or receive (swap_flag=1)
00989  *   Largely copied from set_bvals/send_ix1 and set_bvals/receive_ix1
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 /* Pack single row k=ke of interates into send buffer */
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     /* send contents of buffer to the neighboring grid on R-x3 */
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 /* Receive message and unpack single row k=ke+1 of iterates into ghost zone */
01019 
01020   if (swap_flag == 1) {
01021 
01022     /* Wait to receive the input data from the right grid */
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 /* MPI_PARALLEL */
01036 
01037 /*----------------------------------------------------------------------------*/
01038 /*! \fn void selfg_multig_3d_init(MeshS *pM)
01039  *  \brief Initialize send/receive buffers needed to swap
01040  *   iterates during Jacobi iterations.
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 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
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 /* Allocate memory for send and receive buffers for Phi in MultiGrid
01078  * structure for MPI parallel.
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 /* MPI_PARALLEL */
01129   return;
01130 }
01131 
01132 #endif /* SELF_GRAVITY_USING_MULTIGRID */