Commits

Sam Skillman committed 2fa5cc2 Merge

Merged in ngoldbaum/enzo-3.0 (pull request #33)

  • Participants
  • Parent commits b8c17ce, bb0ed24

Comments (0)

Files changed (21)

run/Hydro/Hydro-3D/CollapseTestAccretingParticleTwoSphere/CollapseTestAccretingParticleTwoSphere.enzo

+#
+#  AMR PROBLEM DEFINITION FILE: Non-cosmological Collapse test
+#  Description: a sphere collapses until becoming pressure supported.
+#
+#  define problem
+#
+ProblemType                = 27         // Collapse test
+TopGridRank                = 3
+TopGridDimensions          = 64 64 64
+SelfGravity                = 1          // gravity on
+TopGridGravityBoundary     = 0          // periodic
+LeftFaceBoundaryCondition  = 3 3 3      // periodic
+RightFaceBoundaryCondition = 3 3 3
+#ExtraOutputs               = 1 2 3 4 5 6 7 8 9
+#
+# problem parameters
+#
+
+CollapseTestRefineAtStart   = 1         // check refinement before running
+CollapseTestNumberOfSpheres = 2
+CollapseTestUseParticles    = 0
+CollapseTestInitialTemperature = 500    // temperature of the background gas
+CollapseTestUniformVelocity = 0.0 0.1 0.0
+CollapseTestSpherePosition[0]   = 0.7 0.2 0.5
+CollapseTestSphereVelocity[0]   = 0.0 0.0 0.0
+CollapseTestSphereRadius[0]     = 0.1
+CollapseTestSphereCoreRadius[0] = 0.05  // only used with sphere type 5
+CollapseTestSphereDensity[0]    = 50   // sphere density, the background density is 1
+CollapseTestSphereTemperature[0] = 10    // put sphere in pressure equilibrium (rho * T is constant)
+CollapseTestSphereRotationPeriod[0] = 0
+CollapseTestSphereType[0]       = 5     // constant density
+                                        // 1: uniform
+					// 2: r^-2 power-law
+					// 3: NFW
+					// 4: Gaussian
+					// 5: r^-2 power-law with a core
+					// 11: Burkert & Bodenheimer setup
+CollapseTestSpherePosition[1]   = 0.3 0.2 0.5
+CollapseTestSphereVelocity[1]   = 0.0 0.0 0.0
+CollapseTestSphereRadius[1]     = 0.1
+CollapseTestSphereCoreRadius[1] = 0.05  // only used with sphere type 5
+CollapseTestSphereDensity[1]    = 50   // sphere density, the background density is 1
+CollapseTestSphereTemperature[1] = 10    // put sphere in pressure equilibrium (rho * T is constant)
+CollapseTestSphereRotationPeriod[1] = 0
+CollapseTestSphereType[1]       = 5     // constant density
+                                        // 1: uniform
+					// 2: r^-2 power-law
+					// 3: NFW
+					// 4: Gaussian
+					// 5: r^-2 power-law with a core
+					// 11: Burkert & Bodenheimer setup
+
+
+
+#  no cosmology for this run
+#
+ComovingCoordinates   = 0              // Expansion OFF
+#
+#  units
+#
+DensityUnits          = 3.82e-20      // 10^4 cm^-3
+LengthUnits           = 5e17          // 1 pc in cm
+TimeUnits             = 1.0e12        // 10^4 yrs
+GravitationalConstant = 0.0320327594   // 4*pi*G_{cgs}*DensityUnits*TimeUnits^2
+#
+#  set I/O and stop/start parameters
+#
+StopTime          = 1.0
+dtDataDump        = 0.01
+#CycleSkipDataDump = 1
+DataDumpDir       = DD
+DataDumpName      = DD
+OutputTemperature = 1                  // Output temperature field.
+#
+#  set hydro parameters
+#
+Gamma                       = 1.0000001
+PPMDiffusionParameter       = 0        // diffusion off
+DualEnergyFormalism         = 1        // use total & internal energy
+InterpolationMethod         = 1        // SecondOrderA
+CourantSafetyNumber         = 0.3
+FluxCorrection              = 1
+ConservativeInterpolation   = 1
+RiemannSolver = 4
+HydroMethod                 = 0        // PPM
+SmallRho                    = 1e-8
+SmallE			    = 1e-8
+#
+#  chemistry/cooling
+#
+MultiSpecies                = 0        // chemistry off
+RadiativeCooling            = 0        // cooling off
+Mu  			    = 3.0
+#
+#  set grid refinement parameters
+#
+StaticHierarchy           = 0          // dynamic hierarchy
+MaximumRefinementLevel    = 1         // use up to 10 levels
+RefineBy                  = 2          // refinement factor
+CellFlaggingMethod        = 6         // use Truelove criterion for refinement 
+MinimumEfficiency         = 0.3
+#OutputFirstTimeAtLevel    = 4         // output when level 4, 5, 6, etc reached (commented out for now)
+#StopFirstTimeAtLevel      = 10         // stop if/when level 10 reached
+RefineByJeansLengthSafetyFactor = 8    // resolve Jeans length by 4 cells (used with CellFlaggingMethod 6)
+
+
+#
+#  set some global parameters
+#
+GreensFunctionMaxNumber   = 10         // # of greens function at any one time
+PotentialIterations	  = 100
+#
+#  active particle parameters
+#
+AppendActiveParticleType = AccretingParticle
+

src/enzo/communication/CommunicationCollectParticles.C

 /         bit inefficient but reduced the number of if-statements.
 /
 ************************************************************************/
-#define DEBUG_AP
- 
 #ifdef USE_MPI
 #include "communicators.h"
 #endif /* USE_MPI */
 #include "ActiveParticle.h"
 #include "SortCompareFunctions.h"
 
+#define NO_DEBUG_AP
+
 void my_exit(int status);
  
 // function prototypes

src/enzo/control/EvolveLevel.C

File contents unchanged.

src/enzo/grid/Make.module.objects

 	Grid_SetSubgridMarkerFromSubgrid.o \
 	Grid_SubgridMarkerPostParallel.o \
 	Grid_SubtractAccretedMassFromSphere.o \
-	Grid_SumGasMass.o \
-	Grid_SumGasMassGZ.o \
 	Grid_WriteCube.o \
 	Grid_WriteCubeInterpolate.o \
 	Grid_WriteGrid.o \

src/enzo/grid/particles/Grid_AccreteOntoAccretingParticle.C

 #include "phys_constants.h"
 
 
-#define NO_DEBUG
+#define NO_DEBUG_AP
 
 float bondi_alpha(float x);
 
   /* Return if this doesn't involve us */
   if (MyProcessorNumber != ProcessorNumber) 
     return SUCCESS;
-  
+
   /* Check whether the cube that circumscribes the accretion zone intersects with this grid */
 
   FLOAT *ParticlePosition = (*ThisParticle)->ReturnPosition();
   float velz = BaryonField[Vel3Num][cgindex];
   FLOAT BondiHoyleRadius;
   float *Temperature = new float[size]();
-  float msink = (*ThisParticle)->ReturnMass()*POW(CellSize,3);
+  float msink = (*ThisParticle)->ReturnMass()*CellVolume;
 
   this->ComputeTemperatureField(Temperature);
 
 	    for (jsub = 0; jsub < NDIV-1; jsub++) {
 	      ydist = CellLeftEdge[1][j] + CellWidth[1][j]*(float(jsub)+0.5)/NDIV - ParticlePosition[1];
 	      for (isub = 0; isub < NDIV-1; isub++) {
-		xdist = CellLeftEdge[0][i] + CellWidth[0][i]*(float(jsub)+0.5)/NDIV - ParticlePosition[1];
-
+		xdist = CellLeftEdge[0][i] + CellWidth[0][i]*(float(jsub)+0.5)/NDIV - ParticlePosition[0];
+		
 		dist = sqrt(xdist*xdist+ydist*ydist+zdist*zdist);
 		if (dist == 0.0)
 		  dist = CellWidth[0][0]/huge;
-
+		
 		// Compute specific angular momentum
 		jsp[0] = ydist*(vgas[2] - vsink[2]) -
 		  zdist*(vgas[1] - vsink[1]);
 		  ydist*(vgas[0] - vsink[0]);
 		
 		jspsqr = jsp[0]*jsp[0]+jsp[1]*jsp[1]+jsp[2]*jsp[2];
-
+		
 		// Compute specific kinetic + gravitational energy
 		esp = (POW((vgas[0] - vsink[0]),2) + 
 		       POW((vgas[1] - vsink[1]),2) +
 		else
 		  rmin = -GravitationalConstant*msink/(2.0*esp) *
 		    (1.0 - sqrt(1.0 + 2.0*jspsqr*esp/POW(GravitationalConstant*msink,2)));
-
+		
 		dxmin = rmin / CellWidth[0][0];
 		if (dxmin >= 0.25)
 		  nexcluded[index]+=1;
-
+		
 	      } // ksub
 	    } // jsub
 	  } // ksub
 	  
 	  if (abs(i-cindex) <= 1 && abs(j-cindex) <= 1 && abs(k-cindex) <= 1)
 	    maxexcluded = max(nexcluded[index],maxexcluded);
-
+	  
 	}
       }
     }
   }
