Commits

John Wise committed 58c31a2

Porting all of the relevant Star routines to active particles.

  • Participants
  • Parent commits b22df5a

Comments (0)

Files changed (6)

File src/enzo/ActiveParticle.h

   
   /* This should return the number of new star particles created, and should
    * create them. */
-  ActiveParticleType(){};
-  ~ActiveParticleType(){};
-  ActiveParticleType(ActiveParticleType*& part){};
-  virtual int GetEnabledParticleID(int id = -1) = 0;
-  star_type ReturnType(void) {return type; };
-#ifdef ACTIVE_PARTICLE_IMPLEMENTED
+//  ActiveParticleType(){};
+//  ~ActiveParticleType(){};
+//  ActiveParticleType(ActiveParticleType*& part){};
+
+  ActiveParticleType(void);
   ActiveParticleType(grid *_grid, int _id, int _level);
-  ActiveParticleType(StarBuffer *buffer, int n);
-  ActiveParticleType(StarBuffer buffer) ;
-  ActiveParticleType* copy(void);
+  ActiveParticleType(ActiveParticleType*& part);
+  ~ActiveParticleType(void);
+
+  void operator=(ActiveParticleType *a);
+
+  template <class active_particle_class> active_particle_class *copy(void);
 
   int   ReturnID(void) { return Identifier; };
   double ReturnMass(void) { return Mass; };
   float ReturnBirthTime(void) { return BirthTime; };
   float ReturnDynamicalTime(void) { return DynamicalTime; };
   star_type ReturnType(void) {return type; };
+  int   ReturnLevel(void) { return level; };
 
-  int   ReturnLevel(void) { return level; };
   void  ReduceLevel(void) { level--; };
   void  IncreaseLevel(void) { level++; };
   void  SetLevel(int i) { level = i; };
 
   FLOAT *ReturnPosition(void) { return pos; }
   float *ReturnVelocity(void) { return vel; }
-  void    ConvertMassToSolar(void);
-  void    Merge(Star a);
-  void    Merge(Star *a);
-  bool    Mergable(Star a);
-  bool  Mergable(Star *a);
-  float Separation(Star a);
-  float Separation(Star *a);
-  float Separation2(Star a);
-  float Separation2(Star *a);
-  float RelativeVelocity2(Star a);
-  float RelativeVelocity2(Star *a);
+  void   ConvertAllMassesToSolar(void);
+  void   ConvertMassToSolar(void);
+  void   Merge(ActiveParticleType *a);
+  bool  Mergable(ActiveParticleType *a);
+  bool  MergableMBH(ActiveParticleType *a);
+  float Separation(ActiveParticleType *a);
+  float Separation2(ActiveParticleType *a);
+  float RelativeVelocity2(ActiveParticleType *a);
   void  UpdatePositionVelocity(void);
-  void    DeleteCopyInGrid(void);
-  int   DeleteCopyInGridGlobal(LevelHierarchyEntry *LevelArray[]);
-  void    CopyToGrid(void);
   void  MirrorToParticle(void);
-  virtual bool  IsARadiationSource(FLOAT Time) { return FALSE };
-  int   DeleteParticle(LevelHierarchyEntry *LevelArray[]);
+  void  CopyFromParticle(grid *_grid, int _id, int _level);
   int   DisableParticle(LevelHierarchyEntry *LevelArray[]);
-  void  ActivateNewStar(FLOAT Time, float Timestep);
+  int   SphereContained(LevelHierarchyEntry *LevelArray[], int level, 
+			float Radius);
+  void  PrintInfo(void);
+  //void  ActivateNewStar(FLOAT Time, float Timestep);
+  //void  DeleteCopyInGrid(void);
+  //int   DeleteCopyInGridGlobal(LevelHierarchyEntry *LevelArray[]);
 
-  int SphereContained(LevelHierarchyEntry *LevelArray[], int level, 
-              float Radius);
+  /* Virtual and pure virtual functions in this base class */
 
-  ActiveParticle* StarBufferToList(StarBuffer *buffer, int n);
-  StarBuffer* StarListToBuffer(int n);
+  virtual bool IsARadiationSource(FLOAT Time) { return FALSE; };
+  virtual int GetEnabledParticleID(int id = -1) = 0;
 
