Source

enzo-3.0 / src / enzo / particles / active_particles / ActiveParticle_AccretingParticle.C

The active_particles branch has multiple heads

Full commit
/***********************************************************************
/
/ Accreting Particle
/
************************************************************************/

#include "ActiveParticle_AccretingParticle.h"

#ifdef NEW_CONFIG

#include "ParameterControl/ParameterControl.h"
extern Configuration Param;

/* Set default parameter values. */

const char config_accreting_particle_defaults[] =
"### ACCRETING PARTICLE DEFAULTS ###\n"
"\n"
"Physics: {\n"
"    ActiveParticles: {\n"
"        AccretingParticle: {\n"
"            OverflowFactor       = 1.01;\n"
"            LinkingLength        = 4;\n   "
"            AccretionRadius      = 4;\n   "
"        };\n"
"    };\n"
"};\n";

#endif

/* 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. */

class AccretingParticleGrid : private grid {
  friend class ActiveParticleType_AccretingParticle;
};

/* Note that we only refer to AccretingParticleGrid here. 
 * Given a grid object, we static cast to get this:
 *
 *    AccretingParticleGrid *thisgrid =
 *      static_cast<AccretingParticleGrid *>(thisgrid_orig); */

float ActiveParticleType_AccretingParticle::OverflowFactor = FLOAT_UNDEFINED;
int ActiveParticleType_AccretingParticle::AccretionRadius = INT_UNDEFINED;
int ActiveParticleType_AccretingParticle::LinkingLength = INT_UNDEFINED;

int ActiveParticleType_AccretingParticle::InitializeParticleType()
{
#ifdef NEW_CONFIG

  // Update the parameter config to include the local defaults.  Note
  // that this does not overwrite any values previously specified.
  Param.Update(config_sink_particle_defaults);

  // Retrieve parameters from Param structure
  Param.GetScalar(OverflowFactor, "Physics.ActiveParticles.AccretingParticle.OverflowFactor");
  Param.GetScalar(LinkingLength, "Physics.ActiveParticles.AccretingParticle.LinkingLength");
  Param.GetScalar(AccretionRadius, "Physics.ActiveParticles.AccretingParticle.AccretionRadius");

#else

  // Leaving these defaults hardcoded for testing. NJG

  OverflowFactor = 1.01;
  LinkingLength = 4;
  AccretionRadius = 4;

#endif

  // Need to turn on particle mass flagging if it isn't already turned on.

  bool TurnOnParticleMassRefinement = true;
  int method;
  for (method = 0; method < MAX_FLAGGING_METHODS; method++)
    if (CellFlaggingMethod[method] == 8 || CellFlaggingMethod[method] == 4) {
      TurnOnParticleMassRefinement = false;
      break;
    }
	
  if (TurnOnParticleMassRefinement) {
    method = 0;
    while(CellFlaggingMethod[method] != INT_UNDEFINED)
      method++;
    CellFlaggingMethod[method] = 4;
  }

  /* Add on the Particle Array Handlers */
  typedef ActiveParticleType_AccretingParticle ap;
  AttributeVector &ah = ap::AttributeHandlers;
  ActiveParticleType::SetupBaseParticleAttributes(ah);

  ah.push_back(new Handler<ap, float, &ap::AccretionRate>("AccretionRate"));

  return SUCCESS;
}

int ActiveParticleType_AccretingParticle::EvaluateFormation
(grid *thisgrid_orig, ActiveParticleFormationData &data)
{
  // No need to do the rest if we're not on the maximum refinement level.
  if (data.level != MaximumRefinementLevel)
    return SUCCESS;

  AccretingParticleGrid *thisGrid =
    static_cast<AccretingParticleGrid *>(thisgrid_orig);
  
  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 JeansDensityUnitConversion = (Gamma*pi*kboltz) / (Mu*mh*GravConst);
  float CellTemperature = 0;
  float JeansDensity = 0;
  float MassRefinementDensity = 0;
  float DensityThreshold = huge_number;
  float ExtraDensity = 0;

  int GridDimension[3] = {thisGrid->GridDimension[0],
                          thisGrid->GridDimension[1],
                          thisGrid->GridDimension[2]};

  FLOAT dx = thisGrid->CellWidth[0][0];
  
  bool HasMetalField = (data.MetalNum != -1 || data.ColourNum != -1);
  bool JeansRefinement = false;
  bool MassRefinement = false;

  // determine refinement criteria
  for (method = 0; method < MAX_FLAGGING_METHODS; method++) {
    if (CellFlaggingMethod[method] == 2) {
      MassRefinement = true;
      MassRefinementMethod = method;
    }
    if (CellFlaggingMethod[method] == 6) 
      JeansRefinement = true;
  }

  for (k = thisGrid->GridStartIndex[2]; k <= thisGrid->GridEndIndex[2]; k++) {
    for (j = thisGrid->GridStartIndex[1]; j <= thisGrid->GridEndIndex[1]; j++) {
      for (i = thisGrid->GridStartIndex[0]; i <= thisGrid->GridEndIndex[0]; i++) {
	index = GRIDINDEX_NOGHOST(i, j, k);

	// If no more room for particles, throw an ENZO_FAIL
	if (data.NumberOfNewParticles >= data.MaxNumberOfNewParticles)
	  return FAIL;

	// Does this cell violate the Jeans condition?
	DensityThreshold = huge_number;
	if (JeansRefinement) {
	  CellTemperature = (JeansRefinementColdTemperature > 0) ? JeansRefinementColdTemperature : data.Temperature[index];
	  JeansDensity = JeansDensityUnitConversion * OverflowFactor * CellTemperature / 
	    POW(data.LengthUnits*dx*4.0,2);
	  JeansDensity /= data.DensityUnits;
	  DensityThreshold = min(DensityThreshold,JeansDensity);
	}
	if (DensityThreshold == huge_number)
	  ENZO_FAIL("Error in Accreting Particles: Must refine by jeans length or overdensity!");
	
	if (density[index] <= DensityThreshold) 
	  continue;

	// Passed creation tests, create sink particle

	ActiveParticleType_AccretingParticle *np = new ActiveParticleType_AccretingParticle();
	data.NewParticles[data.NumberOfNewParticles++] = np;

	ExtraDensity = density[index] - DensityThreshold;
	np->Mass = ExtraDensity;   // Particle 'masses' are actually densities
	np->type = np->GetEnabledParticleID();
	np->BirthTime = thisGrid->ReturnTime();

	np->level = data.level;
	np->GridID = data.GridID;
	np->CurrentGrid = thisGrid;
	
	np->pos[0] = thisGrid->CellLeftEdge[0][i] + 0.5*thisGrid->CellWidth[0][i];
	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];
	
	np->vel[0] = velx[index];
	np->vel[1] = vely[index];
	np->vel[2] = velz[index];
	
	if (HasMetalField)
	  np->Metallicity = data.TotalMetals[index];
	else
	  np->Metallicity = 0.0;

	np->AccretionRate = 0.0;
	
	// Remove mass from grid

	density[index] = DensityThreshold;

      } // i
    } // j
  } // k
 
  return SUCCESS;
}  