-
+  
   // Correct the central cell
   if (nexcluded[cgindex] > 0)
     if (KernelRadius / CellWidth[0][0] >= 0.25)
       nexcluded[cgindex] = maxexcluded;
     else
       nexcluded[cgindex] = 0;
-  
+    
   AverageDensity = WeightedSum/SumOfWeights;
   
   // Eqn 12
 	  if (maccreted > 0.25*mcell) 
 	    maccreted = 0.25*mcell;
 	  
-	  // Scale down maccrete
+	  // Scale down maccreted
 	  maccreted = maccreted/POW(NDIV,3) *
 	    (POW(NDIV,3)-nexcluded[index]);
 
 	    AccretedMomentum[1] += paccrete[1];
 	    AccretedMomentum[2] += paccrete[2];
 
-#ifdef DEBUG
+#ifdef DEBUG_AP
 	    if (index == cgindex)
-	      printf("Sink Density: %"FSYM", Cell Density: %"FSYM", New Density: %"FSYM"\n",
+	      printf("Sink Density: %"GOUTSYM", Cell Density: %"GOUTSYM", New Density: %"GOUTSYM"\n",
 		     maccreted/CellVolume,
 		     BaryonField[DensNum][index],
 		     BaryonField[DensNum][index]-maccreted/CellVolume);
 	    } else // Zeus
 	      BaryonField[TENum][index] = eintnew/mnew;
 
-	    BaryonField[Vel1Num][index] = pnew[0]/mcell;
-	    BaryonField[Vel2Num][index] = pnew[1]/mcell;
-	    BaryonField[Vel3Num][index] = pnew[2]/mcell;
+	    BaryonField[Vel1Num][index] = pnew[0]/mnew;
+	    BaryonField[Vel2Num][index] = pnew[1]/mnew;
+	    BaryonField[Vel3Num][index] = pnew[2]/mnew;
 
 	    // Check if mass or energy is too small, correct if necessary
 	    if (BaryonField[DensNum][index] < SmallRhoFac*SmallRho) {
   return SUCCESS;
 }
 
-#undef DEBUG
+#undef DEBUG_AP

src/enzo/grid/utilities/Grid_CopyActiveZonesFromGrid.C

      The loop breaks when check_overlap returns -1.
    */
   int overlap = 0;
-  while (true) {
-
+  while (true) 
+    {
+   
       for (dim = 0; dim < GridRank; dim++) {
         ActiveLeft[dim]  = GridLeftEdge[dim]  + EdgeOffset[dim];
         ActiveRight[dim] = GridRightEdge[dim] + EdgeOffset[dim];
 
   delete [] shift;
 
-  this->DebugCheck("CopyZonesFromGrid (after)");
+  this->DebugCheck("CopyActiveZonesFromGrid (after)");
 
   return SUCCESS;
 }

src/enzo/grid/utilities/Grid_DepositRefinementZone.C

     }
   } // dim
 