-#endif /* IMPLEMENTED */
-    
+#ifdef TRANSFER
+  RadiationSourceEntry* RadiationSourceInitialize(void);
+#endif
+
 protected:
   grid *CurrentGrid;
   FLOAT	pos[MAX_DIMENSION];

File src/enzo/ActiveParticleRoutines.C

 /           to them.
 /
 ************************************************************************/
+#include <map>
+#include <iostream>
+#include <stdexcept>
 #include <assert.h>
 #include <stdlib.h>
 #include <stdio.h>
 #include "Hierarchy.h"
 #include "TopGridData.h"
 #include "LevelHierarchy.h"
-
-void DeleteStar(ActiveParticleType * &Node);
-ActiveParticleType *PopStar(ActiveParticleType * &Node);
-void InsertStarAfter(ActiveParticleType * &Node, 
-		     ActiveParticleType * &NewNode);
+#include "ActiveParticle.h"
 
 int GetUnits(float *DensityUnits, float *LengthUnits,
 	     float *TemperatureUnits, float *TimeUnits,
 ActiveParticleType::ActiveParticleType(void)
 {
   int dim;
+  CurrentGrid = NULL;
   for (dim = 0; dim < MAX_DIMENSION; dim++)
-    pos[dim] = vel[dim] = delta_vel[dim] = accreted_angmom[dim] = 0.0;
-  accretion_rate = NULL;
-  accretion_time = NULL;
-  NextStar = NULL;
-  PrevStar = NULL;
+    pos[dim] = vel[dim] = 0.0;
+  Mass = BirthTime = DynamicalTime = 0.0;
+  Identifier = level = GridID = type = 0;
+}
+
+ActiveParticleType::ActiveParticleType(ActiveParticleType*& part)
+{
+  int dim;
   CurrentGrid = NULL;
-  Mass = FinalMass = DeltaMass = BirthTime = LifeTime = 
-    last_accretion_rate = NotEjectedMass = Metallicity = deltaZ = 0.0;
-  FeedbackFlag = Identifier = level = GridID = type = naccretions = 0;
-  AddedEmissivity = false;
+  for (dim = 0; dim < MAX_DIMENSION; dim++)
+    pos[dim] = vel[dim] = 0.0;
+  Mass = BirthTime = DynamicalTime = 0.0;
+  Identifier = level = GridID = type = 0;
 }
 
 ActiveParticleType::ActiveParticleType(grid *_grid, int _id, int _level)
   for (dim = 0; dim < MAX_DIMENSION; dim++) {
     pos[dim] = _grid->ParticlePosition[dim][_id];
     vel[dim] = _grid->ParticleVelocity[dim][_id];
-    delta_vel[dim] = 0.0;
-    accreted_angmom[dim] = 0.0;
   }
-  naccretions = 0;
-  accretion_rate = NULL;
-  accretion_time = NULL;
-  NextStar = NULL;
-  PrevStar = NULL;
   CurrentGrid = _grid;
-  DeltaMass = 0.0;
-  AddedEmissivity = false;
-  last_accretion_rate = 0.0;
-  NotEjectedMass = 0.0;
-  deltaZ = 0.0;
   level = _level;
-  FeedbackFlag = NO_FEEDBACK;
 
   GridID = _grid->ID;
   type = _grid->ParticleType[_id];
   Identifier = _grid->ParticleNumber[_id];
-  Mass = FinalMass = (double)(_grid->ParticleMass[_id]);
+  Mass = (double)(_grid->ParticleMass[_id]);
   BirthTime = _grid->ParticleAttribute[0][_id];
-  LifeTime = _grid->ParticleAttribute[1][_id];
+  DynamicalTime = _grid->ParticleAttribute[1][_id];
   Metallicity = _grid->ParticleAttribute[2][_id];
-  this->ConvertAllMassesToSolar();
-}
-
-ActiveParticleType::ActiveParticleType(StarBuffer *buffer, int n) 
-{
-  int i;
-  CurrentGrid = NULL;
-  for (i = 0; i < MAX_DIMENSION; i++) {
-    pos[i] = buffer[n].pos[i];
-    vel[i] = buffer[n].vel[i];
-    delta_vel[i] = buffer[n].delta_vel[i];
-    accreted_angmom[i] = buffer[n].accreted_angmom[i];
-  }
-  naccretions = min(buffer[n].naccretions, MAX_ACCR);
-  if (naccretions > 0) {
-    accretion_time = new FLOAT[naccretions];
-    accretion_rate = new float[naccretions];
-    for (i = 0; i < naccretions; i++) {
-      accretion_time[i] = buffer[n].accretion_time[i];
-      accretion_rate[i] = buffer[n].accretion_rate[i];
-    }
-  } else {
-    accretion_time = NULL;
-    accretion_rate = NULL;
-  }
-  Mass = buffer[n].Mass;
-  FinalMass = buffer[n].FinalMass;
-  DeltaMass = buffer[n].DeltaMass;
-  BirthTime = buffer[n].BirthTime;
-  LifeTime = buffer[n].LifeTime;
-  Metallicity = buffer[n].Metallicity;
-  deltaZ = buffer[n].deltaZ;
-  last_accretion_rate = buffer[n].last_accretion_rate;
-  NotEjectedMass = buffer[n].NotEjectedMass;
-  FeedbackFlag = buffer[n].FeedbackFlag;
-  Identifier = buffer[n].Identifier;
-  level = buffer[n].level;
-  GridID = buffer[n].GridID;
-  type = buffer[n].type;
-  AddedEmissivity = buffer[n].AddedEmissivity;
-  NextStar = NULL;
-  PrevStar = NULL;
-}
-
-ActiveParticleType::ActiveParticleType(StarBuffer buffer) 
-{
-  int i;
-  CurrentGrid = NULL;
-  for (i = 0; i < MAX_DIMENSION; i++) {
-    pos[i] = buffer.pos[i];
-    vel[i] = buffer.vel[i];
-    delta_vel[i] = buffer.delta_vel[i];
-    accreted_angmom[i] = buffer.accreted_angmom[i];
-  }
-  naccretions = min(buffer.naccretions, MAX_ACCR);
-  if (naccretions > 0) {
-    accretion_time = new FLOAT[naccretions];
-    accretion_rate = new float[naccretions];
-    for (i = 0; i < naccretions; i++) {
-      accretion_time[i] = buffer.accretion_time[i];
-      accretion_rate[i] = buffer.accretion_rate[i];
-    }
-  } else {
-    accretion_time = NULL;
-    accretion_rate = NULL;
-  }
-  Mass = buffer.Mass;
-  FinalMass = buffer.FinalMass;
-  DeltaMass = buffer.DeltaMass;
-  BirthTime = buffer.BirthTime;
-  LifeTime = buffer.LifeTime;
-  Metallicity = buffer.Metallicity;
-  deltaZ = buffer.deltaZ;
-  last_accretion_rate = buffer.last_accretion_rate;
-  NotEjectedMass = buffer.NotEjectedMass;
-  FeedbackFlag = buffer.FeedbackFlag;
-  Identifier = buffer.Identifier;
-  level = buffer.level;
-  GridID = buffer.GridID;
-  type = buffer.type;
-  NextStar = NULL;
-  PrevStar = NULL;
+  this->ConvertMassToSolar();
 }
 
 /* No need to delete the accretion arrays because the pointers are
 
 ActiveParticleType::~ActiveParticleType(void)
 {
-  if (accretion_rate != NULL)
-    delete [] accretion_rate;
-  if (accretion_time != NULL)
-    delete [] accretion_time;
-  NextStar = NULL;
-  PrevStar = NULL;
-  CurrentGrid = NULL;
 }
 
 /***************
 
  ***************/
 
