Athena: init_grid.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file init_grid.c 
00004  *  \brief Initializes most variables in the Grid structure.
00005  *
00006  * PURPOSE: Initializes most variables in the Grid structure.  Allocates memory
00007  *   for 3D arrays of Cons, interface B, etc.  With SMR, finds all overlaps
00008  *   between child and parent Grids, and initializes data needed for restriction
00009  *   flux-correction, and prolongation steps.
00010  *
00011  * CONTAINS PUBLIC FUNCTIONS: 
00012  * - init_grid()
00013  *
00014  * PRIVATE FUNCTION PROTOTYPES:
00015  * - checkOverlap() - checks for overlap of cubes, and returns overlap coords
00016  * - checkOverlapTouch() - same as above, but checks for overlap and/or touch */
00017 /*============================================================================*/
00018 
00019 #include <math.h>
00020 #include <stdio.h>
00021 #include <stdlib.h>
00022 #include "defs.h"
00023 #include "athena.h"
00024 #include "globals.h"
00025 #include "prototypes.h"
00026 
00027 /*==============================================================================
00028  * PRIVATE FUNCTION PROTOTYPES:
00029  *  checkOverlap() - checks for overlap of cubes, and returns overlap coords
00030  *  checkOverlapTouch() - same as above, but checks for overlap and/or touch
00031  *============================================================================*/
00032 #ifdef STATIC_MESH_REFINEMENT
00033 /*! \fn int checkOverlap(SideS *pC1, SideS *pC2, SideS *pC3);
00034  *  \brief Checks for overlap of cubes, and returns overlap coords */
00035 int checkOverlap(SideS *pC1, SideS *pC2, SideS *pC3);
00036 /*! \fn int checkOverlapTouch(SideS *pC1, SideS *pC2, SideS *pC3);
00037  *  \brief Same as above, but checks for overlap and/or touch */
00038 int checkOverlapTouch(SideS *pC1, SideS *pC2, SideS *pC3);
00039 #endif
00040 
00041 /*----------------------------------------------------------------------------*/
00042 /*! \fn void init_grid(MeshS *pM)
00043  *  \brief Initializes most variables in the Grid structure.
00044  *
00045  * PURPOSE: Initializes most variables in the Grid structure.  Allocates memory
00046  *   for 3D arrays of Cons, interface B, etc.  With SMR, finds all overlaps
00047  *   between child and parent Grids, and initializes data needed for restriction
00048  *   flux-correction, and prolongation steps.                                 */
00049 
00050 void init_grid(MeshS *pM)
00051 {
00052   DomainS *pD;
00053   GridS *pG;
00054   int nDim,nl,nd,myL,myM,myN;
00055   int i,l,m,n,n1z,n2z,n3z,n1p,n2p,n3p;
00056 #ifdef STATIC_MESH_REFINEMENT
00057   DomainS *pCD,*pPD;
00058   SideS D1,D2,D3,G1,G2,G3;
00059   int isDOverlap,isGOverlap,irefine,ncd,npd,dim,iGrid;
00060   int ncg,nCG,nMyCG,nCB[6],nMyCB[6],nb;
00061   int npg,nPG,nMyPG,nPB[6],nMyPB[6];
00062 #endif
00063 
00064 /* number of dimensions in Grid. */
00065   nDim=1;
00066   for (i=1; i<3; i++) if (pM->Nx[i]>1) nDim++;
00067 
00068 /* Loop over all levels and domains per level */
00069 
00070   for (nl=0; nl<pM->NLevels; nl++){
00071   for (nd=0; nd<pM->DomainsPerLevel[nl]; nd++){
00072     if (pM->Domain[nl][nd].Grid != NULL) {
00073       pD = (DomainS*)&(pM->Domain[nl][nd]);  /* set ptr to Domain */
00074       pG = pM->Domain[nl][nd].Grid;          /* set ptr to Grid */
00075 
00076       pG->time = pM->time;
00077 
00078 /* get (l,m,n) coordinates of Grid being updated on this processor */
00079 
00080       get_myGridIndex(pD, myID_Comm_world, &myL, &myM, &myN);
00081 
00082 /* ---------------------  Intialize grid in 1-direction --------------------- */
00083 /* Initialize is,ie,dx1
00084  * Compute Disp, MinX[0], and MaxX[0] using displacement of Domain and Grid
00085  * location within Domain */
00086 
00087       pG->Nx[0] = pD->GData[myN][myM][myL].Nx[0];
00088 
00089       if(pG->Nx[0] > 1) {
00090         pG->is = nghost;
00091         pG->ie = pG->Nx[0] + nghost - 1;
00092       }
00093       else
00094         pG->is = pG->ie = 0;
00095     
00096       pG->dx1 = pD->dx[0];
00097     
00098       pG->Disp[0] = pD->Disp[0];
00099       pG->MinX[0] = pD->MinX[0];
00100       for (l=1; l<=myL; l++) {
00101         pG->Disp[0] +=        pD->GData[myN][myM][l-1].Nx[0];
00102         pG->MinX[0] += (Real)(pD->GData[myN][myM][l-1].Nx[0])*pG->dx1;
00103       }
00104       pG->MaxX[0] = pG->MinX[0] + (Real)(pG->Nx[0])*pG->dx1;
00105     
00106 /* ---------------------  Intialize grid in 2-direction --------------------- */
00107 /* Initialize js,je,dx2
00108  * Compute Disp, MinX[1], and MaxX[1] using displacement of Domain and Grid
00109  * location within Domain */
00110 
00111       pG->Nx[1] = pD->GData[myN][myM][myL].Nx[1];
00112     
00113       if(pG->Nx[1] > 1) {
00114         pG->js = nghost;
00115         pG->je = pG->Nx[1] + nghost - 1;
00116       }
00117       else
00118         pG->js = pG->je = 0;
00119     
00120       pG->dx2 = pD->dx[1];
00121 
00122       pG->Disp[1] = pD->Disp[1];
00123       pG->MinX[1] = pD->MinX[1];
00124       for (m=1; m<=myM; m++) {
00125         pG->Disp[1] +=        pD->GData[myN][m-1][myL].Nx[1];
00126         pG->MinX[1] += (Real)(pD->GData[myN][m-1][myL].Nx[1])*pG->dx2;
00127       }
00128       pG->MaxX[1] = pG->MinX[1] + (Real)(pG->Nx[1])*pG->dx2;
00129 
00130 /* ---------------------  Intialize grid in 3-direction --------------------- */
00131 /* Initialize ks,ke,dx3
00132  * Compute Disp, MinX[2], and MaxX[2] using displacement of Domain and Grid
00133  * location within Domain */
00134 
00135       pG->Nx[2] = pD->GData[myN][myM][myL].Nx[2];
00136 
00137       if(pG->Nx[2] > 1) {
00138         pG->ks = nghost;
00139         pG->ke = pG->Nx[2] + nghost - 1;
00140       }
00141       else
00142         pG->ks = pG->ke = 0;
00143 
00144       pG->dx3 = pD->dx[2];
00145 
00146       pG->Disp[2] = pD->Disp[2];
00147       pG->MinX[2] = pD->MinX[2];
00148       for (n=1; n<=myN; n++) {
00149         pG->Disp[2] +=        pD->GData[n-1][myM][myL].Nx[2];
00150         pG->MinX[2] += (Real)(pD->GData[n-1][myM][myL].Nx[2])*pG->dx3;
00151       }
00152       pG->MaxX[2] = pG->MinX[2] + (Real)(pG->Nx[2])*pG->dx3;
00153 
00154 /* ---------  Allocate 3D arrays to hold Cons based on size of grid --------- */
00155 
00156       if (pG->Nx[0] > 1)
00157         n1z = pG->Nx[0] + 2*nghost;
00158       else
00159         n1z = 1;
00160 
00161       if (pG->Nx[1] > 1)
00162         n2z = pG->Nx[1] + 2*nghost;
00163       else
00164         n2z = 1;
00165 
00166       if (pG->Nx[2] > 1)
00167         n3z = pG->Nx[2] + 2*nghost;
00168       else
00169         n3z = 1;
00170 
00171 /* Build a 3D array of type ConsS */
00172 
00173       pG->U = (ConsS***)calloc_3d_array(n3z, n2z, n1z, sizeof(ConsS));
00174       if (pG->U == NULL) goto on_error1;
00175     
00176 /* Build 3D arrays to hold interface field */
00177 
00178 #ifdef MHD
00179       pG->B1i = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00180       if (pG->B1i == NULL) goto on_error2;
00181 
00182       pG->B2i = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00183       if (pG->B2i == NULL) goto on_error3;
00184 
00185       pG->B3i = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00186       if (pG->B3i == NULL) goto on_error4;
00187 #endif /* MHD */
00188 
00189 /* Build 3D arrays to magnetic diffusivities */
00190 
00191 #ifdef RESISTIVITY
00192       pG->eta_Ohm = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00193       if (pG->eta_Ohm == NULL) goto on_error5;
00194 
00195       pG->eta_Hall = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00196       if (pG->eta_Hall == NULL) goto on_error6;
00197 
00198       pG->eta_AD = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00199       if (pG->eta_AD == NULL) goto on_error7;
00200 #endif /* RESISTIVITY */
00201 
00202 /* Build 3D arrays to gravitational potential and mass fluxes */
00203 
00204 #ifdef SELF_GRAVITY
00205       pG->Phi = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00206       if (pG->Phi == NULL) goto on_error9;
00207 
00208       pG->Phi_old = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00209       if (pG->Phi_old == NULL) goto on_error10;
00210 
00211       pG->x1MassFlux = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00212       if (pG->x1MassFlux == NULL) goto on_error11;
00213 
00214       pG->x2MassFlux = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00215       if (pG->x2MassFlux == NULL) goto on_error12;
00216 
00217       pG->x3MassFlux = (Real***)calloc_3d_array(n3z, n2z, n1z, sizeof(Real));
00218       if (pG->x3MassFlux == NULL) goto on_error13;
00219 #endif /* SELF_GRAVITY */
00220 
00221 /* Allocate and initialize cylindrical scaling factors */
00222 #ifdef CYLINDRICAL
00223       pG->r = (Real*)calloc_1d_array(n1z, sizeof(Real));
00224       if (pG->r == NULL) goto on_error14;
00225 
00226       pG->ri = (Real*)calloc_1d_array(n1z, sizeof(Real));
00227       if (pG->ri == NULL) goto on_error15;
00228       for (i=pG->is-nghost; i<=pG->ie+nghost; i++) {
00229         pG->ri[i] = pG->MinX[0] + ((Real)(i - pG->is))*pG->dx1;
00230         pG->r[i]  = pG->ri[i] + 0.5*pG->dx1;
00231       }
00232 #endif /* CYLINDRICAL */
00233 
00234 
00235 /*-- Get IDs of neighboring Grids in Domain communicator ---------------------*/
00236 /* If Grid is at the edge of the Domain (so it is either a physical boundary,
00237  * or an internal boundary between fine/coarse grids), then ID is set to -1
00238  */
00239 
00240 /* Left-x1 */
00241       if(myL > 0) pG->lx1_id = pD->GData[myN][myM][myL-1].ID_Comm_Domain;
00242       else pG->lx1_id = -1;
00243 
00244 /* Right-x1 */
00245       if(myL <(pD->NGrid[0])-1)
00246         pG->rx1_id = pD->GData[myN][myM][myL+1].ID_Comm_Domain;
00247       else pG->rx1_id = -1;
00248 
00249 /* Left-x2 */
00250       if(myM > 0) pG->lx2_id = pD->GData[myN][myM-1][myL].ID_Comm_Domain;
00251       else pG->lx2_id = -1;
00252 
00253 /* Right-x2 */
00254       if(myM <(pD->NGrid[1])-1)
00255         pG->rx2_id = pD->GData[myN][myM+1][myL].ID_Comm_Domain;
00256       else pG->rx2_id = -1;
00257 
00258 /* Left-x3 */
00259       if(myN > 0) pG->lx3_id = pD->GData[myN-1][myM][myL].ID_Comm_Domain;
00260       else pG->lx3_id = -1;
00261 
00262 /* Right-x3 */
00263       if(myN <(pD->NGrid[2])-1)
00264         pG->rx3_id = pD->GData[myN+1][myM][myL].ID_Comm_Domain;
00265       else pG->rx3_id = -1;
00266    
00267 #ifdef STATIC_MESH_REFINEMENT
00268 /*---------------------- Initialize variables for SMR ------------------------*/
00269 /* Number of child/parent grids, and data about overlap regions. */
00270 
00271       pG->NCGrid = 0;
00272       pG->NPGrid = 0;
00273       pG->NmyCGrid = 0;  /* can be as large as # of child Domains */
00274       pG->NmyPGrid = 0;  /* must be 0 or 1 */
00275       pG->CGrid = NULL;
00276       pG->PGrid = NULL;
00277 #endif
00278 
00279     }
00280   }}
00281 
00282 #ifdef STATIC_MESH_REFINEMENT
00283 /*-------------- Count number of child Grids, and boundaries -----------------*/
00284 /* First we have to count the number of child Grids, and fine/coarse boundaries
00285  * with child Grids, including child Grids on this and other processors, before
00286  * we can allocate the CGrid array.  This is a shortcoming of making these
00287  * structures 1D arrays. */ 
00288 /* Loop over levels (up to next to last level: maxlevel-1), and domains/level */
00289 
00290   for (nl=0; nl<(pM->NLevels)-1; nl++){
00291   for (nd=0; nd<pM->DomainsPerLevel[nl]; nd++){
00292     if (pM->Domain[nl][nd].Grid != NULL) {
00293       pD = (DomainS*)&(pM->Domain[nl][nd]);  /* set ptr to this Domain */
00294       pG = pM->Domain[nl][nd].Grid;          /* set ptr to this Grid */
00295 
00296 /* edges of this Domain */
00297       for (i=0; i<3; i++) {
00298         D1.ijkl[i] = pD->Disp[i];
00299         D1.ijkr[i] = pD->Disp[i] + pD->Nx[i];
00300       }
00301 
00302 /* edges of this Grid */
00303       for (i=0; i<3; i++) {
00304         G1.ijkl[i] = pG->Disp[i];
00305         G1.ijkr[i] = pG->Disp[i] + pG->Nx[i];
00306       }
00307 
00308 /* For this Domain, check if there is a Domain at next level that overlaps.
00309  * Divide by two to check in units of this (not the child) domain coordinates */
00310 
00311       for (ncd=0; ncd<pM->DomainsPerLevel[nl+1]; ncd++){
00312         pCD = (DomainS*)&(pM->Domain[nl+1][ncd]);  /* ptr to potential child  */
00313 
00314 /* edges of potential child Domain */
00315         for (i=0; i<3; i++) {
00316           D2.ijkl[i] = pCD->Disp[i]/2;
00317           D2.ijkr[i] = 1;
00318           if (pCD->Nx[i] > 1) D2.ijkr[i] = (pCD->Disp[i] + pCD->Nx[i])/2;
00319         }
00320 
00321         isDOverlap = checkOverlap(&D1, &D2, &D3);
00322         if (isDOverlap == 1){
00323 
00324 /*----------------------------------------------------------------------------*/
00325 /* There is a child Domain that overlaps. So on the child Domain, find all the
00326  * Grids that overlap this Grid. */
00327 
00328           for (n=0; n<pCD->NGrid[2]; n++){
00329           for (m=0; m<pCD->NGrid[1]; m++){
00330           for (l=0; l<pCD->NGrid[0]; l++){
00331 
00332 /* edges of child Grid */
00333 /* Divide by two to check in units of this (not the child) grid coordinates */
00334             for (i=0; i<3; i++) {
00335               G2.ijkl[i] = pCD->GData[n][m][l].Disp[i]/2;
00336               G2.ijkr[i] = 1;
00337               if (pCD->Nx[i] > 1) 
00338                 G2.ijkr[i] = (pCD->GData[n][m][l].Disp[i]
00339                             + pCD->GData[n][m][l].Nx[i])/2;
00340             }
00341 
00342             isGOverlap = checkOverlap(&G1, &G2, &G3);
00343 
00344 /* If Grid overlaps, increment Child counters */
00345 
00346             if (isGOverlap == 1){
00347               pG->NCGrid++;
00348               if (pCD->GData[n][m][l].ID_Comm_world == myID_Comm_world)
00349                 pG->NmyCGrid++;
00350             }
00351 
00352           }}} /* end loops over [n,m,l]: all Grids in Domain[nl+1][ncd] */
00353         }
00354       } /* end loop over all child Domains at level [nl+1] */
00355 
00356 /*------------------------- Allocate CGrid array -----------------------------*/
00357 /* Now we know how many child Grids there are for the Grid in Domain[nd] at
00358  * level nl on both this and other processors.    */
00359 
00360       if (pG->NCGrid > 0) {
00361         pG->CGrid =(GridOvrlpS*)calloc_1d_array(pG->NCGrid, sizeof(GridOvrlpS));
00362         if(pG->CGrid==NULL) ath_error("[init_grid]:failed to allocate CGrid\n");
00363 
00364         for (ncg=0; ncg<pG->NCGrid; ncg++){
00365           for (dim=0; dim<6; dim++) {
00366             pG->CGrid[ncg].myFlx[dim] = NULL;
00367 #ifdef MHD
00368             pG->CGrid[ncg].myEMF1[dim] = NULL;
00369             pG->CGrid[ncg].myEMF2[dim] = NULL;
00370             pG->CGrid[ncg].myEMF3[dim] = NULL;
00371 #endif /* MHD */
00372           }
00373         }
00374       }
00375 
00376 /*-------------------------- Fill in CGrid array -----------------------------*/
00377 /* Repeat loop over all domains at next level, and all the logic to find
00378  * overlapping Grids, to fill in data about overlap regions in CGrid array.
00379  * This isn't a particularly pretty way to do it. */
00380 
00381       nMyCG = 0;
00382       nCG = pG->NmyCGrid;
00383 
00384       for (ncd=0; ncd<pM->DomainsPerLevel[nl+1]; ncd++){
00385         pCD = (DomainS*)&(pM->Domain[nl+1][ncd]);   /* ptr to potential child */
00386 
00387 /* edges of potential child Domain */
00388         for (i=0; i<3; i++) {
00389           D2.ijkl[i] = pCD->Disp[i]/2;
00390           D2.ijkr[i] = 1;
00391           if (pCD->Nx[i] > 1) D2.ijkr[i] = (pCD->Disp[i] + pCD->Nx[i])/2;
00392         }
00393 
00394         isDOverlap = checkOverlap(&D1, &D2, &D3);
00395         if (isDOverlap == 1){
00396 
00397 /*----------------------------------------------------------------------------*/
00398 /* Found the Domain that overlaps, so on the child Domain check if there
00399  * is a Grid that overlaps */
00400 
00401           for (n=0; n<pCD->NGrid[2]; n++){
00402           for (m=0; m<pCD->NGrid[1]; m++){
00403           for (l=0; l<pCD->NGrid[0]; l++){
00404 
00405 /* edges of child Grid */
00406 /* Divide by two to check in units of this (not the child) grid coordinates */
00407             for (i=0; i<3; i++) {
00408               G2.ijkl[i] = pCD->GData[n][m][l].Disp[i]/2;
00409               G2.ijkr[i] = 1;
00410               if (pCD->Nx[i] > 1)
00411                 G2.ijkr[i] = (pCD->GData[n][m][l].Disp[i]
00412                             + pCD->GData[n][m][l].Nx[i])/2;
00413             }
00414 
00415             isGOverlap = checkOverlap(&G1, &G2, &G3);
00416             if (isGOverlap == 1){
00417 
00418 /* If Grid OVERLAPS, then:
00419  * (1) fill-in data in CGrid array */
00420 /* Index CGrid array so that child Grids on this processor come first */
00421 
00422               if (pCD->GData[n][m][l].ID_Comm_world == myID_Comm_world) {
00423                 ncg=nMyCG;
00424                 nMyCG++;
00425               } else {
00426                 ncg=nCG;
00427                 nCG++;
00428               }
00429 
00430               pG->CGrid[ncg].ijks[0] = G3.ijkl[0] - pG->Disp[0] + pG->is;
00431               pG->CGrid[ncg].ijke[0] = G3.ijkr[0] - pG->Disp[0] + pG->is - 1;
00432               pG->CGrid[ncg].ijks[1] = G3.ijkl[1] - pG->Disp[1] + pG->js;
00433               pG->CGrid[ncg].ijke[1] = G3.ijkr[1] - pG->Disp[1] + pG->js - 1;
00434               pG->CGrid[ncg].ijks[2] = G3.ijkl[2] - pG->Disp[2] + pG->ks;
00435               pG->CGrid[ncg].ijke[2] = G3.ijkr[2] - pG->Disp[2] + pG->ks - 1;
00436   
00437               pG->CGrid[ncg].DomN = ncd;
00438               pG->CGrid[ncg].ID = pCD->GData[n][m][l].ID_Comm_Parent;
00439 
00440               n1z = pG->CGrid[ncg].ijke[0] - pG->CGrid[ncg].ijks[0] + 1;
00441               n2z = pG->CGrid[ncg].ijke[1] - pG->CGrid[ncg].ijks[1] + 1;
00442               n3z = pG->CGrid[ncg].ijke[2] - pG->CGrid[ncg].ijks[2] + 1;
00443               pG->CGrid[ncg].nWordsRC = n1z*n2z*n3z*(NVAR);
00444               pG->CGrid[ncg].nWordsP  = 0;
00445 if(myID_Comm_world==0){
00446 printf("\nCGrid overlap is %d x %d x %d\n",n1z,n2z,n3z);
00447 }
00448 #ifdef MHD
00449               if (nDim==3) {
00450                 pG->CGrid[ncg].nWordsRC += 
00451                   (n1z+1)*n2z*n3z + n1z*(n2z+1)*n3z + n1z*n2z*(n3z+1);
00452               } else {
00453                 if (nDim==2) {
00454                   pG->CGrid[ncg].nWordsRC += (n1z+1)*n2z + n1z*(n2z+1);
00455                 }
00456               }
00457 #endif /* MHD */
00458 
00459 /* (2) if edge of the child Grid is at edge of the child Domain, then allocate
00460  * memory for fluxes and EMFs for Correction, and count GZ for Prolongation */
00461 
00462               for (dim=0; dim<nDim; dim++){    /* only checks nDim directions */
00463                 if (dim == 0) iGrid=l;
00464                 if (dim == 1) iGrid=m;
00465                 if (dim == 2) iGrid=n;
00466 
00467 /* inner x1/x2/x3 boundary */
00468 /* First check that edge of child Grid is at edge of overlap (otherwise boundary
00469  * is between MPI blocks in parent, and is internal to child Grid).
00470  * Then check that edge of child Grid is not at l-edge of root (so that physical
00471  * BCs are applied), but is at l-edge of own Domain (so it is not an internal
00472  * MPI boundary on the child Domain). */
00473 
00474                 if ((G2.ijkl[dim] == G3.ijkl[dim]) &&
00475                     (pCD->Disp[dim] != 0) &&
00476                     (iGrid == 0)) {
00477 
00478                   if (dim == 0) {
00479                     n1z = G3.ijkr[1] - G3.ijkl[1];
00480                     n2z = G3.ijkr[2] - G3.ijkl[2];
00481                     n1p = n1z;
00482                     n2p = n2z;
00483                     if (pG->Nx[1] > 1) n1p += nghost + 2;
00484                     if (pG->Nx[2] > 1) n2p += nghost + 2;
00485                   }
00486                   if (dim == 1) {
00487                     n1z = G3.ijkr[0] - G3.ijkl[0];
00488                     n2z = G3.ijkr[2] - G3.ijkl[2];
00489                     n1p = n1z + nghost + 2;
00490                     n2p = n2z;
00491                     if (pG->Nx[2] > 1) n2p += nghost + 2;
00492                   }
00493                   if (dim == 2) {
00494                     n1z = G3.ijkr[0] - G3.ijkl[0];
00495                     n2z = G3.ijkr[1] - G3.ijkl[1];
00496                     n1p = n1z + nghost + 2;
00497                     n2p = n2z + nghost + 2;
00498                   }
00499 
00500                   pG->CGrid[ncg].nWordsRC += n1z*n2z*(NVAR); 
00501                   pG->CGrid[ncg].nWordsP  += ((nghost/2)+2)*n1p*n2p*(NVAR); 
00502 
00503 /* Allocate memory for myFlx and my EMFs */
00504 
00505                   pG->CGrid[ncg].myFlx[2*dim] = (ConsS**)calloc_2d_array(
00506                     n2z,n1z, sizeof(ConsS));
00507                   if(pG->CGrid[ncg].myFlx[2*dim] == NULL) ath_error(
00508                    "[init_grid]:failed to allocate CGrid ixb myFlx\n");
00509 if(myID_Comm_world==0){
00510 printf("Allocated %d x %d array for ixb CGrid.myFlx\n",n2z,n1z);
00511 }
00512 #ifdef MHD
00513                   pG->CGrid[ncg].nWordsP += 6*((nghost/2)+2)*n1p*n2p;
00514 
00515                   if (pG->Nx[1] > 1 && dim != 2) {
00516                     pG->CGrid[ncg].nWordsRC += (n1z+1)*n2z; 
00517                     pG->CGrid[ncg].myEMF3[2*dim] = (Real**)calloc_2d_array(
00518                       n2z,n1z+1, sizeof(Real));
00519                     if(pG->CGrid[ncg].myEMF3[2*dim] == NULL) ath_error(
00520                       "[init_grid]:failed to allocate CGrid ixb myEMF3\n");
00521 if(myID_Comm_world==0){
00522 printf("Allocated %d x %d array for ixb CGrid.myEMF3\n",n2z,n1z+1);
00523 }
00524                   }
00525 
00526                   if (pG->Nx[2] > 1  && dim == 0) {
00527                     pG->CGrid[ncg].nWordsRC += n1z*(n2z+1); 
00528                     pG->CGrid[ncg].myEMF2[2*dim] = (Real**)calloc_2d_array(
00529                       n2z+1,n1z, sizeof(Real));
00530                     if(pG->CGrid[ncg].myEMF2[2*dim] == NULL) ath_error(
00531                       "[init_grid]:failed to allocate CGrid ixb myEMF2\n");
00532 if(myID_Comm_world==0){
00533 printf("Allocated %d x %d array for ixb CGrid.myEMF2\n",n2z+1,n1z);
00534 }
00535                   }
00536 
00537                   if (pG->Nx[2] > 1  && dim == 1) {
00538                     pG->CGrid[ncg].nWordsRC += n1z*(n2z+1); 
00539                     pG->CGrid[ncg].myEMF1[2*dim] = (Real**)calloc_2d_array(
00540                       n2z+1,n1z, sizeof(Real));
00541                     if(pG->CGrid[ncg].myEMF1[2*dim] == NULL) ath_error(
00542                       "[init_grid]:failed to allocate CGrid ixb myEMF1\n");
00543 if(myID_Comm_world==0){
00544 printf("Allocated %d x %d array for ixb CGrid.myEMF1\n",n2z+1,n1z);
00545 }
00546                   }
00547 
00548                   if (pG->Nx[2] > 1  && dim == 2) {
00549                     pG->CGrid[ncg].nWordsRC += n1z*(n2z+1) + (n1z+1)*n2z; 
00550                     pG->CGrid[ncg].myEMF1[2*dim] = (Real**)calloc_2d_array(
00551                       n2z+1,n1z, sizeof(Real));
00552                     if(pG->CGrid[ncg].myEMF1[2*dim] == NULL) ath_error(
00553                       "[init_grid]:failed to allocate CGrid ixb myEMF1\n");
00554 if(myID_Comm_world==0){
00555 printf("Allocated %d x %d array for ixb CGrid.myEMF1\n",n2z+1,n1z);
00556 }
00557                     pG->CGrid[ncg].myEMF2[2*dim] = (Real**)calloc_2d_array(
00558                       n2z,n1z+1, sizeof(Real));
00559                     if(pG->CGrid[ncg].myEMF2[2*dim] == NULL) ath_error(
00560                       "[init_grid]:failed to allocate CGrid ixb myEMF2\n");
00561 if(myID_Comm_world==0){
00562 printf("Allocated %d x %d array for ixb CGrid.myEMF2\n",n2z,n1z+1);
00563 }
00564                   }
00565 #endif /* MHD */
00566 
00567                 }
00568 
00569 /* outer x1/x2/x3 boundary */
00570 /* First check that edge of child Grid is at edge of overlap (otherwise boundary
00571  * is between MPI blocks in parent, and is internal to child Grid).
00572  * Then check that edge of child Grid is not at r-edge of root (so that physical
00573  * BCs are applied), but is at r-edge of own Domain (so it is not an internal
00574  * MPI boundary on the child Domain). */
00575 
00576                 irefine = 1;
00577                 for (i=1;i<=(nl+1);i++) irefine *= 2; /* child refinement lev */
00578                 if ( (G2.ijkr[dim] == G3.ijkr[dim]) &&
00579                     ((pCD->Disp[dim] + pCD->Nx[dim])/irefine != pM->Nx[dim]) &&
00580                      (iGrid == (pCD->NGrid[dim]-1)) ) {
00581 
00582                   if (dim == 0) {
00583                     n1z = G3.ijkr[1] - G3.ijkl[1];
00584                     n2z = G3.ijkr[2] - G3.ijkl[2];
00585                     n1p = n1z;
00586                     n2p = n2z;
00587                     if (pG->Nx[1] > 1) n1p += nghost + 2;
00588                     if (pG->Nx[2] > 1) n2p += nghost + 2;
00589                   }
00590                   if (dim == 1) {
00591                     n1z = G3.ijkr[0] - G3.ijkl[0];
00592                     n2z = G3.ijkr[2] - G3.ijkl[2];
00593                     n1p = n1z + nghost + 2;
00594                     n2p = n2z;
00595                     if (pG->Nx[2] > 1) n2p += nghost + 2;
00596                   }
00597                   if (dim == 2) {
00598                     n1z = G3.ijkr[0] - G3.ijkl[0];
00599                     n2z = G3.ijkr[1] - G3.ijkl[1];
00600                     n1p = n1z + nghost + 2;
00601                     n2p = n2z + nghost + 2;
00602                   }
00603 
00604                   pG->CGrid[ncg].nWordsRC += n1z*n2z*(NVAR); 
00605                   pG->CGrid[ncg].nWordsP  += ((nghost/2)+2)*n1p*n2p*(NVAR); 
00606 
00607 /* Allocate memory for myFlx and myEMFs*/
00608 
00609                   pG->CGrid[ncg].myFlx[(2*dim)+1] = (ConsS**)calloc_2d_array(
00610                     n2z,n1z, sizeof(ConsS));
00611                   if(pG->CGrid[ncg].myFlx[(2*dim)+1] == NULL) ath_error(
00612                     "[init_grid]:failed to allocate CGrid oxb myFlx\n");
00613 if(myID_Comm_world==0){
00614 printf("Allocated %d x %d array for oxb CGrid.myFlx\n",n2z,n1z);
00615 }
00616 #ifdef MHD
00617                   pG->CGrid[ncg].nWordsP += 6*((nghost/2)+2)*n1p*n2p;
00618 
00619                   if (pG->Nx[1] > 1 && dim != 2) {
00620                     pG->CGrid[ncg].nWordsRC += (n1z+1)*n2z;
00621                     pG->CGrid[ncg].myEMF3[(2*dim)+1] =(Real**)calloc_2d_array(
00622                       n2z,n1z+1, sizeof(Real));
00623                     if(pG->CGrid[ncg].myEMF3[(2*dim)+1] == NULL) ath_error(
00624                       "[init_grid]:failed to allocate CGrid oxb myEMF3\n");
00625 if(myID_Comm_world==0){
00626 printf("Allocated %d x %d array for oxb CGrid.myEMF3\n",n2z,n1z+1);
00627 }
00628                   }
00629 
00630                   if (pG->Nx[2] > 1  && dim == 0) {
00631                     pG->CGrid[ncg].nWordsRC += n1z*(n2z+1);
00632                     pG->CGrid[ncg].myEMF2[(2*dim)+1] =(Real**)calloc_2d_array(
00633                       n2z+1,n1z, sizeof(Real));
00634                     if(pG->CGrid[ncg].myEMF2[(2*dim)+1] == NULL) ath_error(
00635                       "[init_grid]:failed to allocate CGrid oxb myEMF2\n");
00636 if(myID_Comm_world==0){
00637 printf("Allocated %d x %d array for oxb CGrid.myEMF2\n",n2z+1,n1z);
00638 }
00639                   }
00640 
00641                   if (pG->Nx[2] > 1  && dim == 1) {
00642                     pG->CGrid[ncg].nWordsRC += n1z*(n2z+1);
00643                     pG->CGrid[ncg].myEMF1[(2*dim)+1] =(Real**)calloc_2d_array(
00644                       n2z+1,n1z, sizeof(Real));
00645                     if(pG->CGrid[ncg].myEMF1[(2*dim)+1] == NULL) ath_error(
00646                       "[init_grid]:failed to allocate CGrid oxb myEMF1\n");
00647 if(myID_Comm_world==0){
00648 printf("Allocated %d x %d array for oxb CGrid.myEMF1\n",n2z+1,n1z);
00649 }
00650                   }
00651 
00652                   if (pG->Nx[2] > 1  && dim == 2) {
00653                     pG->CGrid[ncg].nWordsRC += n1z*(n2z+1) + (n1z+1)*n2z;
00654                     pG->CGrid[ncg].myEMF1[(2*dim)+1] =(Real**)calloc_2d_array(
00655                       n2z+1,n1z, sizeof(Real));
00656                     if(pG->CGrid[ncg].myEMF1[(2*dim)+1] == NULL) ath_error(
00657                       "[init_grid]:failed to allocate CGrid oxb myEMF1\n");
00658 if(myID_Comm_world==0){
00659 printf("Allocated %d x %d array for oxb CGrid.myEMF1\n",n2z+1,n1z);
00660 }
00661                     pG->CGrid[ncg].myEMF2[(2*dim)+1] =(Real**)calloc_2d_array(
00662                       n2z,n1z+1, sizeof(Real));
00663                     if(pG->CGrid[ncg].myEMF2[(2*dim)+1] == NULL) ath_error(
00664                       "[init_grid]:failed to allocate CGrid oxb myEMF2\n");
00665 if(myID_Comm_world==0){
00666 printf("Allocated %d x %d array for oxb CGrid.myEMF2\n",n2z,n1z+1);
00667 }
00668                   }
00669 #endif /* MHD */
00670 
00671                 }
00672 
00673               } /* end loop over boundary directions */
00674             } /* end if child Grid overlaps or touches */
00675           }}} /* end loops over all Grids in Domain[nl+1][ncd] */
00676         } /* end if child Domain */
00677       } /* end loop over all Domains at level [nl+1] */
00678     } /* end if Grid on this processor */
00679   }} /* end loops over all levels and domains per level */
00680 
00681 
00682 /************************************/
00683 for (nl=0; nl<(pM->NLevels); nl++){
00684 for (nd=0; nd<pM->DomainsPerLevel[nl]; nd++){
00685   if (pM->Domain[nl][nd].Grid != NULL) {
00686   pG = pM->Domain[nl][nd].Grid;          /* set ptr to this Grid */
00687 
00688 printf("\nProcID=%d level=%d Domain=%d NCgrid=%d NmyCGrid=%d\n",
00689   myID_Comm_world,nl,nd,pG->NCGrid,pG->NmyCGrid);
00690 for (i=0;i<pG->NCGrid; i++){
00691 printf("CGrid=%d, [is,ie,js,je,ks,ke]=%d %d %d %d %d %d\n",i,
00692   pG->CGrid[i].ijks[0],pG->CGrid[i].ijke[0],
00693   pG->CGrid[i].ijks[1],pG->CGrid[i].ijke[1],
00694   pG->CGrid[i].ijks[2],pG->CGrid[i].ijke[2]);
00695 printf("Child_ID=%d DomN=%d nWordsRC=%d nWordsP=%d\n",
00696   pG->CGrid[i].ID,pG->CGrid[i].DomN,pG->CGrid[i].nWordsRC,
00697   pG->CGrid[i].nWordsP);
00698 }
00699 
00700 }}}
00701 /************************************/
00702 
00703 
00704 /*--------------- Count number of parent Grids, and boundaries ---------------*/
00705 /* Now we have to count the number of parent Grids, and fine/coarse boundaries
00706  * with parent Grids, including parent Grids on this and other processors,
00707  * before we can allocate the PGrid array.  This is a shortcoming of making
00708  * these structures 1D arrays. */
00709 /* Loop over levels (except root level=0), and domains per level */
00710 
00711   for (nl=1; nl<(pM->NLevels); nl++){
00712   for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00713     if (pM->Domain[nl][nd].Grid != NULL) {
00714       pD = (DomainS*)&(pM->Domain[nl][nd]);  /* set ptr to Domain */
00715       pG = pM->Domain[nl][nd].Grid;          /* set ptr to Grid */
00716       get_myGridIndex(pD, myID_Comm_world, &myL, &myM, &myN);
00717 
00718 /* edges of this Domain */
00719       for (i=0; i<3; i++) {
00720         D1.ijkl[i] = pD->Disp[i];
00721         D1.ijkr[i] = pD->Disp[i] + pD->Nx[i];
00722       }
00723 
00724 /* edges of this Grid */
00725       for (i=0; i<3; i++) {
00726         G1.ijkl[i] = pG->Disp[i];
00727         G1.ijkr[i] = pG->Disp[i] + pG->Nx[i];
00728       }
00729 
00730 /* For this domain, find domain at last level that overlaps (parent Domain).
00731  * Multiply by two to check in units of this (not the parent) domain coordinates
00732  */
00733 
00734       for (npd=0; npd<(pM->DomainsPerLevel[nl-1]); npd++){
00735         pPD = (DomainS*)&(pM->Domain[nl-1][npd]); /* ptr to potential parent */
00736 
00737 /* edges of potential parent Domain */
00738         for (i=0; i<3; i++) {
00739           D2.ijkl[i] = 2*(pPD->Disp[i]);
00740           D2.ijkr[i] = 2*(pPD->Disp[i] + pPD->Nx[i]);
00741         }
00742 
00743         isDOverlap = checkOverlap(&D1, &D2, &D3);
00744         if (isDOverlap == 1){
00745 
00746 /*----------------------------------------------------------------------------*/
00747 /* Found parent Domain that overlaps.  So on the parent Domain, find all the
00748  * Grids that overlap this Grid. */
00749 
00750           for (n=0; n<pPD->NGrid[2]; n++){
00751           for (m=0; m<pPD->NGrid[1]; m++){
00752           for (l=0; l<pPD->NGrid[0]; l++){
00753 
00754 /* edges of parent Grid */
00755 /* Multiply by two to check in units of this (not the parent) grid coord */
00756             for (i=0; i<3; i++) {
00757               G2.ijkl[i] = 2*(pPD->GData[n][m][l].Disp[i]);
00758               G2.ijkr[i] = 2*(pPD->GData[n][m][l].Disp[i] + 
00759                               pPD->GData[n][m][l].Nx[i]);
00760             }
00761 
00762             isGOverlap = checkOverlap(&G1, &G2, &G3);
00763 
00764 /* If Grid overlaps, increment Parent counters */
00765 
00766             if (isGOverlap == 1){
00767               pG->NPGrid++;
00768               if (pPD->GData[n][m][l].ID_Comm_world == myID_Comm_world)
00769                 pG->NmyPGrid++;
00770             }
00771 
00772           }}} /* end loops over [n,m,l]: all Grids in Domain[nl-1][npd] */
00773         }
00774       } /* end loop over all parent Domains at level [nl-1] */
00775 
00776 /*-------------------------- Allocate PGrid array ----------------------------*/
00777 /* Now we know how many parent Grids there are for the Grid in Domain[nd] at
00778  * level nl on both this and other processors. */
00779 
00780       if (pG->NPGrid > 0) {
00781         pG->PGrid =(GridOvrlpS*)calloc_1d_array(pG->NPGrid, sizeof(GridOvrlpS));
00782         if(pG->PGrid==NULL) ath_error("[init_grid]:failed to allocate PGrid\n");
00783 
00784         for (npg=0; npg<pG->NPGrid; npg++){
00785           for (dim=0; dim<6; dim++) {
00786             pG->PGrid[npg].myFlx[dim] = NULL;
00787 #ifdef MHD
00788             pG->PGrid[npg].myEMF1[dim] = NULL;
00789             pG->PGrid[npg].myEMF2[dim] = NULL;
00790             pG->PGrid[npg].myEMF3[dim] = NULL;
00791 #endif /* MHD */
00792           }
00793         }
00794       }
00795 
00796 /*--------------------------- Fill in PGrid array ----------------------------*/
00797 /* Repeat loop over all domains at last level, and all the logic to find
00798  * overlapping Grids, to fill in data about overlap regions in PGrid array */
00799 
00800       nMyPG = 0;
00801       nPG = pG->NmyPGrid;
00802 
00803       for (npd=0; npd<pM->DomainsPerLevel[nl-1]; npd++){
00804         pPD = (DomainS*)&(pM->Domain[nl-1][npd]);   /* ptr to potential parent*/
00805 
00806 /* edges of potential parent Domain */
00807         for (i=0; i<3; i++) {
00808           D2.ijkl[i] = 2*(pPD->Disp[i]);
00809           D2.ijkr[i] = 2*(pPD->Disp[i] + pPD->Nx[i]);
00810         }
00811 
00812         isDOverlap = checkOverlap(&D1, &D2, &D3);
00813         if (isDOverlap == 1){
00814 
00815 /*----------------------------------------------------------------------------*/
00816 /* Found the Domain that overlaps, so on the parent Domain check if there is a
00817  * Grid that overlaps */
00818 
00819           for (n=0; n<pPD->NGrid[2]; n++){
00820           for (m=0; m<pPD->NGrid[1]; m++){
00821           for (l=0; l<pPD->NGrid[0]; l++){
00822 
00823 /* edges of parent Grid */
00824 /* Multiply by two to check in units of this (not the parent) grid coord */
00825             for (i=0; i<3; i++) {
00826               G2.ijkl[i] = 2*(pPD->GData[n][m][l].Disp[i]);
00827               G2.ijkr[i] = 2*(pPD->GData[n][m][l].Disp[i] +
00828                               pPD->GData[n][m][l].Nx[i]);
00829             }
00830 
00831             isGOverlap = checkOverlap(&G1, &G2, &G3);
00832 
00833             if (isGOverlap == 1){
00834 
00835 /* If Grid OVERLAPS, then:
00836  * (1) fill-in data in PGrid array */
00837 /* Index PGrid array so that parent Grids on this processor come first */
00838 
00839               if (pPD->GData[n][m][l].ID_Comm_world == myID_Comm_world) {
00840                 npg=nMyPG;
00841                 nMyPG++;
00842               } else {
00843                 npg=nPG;
00844                 nPG++;
00845               }
00846 
00847               pG->PGrid[npg].ijks[0] = G3.ijkl[0] - pG->Disp[0] + pG->is;
00848               pG->PGrid[npg].ijke[0] = G3.ijkr[0] - pG->Disp[0] + pG->is - 1;
00849               pG->PGrid[npg].ijks[1] = G3.ijkl[1] - pG->Disp[1] + pG->js;
00850               pG->PGrid[npg].ijke[1] = G3.ijkr[1] - pG->Disp[1] + pG->js - 1;
00851               pG->PGrid[npg].ijks[2] = G3.ijkl[2] - pG->Disp[2] + pG->ks;
00852               pG->PGrid[npg].ijke[2] = G3.ijkr[2] - pG->Disp[2] + pG->ks - 1;
00853   
00854               pG->PGrid[npg].DomN = npd;
00855               pG->PGrid[npg].ID = pPD->GData[n][m][l].ID_Comm_Children;
00856 
00857               n1z =    (pG->PGrid[npg].ijke[0]-pG->PGrid[npg].ijks[0] + 1)/2;
00858               n2z =MAX((pG->PGrid[npg].ijke[1]-pG->PGrid[npg].ijks[1] + 1)/2,1);
00859               n3z =MAX((pG->PGrid[npg].ijke[2]-pG->PGrid[npg].ijks[2] + 1)/2,1);
00860 /*
00861               if (pG->Nx[1]>1) n2z /= 2;
00862               if (pG->Nx[2]>1) n3z /= 2;
00863 */
00864               pG->PGrid[npg].nWordsRC = n1z*n2z*n3z*(NVAR);
00865               pG->PGrid[npg].nWordsP  = 0;
00866 #ifdef MHD
00867               if (pG->Nx[2]>1) {
00868                 pG->PGrid[npg].nWordsRC += 
00869                   (n1z+1)*n2z*n3z + n1z*(n2z+1)*n3z + n1z*n2z*(n3z+1);
00870               } else {
00871                 if (pG->Nx[1]>1) {
00872                   pG->PGrid[npg].nWordsRC += (n1z+1)*n2z + n1z*(n2z+1);
00873                 }
00874               }
00875 #endif /* MHD */
00876 
00877 /* (2) If any edge of this Grid is at the edge of this Domain, then allocate
00878  * memory for fluxes and EMFs for Correction step, and count GZ data passed in
00879  * Prolongation step. */
00880 
00881               for (dim=0; dim<nDim; dim++){    /* only checks nDim directions */
00882                 if (dim == 0) iGrid=myL;
00883                 if (dim == 1) iGrid=myM;
00884                 if (dim == 2) iGrid=myN;
00885 
00886 /* inner x1/x2/x3 boundary */
00887 /* First check that edge of this Grid is at edge of overlap (otherwise boundary
00888  * is between MPI blocks in parent, and is internal to child Grid).
00889  * Then check that edge of this Grid is not at l-edge of root (so that physical
00890  * BCs are applied), but is at l-edge of own Domain (so it is not an internal
00891  * MPI boundary on this Domain). */
00892 
00893                   if ((G1.ijkl[dim] == G3.ijkl[dim]) &&
00894                       (pD->Disp[dim] != 0) &&
00895                       (iGrid == 0)) {
00896 
00897                     if (dim == 0) {
00898                       n1z = G3.ijkr[1] - G3.ijkl[1];
00899                       n2z = G3.ijkr[2] - G3.ijkl[2];
00900                       n1p = MAX((n1z/2),1);
00901                       n2p = MAX((n2z/2),1);
00902                       if (pG->Nx[1] > 1) n1p += nghost + 2;
00903                       if (pG->Nx[2] > 1) n2p += nghost + 2;
00904                     }
00905                     if (dim == 1) {
00906                       n1z = G3.ijkr[0] - G3.ijkl[0];
00907                       n2z = G3.ijkr[2] - G3.ijkl[2];
00908                       n1p = (n1z/2) + nghost + 2;
00909                       n2p = MAX((n2z/2),1);
00910                       if (pG->Nx[2] > 1) n2p += nghost + 2;
00911                     }
00912                     if (dim == 2) {
00913                       n1z = G3.ijkr[0] - G3.ijkl[0];
00914                       n2z = G3.ijkr[1] - G3.ijkl[1];
00915                       n1p = (n1z/2) + nghost + 2;
00916                       n2p = (n2z/2) + nghost + 2;
00917                     }
00918 
00919                     pG->PGrid[npg].nWordsRC += 
00920                       MAX((n1z/2),1)*MAX((n2z/2),1)*(NVAR);
00921                     pG->PGrid[npg].nWordsP  += ((nghost/2)+2)*n1p*n2p*(NVAR);
00922 
00923 /* Allocate memory for myFlx and my EMFS.  Note they have dimension of the
00924  * parent overlap on THIS Grid, which is 2x the transverse dimension of the
00925  * overlap on the parent Grid (the actual number of words sent). */
00926 
00927                     pG->PGrid[npg].myFlx[2*dim] = (ConsS**)calloc_2d_array(
00928                       n2z,n1z, sizeof(ConsS));
00929                     if(pG->PGrid[npg].myFlx[2*dim] == NULL) ath_error(
00930                       "[init_grid]:failed to allocate PGrid ixb myFlx\n");
00931 if(myID_Comm_world==0){
00932 printf("Allocated %d x %d array for ixb PGrid.myFlx\n",n2z,n1z);
00933 }
00934 #ifdef MHD
00935                     pG->PGrid[npg].nWordsP += 6*((nghost/2)+2)*n1p*n2p;
00936 
00937                     if (pG->Nx[1] > 1 && dim != 2) {
00938                       pG->PGrid[npg].nWordsRC +=(n1z/2+1)*MAX((n2z/2),1);
00939                       pG->PGrid[npg].myEMF3[2*dim] = (Real**)calloc_2d_array(
00940                         n2z,n1z+1, sizeof(Real));
00941                       if(pG->PGrid[npg].myEMF3[2*dim]==NULL) ath_error(
00942                         "[init_grid]:failed to allocate PGrid ixb myEMF3\n");
00943 if(myID_Comm_world==0){
00944 printf("Allocated %d x %d array for ixb PGrid.myEMF3\n",n2z,n1z+1);
00945 }
00946                     }
00947 
00948                     if (pG->Nx[2] > 1  && dim == 0) {
00949                       pG->PGrid[npg].nWordsRC += (n1z/2)*(n2z/2+1);
00950                       pG->PGrid[npg].myEMF2[2*dim] = (Real**)calloc_2d_array(
00951                         n2z+1,n1z, sizeof(Real));
00952                       if(pG->PGrid[npg].myEMF2[2*dim]==NULL) ath_error(
00953                         "[init_grid]:failed to allocate PGrid ixb myEMF2\n");
00954 if(myID_Comm_world==0){
00955 printf("Allocated %d x %d array for ixb PGrid.myEMF2\n",n2z+1,n1z);
00956 }
00957                     }
00958 
00959                     if (pG->Nx[2] > 1  && dim == 1) {
00960                       pG->PGrid[npg].nWordsRC += (n1z/2)*(n2z/2+1);
00961                       pG->PGrid[npg].myEMF1[2*dim] = (Real**)calloc_2d_array(
00962                         n2z+1,n1z, sizeof(Real));
00963                       if(pG->PGrid[npg].myEMF1[2*dim]==NULL) ath_error(
00964                         "[init_grid]:failed to allocate PGrid ixb myEMF1\n");
00965 if(myID_Comm_world==0){
00966 printf("Allocated %d x %d array for ixb PGrid.myEMF1\n",n2z+1,n1z);
00967 }
00968                     }
00969 
00970                     if (pG->Nx[2] > 1  && dim == 2) {
00971                       pG->PGrid[npg].nWordsRC += (n1z/2)*(n2z/2+1) 
00972                                              + (n1z/2+1)*(n2z/2);
00973                       pG->PGrid[npg].myEMF1[2*dim] = (Real**)calloc_2d_array(
00974                         n2z+1,n1z, sizeof(Real));
00975                       if(pG->PGrid[npg].myEMF1[2*dim]==NULL) ath_error(
00976                         "[init_grid]:failed to allocate PGrid ixb myEMF1\n");
00977 if(myID_Comm_world==0){
00978 printf("Allocated %d x %d array for ixb PGrid.myEMF1\n",n2z+1,n1z);
00979 }
00980                       pG->PGrid[npg].myEMF2[2*dim] = (Real**)calloc_2d_array(
00981                         n2z,n1z+1, sizeof(Real));
00982                       if(pG->PGrid[npg].myEMF2[2*dim]==NULL) ath_error(
00983                         "[init_grid]:failed to allocate PGrid ixb myEMF2\n");
00984 if(myID_Comm_world==0){
00985 printf("Allocated %d x %d array for ixb PGrid.myEMF2\n",n2z,n1z+1);
00986 }
00987                     }
00988 #endif /* MHD */
00989 
00990                 }
00991 
00992 /* outer x1/x2/x3 boundary */
00993 /* First check that edge of this Grid is at edge of overlap (otherwise boundary
00994  * is between MPI blocks in parent, and is internal to child Grid).
00995  * Then check that edge of this Grid is not at r-edge of root (so that physical
00996  * BCs are applied), but is at r-edge of own Domain (so it is not an internal
00997  * MPI boundary on this Domain). */
00998 
00999                 irefine = 1;
01000                 for (i=1;i<=nl;i++) irefine *= 2; /* this level refinement */
01001                 if ( (G1.ijkr[dim] == G3.ijkr[dim]) &&
01002                     ((pD->Disp[dim] + pD->Nx[dim])/irefine != pM->Nx[dim]) &&
01003                      (iGrid == (pD->NGrid[dim]-1)) ) {
01004   
01005                     if (dim == 0) {
01006                       n1z = G3.ijkr[1] - G3.ijkl[1];
01007                       n2z = G3.ijkr[2] - G3.ijkl[2];
01008                       n1p = MAX((n1z/2),1);
01009                       n2p = MAX((n2z/2),1);
01010                       if (pG->Nx[1] > 1) n1p += nghost + 2;
01011                       if (pG->Nx[2] > 1) n2p += nghost + 2;
01012                     }
01013                     if (dim == 1) {
01014                       n1z = G3.ijkr[0] - G3.ijkl[0];
01015                       n2z = G3.ijkr[2] - G3.ijkl[2];
01016                       n1p = (n1z/2) + nghost + 2;
01017                       n2p = MAX((n2z/2),1);
01018                       if (pG->Nx[2] > 1) n2p += nghost + 2;
01019                     }
01020                     if (dim == 2) {
01021                       n1z = G3.ijkr[0] - G3.ijkl[0];
01022                       n2z = G3.ijkr[1] - G3.ijkl[1];
01023                       n1p = (n1z/2) + nghost + 2;
01024                       n2p = (n2z/2) + nghost + 2;
01025                     }
01026 
01027                     pG->PGrid[npg].nWordsRC += 
01028                       MAX((n1z/2),1)*MAX((n2z/2),1)*(NVAR);
01029                     pG->PGrid[npg].nWordsP  += ((nghost/2)+2)*n1p*n2p*(NVAR);
01030 
01031 /* Allocate memory for myFlx and my EMFS.  Note they have dimension of the
01032  * parent overlap on THIS Grid, which is 2x the transverse dimension of the
01033  * overlap on the parent Grid (the actual number of words sent). */
01034 
01035                     pG->PGrid[npg].myFlx[(2*dim)+1] =
01036                       (ConsS**)calloc_2d_array(n2z,n1z, sizeof(ConsS));
01037                     if(pG->PGrid[npg].myFlx[(2*dim)+1] == NULL) ath_error(
01038                       "[init_grid]:failed to allocate PGrid oxb myFlx\n");
01039 if(myID_Comm_world==0){
01040 printf("Allocated %d x %d array for oxb PGrid.myFlx\n",n2z,n1z);
01041 }
01042 #ifdef MHD
01043                     pG->PGrid[npg].nWordsP += 6*((nghost/2)+2)*n1p*n2p;
01044 
01045                     if (pG->Nx[1] > 1 && dim != 2) {
01046                       pG->PGrid[npg].nWordsRC += (n1z/2+1)*MAX(n2z/2,1);
01047                       pG->PGrid[npg].myEMF3[(2*dim)+1] =(Real**)calloc_2d_array(
01048                         n2z,n1z+1, sizeof(Real));
01049                       if(pG->PGrid[npg].myEMF3[(2*dim)+1]==NULL) ath_error(
01050                         "[init_grid]:failed to allocate PGrid oxb myEMF3\n");
01051 if(myID_Comm_world==0){
01052 printf("Allocated %d x %d array for oxb PGrid.myEMF3\n",n2z,n1z+1);
01053 }
01054                     }
01055 
01056                     if (pG->Nx[2] > 1  && dim == 0) {
01057                       pG->PGrid[npg].nWordsRC += (n1z/2)*(n2z/2+1);
01058                       pG->PGrid[npg].myEMF2[(2*dim)+1] =(Real**)calloc_2d_array(
01059                         n2z+1,n1z, sizeof(Real));
01060                       if(pG->PGrid[npg].myEMF2[(2*dim)+1]==NULL) ath_error(
01061                         "[init_grid]:failed to allocate PGrid oxb myEMF2\n");
01062 if(myID_Comm_world==0){
01063 printf("Allocated %d x %d array for oxb PGrid.myEMF2\n",n2z+1,n1z);
01064 }
01065                     }
01066 
01067                     if (pG->Nx[2] > 1  && dim == 1) {
01068                       pG->PGrid[npg].nWordsRC += (n1z/2)*(n2z/2+1);
01069                       pG->PGrid[npg].myEMF1[(2*dim)+1] =(Real**)calloc_2d_array(
01070                         n2z+1,n1z, sizeof(Real));
01071                       if(pG->PGrid[npg].myEMF1[(2*dim)+1]==NULL) ath_error(
01072                         "[init_grid]:failed to allocate PGrid oxb myEMF1\n");
01073 if(myID_Comm_world==0){
01074 printf("Allocated %d x %d array for oxb PGrid.myEMF1\n",n2z+1,n1z);
01075 }
01076                     }
01077 
01078                     if (pG->Nx[2] > 1  && dim == 2) {
01079                       pG->PGrid[npg].nWordsRC += (n1z/2)*(n2z/2+1) 
01080                                                + (n1z/2+1)*(n2z/2);
01081                       pG->PGrid[npg].myEMF1[(2*dim)+1] =(Real**)calloc_2d_array(
01082                         n2z+1,n1z, sizeof(Real));
01083                       if(pG->PGrid[npg].myEMF1[(2*dim)+1]==NULL) ath_error(
01084                         "[init_grid]:failed to allocate PGrid oxb myEMF1\n");
01085 if(myID_Comm_world==0){
01086 printf("Allocated %d x %d array for oxb PGrid.myEMF1\n",n2z+1,n1z);
01087 }
01088                       pG->PGrid[npg].myEMF2[(2*dim)+1] =(Real**)calloc_2d_array(
01089                         n2z,n1z+1, sizeof(Real));
01090                       if(pG->PGrid[npg].myEMF2[(2*dim)+1]==NULL) ath_error(
01091                         "[init_grid]:failed to allocate PGrid oxb myEMF2\n");
01092 if(myID_Comm_world==0){
01093 printf("Allocated %d x %d array for oxb PGrid.myEMF2\n",n2z,n1z+1);
01094 }
01095                     }
01096 #endif /* MHD */
01097 
01098                 }
01099 
01100               } /* end loop over boundary directions */
01101             } /* end if parent Grid overlaps or touches */
01102           }}} /* end loops over all Grids in Domain[nl-1][npd] */
01103         } /* end if parent Domain */
01104       } /* end loop over all Domains at level [nl-1] */
01105     } /* end if Grid on this processor */
01106   }} /* end loops over all levels and domains per level */
01107 
01108 /************************************/
01109 for (nl=0; nl<(pM->NLevels); nl++){
01110 for (nd=0; nd<pM->DomainsPerLevel[nl]; nd++){
01111   if (pM->Domain[nl][nd].Grid != NULL) {
01112   pG = pM->Domain[nl][nd].Grid;          /* set ptr to this Grid */
01113 
01114 printf("\nProcID=%d level=%d Domain=%d NPgrid=%d NmyPGrid=%d\n",
01115   myID_Comm_world,nl,nd,pG->NPGrid,pG->NmyPGrid);
01116 for (i=0;i<pG->NPGrid; i++){
01117 printf("PGrid=%d, [is,ie,js,je,ks,ke]=%d %d %d %d %d %d\n",i,
01118   pG->PGrid[i].ijks[0],pG->PGrid[i].ijke[0],
01119   pG->PGrid[i].ijks[1],pG->PGrid[i].ijke[1],
01120   pG->PGrid[i].ijks[2],pG->PGrid[i].ijke[2]);
01121 printf("Parent_ID=%d DomN=%d nWordsRC=%d nWordsP=%d\n",
01122   pG->PGrid[i].ID,pG->PGrid[i].DomN,pG->PGrid[i].nWordsRC,
01123   pG->PGrid[i].nWordsP);
01124 }
01125 
01126 }}}
01127 /************************************/
01128 
01129 
01130 #endif /* STATIC_MESH_REFINEMENT */
01131 
01132   return;
01133 
01134 /*--- Error messages ---------------------------------------------------------*/
01135 
01136 #ifdef CYLINDRICAL
01137   on_error15:
01138     free_1d_array(pG->ri);
01139   on_error14:
01140     free_1d_array(pG->r);
01141 #endif
01142 #ifdef SELF_GRAVITY
01143   on_error13:
01144     free_3d_array(pG->x3MassFlux);
01145   on_error12:
01146     free_3d_array(pG->x2MassFlux);
01147   on_error11:
01148     free_3d_array(pG->x1MassFlux);
01149   on_error10:
01150     free_3d_array(pG->Phi_old);
01151   on_error9:
01152     free_3d_array(pG->Phi);
01153 #endif
01154 #ifdef RESISTIVITY
01155   on_error7:
01156     free_3d_array(pG->eta_AD);
01157   on_error6:
01158     free_3d_array(pG->eta_Hall);
01159   on_error5:
01160     free_3d_array(pG->eta_Ohm);
01161 #endif
01162 #ifdef MHD
01163   on_error4:
01164     free_3d_array(pG->B3i);
01165   on_error3:
01166     free_3d_array(pG->B2i);
01167   on_error2:
01168     free_3d_array(pG->B1i);
01169 #endif
01170   on_error1:
01171     free_3d_array(pG->U);
01172     ath_error("[init_grid]: Error allocating memory\n");
01173 }
01174 
01175 #ifdef STATIC_MESH_REFINEMENT
01176 /*=========================== PRIVATE FUNCTIONS ==============================*/
01177 /*----------------------------------------------------------------------------*/
01178 /*! \fn int checkOverlap(SideS *pC1, SideS *pC2, SideS *pC3)
01179  *  \brief Checks if two cubes are overlapping.
01180  *
01181  *  - If yes returns true and sides of overlap region in Cube3
01182  *  - If no  returns false and -1 in Cube3
01183  *
01184  * Arguments are Side structures, containing indices of the 6 sides of cube */
01185 
01186 int checkOverlap(SideS *pC1, SideS *pC2, SideS *pC3)
01187 {
01188   int isOverlap=0;
01189 
01190   if (pC1->ijkl[0] < pC2->ijkr[0] && pC1->ijkr[0] > pC2->ijkl[0] &&
01191       pC1->ijkl[1] < pC2->ijkr[1] && pC1->ijkr[1] > pC2->ijkl[1] &&
01192       pC1->ijkl[2] < pC2->ijkr[2] && pC1->ijkr[2] > pC2->ijkl[2]) isOverlap=1;
01193 
01194   if (isOverlap==1) {
01195     pC3->ijkl[0] = MAX(pC1->ijkl[0], pC2->ijkl[0]);
01196     pC3->ijkr[0] = MIN(pC1->ijkr[0], pC2->ijkr[0]);
01197     pC3->ijkl[1] = MAX(pC1->ijkl[1], pC2->ijkl[1]);
01198     pC3->ijkr[1] = MIN(pC1->ijkr[1], pC2->ijkr[1]);
01199     pC3->ijkl[2] = MAX(pC1->ijkl[2], pC2->ijkl[2]);
01200     pC3->ijkr[2] = MIN(pC1->ijkr[2], pC2->ijkr[2]);
01201   } else {
01202     pC3->ijkl[0] = -1;
01203     pC3->ijkr[0] = -1;
01204     pC3->ijkl[1] = -1;
01205     pC3->ijkr[1] = -1;
01206     pC3->ijkl[2] = -1;
01207     pC3->ijkr[2] = -1;
01208   }
01209 
01210   return isOverlap;
01211 }
01212 
01213 /*----------------------------------------------------------------------------*/
01214 /*! \fn int checkOverlapTouch(SideS *pC1, SideS *pC2, SideS *pC3)
01215  *  \brief Checks if two cubes are overlapping or touching.
01216  *
01217  *  - If yes returns true and sides of overlap region in Cube3
01218  *  - If no  returns false and -1 in Cube3
01219  *
01220  * Arguments are Side structures, containing indices of the 6 sides of cube */
01221 
01222 int checkOverlapTouch(SideS *pC1, SideS *pC2, SideS *pC3)
01223 {
01224   int isOverlap=0;
01225 
01226   if (pC1->ijkl[0] <= pC2->ijkr[0] && pC1->ijkr[0] >= pC2->ijkl[0] &&
01227       pC1->ijkl[1] <= pC2->ijkr[1] && pC1->ijkr[1] >= pC2->ijkl[1] &&
01228       pC1->ijkl[2] <= pC2->ijkr[2] && pC1->ijkr[2] >= pC2->ijkl[2]) isOverlap=1;
01229 
01230   if (isOverlap==1) {
01231     pC3->ijkl[0] = MAX(pC1->ijkl[0], pC2->ijkl[0]);
01232     pC3->ijkr[0] = MIN(pC1->ijkr[0], pC2->ijkr[0]);
01233     pC3->ijkl[1] = MAX(pC1->ijkl[1], pC2->ijkl[1]);
01234     pC3->ijkr[1] = MIN(pC1->ijkr[1], pC2->ijkr[1]);
01235     pC3->ijkl[2] = MAX(pC1->ijkl[2], pC2->ijkl[2]);
01236     pC3->ijkr[2] = MIN(pC1->ijkr[2], pC2->ijkr[2]);
01237   } else {
01238     pC3->ijkl[0] = -1;
01239     pC3->ijkr[0] = -1;
01240     pC3->ijkl[1] = -1;
01241     pC3->ijkr[1] = -1;
01242     pC3->ijkl[2] = -1;
01243     pC3->ijkr[2] = -1;
01244   }
01245 
01246   return isOverlap;
01247 }
01248 #endif /* STATIC_MESH_REFINEMENT */