00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
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
00029
00030
00031
00032 #ifdef STATIC_MESH_REFINEMENT
00033
00034
00035 int checkOverlap(SideS *pC1, SideS *pC2, SideS *pC3);
00036
00037
00038 int checkOverlapTouch(SideS *pC1, SideS *pC2, SideS *pC3);
00039 #endif
00040
00041
00042
00043
00044
00045
00046
00047
00048
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
00065 nDim=1;
00066 for (i=1; i<3; i++) if (pM->Nx[i]>1) nDim++;
00067
00068
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]);
00074 pG = pM->Domain[nl][nd].Grid;
00075
00076 pG->time = pM->time;
00077
00078
00079
00080 get_myGridIndex(pD, myID_Comm_world, &myL, &myM, &myN);
00081
00082
00083
00084
00085
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
00107
00108
00109
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
00131
00132
00133
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
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
00172
00173 pG->U = (ConsS***)calloc_3d_array(n3z, n2z, n1z, sizeof(ConsS));
00174 if (pG->U == NULL) goto on_error1;
00175
00176
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
00188
00189
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
00201
00202
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
00220
00221
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
00233
00234
00235
00236
00237
00238
00239
00240
00241 if(myL > 0) pG->lx1_id = pD->GData[myN][myM][myL-1].ID_Comm_Domain;
00242 else pG->lx1_id = -1;
00243
00244
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
00250 if(myM > 0) pG->lx2_id = pD->GData[myN][myM-1][myL].ID_Comm_Domain;
00251 else pG->lx2_id = -1;
00252
00253
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
00259 if(myN > 0) pG->lx3_id = pD->GData[myN-1][myM][myL].ID_Comm_Domain;
00260 else pG->lx3_id = -1;
00261
00262
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
00269
00270
00271 pG->NCGrid = 0;
00272 pG->NPGrid = 0;
00273 pG->NmyCGrid = 0;
00274 pG->NmyPGrid = 0;
00275 pG->CGrid = NULL;
00276 pG->PGrid = NULL;
00277 #endif
00278
00279 }
00280 }}
00281
00282 #ifdef STATIC_MESH_REFINEMENT
00283
00284
00285
00286
00287
00288
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]);
00294 pG = pM->Domain[nl][nd].Grid;
00295
00296
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
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
00309
00310
00311 for (ncd=0; ncd<pM->DomainsPerLevel[nl+1]; ncd++){
00312 pCD = (DomainS*)&(pM->Domain[nl+1][ncd]);
00313
00314
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
00326
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
00333
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
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 }}}
00353 }
00354 }
00355
00356
00357
00358
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
00372 }
00373 }
00374 }
00375
00376
00377
00378
00379
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]);
00386
00387
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
00399
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
00406
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
00419
00420
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
00458
00459
00460
00461
00462 for (dim=0; dim<nDim; dim++){
00463 if (dim == 0) iGrid=l;
00464 if (dim == 1) iGrid=m;
00465 if (dim == 2) iGrid=n;
00466
00467
00468
00469
00470
00471
00472
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
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
00566
00567 }
00568
00569
00570
00571
00572
00573
00574
00575
00576 irefine = 1;
00577 for (i=1;i<=(nl+1);i++) irefine *= 2;
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
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
00670
00671 }
00672
00673 }
00674 }
00675 }}}
00676 }
00677 }
00678 }
00679 }}
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;
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
00705
00706
00707
00708
00709
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]);
00715 pG = pM->Domain[nl][nd].Grid;
00716 get_myGridIndex(pD, myID_Comm_world, &myL, &myM, &myN);
00717
00718
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
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
00731
00732
00733
00734 for (npd=0; npd<(pM->DomainsPerLevel[nl-1]); npd++){
00735 pPD = (DomainS*)&(pM->Domain[nl-1][npd]);
00736
00737
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
00748
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
00755
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
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 }}}
00773 }
00774 }
00775
00776
00777
00778
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
00792 }
00793 }
00794 }
00795
00796
00797
00798
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]);
00805
00806
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
00817
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
00824
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
00836
00837
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
00862
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
00876
00877
00878
00879
00880
00881 for (dim=0; dim<nDim; dim++){
00882 if (dim == 0) iGrid=myL;
00883 if (dim == 1) iGrid=myM;
00884 if (dim == 2) iGrid=myN;
00885
00886
00887
00888
00889
00890
00891
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
00924
00925
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
00989
00990 }
00991
00992
00993
00994
00995
00996
00997
00998
00999 irefine = 1;
01000 for (i=1;i<=nl;i++) irefine *= 2;
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
01032
01033
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
01097
01098 }
01099
01100 }
01101 }
01102 }}}
01103 }
01104 }
01105 }
01106 }}
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;
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
01131
01132 return;
01133
01134
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
01177
01178
01179
01180
01181
01182
01183
01184
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
01215
01216
01217
01218
01219
01220
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