-  //if ((LeftCorner[0] > ParticlePosition[0]+RefinementRadius) || (RightCorner[0] < ParticlePosition[0]-RefinementRadius) ||
-  //    (LeftCorner[1] > ParticlePosition[1]+RefinementRadius) || (RightCorner[1] < ParticlePosition[1]-RefinementRadius) ||
-  //   (LeftCorner[2] > ParticlePosition[2]+RefinementRadius) || (RightCorner[2] < ParticlePosition[2]-RefinementRadius))
-  //  return SUCCESS;
-
   for (dim = 0; dim < GridRank; dim++) {
     if (overlaps[dim] == false) return SUCCESS;
   }
     for (j = 0; j < GridDimension[1]; j++) {
       for (k = 0; k < GridDimension[2]; k++) {
 	index = (k*GridDimension[1]+j)*GridDimension[0]+i;
-    dist2 = calc_dist2(CellLeftEdge[0][i] + CellSize/2.,
-      CellLeftEdge[1][j] + CellSize/2.,
-      CellLeftEdge[2][k] + CellSize/2.,
-      ParticlePosition[0], ParticlePosition[1], ParticlePosition[2],
-      period);
-    if (dist2 <= rad2) {
-      FlaggingField[index] = 1;
-      NumberOfFlaggedCells++;
-    }
-    
-
-// 	for (dim = 0; dim < GridRank; dim++) overlaps[dim] = false;
-// 	// Need to do this sort of expensive check since a derefined
-// 	// cell could enclose the accretion zone yet still be centered
-// 	// outside the accretion radius
-// 	for (dim = 0; dim < GridRank; dim++) {
-// 	  left = ParticlePosition[dim]-RefinementRadius;
-// 	  right = ParticlePosition[dim]+RefinementRadius;
-// 	  if (dim == 0) {a = i;}
-// 	  else if (dim == 1) {a = j;}
-// 	  else {a = k;}
-// 	  if (!((CellLeftEdge[dim][a] > right) ||
-// 	        (CellLeftEdge[dim][a] + CellSize < left))) {
-// 	    overlaps[dim] = true;
-// 	  }
-// 	  // Periodicity.
-// 	  if (overlaps[dim] == false) {
-//         pleft = max(period[dim] - fabs(left), left);
-//         pright = min(fabs(right - period[dim]), right);
-//         if ((pleft != left) && (CellLeftEdge[dim][a] + CellSize > pleft)) {
-//           overlaps[dim] = true;
-//         }
-//         if ((pright != right) && (CellLeftEdge[dim][a] < pright)) {
-//           overlaps[dim] = true;
-//         }
-// 	  } // overlaps
-// 	} // dim
-//     
-//     if ((overlaps[0] == true) && (overlaps[1] == true) && (overlaps[2])) {
-//       //printf("%f %f %f\n", CellLeftEdge[0][i], CellLeftEdge[1][j], CellLeftEdge[2][k]);
-//       FlaggingField[index] = 1;
-//       NumberOfFlaggedCells++;
-//     }
+	dist2 = calc_dist2(CellLeftEdge[0][i] + CellSize/2.,
+		  CellLeftEdge[1][j] + CellSize/2.,
+		  CellLeftEdge[2][k] + CellSize/2.,
+		  ParticlePosition[0], ParticlePosition[1], ParticlePosition[2],
+		  period);
+	if (dist2 <= rad2) {
+	  FlaggingField[index] = 1;
+	  NumberOfFlaggedCells++;
+	}
       } // k
     } // j
   } // i

src/enzo/grid/utilities/Grid_SumGasMass.C

-/***********************************************************************
-/
-/  GET TOTAL GAS MASS ON THIS GRID
-/
-/  written by: Nathan Goldbaum
-/  date:       June, 2012
-/  modified1:
-/
-/
-************************************************************************/
-#include "preincludes.h"
-
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "phys_constants.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "TopGridData.h"
-#include "LevelHierarchy.h"
-#include "CosmologyParameters.h"
-#include "ActiveParticle.h"
-
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, FLOAT Time);
-
-int grid::SumGasMass(float *mass)
-{
-
-  if (MyProcessorNumber != ProcessorNumber)
-    return SUCCESS;
-
-  /* Get indices in BaryonField for density, internal energy, thermal energy, velocity */
-  int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num;
-  if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
-				       Vel3Num, TENum) == FAIL) {
-    ENZO_FAIL("Error in IdentifyPhysicalQuantities.");
-  }
-
-  float CellVolume = POW(CellWidth[0][0],3);
-  float MassOnGrid = 0.0;
-  int i,j,k,index;
-  
-  for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
-    for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
-      index = GRIDINDEX_NOGHOST(GridStartIndex[0],j,k);
-      for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) {
-	MassOnGrid+=BaryonField[DensNum][index];
-      }
-    }
-  }
-
-  for (i = 0; i < NumberOfActiveParticles; i++)
-    MassOnGrid+=this->ActiveParticles[i]->ReturnMass();
-  
-  *mass += MassOnGrid;//*CellVolume;
-
-  return SUCCESS;
-}

src/enzo/grid/utilities/Grid_SumGasMassGZ.C

-/***********************************************************************
-/
-/  GET TOTAL GAS MASS ON THIS GRID
-/
-/  written by: Nathan Goldbaum
-/  date:       June, 2012
-/  modified1:
-/
-/
-************************************************************************/
-#include "preincludes.h"
-
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "phys_constants.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "TopGridData.h"
-#include "LevelHierarchy.h"
-#include "CosmologyParameters.h"
-#include "ActiveParticle.h"
-
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, FLOAT Time);
-
-float grid::SumGasMassGZ()
-{
-
-  if (MyProcessorNumber != ProcessorNumber)
-    return 0;
-
-  /* Get indices in BaryonField for density, internal energy, thermal energy, velocity */
-  int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num;
-  if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
-				       Vel3Num, TENum) == FAIL) {
-    ENZO_FAIL("Error in IdentifyPhysicalQuantities.");
-  }
-
-  float CellVolume = POW(CellWidth[0][0],3);
-  float MassOnGrid = 0.0;
-  int size=1,index=0, i;
-  
-  for (i = 0; i < 3; i++)
-    size *= GridDimension[i];
-
-  for (index = 0; index < size; index++)
-    MassOnGrid+=BaryonField[DensNum][index];
-
-  for (i = 0; i < NumberOfActiveParticles; i++)
-    MassOnGrid+=this->ActiveParticles[i]->ReturnMass();
-  
-  return MassOnGrid;//*CellVolume;
-
-}

