Commits

Nathan Goldbaum committed 1beb25e

Intermediate commit: fixing a bug in sink particle accretion.

Comments (0)

Files changed (12)

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 DEBUG_AP
+
 void my_exit(int status);
  
 // function prototypes

src/enzo/control/EvolveLevel.C

  
     /* Finalize (accretion, feedback, etc.) star particles */
 
+    float mass1 = 0;
+
+    for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
+      Grids[grid1]->GridData->SumGasMass(&mass1);
+      printf("mass = %"FSYM" \n", mass1);
+    }
+
     ActiveParticleFinalize(Grids, MetaData, NumberOfGrids, LevelArray,
 			   level, NumberOfNewActiveParticles);
 
+    float mass2 = 0;
+
+    for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
+      Grids[grid1]->GridData->SumGasMass(&mass2);
+      printf("mass = %"FSYM" \n", mass2);
+    }
+
+    printf("mass1 - mass2 = %"FSYM" \n", mass1 - mass2);
+
+    for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
+      for (int pn = 0; pn < Grids[grid1]->GridData->ReturnNumberOfActiveParticles(); pn++) {
+	float mass = Grids[grid1]->GridData->ReturnActiveParticles()[pn]->ReturnMass();
+	printf("ActiveParticles[%"ISYM"]->mass = %"FSYM" \n", pn, mass);
+      }
+    }
+
     /* For each grid: a) interpolate boundaries from the parent grid.
                       b) copy any overlapping zones from siblings. */
  

src/enzo/grid/particles/Grid_AccreteOntoAccretingParticle.C

   /* Return if this doesn't involve us */
   if (MyProcessorNumber != ProcessorNumber) 
     return SUCCESS;
-  
+
+  printf("InitialMass    = %"FSYM" \n", (*ThisParticle)->ReturnMass());
+
   /* Check whether the cube that circumscribes the accretion zone intersects with this grid */
 
   FLOAT *ParticlePosition = (*ThisParticle)->ReturnPosition();
   (*ThisParticle)->AddMass(AccretedMass);
   (*ThisParticle)->SetVelocity(NewVelocity);
 
+  printf("AccretedMass = %"FSYM" \n", AccretedMass);
+  printf("FinalMass    = %"FSYM" \n", (*ThisParticle)->ReturnMass());
+
   delete [] nexcluded;
 
   return SUCCESS;

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/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);
 

src/enzo/hierarchy/ConstructFeedbackZones.C

   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];
+    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[i]+DEFAULT_GHOST_ZONES)+1);

src/enzo/hierarchy/DistributeFeedbackZones.C

 #include "CommunicationUtilities.h"
 #include "communication.h"
 
+#define DEBUG_AP
+
 int CommunicationBufferPurge(void);
 int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
 				int NumberOfSubgrids[] = NULL,

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

 int DistributeFeedbackZones(grid** FeedbackZones, int NumberOfFeedbackZones,
 			    HierarchyEntry** Grids, int NumberOfGrids, int SendField);
 
+void CheckFeedbackZones(grid** FeedbackZones, int NumberOfFeedbackZones) {
+  for (int i = 0; i < NumberOfFeedbackZones; i++)
+    for (int j = 0; j < NumberOfFeedbackZones; j++) {
+      if (i == j)
+	continue;
+      grid* ThisFeedbackZone = FeedbackZones[i];
+      grid* OtherFeedbackZone = FeedbackZones[j];
+      if (!(ThisFeedbackZone->GridLeftEdge[0] > OtherFeedbackZone->GridRightEdge[0] ||
+	    ThisFeedbackZone->GridLeftEdge[0] > OtherFeedbackZone->GridRightEdge[1] ||
+	    ThisFeedbackZone->GridLeftEdge[0] > OtherFeedbackZone->GridRightEdge[2] ||
+	    ThisFeedbackZone->GridRightEdge[0] < OtherFeedbackZone->GridLeftEdge[0] ||
+	    ThisFeedbackZone->GridRightEdge[1] < OtherFeedbackZone->GridLeftEdge[1] ||
+	    ThisFeedbackZone->GridRightEdge[2] < OtherFeedbackZone->GridLeftEdge[2]))
+	ENZO_FAIL("Overlapping feedback zones!");
+    }
+}
+
+#define DEBUG_AP
+
 int ActiveParticleType_AccretingParticle::Accrete(int nParticles, ActiveParticleType** ParticleList,
 						  int AccretionRadius, FLOAT dx, 
 						  LevelHierarchyEntry *LevelArray[], int ThisLevel)
     FeedbackRadius[i] = AccretionRadius;
   }
   
+  for (i = 0; i < nParticles; i++)
+    printf("ActiveParticles[%"ISYM"]->mass = %"FSYM" \n", i, Grids[0]->GridData->ReturnActiveParticles()[i]->ReturnMass());
+
   grid** FeedbackZones = ConstructFeedbackZones(ParticleList, nParticles,
     FeedbackRadius, dx, Grids, NumberOfGrids, ALL_FIELDS);
 
+  CheckFeedbackZones(FeedbackZones, nParticles);
+
   delete FeedbackRadius;
 
   for (i = 0; i < nParticles; i++) {
     if (MyProcessorNumber == FeedbackZone->ReturnProcessorNumber()) {
     
       float AccretionRate = 0;
-        
-      if (FeedbackZone->AccreteOntoAccretingParticle(&ParticleList[i],AccretionRadius*dx,&AccretionRate) == FAIL)
+
+      float mass = 0;
+      FeedbackZone->SumGasMass(&mass);
+      printf("Feedbackzone %"ISYM" init mass = %"FSYM" \n", i, mass);
+
+      float initmass = ParticleList[i]->ReturnMass();
+
+      if (FeedbackZone->AccreteOntoAccretingParticle(&ParticleList[i],
+				    AccretionRadius*dx,&AccretionRate) == FAIL)
 	return FAIL;
-  
+      
+      mass = 0;
+      FeedbackZone->SumGasMass(&mass);
+      printf("Feedbackzone %"ISYM" final mass = %"FSYM" \n", i, mass+(ParticleList[i]->ReturnMass() - initmass));
+
       // 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;
     }
   if (AssignActiveParticlesToGrids(ParticleList, nParticles, LevelArray) == FAIL)
     return FAIL;
 
+  for (i = 0; i < nParticles; i++)
+    printf("ActiveParticles[%"ISYM"]->mass = %"FSYM" \n", i, Grids[0]->GridData->ReturnActiveParticles()[i]->ReturnMass());
+
   delete [] Grids;
   return SUCCESS;
 }

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());
 	    }