Commits

Nathan Goldbaum committed 480a566

Some updates. I now mandate that the feedback zone has ghost zones
(although they probably shouldn't be used).

  • Participants
  • Parent commits a9c6772

Comments (0)

Files changed (5)

File src/enzo/build/machines/Make.mach.nasa-pleiades

 #-----------------------------------------------------------------------
 
 LOCAL_MPI_INSTALL    = /nasa/sgi/mpt/2.04.10789/
-LOCAL_HDF5_INSTALL   = /u/ngoldbau/yt-x86_64/
-LOCAL_PYTHON_INSTALL = /u/ngoldbau/yt-x86_64/
+LOCAL_HDF5_INSTALL   = /u/ngoldbau/Tools/yt-x86_64
+LOCAL_PYTHON_INSTALL = /u/ngoldbau/Tools/yt-x86_64
 LOCAL_COMPILER_DIR    = /nasa/intel/Compiler/2011.7.256/
 
 #-----------------------------------------------------------------------

File src/enzo/grid/particles/Grid_AccreteOntoAccretingParticle.C

 
   for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
     for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
-      for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) {
-	index = GRIDINDEX_NOGHOST(i,j,k);
+      index = GRIDINDEX_NOGHOST(GridStartIndex[0],j,k);
+      for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) {
 	radius2 = 
 	  POW((CellLeftEdge[0][i] + 0.5*CellWidth[0][i]) - ParticlePosition[0],2) +
 	  POW((CellLeftEdge[1][j] + 0.5*CellWidth[1][j]) - ParticlePosition[1],2) +
   vsink[0] = ThisParticle->ReturnVelocity()[0];
   vsink[1] = ThisParticle->ReturnVelocity()[1];
   vsink[2] = ThisParticle->ReturnVelocity()[2];
-
+  
   for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
     for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
       index = GRIDINDEX_NOGHOST(GridStartIndex[0],j,k);
 	  POW((CellLeftEdge[1][j] + 0.5*CellWidth[1][j]) - ParticlePosition[1],2) +
 	  POW((CellLeftEdge[2][k] + 0.5*CellWidth[2][k]) - ParticlePosition[2],2);   
 	if ((AccretionRadius*AccretionRadius) > radius2) {
-#ifdef DEBUG
-	  fprintf(stderr,
-		  "CellLeftEdge[0][i] = %"GSYM", CellRightEdge[0][i] =%"GSYM"\n"
-		  "CellLeftEdge[1][j] = %"GSYM", CellRightEdge[1][j] =%"GSYM"\n"
-		  "CellLeftEdge[2][k] = %"GSYM", CellRightEdge[2][k] =%"GSYM"\n",
-		  CellLeftEdge[0][i],CellLeftEdge[0][i]+CellWidth[0][i],CellLeftEdge[1][j],CellLeftEdge[1][j]+CellWidth[1][j],
-		  CellLeftEdge[2][k],CellLeftEdge[2][k]+CellWidth[2][k]);
-#endif
-	  if ((CellLeftEdge[0][i] < ParticlePosition[0]) && (CellLeftEdge[0][i]+CellWidth[0][i] > ParticlePosition[0]) &&
-	      (CellLeftEdge[1][j] < ParticlePosition[1]) && (CellLeftEdge[1][j]+CellWidth[1][j] > ParticlePosition[1]) &&
-	      (CellLeftEdge[2][k] < ParticlePosition[2]) && (CellLeftEdge[2][k]+CellWidth[2][k] > ParticlePosition[2]))
-	    fprintf(stderr,"Particle in this cell!\n");
-
-
 	  // useful shorthand
 	  vgas[0] = BaryonField[Vel1Num][index];
 	  vgas[1] = BaryonField[Vel2Num][index];

File src/enzo/grid/particles/Grid_ConstructFeedbackZone.C

 #include "ActiveParticle.h"
 #include "phys_constants.h"
 
-int grid::ConstructFeedbackZone(ActiveParticleType* ThisParticle, int FeedbackRadius, FLOAT dx, grid* FeedbackZone)
+grid* grid::ConstructFeedbackZone(ActiveParticleType* ThisParticle, int FeedbackRadius, FLOAT dx)
 {
   FLOAT* ParticlePosition = ThisParticle->ReturnPosition();
 
   if ((GridLeftEdge[0] > ParticlePosition[0]+FeedbackRadius) || (GridRightEdge[0] < ParticlePosition[0]-FeedbackRadius) ||
       (GridLeftEdge[1] > ParticlePosition[1]+FeedbackRadius) || (GridRightEdge[1] < ParticlePosition[1]-FeedbackRadius) ||
       (GridLeftEdge[2] > ParticlePosition[2]+FeedbackRadius) || (GridRightEdge[2] < ParticlePosition[2]-FeedbackRadius))
-    return FAIL;
+    ENZO_FAIL("Particle outside own grid!");
 
   /* Setup grid properties */
 
-  int FeedbackZoneRank = FeedbackZone->GetGridRank();
+  int FeedbackZoneRank = this->GetGridRank();
 
-  // Since the grid creation machinery assume we want ghost zones, we need to
-  // trick it into giving us a fake grid without ghost zones.  We therefore ask
-  // for a cubic grid of dimension (2*(FeedbackRadius-3)+1)^3
-  
-  int FeedbackZoneDimension[FeedbackZoneRank], size;
-  FLOAT LeftCellOffset[FeedbackZoneRank],RightCellOffset[FeedbackZoneRank],
-    FeedbackZoneLeftEdge[FeedbackZoneRank], FeedbackZoneRightEdge[FeedbackZoneRank];
+  int FeedbackZoneDimension[FeedbackZoneRank], size=1;
+  FLOAT LeftCellOffset[FeedbackZoneRank], FeedbackZoneLeftEdge[FeedbackZoneRank], FeedbackZoneRightEdge[FeedbackZoneRank];
+  FLOAT CellSize = CellWidth[0][0], ncells[FeedbackZoneRank];
 
   for (int i = 0; i < FeedbackZoneRank; i++) {
-    if (FeedbackRadius > DEFAULT_GHOST_ZONES)
-      FeedbackZoneDimension[i] = 2*FeedbackRadius+1;
-    else
-      FeedbackZoneDimension[i] = 2*DEFAULT_GHOST_ZONES+1;
+    FeedbackZoneDimension[i] = (2*(FeedbackRadius+DEFAULT_GHOST_ZONES)+1);
     size *= FeedbackZoneDimension[i];
 
-    LeftCellOffset[i]        = fmod(ParticlePosition[i],dx);
-    RightCellOffset[i]       = dx-LeftCellOffset[i];
-    
-    FeedbackZoneLeftEdge[i]  = ParticlePosition[i]-FeedbackRadius*dx-LeftCellOffset[i];
-    FeedbackZoneRightEdge[i] = ParticlePosition[i]+FeedbackRadius*dx+RightCellOffset[i];
+    LeftCellOffset[i]        = modf((ParticlePosition[i]-CellLeftEdge[i][0])/CellSize,&ncells[i]);
+        
+    FeedbackZoneLeftEdge[i]  = CellLeftEdge[0][0]+CellSize*(ncells[i]-FeedbackRadius);
+    FeedbackZoneRightEdge[i] = CellLeftEdge[0][0]+CellSize*(ncells[i]+FeedbackRadius+1);
   }
 
-  /* Intialize the fake grid */
-  FeedbackZone = new grid;
+  grid *FeedbackZone = new grid;
   
   FeedbackZone->InheritProperties(this);
   
   
   FeedbackZone->SetProcessorNumber(MyProcessorNumber);
 
-  FeedbackZone->AllocateAndZeroBaryonField();
-    
+  FeedbackZone->SetTimeStep(this->ReturnTimeStep());
+
+  if (FeedbackZone->AllocateAndZeroBaryonField() == FAIL)
+    ENZO_FAIL("FeedbackZone BaryonField allocation failed");
+
   // Allocate flagging field of the same size as BaryonField. If FlaggingField =
   // 0, the corresponding zone in the BaryonField has not been copied yet.
     
   // Note, using ZeroVector here will break if a FeedbackZone overlaps with a
   // domain boundary
   float ZeroVector[] = {0,0,0};
-  FeedbackZone->CopyZonesFromGrid(this,ZeroVector);
+  if (FeedbackZone->CopyZonesFromGrid(this,ZeroVector) == FAIL)
+    ENZO_FAIL("FeedbackZone copy failed!");
 
   // if the grid is filled, return
   
 
   delete [] FlaggingField;
 
-  return SUCCESS;
+  return FeedbackZone;
   
 }

File src/enzo/headers/Grid.h

     if (MyProcessorNumber != ProcessorNumber)
       return SUCCESS;
 
-    if (BaryonField != NULL)
+    if (BaryonField[0] != NULL)
       return FAIL;
 
     int size = this->GetGridSize();
 
     for (int field = 0; field < NumberOfBaryonFields; field++) {
-      BaryonField[field] = new float[size]();
+      BaryonField[field] = new float[size];
     }
 
     return SUCCESS;
   int DetachActiveParticles(void);
   int MirrorActiveParticles(void);
   int DebugActiveParticles(int level);
-  int ConstructFeedbackZone(ActiveParticleType* ThisParticle, int FeedbackRadius, FLOAT dx, grid* FeedbackZone);
+  grid* ConstructFeedbackZone(ActiveParticleType* ThisParticle, int FeedbackRadius, FLOAT dx);
   int DistributeFeedbackZone(ActiveParticleType* ThisParticle, FLOAT FeedbackRadius);
 
   /* Returns averaged velocity from the 6 neighbor cells and itself */

File src/enzo/particles/active_particles/ActiveParticle_AccretingParticle.C

       return FAIL;
     }
     
-    if (sinkGrid->ConstructFeedbackZone(ParticleList[i],AccretionRadius, dx, FeedbackZone) == FAIL) 
-      return FAIL;
-    
+    FeedbackZone = sinkGrid->ConstructFeedbackZone(ParticleList[i],AccretionRadius, dx);
+						       
     float AccretionRate = 0;
     ActiveParticleType_AccretingParticle* temp = static_cast<ActiveParticleType_AccretingParticle*>(ParticleList[i]);