src/enzo/headers/ActiveParticle.h

 
   FLOAT *ReturnPosition(void) { return pos; };
   float *ReturnVelocity(void) { return vel; };
+  float ReturnMomentum(int dim) { return Mass*vel[dim]; };
   void   ConvertAllMassesToSolar(void);
   void   ConvertMassToSolar(void);
   void   Merge(ActiveParticleType *a);

src/enzo/headers/ActiveParticle_AccretingParticle.h

       /* Do accretion */
       if (Accrete(nParticles,ParticleList,AccretionRadius,dx,LevelArray,ThisLevel) == FAIL)
 	ENZO_FAIL("Accreting Particle accretion failed. \n");
-     
+      
       for (i = 0; i < nParticles; i++)
 	if (ParticleList[i]->ReturnCurrentGrid()->ReturnProcessorNumber() != MyProcessorNumber) {
 	  delete ParticleList[i];

src/enzo/headers/Grid.h

    int ReturnNumberOfParticles() {return NumberOfParticles;};
    int ReturnNumberOfActiveParticles() {return NumberOfActiveParticles;};
    int ReturnNumberOfActiveParticlesOfThisType(int ActiveParticleIDToFind);
+   ActiveParticleType** ReturnActiveParticles() {return ActiveParticles;};
 
    int ReturnNumberOfStarParticles(void);
 
 			     float &metallicity3,
 			     float &coldgas_mass, float AvgVelocity[]);
 
-  /* Sum mass on this grid */
-
-  int SumGasMass(float *mass);
-
-  float SumGasMassGZ();
-
   int RemoveParticle(int ID, bool disable=false);
 
   int RemoveActiveParticle(PINT ID, int NewProcessorNumber);

src/enzo/hierarchy/ConstructFeedbackZone.C

+/***********************************************************************
+/
+/  Construct a fake grid for feedback algorithms based on a 
+/  list of active particles
+/
+/  written by: Nathan Goldbaum
+/  date:       June 2012
+/
+************************************************************************/
+
+#ifdef USE_MPI
+#include "communicators.h"
+#endif
+
+#include "preincludes.h"
+
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "phys_constants.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "ActiveParticle.h"
+#include "phys_constants.h"
+#include "CommunicationUtilities.h"
+#include "communication.h"
+
+int CommunicationBufferPurge(void);
+int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
+				int NumberOfSubgrids[] = NULL,
+				int FluxFlag = FALSE,
+				TopGridData* MetaData = NULL);
+
+
+
+grid* ConstructFeedbackZone(ActiveParticleType* ThisParticle,int FeedbackRadius, 
+			    FLOAT dx, HierarchyEntry** Grids, int NumberOfGrids,
+			    int SendField)
+{
+  int i,j,dim,size;
+  int FeedbackZoneRank;
+  FLOAT FBRdx;
+  FLOAT ParticlePosition[3] = {ThisParticle->ReturnPosition()[0],
+			       ThisParticle->ReturnPosition()[1],
+			       ThisParticle->ReturnPosition()[2]};
+
+  /* Build array of AP grids and check for errors */
+  grid* APGrid = new grid;
+
+  APGrid = ThisParticle->ReturnCurrentGrid();
+  FBRdx = dx * FLOAT(FeedbackRadius);
+    
+  if (APGrid == NULL)
+    ENZO_FAIL("Particle CurrentGrid is invalid!\n"); 
+    
+  // This should only happen if the grid pointer is invalid
+  if ((APGrid->GetGridLeftEdge(0) > ParticlePosition[0]+FBRdx) || 
+      (APGrid->GetGridLeftEdge(1) > ParticlePosition[1]+FBRdx) || 
+      (APGrid->GetGridLeftEdge(2) > ParticlePosition[2]+FBRdx) || 
+      (APGrid->GetGridRightEdge(0) < ParticlePosition[0]-FBRdx) ||
+      (APGrid->GetGridRightEdge(1) < ParticlePosition[1]-FBRdx) ||
+      (APGrid->GetGridRightEdge(2) < ParticlePosition[2]-FBRdx))
+    ENZO_FAIL("Particle outside own grid!\n");
+  
+  
+  /* Setup Feedback Zones before copying data */
+  
+  FeedbackZoneRank = APGrid->GetGridRank();
+
+  int FeedbackZoneDimension[MAX_DIMENSION];
+  FLOAT LeftCellOffset[MAX_DIMENSION],FeedbackZoneLeftEdge[MAX_DIMENSION], 
+    FeedbackZoneRightEdge[MAX_DIMENSION], ncells[MAX_DIMENSION];
+  FLOAT CellSize, GridGZLeftEdge;
+
+  for (dim = 0; dim < FeedbackZoneRank; dim++) {
+    FeedbackZoneDimension[dim] = (2*(FeedbackRadius+DEFAULT_GHOST_ZONES)+1);
+    CellSize = APGrid->GetCellWidth(dim,0);
+    GridGZLeftEdge = APGrid->GetCellLeftEdge(dim,0);
+      
+    LeftCellOffset[dim] = MODF((ParticlePosition[dim]-GridGZLeftEdge)/CellSize,&ncells[dim]);
+
+    FeedbackZoneLeftEdge[dim]  = GridGZLeftEdge + CellSize*(ncells[dim]-FeedbackRadius);
+    FeedbackZoneRightEdge[dim] = GridGZLeftEdge + CellSize*(ncells[dim]+FeedbackRadius+1);
+  }
+
+  grid* FeedbackZone = new grid;
+      
+  FeedbackZone->InheritProperties(APGrid);
+    
+  FeedbackZone->PrepareGrid(FeedbackZoneRank, FeedbackZoneDimension, 
+			    FeedbackZoneLeftEdge,FeedbackZoneRightEdge,0);
+    
+  FeedbackZone->SetProcessorNumber(APGrid->ReturnProcessorNumber());
+  
+  FeedbackZone->SetTimeStep(APGrid->ReturnTimeStep());
+        
+  // This will only allocate the BaryonField on the host processor
+  if (FeedbackZone->AllocateAndZeroBaryonField() == FAIL)
+    ENZO_FAIL("FeedbackZone BaryonField allocation failed\n");
+
+  if (SendField == GRAVITATING_MASS_FIELD) {
+    size = 1;
+    // If we're doing gravity, we need to init that field.
+    FeedbackZone->InitializeGravitatingMassField(RefineBy);
+    for (dim = 0; dim < FeedbackZoneRank; dim++) {
+      size *= FeedbackZone->ReturnGravitatingMassFieldDimension(dim);
+    }
+    FeedbackZone->InitGravitatingMassField(size);
+  }
+
+  // Copy zones from this grid (which must overlap the position of the AP).
+  FLOAT ZeroVector[] = {0,0,0};
+
+  /* Post receives */
+
+  CommunicationReceiveIndex = 0;
+  CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE;
+  CommunicationDirection = COMMUNICATION_POST_RECEIVE;
+
+  for (j = 0; j < NumberOfGrids; j++) 
+    if (FeedbackZone->CopyActiveZonesFromGrid(Grids[j]->GridData,ZeroVector,SendField) == FAIL)
+      ENZO_FAIL("FeedbackZone copy failed!\n");
+    
+  /* Send data */
+
+  CommunicationDirection = COMMUNICATION_SEND;
+
+  for (j = 0; j < NumberOfGrids; j++)
+    if (FeedbackZone->CopyActiveZonesFromGrid(Grids[j]->GridData,ZeroVector,SendField) == FAIL)
+      ENZO_FAIL("FeedbackZone copy failed!\n");
+ 
+  /* Receive data */
+
+  if (CommunicationReceiveHandler() == FAIL)
+    ENZO_FAIL("CommunicationReceiveHandler() failed!\n");
+
+#ifdef USE_MPI
+  CommunicationBufferPurge();
+#endif
+
+  return FeedbackZone;
+  
+}