int ActiveParticleType_AccretingParticle::EvaluateFeedback(grid *thisgrid_orig, ActiveParticleFormationData &data)
{
  AccretingParticleGrid *thisGrid =
    static_cast<AccretingParticleGrid *>(thisgrid_orig);
  
  int npart = thisGrid->NumberOfActiveParticles;

  for (int n = 0; n < npart; n++) {
    ActiveParticleType_AccretingParticle *ThisParticle = 
      static_cast<ActiveParticleType_AccretingParticle*>(thisGrid->ActiveParticles[n]);
  
    ThisParticle->level = data.level;
  }

  return SUCCESS;
}

void ActiveParticleType_AccretingParticle::DescribeSupplementalData(ActiveParticleFormationDataFlags &flags)
{
  flags.DarkMatterDensity = true;
  flags.Temperature = true;
  flags.UnitConversions = true;
  flags.DataFieldNumbers = true;
  flags.MetalField = true;
}


grid* ConstructFeedbackZone(ActiveParticleType* ThisParticle, int FeedbackRadius, 
			     FLOAT dx, 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, 
						  LevelHierarchyEntry *LevelArray[], int ThisLevel)
{
  
  /* Skip accretion if we're not on the maximum refinement level. 
     This should only ever happen right after creation and then
     only in pathological cases where sink creation is happening at 
     the edges of two regions at the maximum refinement level */

  if (ThisLevel < MaximumRefinementLevel)
    return SUCCESS;

  /* For each particle, loop over all of the grids and do accretion 
     if the grid overlaps with the accretion zone                   */
  
  int i, NumberOfGrids;
  int *FeedbackRadius = NULL;
  HierarchyEntry **Grids = NULL;
  grid *sinkGrid = NULL;
  
  bool SinkIsOnThisProc, SinkIsOnThisGrid;
  
  float SubtractedMass, SubtractedMomentum[3] = {};
  
  NumberOfGrids = GenerateGridArray(LevelArray, ThisLevel, &Grids);
  
  for (i = 0; i < nParticles; i++) {
    grid* FeedbackZone = ConstructFeedbackZone(ParticleList[i], AccretionRadius, 
					       dx, Grids, NumberOfGrids, ALL_FIELDS);

    if (MyProcessorNumber == FeedbackZone->ReturnProcessorNumber()) {
    
      float AccretionRate = 0;
      
      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;
  }
  
  if (AssignActiveParticlesToGrids(ParticleList, nParticles, LevelArray) == FAIL)
    return FAIL;

  delete [] Grids;
  return SUCCESS;
}

int ActiveParticleType_AccretingParticle::SetFlaggingField(LevelHierarchyEntry *LevelArray[], int level, 
							   int TopGridDims[], int AccretingParticleID)
{
  /* Generate a list of all sink particles in the simulation box */
  int i, nParticles;
  FLOAT *pos = NULL, dx=0;
  ActiveParticleType **AccretingParticleList = NULL ;
  LevelHierarchyEntry *Temp = NULL;
  
  AccretingParticleList = ActiveParticleFindAll(LevelArray, &nParticles, AccretingParticleID);
  
  /* Calculate CellWidth on maximum refinement level */
  
  // this will fail for noncubic boxes or simulations with MinimimMassForRefinementLevelExponent
  dx = (DomainRightEdge[0] - DomainLeftEdge[0]) /
    (TopGridDims[0]*POW(FLOAT(RefineBy),FLOAT(MaximumRefinementLevel)));
  
  for (i=0 ; i<nParticles; i++){
    pos = AccretingParticleList[i]->ReturnPosition();
    for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel)
      if (Temp->GridData->DepositRefinementZone(level,pos,AccretionRadius*dx) == FAIL) {
	ENZO_FAIL("Error in grid->DepositRefinementZone.\n")
	  }
  }

  if (NumberOfProcessors > 1)
    for (i = 0; i < nParticles; i++) {
      delete AccretingParticleList[i];
      AccretingParticleList[i] = NULL;
    }

  delete [] AccretingParticleList;

  return SUCCESS;
}

namespace {
  ActiveParticleType_info *AccretingParticleInfo = 
    register_ptype <ActiveParticleType_AccretingParticle> 
    ("AccretingParticle");
}

std::vector<ParticleAttributeHandler*>
  ActiveParticleType_AccretingParticle::AttributeHandlers;