-void ActiveParticleType::operator=(ActiveParticleType a)
+void ActiveParticleType::operator=(ActiveParticleType *a)
 {
   int i, dim;
-  //NextStar = a.NextStar;
-  CurrentGrid = a.CurrentGrid;
+  CurrentGrid = a->CurrentGrid;
   for (dim = 0; dim < MAX_DIMENSION; dim++) {
-    pos[dim] = a.pos[dim];
-    vel[dim] = a.vel[dim];
-    delta_vel[dim] = a.delta_vel[dim];
-    accreted_angmom[dim] = a.accreted_angmom[dim];
+    pos[dim] = a->pos[dim];
+    vel[dim] = a->vel[dim];
   }
-  naccretions = a.naccretions;
-  Mass = a.Mass;
-  FinalMass = a.FinalMass;
-  DeltaMass = a.DeltaMass;
-  BirthTime = a.BirthTime;
-  LifeTime = a.LifeTime;
-  Metallicity = a.Metallicity;
-  deltaZ = a.deltaZ;
-  last_accretion_rate = a.last_accretion_rate;
-  NotEjectedMass = a.NotEjectedMass;
-  FeedbackFlag = a.FeedbackFlag;
-  Identifier = a.Identifier;
-  level = a.level;
-  GridID = a.GridID;
-  type = a.type;
-  AddedEmissivity = a.AddedEmissivity;
-  if (accretion_rate != NULL)
-    delete [] accretion_rate;
-  if (accretion_time != NULL)
-    delete [] accretion_time;
-  if (naccretions > 0) {
-    accretion_rate = new float[naccretions];
-    accretion_time = new FLOAT[naccretions];
-    for (i = 0; i < naccretions; i++) {
-      accretion_rate[i] = a.accretion_rate[i];
-      accretion_time[i] = a.accretion_time[i];
-    }
-  } else {
-    accretion_rate = NULL;
-    accretion_time = NULL;
-  }
+  Mass = a->Mass;
+  BirthTime = a->BirthTime;
+  DynamicalTime = a->DynamicalTime;
+  Metallicity = a->Metallicity;
+  Identifier = a->Identifier;
+  level = a->level;
+  GridID = a->GridID;
+  type = a->type;
   return;
 }
 
 
  **********************/
 