src/enzo/hierarchy/ConstructFeedbackZones.C

-/***********************************************************************
-/
-/  Construct a fake grid for feedback algorithms based on a 
-/  list of active particles
-/
-/  written by: Nathan Goldbaum
-/  date:       June 2012
-/
-************************************************************************/
-
-#ifdef USE_MPI
-#include "communicators.h"
-#endif
-
-#include "preincludes.h"
-
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "phys_constants.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "ActiveParticle.h"
-#include "phys_constants.h"
-#include "CommunicationUtilities.h"
-#include "communication.h"
-
-int CommunicationBufferPurge(void);
-int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
-				int NumberOfSubgrids[] = NULL,
-				int FluxFlag = FALSE,
-				TopGridData* MetaData = NULL);
-
-
-
-grid** ConstructFeedbackZones(ActiveParticleType** ParticleList, int nParticles,
-    int *FeedbackRadius, FLOAT dx, HierarchyEntry** Grids, int NumberOfGrids,
-    int SendField)
-{
-  int i,j,dim,size;
-  int FeedbackZoneRank;
-  FLOAT FBRdx;
-  FLOAT** ParticlePosition = new FLOAT*[nParticles]();
-
-  /* Build array of AP grids and check for errors */
-  grid** APGrids = new grid*[nParticles];
-
-  for (i = 0; i < nParticles; i++) {
-    APGrids[i] = ParticleList[i]->ReturnCurrentGrid();
-    FBRdx = dx * FLOAT(FeedbackRadius[i]);
-    
-    if (APGrids[i] == NULL)
-      ENZO_FAIL("Particle CurrentGrid is invalid!\n");
-
-    ParticlePosition[i] = ParticleList[i]->ReturnPosition();
-    
-    // This should only happen if the grid pointer is invalid
-    if ((APGrids[i]->GetGridLeftEdge(0) > ParticlePosition[i][0]+FBRdx) || 
-	(APGrids[i]->GetGridLeftEdge(1) > ParticlePosition[i][1]+FBRdx) || 
-	(APGrids[i]->GetGridLeftEdge(2) > ParticlePosition[i][2]+FBRdx) || 
-	(APGrids[i]->GetGridRightEdge(0) < ParticlePosition[i][0]-FBRdx) ||
-	(APGrids[i]->GetGridRightEdge(1) < ParticlePosition[i][1]-FBRdx) ||
-	(APGrids[i]->GetGridRightEdge(2) < ParticlePosition[i][2]-FBRdx))
-      ENZO_FAIL("Particle outside own grid!\n");
-  }
-
-  /* Setup Feedback Zones before copying data */
-
-  grid** FeedbackZones = new grid*[nParticles];
-  
-  for (i = 0; i < nParticles; i++) {
-    FeedbackZoneRank = APGrids[i]->GetGridRank();
-
-    int FeedbackZoneDimension[FeedbackZoneRank];
-    FLOAT LeftCellOffset[FeedbackZoneRank],FeedbackZoneLeftEdge[FeedbackZoneRank], 
-      FeedbackZoneRightEdge[FeedbackZoneRank];
-    FLOAT CellSize, GridGZLeftEdge, ncells[FeedbackZoneRank];
-
-    for (dim = 0; dim < FeedbackZoneRank; dim++) {
-      FeedbackZoneDimension[dim] = (2*(FeedbackRadius[i]+DEFAULT_GHOST_ZONES)+1);
-      CellSize = APGrids[i]->GetCellWidth(dim,0);
-      GridGZLeftEdge = APGrids[i]->GetCellLeftEdge(dim,0);
-      
-      LeftCellOffset[dim] = MODF((ParticlePosition[i][dim]-GridGZLeftEdge)/CellSize,&ncells[dim]);
-
-      FeedbackZoneLeftEdge[dim]  = GridGZLeftEdge + CellSize*(ncells[dim]-FeedbackRadius[i]);
-      FeedbackZoneRightEdge[dim] = GridGZLeftEdge + CellSize*(ncells[dim]+FeedbackRadius[i]+1);
-    }
-    
-    grid *FeedbackZone = new grid;
-    
-    FeedbackZone->InheritProperties(APGrids[i]);
-    
-    FeedbackZone->PrepareGrid(FeedbackZoneRank, FeedbackZoneDimension, 
-			      FeedbackZoneLeftEdge,FeedbackZoneRightEdge,0);
-    
-    FeedbackZone->SetProcessorNumber(APGrids[i]->ReturnProcessorNumber());
-    
-    FeedbackZone->SetTimeStep(APGrids[i]->ReturnTimeStep());
-        
-    // This will only allocate the BaryonField on the host processor
-    if (FeedbackZone->AllocateAndZeroBaryonField() == FAIL)
-      ENZO_FAIL("FeedbackZone BaryonField allocation failed\n");
-
-	if (SendField == GRAVITATING_MASS_FIELD) {
-	    size = 1;
-	    // If we're doing gravity, we need to init that field.
-	    FeedbackZone->InitializeGravitatingMassField(RefineBy);
-	    for (dim = 0; dim < FeedbackZoneRank; dim++) {
-	      size *= FeedbackZone->ReturnGravitatingMassFieldDimension(dim);
-	    }
-	    FeedbackZone->InitGravitatingMassField(size);
-    }
-
-    FeedbackZones[i] = FeedbackZone;
-  }
-
-  delete [] APGrids;
-  delete [] ParticlePosition; 
-
-  // Copy zones from this grid (which must overlap the position of the AP).
-  FLOAT ZeroVector[] = {0,0,0};
-
-
-  /* Post receives */
-
-  CommunicationReceiveIndex = 0;
-  CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE;
-  CommunicationDirection = COMMUNICATION_POST_RECEIVE;
-
-  for (i = 0; i < nParticles; i++) 
-    for (j = 0; j < NumberOfGrids; j++) 
-      if (FeedbackZones[i]->CopyActiveZonesFromGrid(Grids[j]->GridData,ZeroVector,SendField) == FAIL)
-	ENZO_FAIL("FeedbackZone copy failed!\n");
-    
-  /* Send data */
-
-  CommunicationDirection = COMMUNICATION_SEND;
-
-  for (i = 0; i < nParticles; i++) {
-    for (j = 0; j < NumberOfGrids; j++) {
-      if (FeedbackZones[i]->CopyActiveZonesFromGrid(Grids[j]->GridData,ZeroVector,SendField) == FAIL)
-	ENZO_FAIL("FeedbackZone copy failed!\n");
-    }
-  }
-
-  /* Receive data */
-
-  if (CommunicationReceiveHandler() == FAIL)
-    ENZO_FAIL("CommunicationReceiveHandler() failed!\n");
-
-#ifdef USE_MPI
-  CommunicationBufferPurge();
-#endif
-
-  return FeedbackZones;
-  
-}

