Commits

Nathan Goldbaum committed 148fe9e

Fixing accretion for the AccretingParticle. Removing some debug code.

Comments (0)

Files changed (9)

src/enzo/control/EvolveLevel.C

  
     /* Finalize (accretion, feedback, etc.) star particles */
 
-    float mymass1 = 0, mass1 = 0;
-
-    for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
-      if (Grids[grid1]->GridData->ReturnProcessorNumber() == MyProcessorNumber) { 
-	Grids[grid1]->GridData->SumGasMass(&mymass1);
-	printf("init  mass on grid %"ISYM" = %"FSYM" \n", grid1, mymass1);
-      }
-    }
-
-    MPI_Allreduce(&mymass1, &mass1, 1, MPI_DOUBLE, MPI_SUM, EnzoTopComm);
-
     ActiveParticleFinalize(Grids, MetaData, NumberOfGrids, LevelArray,
 			   level, NumberOfNewActiveParticles);
 
-    float mymass2 = 0, mass2 = 0;
-
-    for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
-      if (Grids[grid1]->GridData->ReturnProcessorNumber() == MyProcessorNumber) {
-	Grids[grid1]->GridData->SumGasMass(&mymass2);
-	printf("final mass on grid %"ISYM" = %"FSYM" \n", grid1, mymass2);
-      }
-    }
-
-    MPI_Allreduce(&mymass2, &mass2, 1, MPI_DOUBLE, MPI_SUM, EnzoTopComm);
-
-    if (MyProcessorNumber == ROOT_PROCESSOR)
-      printf("mass1 - mass2 = %"FSYM" \n", mass1 - mass2);
-
     /* For each grid: a) interpolate boundaries from the parent grid.
                       b) copy any overlapping zones from siblings. */
  

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);
 
   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();
   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);
 
     KernelRadius = AccretionRadius/2.0;
 
   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++) {
-	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) +
-	  POW((CellLeftEdge[2][k] + 0.5*CellWidth[2][k]) - ParticlePosition[2],2);   
-	if ((AccretionRadius*AccretionRadius) > radius2) {
-	  WeightedSum += BaryonField[DensNum][index]*exp(-radius2/(KernelRadius*KernelRadius)); 
-	  SumOfWeights += exp(-radius2/(KernelRadius*KernelRadius));
-	  NumberOfCells++;
-	  vgas[0] = BaryonField[Vel1Num][index];
-	  vgas[1] = BaryonField[Vel2Num][index];
-	  vgas[2] = BaryonField[Vel3Num][index];
-	  
-	  /* The true accretion rate is somewhat less than this due to
-	     angular momentum conservation.  Subdivide the cell into
-	     NDIV^2 subcells and estimate the reduction assuming
-	     ballistic orbits. See the discussion near Eqn 15. */	  
-	  for (ksub = 0; ksub < NDIV-1; ksub++) {
-	    zdist = CellLeftEdge[2][k] + CellWidth[2][k]*(float(ksub)+0.5)/NDIV - ParticlePosition[2];
-	    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];
-
-		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]);
-		jsp[1] = zdist*(vgas[0] - vsink[0]) - 
-		  xdist*(vgas[2] - vsink[2]);
-		jsp[2] = xdist*(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) +
-		       POW((vgas[2] - vsink[2]),2)) / 2.0 - GravitationalConstant * msink/dist;
-		
-		// Compute distance of closest approach
-		if (esp > 0.0)
-		  rmin = huge*CellWidth[0][0];
-		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);
-
-	}
-      }
-    }
+     for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
+       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) +
+	   POW((CellLeftEdge[2][k] + 0.5*CellWidth[2][k]) - ParticlePosition[2],2);   
+	 if ((AccretionRadius*AccretionRadius) > radius2) {
+	   WeightedSum += BaryonField[DensNum][index]*exp(-radius2/(KernelRadius*KernelRadius)); 
+	   SumOfWeights += exp(-radius2/(KernelRadius*KernelRadius));
+	   NumberOfCells++;
+	   vgas[0] = BaryonField[Vel1Num][index];
+	   vgas[1] = BaryonField[Vel2Num][index];
+	   vgas[2] = BaryonField[Vel3Num][index];
+	   
+	   /* The true accretion rate is somewhat less than this due to
+	      angular momentum conservation.  Subdivide the cell into
+	      NDIV^2 subcells and estimate the reduction assuming
+	      ballistic orbits. See the discussion near Eqn 15. */	  
+	   for (ksub = 0; ksub < NDIV-1; ksub++) {
+	     zdist = CellLeftEdge[2][k] + CellWidth[2][k]*(float(ksub)+0.5)/NDIV - ParticlePosition[2];
+	     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[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]);
+		 jsp[1] = zdist*(vgas[0] - vsink[0]) - 
+		   xdist*(vgas[2] - vsink[2]);
+		 jsp[2] = xdist*(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) +
+			POW((vgas[2] - vsink[2]),2)) / 2.0 - GravitationalConstant * msink/dist;
+		 
+		 // Compute distance of closest approach
+		 if (esp > 0.0)
+		   rmin = huge*CellWidth[0][0];
+		 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) {
   (*ThisParticle)->AddMass(AccretedMass);
   (*ThisParticle)->SetVelocity(NewVelocity);
 
-  printf("AccretedMass = %"FSYM" \n", AccretedMass);
-  printf("FinalMass    = %"FSYM" \n", (*ThisParticle)->ReturnMass());
-
   delete [] nexcluded;
 
   return SUCCESS;
 }
 
-#undef DEBUG
+#undef DEBUG_AP

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

       ParticleList = ActiveParticleFindAll(LevelArray, &nParticles, AccretingParticleID);
 
       /* Do accretion */
-      /*if (Accrete(nParticles,ParticleList,AccretionRadius,dx,LevelArray,ThisLevel) == FAIL)
+      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) {

src/enzo/headers/Grid.h

 			     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/particles/active_particles/ActiveParticle_AccretingParticle.C

   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
     
       float AccretionRate = 0;
       
-      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;
     }