-ActiveParticleType *ActiveParticleType::copy(void)
+template<class active_particle_class>
+active_particle_class *ActiveParticleType::copy(void)
 {
   int i, dim;
-  ActiveParticleType *a = new ActiveParticleType;
-  a->NextStar = NULL;
-  a->PrevStar = NULL;
+  active_particle_class *a = new active_particle_class();
   a->CurrentGrid = CurrentGrid;
   for (dim = 0; dim < MAX_DIMENSION; dim++) {
     a->pos[dim] = pos[dim];
     a->vel[dim] = vel[dim];
-    a->delta_vel[dim] = delta_vel[dim];
-    a->accreted_angmom[dim] = accreted_angmom[dim];
   }
-  a->naccretions = naccretions;
   a->Mass = Mass;
-  a->FinalMass = FinalMass;
-  a->DeltaMass = DeltaMass;
   a->BirthTime = BirthTime;
-  a->LifeTime = LifeTime;
+  a->DynamicalTime = DynamicalTime;
   a->Metallicity = Metallicity;
-  a->deltaZ = deltaZ;
-  a->last_accretion_rate = last_accretion_rate;
-  a->NotEjectedMass = NotEjectedMass;
-  a->FeedbackFlag = FeedbackFlag;
   a->Identifier = Identifier;
   a->level = level;
   a->GridID = GridID;
   a->type = type;
-  a->AddedEmissivity = AddedEmissivity;
-  if (naccretions > 0) {
-    a->accretion_rate = new float[naccretions];
-    a->accretion_time = new FLOAT[naccretions];
-    for (i = 0; i < naccretions; i++) {
-      a->accretion_rate[i] = accretion_rate[i];
-      a->accretion_time[i] = accretion_time[i];
-    }
-  } else {
-    a->accretion_rate = NULL;
-    a->accretion_time = NULL;
-  }
   return a;
 }
 
   dx = LengthUnits * CurrentGrid->CellWidth[0][0];
   MassConversion = (float) (dx*dx*dx * double(DensityUnits) / Msun);
   this->Mass *= MassConversion;
-  this->FinalMass *= MassConversion;
   return;
 }
 
   return;
 }
 