src/enzo/hierarchy/DistributeFeedbackZone.C

+/***********************************************************************
+/
+/  Copy data from a 'fake' feedback zone grid back to 
+/  the real grids
+/
+/  written by: Nathan Goldbaum
+/  date:       June 2012
+/
+************************************************************************/
+
+#ifdef USE_MPI
+#include "communicators.h"
+#endif
+
+#include "preincludes.h"
+
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "phys_constants.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "ActiveParticle.h"
+#include "phys_constants.h"
+#include "CommunicationUtilities.h"
+#include "communication.h"
+
+int CommunicationBufferPurge(void);
+int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
+				int NumberOfSubgrids[] = NULL,
+				int FluxFlag = FALSE,
+				TopGridData* MetaData = NULL);
+
+int DistributeFeedbackZone(grid* FeedbackZone, HierarchyEntry** Grids, 
+			   int NumberOfGrids, int SendField)
+{
+  int i,j;
+
+  FLOAT ZeroVector[] = {0, 0, 0};
+
+  /* Post receives */
+  CommunicationReceiveIndex = 0;
+  CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE;
+  CommunicationDirection = COMMUNICATION_POST_RECEIVE;
+
+  for (i = 0; i < NumberOfGrids; i++) 
+    if (Grids[i]->GridData->CopyActiveZonesFromGrid(FeedbackZone,ZeroVector,SendField) == FAIL)
+      ENZO_FAIL("FeedbackZone copy failed!\n");
+    
+  /* Send data */
+    
+  CommunicationDirection = COMMUNICATION_SEND;
+
+  for (i = 0; i < NumberOfGrids; i++) 
+    if (Grids[i]->GridData->CopyActiveZonesFromGrid(FeedbackZone,ZeroVector,SendField) == FAIL)
+      ENZO_FAIL("FeedbackZone copy failed!\n");
+  
+  /* Receive data */
+  
+  if (CommunicationReceiveHandler() == FAIL)
+    ENZO_FAIL("CommunicationReceiveHandler() failed!\n");
+
+#ifdef USE_MPI
+  CommunicationBufferPurge();
+#endif
+
+  return SUCCESS;
+}

src/enzo/hierarchy/DistributeFeedbackZones.C

-/***********************************************************************
-/
-/  Copy data from a 'fake' feedback zone grid back to 
-/  the real grids
-/
-/  written by: Nathan Goldbaum
-/  date:       June 2012
-/
-************************************************************************/
-
-#ifdef USE_MPI
-#include "communicators.h"
-#endif
-
-#include "preincludes.h"
-
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "phys_constants.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "ActiveParticle.h"
-#include "phys_constants.h"
-#include "CommunicationUtilities.h"
-#include "communication.h"
-
-int CommunicationBufferPurge(void);
-int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
-				int NumberOfSubgrids[] = NULL,
-				int FluxFlag = FALSE,
-				TopGridData* MetaData = NULL);
-
-int DistributeFeedbackZones(grid** FeedbackZones, int NumberOfFeedbackZones,
-			    HierarchyEntry** Grids, int NumberOfGrids, int SendField)
-{
-  int i,j;
-
-  FLOAT ZeroVector[] = {0, 0, 0};
-
-  /* Post receives */
-  CommunicationReceiveIndex = 0;
-  CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE;
-  CommunicationDirection = COMMUNICATION_POST_RECEIVE;
-
-  for (i = 0; i < NumberOfGrids; i++) 
-    for (j = 0; j < NumberOfFeedbackZones; j++) 
-      if (Grids[i]->GridData->CopyActiveZonesFromGrid(FeedbackZones[j],ZeroVector,SendField) == FAIL)
-	ENZO_FAIL("FeedbackZone copy failed!\n");
-    
-  /* Send data */
-    
-  CommunicationDirection = COMMUNICATION_SEND;
-
-  for (i = 0; i < NumberOfGrids; i++) 
-    for (j = 0; j < NumberOfFeedbackZones; j++) 
-      if (Grids[i]->GridData->CopyActiveZonesFromGrid(FeedbackZones[j],ZeroVector,SendField) == FAIL)
-	ENZO_FAIL("FeedbackZone copy failed!\n");
-
-  /* Receive data */
-  
-  if (CommunicationReceiveHandler() == FAIL)
-    ENZO_FAIL("CommunicationReceiveHandler() failed!\n");
-
-#ifdef USE_MPI
-  CommunicationBufferPurge();
-#endif
-
-  return SUCCESS;
-}