-void ActiveParticleType::Merge(ActiveParticleType a)
+void ActiveParticleType::Merge(ActiveParticleType *a)
 {
   int dim;
   double ratio1, ratio2;
-  ratio1 = Mass / (Mass + a.Mass);
+  ratio1 = Mass / (Mass + a->Mass);
   ratio2 = 1.0 - ratio1;
-  Metallicity = ratio1 * Metallicity + ratio2 * a.Metallicity;
+  Metallicity = ratio1 * Metallicity + ratio2 * a->Metallicity;
   for (dim = 0; dim < MAX_DIMENSION; dim++) {
-    pos[dim] = ratio1 * pos[dim] + ratio2 * a.pos[dim];
-    vel[dim] = ratio1 * vel[dim] + ratio2 * a.vel[dim];
-    accreted_angmom[dim] = ratio1 * accreted_angmom[dim] + ratio2 * a.accreted_angmom[dim];
+    pos[dim] = ratio1 * pos[dim] + ratio2 * a->pos[dim];
+    vel[dim] = ratio1 * vel[dim] + ratio2 * a->vel[dim];
+    //accreted_angmom[dim] = ratio1 * accreted_angmom[dim] + ratio2 * a.accreted_angmom[dim];
   }
-  Mass += a.Mass;
+  Mass += a->Mass;
   //FinalMass += a.FinalMass;
-  DeltaMass += a.DeltaMass;
-  last_accretion_rate += a.last_accretion_rate;
-  NotEjectedMass += a.NotEjectedMass;
+  //DeltaMass += a.DeltaMass;
+  //last_accretion_rate += a.last_accretion_rate;
+  //NotEjectedMass += a.NotEjectedMass;
   return;
 }
-void ActiveParticleType::Merge(ActiveParticleType *a) { this->Merge(*a); };
 
-bool ActiveParticleType::Mergable(ActiveParticleType a)
+bool ActiveParticleType::Mergable(ActiveParticleType *a)
 {
   // Only merge yet-to-be born stars
-  return type == a.type && type < 0;
+  return type == a->type && type < 0;
 }
-bool ActiveParticleType::Mergable(ActiveParticleType *a) { return this->Mergable(*a); }
 
-bool ActiveParticleType::MergableMBH(ActiveParticleType a)
+bool ActiveParticleType::MergableMBH(ActiveParticleType *a)
 {
   // Merge MBH particle with another 
-  return type == a.type && type == MBH;
+  return type == a->type && type == MBH;
 }
-bool ActiveParticleType::MergableMBH(ActiveParticleType *a) { return this->MergableMBH(*a); }
 
-float ActiveParticleType::Separation2(ActiveParticleType a)
+float ActiveParticleType::Separation2(ActiveParticleType *a)
 {
   int dim;
   float dr, result = 0;
   for (dim = 0; dim < MAX_DIMENSION; dim++) {
-    dr = pos[dim] - a.pos[dim];
+    dr = pos[dim] - a->pos[dim];
     result += dr*dr;
   }
   return result;
 }
-float ActiveParticleType::Separation2(ActiveParticleType *a) { return this->Separation2(*a); };
 
-float ActiveParticleType::Separation(ActiveParticleType a)  { return sqrt(this->Separation2(a)); }
-float ActiveParticleType::Separation(ActiveParticleType *a) { return this->Separation(*a); };
+float ActiveParticleType::Separation(ActiveParticleType *a)  { return sqrt(this->Separation2(a)); }
 
-float ActiveParticleType::RelativeVelocity2(ActiveParticleType a)
+float ActiveParticleType::RelativeVelocity2(ActiveParticleType *a)
 {
   int dim;
   float dv, result = 0;
   for (dim = 0; dim < MAX_DIMENSION; dim++) {
-    dv = vel[dim] - a.vel[dim];
+    dv = vel[dim] - a->vel[dim];
     result += dv*dv;
   }
   return result;
 }
-float ActiveParticleType::RelativeVelocity2(ActiveParticleType *a) { return this->RelativeVelocity2(*a); };
-
-void ActiveParticleType::CopyToGrid()
-{
-  ActiveParticleType *cstar;
-  if (CurrentGrid != NULL)   // NULL => On another processor
-    for (cstar = CurrentGrid->Stars; cstar; cstar = cstar->NextStar)
-      if (Identifier == cstar->Identifier) {
-	//cstar = this->copy();
-	*cstar = *this;
-	break;
-      } // ENDIF match
-  return;
-}
 
 void ActiveParticleType::UpdatePositionVelocity(void)
 {
   level = _level;
   GridID = _grid->ID;
   BirthTime = _grid->ParticleAttribute[0][_id];
-  LifeTime = _grid->ParticleAttribute[1][_id];
+  DynamicalTime = _grid->ParticleAttribute[1][_id];
   Metallicity = _grid->ParticleAttribute[2][_id];
 
   // below is removed because we want to keep Star->Mass as double 
   return;
 }
 
-void ActiveParticleType::DeleteCopyInGrid(void)
-{
-  ActiveParticleType *cstar, *MoveStar;
-  if (CurrentGrid != NULL) {   // NULL => On another processor
-    cstar = CurrentGrid->Stars;
-    CurrentGrid->Stars = NULL;
-    while (cstar) {
-      MoveStar = PopStar(cstar);
-      if (Identifier == MoveStar->Identifier)
-	delete MoveStar;
-      else
-	InsertStarAfter(CurrentGrid->Stars, MoveStar);
-    } // ENDWHILE stars
-  } // ENDIF grid != NULL
-  return;
-}
-
 void ActiveParticleType::PrintInfo(void)
 {
-  printf("[P%d] Star %"ISYM": pos = %"PSYM" %"PSYM" %"PSYM", vel = %"FSYM" %"FSYM" %"FSYM"\n",
+  printf("[P%d] ActiveParticle %"ISYM": pos = %"PSYM" %"PSYM" %"PSYM", vel = %"FSYM" %"FSYM" %"FSYM"\n",
 	 MyProcessorNumber, Identifier, pos[0], pos[1], pos[2], vel[0], vel[1], vel[2]);
-  printf("\t delta_vel = %"FSYM" %"FSYM" %"FSYM"\n", delta_vel[0], delta_vel[1],
-	 delta_vel[2]);
-  printf("\t naccr = %"ISYM, naccretions);
-  if (naccretions > 0)
-    printf(", accr_rate[0] = %"GSYM", accr_time[0] = %"GSYM"\n", 
-	   accretion_rate[0], accretion_time[0]);
-  else
-    printf("\n");
-  printf("\t birthtime = %"FSYM", lifetime = %"FSYM"\n", BirthTime, LifeTime);
-  printf("\t Z = %"GSYM", deltaZ = %"GSYM"\n", Metallicity, deltaZ);
-  printf("\t mass = %"GSYM", dmass = %"GSYM", fmass = %"GSYM", type = %"ISYM", grid %"ISYM","
-	 " lvl %"ISYM"\n", Mass, DeltaMass, FinalMass, type, GridID, level);
-  printf("\t FeedbackFlag = %"ISYM"\n", FeedbackFlag);
-  printf("\t accreted_angmom = %"FSYM" %"FSYM" %"FSYM"\n", accreted_angmom[0],
-	 accreted_angmom[1], accreted_angmom[2]);
-  printf("\t this = %x, PrevStar = %x, NextStar = %x\n", this, PrevStar, NextStar);
+  printf("\t birthtime = %"FSYM", tdyn = %"FSYM"\n", BirthTime, DynamicalTime);
+  printf("\t Z = %"GSYM"\n", Metallicity);
+  printf("\t mass = %"GSYM", type = %"ISYM", grid %"ISYM", lvl %"ISYM"\n", 
+	 Mass, type, GridID, level);
   return;
 }
 
   source->GridID         = GridID;
   source->GridLevel      = level;
   source->Type           = type;
-  source->LifeTime       = LifeTime;
+  source->LifeTime       = DynamicalTime;
   source->CreationTime   = BirthTime;
   source->AddedEmissivity = false;
   source->Position       = new FLOAT[3];

File src/enzo/ActiveParticle_DisableParticle.C

+/***********************************************************************
+/
+/  DISABLE THE ASSOCIATED NORMAL PARTICLE (set mass to tiny)
+/
+/  written by: John Wise
+/  date:       December, 2009
+/  modified1:
+/
+************************************************************************/
+
+#include <map>
+#include <iostream>
+#include <stdexcept>
+
+#ifdef USE_MPI
+#include "mpi.h"
+#endif
+#include <stdlib.h>
+#include <stdio.h>
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "TopGridData.h"
+#include "LevelHierarchy.h"
+#include "CommunicationUtilities.h"
+#include "ActiveParticle.h"
+
+int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level,
+		      HierarchyEntry **Grids[]);
+
+int ActiveParticleType::DisableParticle(LevelHierarchyEntry *LevelArray[])
+{
+
+  int i, nPart, NumberOfGrids, changedGrid = INT_UNDEFINED, found = FALSE;
+  HierarchyEntry **Grids;
+  
+  NumberOfGrids = GenerateGridArray(LevelArray, this->level, &Grids);
+  for (i = 0; i < NumberOfGrids; i++) {
+    found = Grids[i]->GridData->RemoveParticle(this->Identifier, true);
+    if (found) {
+      changedGrid = i;
+      break;
+    }
+  } // ENDFOR grids
+
+  /* Now clean up deleted particle on the local processor and adjust
+     NumberOfParticles on others */
+
+#ifdef USE_MPI
+  CommunicationAllReduceValues(&changedGrid, 1, MPI_MAX);
+#endif
+
+  if (changedGrid == INT_UNDEFINED) {
+    if (debug)
+      this->PrintInfo();
+    ENZO_VFAIL("DisableParticle: WARNING -- "
+	       "particle %"ISYM" not found...\n", this->Identifier)
+  }
+
+  Grids[changedGrid]->GridData->NumberOfStars--;
+
+  delete [] Grids;
+
+  return SUCCESS;
+
+}