src/enzo/hierarchy/Make.module.objects

 ENZO_OBJS +=\
 	AdjustMustRefineParticlesRefineToLevel.o \
 	AdjustRefineRegion.o \
-	ConstructFeedbackZones.o \
+	ConstructFeedbackZone.o \
 	CopyOverlappingZones.o \
 	CopyZonesFromOldGrids.o \
 	CreateSUBlingList.o \
 	CreateSiblingList.o \
 	DeleteSUBlingList.o \
 	DetermineSubgridSizeExtrema.o \
-	DistributeFeedbackZones.o \
+	DistributeFeedbackZone.o \
 	ExternalBoundary_AddField.o \
 	ExternalBoundary_AppendForcingToBaryonFields.o \
 	ExternalBoundary_DeleteObsoleteFields.o \

src/enzo/particles/active_particles/ActiveParticleRoutines.C

   GridID = part->GridID;
   type = part->type;
   CurrentGrid = part->CurrentGrid;
+  dest_processor = part->dest_processor;
 }
 
 ActiveParticleType::ActiveParticleType(grid *_grid, ActiveParticleFormationData &data)
   /* The correct indices are assigned in CommunicationUpdateActiveParticleCount 
      in ActiveParticleFinalize.*/
   Identifier = INT_UNDEFINED;
+  dest_processor = -1;
 }
 
 

src/enzo/particles/active_particles/ActiveParticle_AccretingParticle.C

 
 #endif
 
-#define NO_DEBUG
-
 /* We need to make sure that we can operate on the grid, so this dance is
  * necessary to make sure that grid is 'friend' to this particle type. */
 
   AccretingParticleGrid *thisGrid =
     static_cast<AccretingParticleGrid *>(thisgrid_orig);
   
-
-  float *tvel;
-
   int i,j,k,index,method,MassRefinementMethod;
 
   float *density = thisGrid->BaryonField[data.DensNum];
-  float *velx = thisGrid->BaryonField[data.Vel1Num];
-  float *vely = thisGrid->BaryonField[data.Vel2Num];
-  float *velz = thisGrid->BaryonField[data.Vel3Num];
+  
+  float* velx = thisGrid->BaryonField[data.Vel1Num]; 
+  float* vely = thisGrid->BaryonField[data.Vel2Num];
+  float* velz = thisGrid->BaryonField[data.Vel3Num];
+  
   float JeansDensityUnitConversion = (Gamma*pi*kboltz) / (Mu*mh*GravConst);
   float CellTemperature = 0;
   float JeansDensity = 0;
 	np->pos[1] = thisGrid->CellLeftEdge[1][j] + 0.5*thisGrid->CellWidth[1][j];
 	np->pos[2] = thisGrid->CellLeftEdge[2][k] + 0.5*thisGrid->CellWidth[2][k];
 	
-	tvel = thisGrid->AveragedVelocityAtCell(index,data.DensNum,data.Vel1Num);
-	  
-	np->vel[0] = tvel[0];
-	np->vel[1] = tvel[1];
-	np->vel[2] = tvel[2];
+	np->vel[0] = velx[index];
+	np->vel[1] = vely[index];
+	np->vel[2] = velz[index];
 	
 	if (HasMetalField)
 	  np->Metallicity = data.TotalMetals[index];
 
 	density[index] = DensityThreshold;
 
-	// Clean up
-	delete [] tvel;
-
       } // i
     } // j
   } // k
 }
 
 
-grid** ConstructFeedbackZones(ActiveParticleType** ParticleList, int nParticles,
-    int *FeedbackRadius, FLOAT dx, HierarchyEntry** Grids, int NumberOfGrids,
-    int SendField);
+grid* ConstructFeedbackZone(ActiveParticleType* ThisParticle, int FeedbackRadius, 
+			     FLOAT dx, HierarchyEntry** Grids, int NumberOfGrids,
+			     int SendField);
 
-int DistributeFeedbackZones(grid** FeedbackZones, int NumberOfFeedbackZones,
-			    HierarchyEntry** Grids, int NumberOfGrids, int SendField);
+int DistributeFeedbackZone(grid* FeedbackZone, HierarchyEntry** Grids, 
+			   int NumberOfGrids, int SendField);
 
 int ActiveParticleType_AccretingParticle::Accrete(int nParticles, ActiveParticleType** ParticleList,
 						  int AccretionRadius, FLOAT dx, 
   
   NumberOfGrids = GenerateGridArray(LevelArray, ThisLevel, &Grids);
   
-  // Some active particles have different feedback radii for each particle,
-  // so we need to set that up here, despite the fact that these accreting
-  // particles don't (right now).
-  FeedbackRadius = new int[nParticles];
   for (i = 0; i < nParticles; i++) {
-    FeedbackRadius[i] = AccretionRadius;
-  }
-  
-  grid** FeedbackZones = ConstructFeedbackZones(ParticleList, nParticles,
-    FeedbackRadius, dx, Grids, NumberOfGrids, ALL_FIELDS);
+    grid* FeedbackZone = ConstructFeedbackZone(ParticleList[i], AccretionRadius, 
+					       dx, Grids, NumberOfGrids, ALL_FIELDS);
 
-  delete FeedbackRadius;
-
-  for (i = 0; i < nParticles; i++) {
-    grid* FeedbackZone = FeedbackZones[i];
     if (MyProcessorNumber == FeedbackZone->ReturnProcessorNumber()) {
     
       float AccretionRate = 0;
-        
-      if (FeedbackZone->AccreteOntoAccretingParticle(&ParticleList[i],AccretionRadius*dx,&AccretionRate) == FAIL)
+      
+      if (FeedbackZone->AccreteOntoAccretingParticle(&ParticleList[i],
+			      AccretionRadius*dx,&AccretionRate) == FAIL)
 	return FAIL;
-  
+      
       // No need to communicate the accretion rate to the other CPUs since this particle is already local.
       static_cast<ActiveParticleType_AccretingParticle*>(ParticleList[i])->AccretionRate = AccretionRate;
     }
+  
+    DistributeFeedbackZone(FeedbackZone, Grids, NumberOfGrids, ALL_FIELDS);
+
+    delete FeedbackZone;
   }
   
-  DistributeFeedbackZones(FeedbackZones, nParticles, Grids, NumberOfGrids, ALL_FIELDS);
-
-  for (i = 0; i < nParticles; i++) {
-    delete FeedbackZones[i];    
-  }
-
-  delete [] FeedbackZones;
-
   if (AssignActiveParticlesToGrids(ParticleList, nParticles, LevelArray) == FAIL)
     return FAIL;
 
 std::vector<ParticleAttributeHandler*>
   ActiveParticleType_AccretingParticle::AttributeHandlers;
 
-#undef DEBUG

src/enzo/particles/active_particles/ActiveParticle_GalaxyParticle.C

   return SUCCESS;
 }
 