File src/enzo/ActiveParticle_MirrorToParticle.C

+/***********************************************************************
+/
+/  REFLECTS CHANGES IN STAR PARTICLE IN NORMAL PARTICLE
+/
+/  written by: John Wise
+/  date:       March, 2009
+/  modified1:  November, 2011 (JHW) -- porting to active particles
+/
+************************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <map>
+#include <iostream>
+#include <stdexcept>
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "TopGridData.h"
+#include "LevelHierarchy.h"
+#include "ActiveParticle.h"
+
+int GetUnits(float *DensityUnits, float *LengthUnits,
+	     float *TemperatureUnits, float *TimeUnits,
+	     float *VelocityUnits, FLOAT Time);
+
+void ActiveParticleType::MirrorToParticle(void)
+{
+
+  if (CurrentGrid == NULL)
+    return;
+
+  const double Msun = 1.989e33;
+  int i, dim, place = -1;
+  float MassConversion;
+
+  float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits,
+    VelocityUnits;
+  GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
+	   &TimeUnits, &VelocityUnits, CurrentGrid->Time);
+
+  double dx = LengthUnits * CurrentGrid->CellWidth[0][0];
+  MassConversion = (float) (dx*dx*dx * double(DensityUnits) / Msun);
+
+  // Find where this star particle is stored in main arrays
+  for (i = 0; i < CurrentGrid->NumberOfParticles; i++) 
+    if (CurrentGrid->ParticleNumber[i] == this->Identifier) {
+      place = i;
+      break;
+    }
+
+  if (place < 0) {
+    printf("star::MTP: CurrentGrid->NumberOfParticles = %d, level = %d, "
+	   "place =%d, Mass = %d, GridID = %d\n", 
+	   CurrentGrid->NumberOfParticles, level, place, Mass, GridID); 
+    printf("star::MTP: LeftEdge // RightEdge = %"PSYM" %"PSYM" %"PSYM
+	   " // %"PSYM" %"PSYM" %"PSYM"\n",
+	   CurrentGrid->GridLeftEdge[0], CurrentGrid->GridLeftEdge[1], 
+	   CurrentGrid->GridLeftEdge[2], 
+	   CurrentGrid->GridRightEdge[0], CurrentGrid->GridRightEdge[1], 
+	   CurrentGrid->GridRightEdge[2]);
+    this->PrintInfo();
+    ENZO_FAIL("Cannot find matching particle for this star.");
+  }
+  //assert(place >= 0);
+
+  // Change all particle data in favor of updated Star particle
+  for (dim = 0; dim < MAX_DIMENSION; dim++) {
+    CurrentGrid->ParticlePosition[dim][place] = this->pos[dim];
+    CurrentGrid->ParticleVelocity[dim][place] = this->vel[dim];
+  }
+  CurrentGrid->ParticleMass[place] = (float)(this->Mass) / MassConversion;
+  CurrentGrid->ParticleType[place] = this->type;
+  CurrentGrid->ParticleAttribute[0][place] = this->BirthTime;
+  CurrentGrid->ParticleAttribute[1][place] = this->DynamicalTime;
+  CurrentGrid->ParticleAttribute[2][place] = this->Metallicity;
+  
+  return;
+}

File src/enzo/ActiveParticle_SphereContained.C

+/***********************************************************************
+/
+/  FIND ACCRETION SPHERE
+/
+/  written by: John Wise
+/  date:       March, 2009
+/  modified1:  November, 2011 (JHW) -- porting to active particles
+/
+/  PURPOSE: When we remove baryons from the grid to add to the star
+/           particle, look for a sphere that contains twice its mass.
+/           Stepping outward by a cell width.
+/
+************************************************************************/
+#ifdef USE_MPI
+#include "mpi.h"
+#endif /* USE_MPI */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <map>
+#include <iostream>
+#include <stdexcept>
+#include "performance.h"
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "TopGridData.h"
+#include "LevelHierarchy.h"
+#include "CommunicationUtilities.h"
+#include "ActiveParticle.h"
+
+int ActiveParticleType::SphereContained(LevelHierarchyEntry *LevelArray[], int level, 
+					float Radius)
+{
+
+  LevelHierarchyEntry *Temp;
+  int i, dim, direction, cornersContained, Rank, result;
+  bool inside;
+  int cornerDone[8], Dims[MAX_DIMENSION];
+  FLOAT corners[MAX_DIMENSION][8];
+  FLOAT LeftEdge[MAX_DIMENSION], RightEdge[MAX_DIMENSION];
+
+  /**************************************************************
+
+     Compute corners of cube that contains a sphere of r=Radius
+
+   **************************************************************/
+
+  LCAPERF_START("star_SphereContained");
+
+  for (i = 0; i < 8; i++) {
+    for (dim = 0; dim < MAX_DIMENSION; dim++) {
+
+      // If the bit is true, forward.  If not, reverse.
+      direction = (i >> dim & 1) ? 1 : -1;
+      corners[dim][i] = pos[dim] + direction * Radius;
+    }
+    cornerDone[i] = 0;
+  }
+
+  /* Check if the influenced sphere is contained within the grids on
+     this level */
+
+  Temp = LevelArray[level];
+  for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel) {
+    if (Temp->GridData->ReturnProcessorNumber() == MyProcessorNumber) {
+      Temp->GridData->ReturnGridInfo(&Rank, Dims, LeftEdge, RightEdge);
+
+      for (i = 0; i < 8; i++) {
+	if (cornerDone[i]) continue;  // Skip if already locally found
+	inside = true;
+	for (dim = 0; dim < MAX_DIMENSION; dim++)
+	  inside &= (corners[dim][i] >= LeftEdge[dim] &&
+		     corners[dim][i] <= RightEdge[dim]);
+	if (inside)
+	  cornerDone[i] = 1;
+      } // ENDFOR corners
+
+    } // ENDIF MyProcessorNumber == ProcessorNumber
+  } // ENDFOR grids
+
+  /* Take the MPI_MAX of cornerDone flags, then sum them to see if
+     they equal 8.  If so, the sphere is contained within grids on
+     this level. */
+
+#ifdef USE_MPI
+  CommunicationAllReduceValues(cornerDone, 8, MPI_MAX);
+#endif
+
+  cornersContained = 0;
+  for (i = 0; i < 8; i++)
+    cornersContained += cornerDone[i];
+
+  result = (cornersContained == 8);
+  LCAPERF_STOP("star_SphereContained");
+  return result;
+
+}

File src/enzo/Make.config.objects

 OBJS_CONFIG_LIB = \
         acml_st1.o \
         ActiveParticle.o \
+	ActiveParticleRoutines.o \
         ActiveParticle_CenOstriker.o \
+	ActiveParticle_DisableParticle.o \
         ActiveParticle_Kravtsov.o \
         ActiveParticle_PopIII.o \
+	ActiveParticle_MirrorToParticle.o \
         ActiveParticle_SampleParticle.o \
+	ActiveParticle_SphereContained.o \
         ActiveParticle_SpringelHernquist.o \
 	ActiveParticle_SinkParticle.o \
         AdiabaticExpansionInitialize.o \