-grid** ConstructFeedbackZones(ActiveParticleType** ParticleList, int nParticles,
-    int *FeedbackRadius, FLOAT dx, HierarchyEntry** Grids, int NumberOfGrids,
-    int SendField);
+grid* ConstructFeedbackZone(ActiveParticleType* ThisParticle, int FeedbackRadius, 
+			    FLOAT dx, HierarchyEntry** Grids, int NumberOfGrids,
+			    int SendField);
 
-int DistributeFeedbackZones(grid** FeedbackZones, int NumberOfFeedbackZones,
-			    HierarchyEntry** Grids, int NumberOfGrids, int SendField);
+int DistributeFeedbackZone(grid* FeedbackZones, HierarchyEntry** Grids, 
+			   int NumberOfGrids, int SendField);
 
 int ActiveParticleType_GalaxyParticle::SubtractMassFromGrid(int nParticles,
     ActiveParticleType** ParticleList, LevelHierarchyEntry *LevelArray[],
     FeedbackRadius[i] = nint(static_cast<ActiveParticleType_GalaxyParticle*>(ParticleList[i])->Radius / dx);
   }
   
-  grid** FeedbackZones = ConstructFeedbackZones(ParticleList, nParticles,
-    FeedbackRadius, dx, Grids, NumberOfGrids, ALL_FIELDS);
+  for (i = 0; i < nParticles; i++) {
+    grid* FeedbackZone = ConstructFeedbackZone(ParticleList[i], FeedbackRadius[i], dx, 
+					       Grids, NumberOfGrids, ALL_FIELDS);
 
-  delete [] FeedbackRadius;
-
-  for (i = 0; i < nParticles; i++) {
-    grid* FeedbackZone = FeedbackZones[i];
     if (MyProcessorNumber == FeedbackZone->ReturnProcessorNumber()) {
         
       if (FeedbackZone->ApplyGalaxyParticleFeedback(&ParticleList[i]) == FAIL)
 	return FAIL;
-  
+
     }
-  }
-  
-  DistributeFeedbackZones(FeedbackZones, nParticles, Grids, NumberOfGrids, ALL_FIELDS);
 
-  for (i = 0; i < nParticles; i++) {
-    delete FeedbackZones[i];    
+    DistributeFeedbackZone(FeedbackZone, Grids, NumberOfGrids, ALL_FIELDS);
+
+    delete FeedbackZone;
   }
 
-  delete [] FeedbackZones;
+  delete [] FeedbackRadius;
+  
   if (AssignActiveParticlesToGrids(ParticleList, nParticles, LevelArray) == FAIL)
     return FAIL;
 
     FeedbackRadius[i] = nint(static_cast<ActiveParticleType_GalaxyParticle*>(ParticleList[i])->Radius / dx);
   }
   
-  grid** FeedbackZones = ConstructFeedbackZones(ParticleList, nParticles,
-    FeedbackRadius, dx, Grids, NumberOfGrids, GRAVITATING_MASS_FIELD);
+  for (i = 0; i < nParticles; i++) {
+    grid* FeedbackZone = ConstructFeedbackZone(ParticleList[i], FeedbackRadius[i], 
+				     dx, Grids, NumberOfGrids, GRAVITATING_MASS_FIELD);
 
-  delete [] FeedbackRadius;
-
-  for (i = 0; i < nParticles; i++) {
-    grid* FeedbackZone = FeedbackZones[i];
     if (MyProcessorNumber == FeedbackZone->ReturnProcessorNumber()) {
         
       if (FeedbackZone->ApplyGalaxyParticleGravity(&ParticleList[i]) == FAIL)
 	return FAIL;
-  
+
     }
-  }
-  
-  DistributeFeedbackZones(FeedbackZones, nParticles, Grids, NumberOfGrids,
-    GRAVITATING_MASS_FIELD);
 
-  for (i = 0; i < nParticles; i++) {
-    delete FeedbackZones[i];    
+    DistributeFeedbackZone(FeedbackZone, Grids, NumberOfGrids, GRAVITATING_MASS_FIELD);
+    
+    delete FeedbackZone;
   }
 
-  delete [] FeedbackZones;
+  delete [] FeedbackRadius;
+  
   delete [] Grids;
   return SUCCESS;
 }
 }
 
 std::vector<ParticleAttributeHandler*>
-  ActiveParticleType_GalaxyParticle::AttributeHandlers;
+  ActiveParticleType_GalaxyParticle::AttributeHandlers;

src/enzo/particles/active_particles/AssignActiveParticlesToGrids.C

 	      OldGrid->CleanUpMovedParticles();
 	    }
       }
-      // If the particle didn't change grids, we still need to mirror the AP data to the particle list.
+      // The particle didn't move, still need to update the AP data and mirror the new data to the particle list.
       else {
 	LevelGrids[SavedGrid]->GridData->UpdateParticleWithActiveParticle(ParticleList[i]->ReturnID());
       }
 	    if (LevelGrids[SavedGrid]->GridData->AddActiveParticle(static_cast<ActiveParticleType*>(ParticleList[i])) == FAIL) {
 	      ENZO_FAIL("Active particle grid assignment failed"); 
 	    } 
-	    // Still need to mirror the AP data to the particle list.
+	    // Still need to update the AP data and mirror the new data to the particle list.
 	    else {
 	      LevelGrids[SavedGrid]->GridData->UpdateParticleWithActiveParticle(ParticleList[i]->ReturnID());
 	    }